#NEWLINE# def main( argv = None ) : #NEWLINE# #INDENT# 'script main.\n\n    parses command line options in sys.argv, unless *argv* is given.\n    ' #NEWLINE# if ( not argv ) : #NEWLINE# #INDENT# argv = sys . argv #NEWLINE# #UNINDENT# parser = E . OptionParser ( version = '%prog version: $Id$' , usage = globals ( ) [ '__doc__' ] ) #NEWLINE# parser . add_option ( '-g' , '--gtf-file' , dest = 'filename_gtf' , type = 'string' , help = 'filename with gene models in gtf format [%default]' ) #NEWLINE# parser . add_option ( '-m' , '--filename-mismapped' , dest = 'filename_mismapped' , type = 'string' , help = 'output bam file for mismapped reads [%default]' ) #NEWLINE# parser . add_option ( '-j' , '--junctions-bed-file' , dest = 'filename_junctions' , type = 'string' , help = 'bam file with reads mapped across junctions [%default]' ) #NEWLINE# parser . add_option ( '-r' , '--filename-regions' , dest = 'filename_regions' , type = 'string' , help = 'filename with regions to remove in bed format [%default]' ) #NEWLINE# parser . add_option ( '-t' , '--transcripts-gtf-file' , dest = 'filename_transcriptome' , type = 'string' , help = 'bam file with reads mapped against transcripts [%default]' ) #NEWLINE# parser . add_option ( '-p' , '--map-tsv-file' , dest = 'filename_map' , type = 'string' , help = 'filename mapping transcript numbers (used by --filename-transciptome) to transcript names (used by --filename-gtf) [%default]' ) #NEWLINE# parser . add_option ( '-s' , '--filename-stats' , dest = 'filename_stats' , type = 'string' , help = 'filename to output stats to [%default]' ) #NEWLINE# parser . add_option ( '-o' , '--colour' , dest = 'colour_mismatches' , action = 'store_true' , help = 'mismatches will use colour differences (CM tag) [%default]' ) #NEWLINE# parser . add_option ( '-i' , '--ignore-mismatches' , dest = 'ignore_mismatches' , action = 'store_true' , help = 'ignore mismatches [%default]' ) #NEWLINE# parser . add_option ( '-c' , '--remove-contigs' , dest = 'remove_contigs' , type = 'string' , help = "','-separated list of contigs to remove [%default]" ) #NEWLINE# parser . add_option ( '-f' , '--force-output' , dest = 'force' , action = 'store_true' , help = 'force overwriting of existing files [%default]' ) #NEWLINE# parser . add_option ( '-u' , '--unique' , dest = 'unique' , action = 'store_true' , help = 'remove reads not matching uniquely [%default]' ) #NEWLINE# parser . add_option ( '--output-sam' , dest = 'output_sam' , action = 'store_true' , help = 'output in sam format [%default]' ) #NEWLINE# parser . set_defaults ( filename_gtf = None , filename_mismapped = None , filename_junctions = None , filename_transcriptome = None , filename_map = None , remove_contigs = None , force = False , unique = False , colour_mismatches = False , ignore_mismatches = False , output_sam = False , filename_table = None ) #NEWLINE# ( options , args ) = E . Start ( parser , argv = argv ) #NEWLINE# if ( len ( args ) != 1 ) : #NEWLINE# #INDENT# raise ValueError ( 'please supply one bam file' ) #NEWLINE# #UNINDENT# bamfile_genome = args [ 0 ] #NEWLINE# genome_samfile = pysam . Samfile ( bamfile_genome , 'rb' ) #NEWLINE# if options . remove_contigs : #NEWLINE# #INDENT# options . remove_contigs = options . remove_contigs . split ( ',' ) #NEWLINE# #UNINDENT# if options . filename_map : #NEWLINE# #INDENT# E . info ( 'reading map' ) #NEWLINE# id_map = IOTools . readMap ( IOTools . openFile ( options . filename_map ) , has_header = True ) #NEWLINE# id_map = dict ( [ ( regions_to_remove , x ) for ( x , y ) in id_map . iteritems ( ) ] ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# id_map = None #NEWLINE# #UNINDENT# transcripts = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# if options . filename_gtf : #NEWLINE# #INDENT# E . info ( 'indexing geneset' ) #NEWLINE# ( mapped , missed ) = ( 0 , 0 ) #NEWLINE# for gtf in GTF . transcript_iterator ( GTF . iterator ( IOTools . openFile ( options . filename_gtf ) ) ) : #NEWLINE# #INDENT# gtf . sort ( key = ( lambda x : x . start ) ) #NEWLINE# transcript_id = gtf [ 0 ] . transcript_id #NEWLINE# if id_map : #NEWLINE# #INDENT# try : #NEWLINE# #INDENT# transcript_id = id_map [ transcript_id ] #NEWLINE# mapped += 1 #NEWLINE# #UNINDENT# except KeyError : #NEWLINE# #INDENT# missed += 1 #NEWLINE# continue #NEWLINE# #UNINDENT# #UNINDENT# transcripts [ transcript_id ] = gtf #NEWLINE# #UNINDENT# E . info ( ( 'read %i transcripts from geneset (%i mapped, %i missed)' % ( len ( transcripts ) , mapped , missed ) ) ) #NEWLINE# #UNINDENT# regions_to_remove = None #NEWLINE# if options . filename_regions : #NEWLINE# #INDENT# E . info ( 'indexing regions' ) #NEWLINE# regions_to_remove = IndexedGenome . Simple ( ) #NEWLINE# for bed in Bed . iterator ( IOTools . openFile ( options . filename_regions ) ) : #NEWLINE# #INDENT# regions_to_remove . add ( bed . contig , bed . start , bed . end ) #NEWLINE# #UNINDENT# E . info ( ( 'read %i regions' % len ( regions_to_remove ) ) ) #NEWLINE# #UNINDENT# if options . filename_transcriptome : #NEWLINE# #INDENT# transcripts_samfile = pysam . Samfile ( options . filename_transcriptome , 'rb' ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# transcripts_samfile = None #NEWLINE# #UNINDENT# if options . output_sam : #NEWLINE# #INDENT# output_samfile = pysam . Samfile ( '-' , 'wh' , template = genome_samfile ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# output_samfile = pysam . Samfile ( '-' , 'wb' , template = genome_samfile ) #NEWLINE# #UNINDENT# if options . filename_mismapped : #NEWLINE# #INDENT# if ( ( not options . force ) and os . path . exists ( options . filename_mismapped ) ) : #NEWLINE# #INDENT# raise IOError ( ( 'output file %s already exists' % options . filename_mismapped ) ) #NEWLINE# #UNINDENT# output_mismapped = pysam . Samfile ( options . filename_mismapped , 'wb' , template = genome_samfile ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# output_mismapped = None #NEWLINE# #UNINDENT# if options . filename_junctions : #NEWLINE# #INDENT# junctions_samfile = pysam . Samfile ( options . filename_junctions , 'rb' ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# junctions_samfile = None #NEWLINE# #UNINDENT# c = _bams2bam . filter ( genome_samfile , output_samfile , output_mismapped , transcripts_samfile , junctions_samfile , transcripts , regions = regions_to_remove , unique = options . unique , remove_contigs = options . remove_contigs , colour_mismatches = options . colour_mismatches , ignore_mismatches = options . ignore_mismatches , ignore_transcripts = ( transcripts_samfile is None ) , ignore_junctions = ( junctions_samfile is None ) ) #NEWLINE# if options . filename_stats : #NEWLINE# #INDENT# outf = IOTools . openFile ( options . filename_stats , 'w' ) #NEWLINE# outf . write ( ( 'category\tcounts\n%s\n' % c . asTable ( ) ) ) #NEWLINE# outf . close ( ) #NEWLINE# #UNINDENT# if options . filename_transcriptome : #NEWLINE# #INDENT# transcripts_samfile . close ( ) #NEWLINE# #UNINDENT# genome_samfile . close ( ) #NEWLINE# output_samfile . close ( ) #NEWLINE# if output_mismapped : #NEWLINE# #INDENT# output_mismapped . close ( ) #NEWLINE# #UNINDENT# E . Stop ( )