Mercurial > repos > fubar > htseq_bams_to_count_matrix
changeset 39:cc4efb5bd07e draft
try iuc dependencies from main toolshed
| author | fubar |
|---|---|
| date | Tue, 19 Nov 2013 20:52:13 -0500 (2013-11-20) |
| parents | 7a4061baef5c |
| children | 7da5a9e8ef58 |
| files | htseq_bams_to_count_matrix/htseqsams2mx.py htseq_bams_to_count_matrix/htseqsams2mx.xml htseq_bams_to_count_matrix/tool_dependencies.xml |
| diffstat | 3 files changed, 117 insertions(+), 71 deletions(-) [+] |
line wrap: on
line diff
--- a/htseq_bams_to_count_matrix/htseqsams2mx.py Tue Nov 19 20:25:26 2013 -0500 +++ b/htseq_bams_to_count_matrix/htseqsams2mx.py Tue Nov 19 20:52:13 2013 -0500 @@ -60,7 +60,7 @@ self.msg = msg -def htseqMX(gff_filename,sam_filenames,colnames,opts): +def htseqMX(gff_filename,sam_filenames,colnames,sam_exts,sam_bais,opts): """ Code taken from count.py in Simon Anders HTSeq distribution Wrapped in a loop to accept multiple bam/sam files and their names from galaxy to @@ -90,11 +90,13 @@ features = HTSeq.GenomicArrayOfSets( "auto", opts.stranded != "no" ) mapqMin = int(opts.mapqMin) counts = {} + nreads = 0 empty = 0 ambiguous = 0 notaligned = 0 lowqual = 0 nonunique = 0 + filtered = 0 # new filter_extras - need a better way to do this - independent filter tool? gff = HTSeq.GFF_Reader( gff_filename ) try: for i,f in enumerate(gff): @@ -102,14 +104,17 @@ try: feature_id = f.attr[ opts.id_attribute ] except KeyError: - sys.exit( "Feature at row %d %s does not contain a '%s' attribute" % - ( (i+1), f.name, opts.id_attribute ) ) + try: + feature_id = f.attr[ 'gene_id' ] + except KeyError: + sys.exit( "Feature at row %d %s does not contain a '%s' attribute OR a gene_id attribute - faulty GFF?" % + ( (i+1), f.name, opts.id_attribute ) ) if opts.stranded != "no" and f.iv.strand == ".": sys.exit( "Feature %s at %s does not have strand information but you are " "running htseq-count in stranded mode. Use '--stranded=no'." % ( f.name, f.iv ) ) features[ f.iv ] += feature_id - counts[ f.attr[ opts.id_attribute ] ] = [0 for x in colnames] # we use sami as an index here to bump counts later + counts[ feature_id ] = [0 for x in colnames] # we use sami as an index here to bump counts later except: sys.stderr.write( "Error occured in %s.\n" % gff.get_line_number_string() ) raise @@ -121,19 +126,29 @@ sys.stdout.write( "Warning: No features of type '%s' found.\n" % opts.feature_type ) for sami,sam_filename in enumerate(sam_filenames): colname = colnames[sami] - isbam = colname.endswith('.bam') + isbam = sam_exts[sami] == 'bam' + hasbai = sam_bais[sami] > '' + if hasbai: + tempname = os.path.splitext(os.path.basename(sam_filename))[0] + tempbam = '%s.bam' % tempname + tempbai = '%s.bai' % tempname + os.link(sam_filename,tempbam) + os.link(sam_bais[sami],tempbai) try: if isbam: - read_seq = HTSeq.BAM_Reader( sam_filename ) + if hasbai: + read_seq = HTSeq.BAM_Reader ( tempbam ) + else: + read_seq = HTSeq.BAM_Reader( sam_filename ) else: read_seq = HTSeq.SAM_Reader( sam_filename ) first_read = iter(read_seq).next() pe_mode = first_read.paired_end except: if isbam: - sys.stderr.write( "Error occured when reading first line of bam file %s colname=%s \n" % (sam_filename,colname) ) + print >> sys.stderr, "Error occured when reading first line of bam file %s colname=%s \n" % (sam_filename,colname ) else: - sys.stderr.write( "Error occured when reading first line of sam file %s colname=%s \n" % (sam_filename,colname )) + print >> sys.stderr, "Error occured when reading first line of sam file %s colname=%s \n" % (sam_filename,colname ) raise try: @@ -141,11 +156,17 @@ read_seq_pe_file = read_seq read_seq = HTSeq.pair_SAM_alignments( read_seq ) for seqi,r in enumerate(read_seq): + nreads += 1 if not pe_mode: if not r.aligned: notaligned += 1 continue try: + if len(opts.filter_extras) > 0: + for extra in opts.filter_extras: + if r.optional_field(extra): + filtered += 1 + continue if r.optional_field( "NH" ) > 1: nonunique += 1 continue @@ -222,10 +243,10 @@ else: rr = r[0] if r[0] is not None else r[1] empty += 1 - #if not quiet: - # sys.stderr.write( ( "Warning: Skipping read '%s', because chromosome " + - # "'%s', to which it has been aligned, did not appear in the GFF file.\n" ) % - # ( rr.read.name, iv.chrom ) ) + if not opts.quiet: + sys.stdout.write( ( "Warning: Skipping read '%s', because chromosome " + + "'%s', to which it has been aligned, did not appear in the GFF file.\n" ) % + ( rr.read.name, iv.chrom ) ) except: if not pe_mode: sys.stderr.write( "Error occured in %s.\n" % read_seq.get_line_number_string() ) @@ -234,13 +255,13 @@ raise if not opts.quiet: - sys.stdout.write( "%d sam %s processed.\n" % ( seqi, "lines " if not pe_mode else "line pairs" ) ) - return counts,empty,ambiguous,lowqual,notaligned,nonunique + sys.stdout.write( "%d sam %s processed for %s.\n" % ( seqi, "lines " if not pe_mode else "line pairs", colname ) ) + return counts,empty,ambiguous,lowqual,notaligned,nonunique,filtered,nreads warnings.showwarning = my_showwarning assert os.path.isfile(gff_filename),'## unable to open supplied gff file %s' % gff_filename try: - counts,empty,ambiguous,lowqual,notaligned,nonunique = count_reads_in_features( sam_filenames, colnames, gff_filename,opts) + counts,empty,ambiguous,lowqual,notaligned,nonunique,filtered,nreads = count_reads_in_features( sam_filenames, colnames, gff_filename,opts) except: sys.stderr.write( "Error: %s\n" % str( sys.exc_info()[1] ) ) sys.stderr.write( "[Exception type: %s, raised in %s:%d]\n" % @@ -248,7 +269,7 @@ os.path.basename(traceback.extract_tb( sys.exc_info()[2] )[-1][0]), traceback.extract_tb( sys.exc_info()[2] )[-1][1] ) ) sys.exit( 1 ) - return counts,empty,ambiguous,lowqual,notaligned,nonunique + return counts,empty,ambiguous,lowqual,notaligned,nonunique,filtered,nreads def usage(): @@ -278,6 +299,7 @@ op.add_option('-o', '--outfname', default=None) op.add_option('-f','--forceName', default="false") op.add_option('--samf', default=[], action="append") + op.add_option('--filter_extras', default=[], action="append") op.add_option('--mapqMin', default='0') op.add_option( "-t", "--type", type="string", dest="featuretype", default = "exon", help = "feature type (3rd column in GFF file) to be used, " + @@ -297,24 +319,43 @@ outfname = opts.outfname sam_filenames = [] colnames = [] - samdat = opts.samf - samf = [x.split(',')[0].replace("'",'').replace('"','') for x in samdat] # get rid of wrapper supplied quotes - assert len(set(samf)) == len(samf),'## ERROR sams2mx: Duplicate input sam file in %s' % ','.join(samf) - scolnames = [x.split(',')[1].replace("'",'').replace('"','') for x in samdat] - assert len(samf) == len(scolnames), '##ERROR sams2mx: Count of sam/cname not consistent - %d/%d' % (len(samf),len(scolname)) - for i,b in enumerate(samf): - if b != 'None': - assert os.path.isfile(b),'## Supplied input sam file "%s" not found' % b - sam_filenames.append(b) - sampName = scolnames[i] # better be unique - if sampName == "": - sampName = b # for test - sampName = sampName.replace('#','') # for R - sampName = sampName.replace('(','') # for R - sampName = sampName.replace(')','') # for R - sampName = sampName.replace(' ','_') # for R - colnames.append(sampName) - counts,empty,ambiguous,lowqual,notaligned,nonunique = htseqMX(gff_file, sam_filenames,colnames,opts) + samf = opts.samf + samfsplit = [x.split(',') for x in samf] # one per samf set + samsets = [] + for samfs in samfsplit: + samset = [x.replace("'","") for x in samfs] + samset = [x.replace('"','') for x in samset] + samsets.append(samset) + samsets = [x for x in samsets if x[0].lower() != 'none'] + # just cannot stop getting these on cl! wtf in cheetah for a repeat group? + samfnames = [x[0] for x in samsets] + if len(set(samfnames)) != len(samfnames): + samnames = [] + delme = [] + for i,s in enumerate(samfnames): + if s in samnames: + delme.append(i) + print sys.stdout,'## WARNING htseqsams2mx: Duplicate input sam file %s in %s - ignoring dupe in 0 based position %s' %\ + (s,','.join(samfnames), str(delme)) + else: + samnames.append(s) # first time + samsets = [x for i,x in enumerate(samsets) if not (i in delme)] + samfnames = [x[0] for x in samsets] + scolnames = [x[1]for x in samsets] + assert len(samfnames) == len(scolnames), '##ERROR sams2mx: Count of sam/cname not consistent - %d/%d' % (len(samfnames),len(scolnames)) + sam_exts = [x[2] for x in samsets] + assert len(samfnames) == len(sam_exts), '##ERROR sams2mx: Count of extensions not consistent - %d/%d' % (len(samfnames),len(sam_exts)) + sam_bais = [x[3] for x in samsets] # these only exist for bams and need to be finessed with a symlink so pysam will just work + for i,b in enumerate(samfnames): + assert os.path.isfile(b),'## Supplied input sam file "%s" not found' % b + sam_filenames.append(b) + sampName = scolnames[i] # better be unique + sampName = sampName.replace('#','') # for R + sampName = sampName.replace('(','') # for R + sampName = sampName.replace(')','') # for R + sampName = sampName.replace(' ','_') # for R + colnames.append(sampName) + counts,empty,ambiguous,lowqual,notaligned,nonunique,filtered,nreads = htseqMX(gff_file, sam_filenames,colnames,sam_exts,sam_bais,opts) heads = '\t'.join(['Contig',] + colnames) res = [heads,] contigs = counts.keys() @@ -334,9 +375,10 @@ outf.write('\n') outf.close() walltime = int(time.time() - starttime) - accumulatornames = ('walltimeseconds','totreadscounted','ncontigs','emptyreads','ambiguousreads','lowqualreads', - 'notalignedreads','nonuniquereads','emptycontigs') - accums = (walltime,totalc,len(contigs),empty,ambiguous,lowqual,notaligned,nonunique,emptycontigs) - notes = ['%s=%d' % (accumulatornames[i],x) for i,x in enumerate(accums)] - print >> sys.stdout, ','.join(notes) + accumulatornames = ('walltime (seconds)','total reads read','total reads counted','number of contigs','total empty reads','total ambiguous reads','total low quality reads', + 'total not aligned reads','total not unique mapping reads','extra filtered reads','empty contigs') + accums = (walltime,nreads,totalc,len(contigs),empty,ambiguous,lowqual,notaligned,nonunique,filtered,emptycontigs) + fracs = (1.0,1.0,float(totalc)/nreads,1.0,float(empty)/nreads,float(ambiguous)/nreads,float(lowqual)/nreads,float(notaligned)/nreads,float(nonunique)/nreads,float(filtered)/nreads,float(emptycontigs)/len(contigs)) + notes = ['%s = %d (%2.3f)' % (accumulatornames[i],x,100.0*fracs[i]) for i,x in enumerate(accums)] + print >> sys.stdout, '\n'.join(notes) sys.exit(0)
--- a/htseq_bams_to_count_matrix/htseqsams2mx.xml Tue Nov 19 20:25:26 2013 -0500 +++ b/htseq_bams_to_count_matrix/htseqsams2mx.xml Tue Nov 19 20:52:13 2013 -0500 @@ -1,5 +1,5 @@ -<tool id="htseqsams2mx" name="SAM/BAM to count matrix" version="0.4"> - <description>using HTSeq code</description> +<tool id="htseqsams2mxlocal" name="SAM/BAM to count matrix" version="0.4"> + <description>using HTSeq code</description> <stdio> <regex match=".*" source="both" level="warning" description="chatter from HTSeq:"/> </stdio> @@ -12,21 +12,18 @@ </requirements> <command interpreter="python"> htseqsams2mx.py -g "$gfffile" -o "$outfile" -m "$model" --id_attribute "$id_attr" --feature_type "$feature_type" - --samf "'$firstsamf','${firstsamf.name}'" - #if $secondsamf: - --samf "'$secondsamf','${secondsamf.name}'" - #end if - #if $thirdsamf: - --samf "'$thirdsamf','${thirdsamf.name}'" - #end if - #if $fourthsamf: - --samf "'$fourthsamf','${fourthsamf.name}'" + --mapqMin $mapqMin --samf "'${firstsamf}','${firstsamf.name}','${firstsamf.ext}','${firstsamf.metadata.bam_index}'" + #if $secondsamf.ext != 'data': + --samf "'${secondsamf}','${secondsamf.name}','${secondsamf.ext}','${secondsamf.metadata.bam_index}'" #end if #for $s in $samfiles: - #if $s.samf: - --samf "'${s.samf}','${s.samf.name}'" + #if $s.samf.ext != 'data': + --samf "'${s.samf}','${s.samf.name}','${s.samf.ext}','${s.samf.metadata.bam_index}'" #end if #end for + #if $filter_extras: + --filter_extras "$filter_extras" + #end if </command> <inputs> <param format="gtf" name="gfffile" type="data" label="Gene model (GFF) file to count reads over from your current history" size="100" /> @@ -56,16 +53,23 @@ <option value="UTR">UTR</option> <option value="transcript">transcript</option> </param> - <param name="firstsamf" type="data" label="First bam/sam file from your history to count reads overlapping gene model regions" - format="sam,bam" optional="false"/> - <param name="secondsamf" type="data" label="Another bam/sam file from your history to count reads overlapping gene model regions" - format="sam,bam" optional="false"/> - <param name="thirdsamf" type="data" label="Another bam/sam file from your history to count reads overlapping gene model regions" - format="sam,bam" optional="true"/> - <param name="fourthsamf" type="data" label="Another bam/sam file from your history to count reads overlapping gene model regions" - format="sam,bam" optional="true"/> - <repeat name="samfiles" title="Use this to add all needed additional bam/sam files from your history to count reads overlapping gene model regions"> - <param name="samf" type="data" label="Additional bam/sam file from your history" format="sam,bam" size="100" optional="true"/> + <param name="filter_extras" type="select" label="Filter any read with one or more flags" + help="eg the XS tag created by bowtie for multiple reads" optional="true" mutliple="true"> + <option value="">None</option> + <option value="XS">XS:i > 0 - More than one mapping position Bowtie</option> + <option value="XS:A">Might be useful for tophat</option> + </param> + + <param name="firstsamf" type="data" label="bam/sam file from your history" format="sam,bam" size="100" + help="Each sam/bam contributes a column of read counts overlapping the specified gene model contigs" + optional="false"/> + <param name="secondsamf" type="data" label="Additional bam/sam file from your history" format="sam,bam" size="100" + help="Each sam/bam contributes a column of read counts overlapping the specified gene model contigs" + optional="false"/> + <repeat name="samfiles" min="16" + title="Specify additional bam/sam file inputs" help="Each sam/bam contributes a column of read counts overlapping the specified gene model contigs"> + <param name="samf" type="data" label="Additional bam/sam file from your history" format="sam,bam" size="100" + optional="true"/> </repeat> </inputs> <outputs> @@ -136,7 +140,7 @@ was written by Ross Lazarus and is licensed to you under the LGPL_ like other rgenetics artefacts -Sorry, don't use so can't be buggered with read groups - contributions welcome - send code +Sorry, I don't use readgroups so had no reason to code read groups. Contributions welcome. Send code .. _LGPL: http://www.gnu.org/copyleft/lesser.html .. _HTSeq: http://www-huber.embl.de/users/anders/HTSeq/doc/index.html
--- a/htseq_bams_to_count_matrix/tool_dependencies.xml Tue Nov 19 20:25:26 2013 -0500 +++ b/htseq_bams_to_count_matrix/tool_dependencies.xml Tue Nov 19 20:52:13 2013 -0500 @@ -1,16 +1,16 @@ <?xml version="1.0"?> <tool_dependency> <package name="numpy" version="1.7"> - <repository changeset_revision="84125ffacb90" name="package_numpy_1_7" owner="iuc" prior_installation_required="True" toolshed="http://testtoolshed.g2.bx.psu.edu" /> + <repository changeset_revision="84125ffacb90" name="package_numpy_1_7" owner="iuc" prior_installation_required="True" toolshed="http://toolshed.g2.bx.psu.edu" /> </package> <package name="pysam" version="0.7.6"> - <repository changeset_revision="247e5e5bee87" name="package_pysam_0_7_6" owner="iuc" prior_installation_required="True" toolshed="http://testtoolshed.g2.bx.psu.edu" /> + <repository changeset_revision="247e5e5bee87" name="package_pysam_0_7_6" owner="iuc" prior_installation_required="True" toolshed="http://toolshed.g2.bx.psu.edu" /> </package> <package name="freetype" version="2.4"> - <repository changeset_revision="fe5cfaf931ff" name="package_freetype_2_4" owner="iuc" prior_installation_required="True" toolshed="http://testtoolshed.g2.bx.psu.edu/" /> + <repository changeset_revision="fe5cfaf931ff" name="package_freetype_2_4" owner="iuc" prior_installation_required="True" toolshed="http://toolshed.g2.bx.psu.edu/" /> </package> <package name="matplotlib" version="1.2"> - <repository changeset_revision="966f29c955b9" name="package_matplotlib_1_2" owner="iuc" prior_installation_required="True" toolshed="http://testtoolshed.g2.bx.psu.edu/" /> + <repository changeset_revision="966f29c955b9" name="package_matplotlib_1_2" owner="iuc" prior_installation_required="True" toolshed="http://toolshed.g2.bx.psu.edu/" /> </package> <package name="htseq" version="0.5.4p3"> <install version="1.0"> @@ -19,16 +19,16 @@ <action type="make_directory">$INSTALL_DIR/lib/python</action> <!-- Not sure why these must be made apriori, but install fails otherwise --> <action type="make_directory">$INSTALL_DIR/lib64/python</action> <!-- Not sure why these must be made apriori, but install fails otherwise --> <action type="set_environment_for_install"> - <repository changeset_revision="84125ffacb90" name="package_numpy_1_7" owner="iuc" prior_installation_required="True" toolshed="http://testtoolshed.g2.bx.psu.edu"> + <repository changeset_revision="84125ffacb90" name="package_numpy_1_7" owner="iuc" prior_installation_required="True" toolshed="http://toolshed.g2.bx.psu.edu"> <package name="numpy" version="1.7" /> </repository> - <repository changeset_revision="247e5e5bee87" name="package_pysam_0_7_6" owner="iuc" prior_installation_required="True" toolshed="http://testtoolshed.g2.bx.psu.edu"> + <repository changeset_revision="247e5e5bee87" name="package_pysam_0_7_6" owner="iuc" prior_installation_required="True" toolshed="http://toolshed.g2.bx.psu.edu"> <package name="pysam" version="0.7.6" /> </repository> - <repository changeset_revision="fe5cfaf931ff" name="package_freetype_2_4" owner="iuc" prior_installation_required="True" toolshed="http://testtoolshed.g2.bx.psu.edu/"> + <repository changeset_revision="fe5cfaf931ff" name="package_freetype_2_4" owner="iuc" prior_installation_required="True" toolshed="http://toolshed.g2.bx.psu.edu/"> <package name="freetype" version="2.4" /> </repository> - <repository changeset_revision="966f29c955b9" name="package_matplotlib_1_2" owner="iuc" prior_installation_required="True" toolshed="http://testtoolshed.g2.bx.psu.edu/"> + <repository changeset_revision="966f29c955b9" name="package_matplotlib_1_2" owner="iuc" prior_installation_required="True" toolshed="http://toolshed.g2.bx.psu.edu/"> <package name="matplotlib_1_2" version="1.2" /> </repository> </action>