# HG changeset patch # User fubar # Date 1384912333 18000 # Node ID cc4efb5bd07ef36c5261028eb69aeedb0fd4ffcf # Parent 7a4061baef5c8fe47c263cbc0c1370351528f496 try iuc dependencies from main toolshed diff -r 7a4061baef5c -r cc4efb5bd07e htseq_bams_to_count_matrix/htseqsams2mx.py --- 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) diff -r 7a4061baef5c -r cc4efb5bd07e htseq_bams_to_count_matrix/htseqsams2mx.xml --- 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 @@ - - using HTSeq code + + using HTSeq code @@ -12,21 +12,18 @@ 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 @@ -56,16 +53,23 @@ - - - - - - + + + + + + + + + + @@ -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 diff -r 7a4061baef5c -r cc4efb5bd07e htseq_bams_to_count_matrix/tool_dependencies.xml --- 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 @@ - + - + - + - + @@ -19,16 +19,16 @@ $INSTALL_DIR/lib/python $INSTALL_DIR/lib64/python - + - + - + - +