#NEWLINE# def main( argv = None ) : #NEWLINE# #INDENT# if ( not argv ) : #NEWLINE# #INDENT# argv = sys . argv #NEWLINE# #UNINDENT# parser = E . OptionParser ( version = '%prog version: $Id$' , usage = globals ( ) [ '__doc__' ] ) #NEWLINE# parser . add_option ( '--merge-exons-distance' , dest = 'merge_exons_distance' , type = 'int' , help = 'distance in nucleotides between exons to be merged [%default].' ) #NEWLINE# parser . add_option ( '--pattern-identifier' , dest = 'pattern' , type = 'string' , help = 'pattern to use for renaming genes/transcripts. The pattern should contain a %i, for example --pattern-identifier=ENSG%010i [%default].' ) #NEWLINE# parser . add_option ( '--sort-order' , dest = 'sort_order' , type = 'choice' , choices = ( 'gene' , 'gene+transcript' , 'transcript' , 'position' , 'contig+gene' , 'position+gene' , 'gene+position' , 'gene+exon' ) , help = 'sort input data [%default].' ) #NEWLINE# parser . add_option ( '--mark-utr' , dest = 'mark_utr' , action = 'store_true' , help = 'mark utr for method --merge-exons. [%default].' ) #NEWLINE# parser . add_option ( '--without-utr' , dest = 'with_utr' , action = 'store_false' , help = 'exclude UTR in methods --merge-exons, merge-transcripts and intersect-transripts. Setting this option will remove non-coding transcripts. [%default].' ) #NEWLINE# parser . add_option ( '--filter-method' , dest = 'filter_method' , type = 'choice' , choices = ( 'gene' , 'transcript' , 'longest-gene' , 'longest-transcript' , 'representative-transcript' , 'proteincoding' , 'lincrna' ) , help = "Filter method to apply. Available filters are: 'gene': filter by gene_id given in ``--map-tsv-file``, 'transcript': filter by transcript_id given in ``--map-tsv-file``, 'longest-gene': output the longest gene for overlapping genes ,'longest-transcript': output the longest transcript per gene,'representative-transcript': output the representative transcript per gene. The representative transcript is the transcript that shares most exons with other transcripts in a gene. The input needs to be sorted by gene. 'proteincoding': only output protein coding features. 'lincrna': only output lincRNA features. [%default]." ) #NEWLINE# parser . add_option ( '-a' , '--map-tsv-file' , dest = 'filename_filter' , type = 'string' , metavar = 'tsv' , help = 'filename of ids to map/filter [%default].' ) #NEWLINE# parser . add_option ( '--gff-file' , dest = 'filename_gff' , type = 'string' , metavar = 'GFF' , help = 'second filename of features (see --remove-overlapping) [%default]' ) #NEWLINE# parser . add_option ( '--invert-filter' , dest = 'invert_filter' , action = 'store_true' , help = 'when using --filter, invert selection (like grep -v). [%default].' ) #NEWLINE# parser . add_option ( '--sample-size' , dest = 'sample_size' , type = 'int' , help = "extract a random sample of size # if the option '--method=filter --filter-method' is set [%default]." ) #NEWLINE# parser . add_option ( '--intron-min-length' , dest = 'intron_min_length' , type = 'int' , help = 'minimum length for introns (for --exons-file2introns) [%default].' ) #NEWLINE# parser . add_option ( '--min-exons-length' , dest = 'min_exons_length' , type = 'int' , help = 'minimum length for gene (sum of exons) (--sam-fileple-size) [%default].' ) #NEWLINE# parser . add_option ( '--intron-border' , dest = 'intron_border' , type = 'int' , help = 'number of residues to exclude at intron at either end (--exons-file2introns) [%default].' ) #NEWLINE# parser . add_option ( '--ignore-strand' , dest = 'ignore_strand' , action = 'store_true' , help = "remove strandedness of features (set to '.') when using ``transcripts2genes`` or ``filter``[%default]." ) #NEWLINE# parser . add_option ( '--permit-duplicates' , dest = 'strict' , action = 'store_false' , help = 'permit duplicate genes. [%default]' ) #NEWLINE# parser . add_option ( '--duplicate-feature' , dest = 'duplicate_feature' , type = 'choice' , choices = ( 'gene' , 'transcript' , 'both' , 'ucsc' , 'coordinates' ) , help = 'remove duplicates by gene/transcript. If ``ucsc`` is chosen, transcripts ending on _dup# are removed. This is necessary to remove duplicate entries that are next to each other in the sort order [%default]' ) #NEWLINE# parser . add_option ( '--use-gene-id' , dest = 'use_geneid' , action = 'store_true' , help = 'when merging transcripts, exons or introns, use the parent gene_id as the transcript id.' ) #NEWLINE# parser . add_option ( '-m' , '--method' , dest = 'method' , type = 'choice' , action = 'append' , choices = ( 'add-protein-id' , 'exons2introns' , 'filter' , 'find-retained-introns' , 'genes-to-unique-chunks' , 'intersect-transcripts' , 'join-exons' , 'merge-exons' , 'merge-transcripts' , 'merge-genes' , 'merge-introns' , 'remove-overlapping' , 'remove-duplicates' , 'rename-genes' , 'rename-transcripts' , 'rename-duplicates' , 'renumber-genes' , 'renumber-transcripts' , 'set-transcript-to-gene' , 'set-gene-to-transcript' , 'set-protein-to-transcript' , 'set-score-to-distance' , 'set-gene_biotype-to-source' , 'set-source-to-transcript_biotype' , 'sort' , 'transcript2genes' , 'unset-genes' ) , help = 'Method to apply [%default].Please only select one.' ) #NEWLINE# parser . set_defaults ( sort_order = 'gene' , filter_method = 'gene' , pattern = '%i' , merge_exons_distance = 0 , filename_filter = None , intron_border = None , intron_min_length = None , sample_size = 0 , min_exons_length = 0 , ignore_strand = False , mark_utr = False , with_utr = True , invert_filter = False , duplicate_feature = None , strict = True , method = None , use_geneid = False ) #NEWLINE# ( options , args ) = E . Start ( parser , argv = argv ) #NEWLINE# ( ninput , noutput , nfeatures , ndiscarded ) = ( 0 , 0 , 0 , 0 ) #NEWLINE# if ( options . method is None ) : #NEWLINE# #INDENT# raise ValueError ( 'please specify a --method' ) #NEWLINE# #UNINDENT# if ( len ( options . method ) > 1 ) : #NEWLINE# #INDENT# raise ValueError ( 'multiple --method arguements specified' ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# options . method = options . method [ 0 ] #NEWLINE# #UNINDENT# if ( options . method == 'set-transcript-to-gene' ) : #NEWLINE# #INDENT# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# gff . setAttribute ( 'transcript_id' , gff . gene_id ) #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# noutput += 1 #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'set-gene_biotype-to-source' ) : #NEWLINE# #INDENT# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# if ( 'gene_biotype' not in gff . attributes ) : #NEWLINE# #INDENT# gff . setAttribute ( 'gene_biotype' , gff . source ) #NEWLINE# #UNINDENT# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# noutput += 1 #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'set-source-to-transcript_biotype' ) : #NEWLINE# #INDENT# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# try : #NEWLINE# #INDENT# gff . source = gff . transcript_biotype #NEWLINE# #UNINDENT# except AttributeError : #NEWLINE# #INDENT# pass #NEWLINE# #UNINDENT# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# noutput += 1 #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'remove-duplicates' ) : #NEWLINE# #INDENT# counts = collections . defaultdict ( int ) #NEWLINE# if ( options . duplicate_feature == 'ucsc' ) : #NEWLINE# #INDENT# store = [ ] #NEWLINE# remove = set ( ) #NEWLINE# f = ( lambda x : x [ 0 ] . transcript_id ) #NEWLINE# gffs = GTF . transcript_iterator ( GTF . iterator ( options . stdin ) , strict = False ) #NEWLINE# outf = ( lambda x : '\n' . join ( [ str ( y ) for y in x ] ) ) #NEWLINE# for entry in gffs : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# store . append ( entry ) #NEWLINE# id = f ( entry ) #NEWLINE# if ( '_dup' in id ) : #NEWLINE# #INDENT# remove . add ( re . sub ( '_dup\\d+' , '' , id ) ) #NEWLINE# remove . add ( id ) #NEWLINE# #UNINDENT# #UNINDENT# for entry in store : #NEWLINE# #INDENT# id = f ( entry ) #NEWLINE# if ( id not in remove ) : #NEWLINE# #INDENT# options . stdout . write ( ( outf ( entry ) + '\n' ) ) #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# E . info ( ( 'discarded duplicates for %s' % id ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . duplicate_feature == 'gene' ) : #NEWLINE# #INDENT# gffs = GTF . gene_iterator ( GTF . iterator ( options . stdin ) , strict = False ) #NEWLINE# f = ( lambda x : x [ 0 ] [ 0 ] . gene_id ) #NEWLINE# outf = ( lambda x : '\n' . join ( [ '\n' . join ( [ str ( y ) for y in xx ] ) for xx in x ] ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . duplicate_feature == 'transcript' ) : #NEWLINE# #INDENT# gffs = GTF . transcript_iterator ( GTF . iterator ( options . stdin ) , strict = False ) #NEWLINE# f = ( lambda x : x [ 0 ] . transcript_id ) #NEWLINE# outf = ( lambda x : '\n' . join ( [ str ( y ) for y in x ] ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . duplicate_feature == 'coordinates' ) : #NEWLINE# #INDENT# gffs = GTF . chunk_iterator ( GTF . iterator ( options . stdin ) ) #NEWLINE# f = ( lambda x : ( ( ( ( x [ 0 ] . contig + '_' ) + str ( x [ 0 ] . start ) ) + '-' ) + str ( x [ 0 ] . end ) ) ) #NEWLINE# outf = ( lambda x : '\n' . join ( [ str ( y ) for y in x ] ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# store = [ ] #NEWLINE# for entry in gffs : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# store . append ( entry ) #NEWLINE# id = f ( entry ) #NEWLINE# counts [ id ] += 1 #NEWLINE# #UNINDENT# last_id = '' #NEWLINE# if ( options . duplicate_feature == 'coordinates' ) : #NEWLINE# #INDENT# for entry in store : #NEWLINE# #INDENT# id = f ( entry ) #NEWLINE# if ( id == last_id ) : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# E . info ( ( 'discarded duplicates for %s: %i' % ( id , counts [ id ] ) ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# options . stdout . write ( ( outf ( entry ) + '\n' ) ) #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# last_id = id #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# for entry in store : #NEWLINE# #INDENT# id = f ( entry ) #NEWLINE# if ( counts [ id ] == 1 ) : #NEWLINE# #INDENT# options . stdout . write ( ( outf ( entry ) + '\n' ) ) #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# E . info ( ( 'discarded duplicates for %s: %i' % ( id , counts [ id ] ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( 'sort' == options . method ) : #NEWLINE# #INDENT# for gff in GTF . iterator_sorted ( GTF . iterator ( options . stdin ) , sort_order = options . sort_order ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# noutput += 1 #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( 'set-gene-to-transcript' == options . method ) : #NEWLINE# #INDENT# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# gff . setAttribute ( 'gene_id' , gff . transcript_id ) #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# noutput += 1 #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( 'set-protein-to-transcript' == options . method ) : #NEWLINE# #INDENT# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# gff . setAttribute ( 'protein_id' , gff . transcript_id ) #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# noutput += 1 #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( 'add-protein-id' == options . method ) : #NEWLINE# #INDENT# transcript2protein = IOTools . readMap ( IOTools . openFile ( options . filename_filter , 'r' ) ) #NEWLINE# missing = set ( ) #NEWLINE# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# if ( gff . transcript_id not in transcript2protein ) : #NEWLINE# #INDENT# if ( gff . transcript_id not in missing ) : #NEWLINE# #INDENT# E . debug ( ( "removing transcript '%s' due to missing protein id" % gff . transcript_id ) ) #NEWLINE# missing . add ( gff . transcript_id ) #NEWLINE# #UNINDENT# ndiscarded += 1 #NEWLINE# continue #NEWLINE# #UNINDENT# gff . setAttribute ( 'protein_id' , transcript2protein [ gff . transcript_id ] ) #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# noutput += 1 #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# E . info ( ( 'transcripts removed due to missing protein ids: %i' % len ( missing ) ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( 'join-exons' == options . method ) : #NEWLINE# #INDENT# for exons in GTF . transcript_iterator ( GTF . iterator ( options . stdin ) ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# strand = Genomics . convertStrand ( exons [ 0 ] . strand ) #NEWLINE# contig = exons [ 0 ] . contig #NEWLINE# transid = exons [ 0 ] . transcript_id #NEWLINE# geneid = exons [ 0 ] . gene_id #NEWLINE# biotype = exons [ 0 ] . source #NEWLINE# ( all_start , all_end ) = ( min ( [ x . start for x in exons ] ) , max ( [ x . end for x in exons ] ) ) #NEWLINE# y = GTF . Entry ( ) #NEWLINE# y . contig = contig #NEWLINE# y . source = biotype #NEWLINE# y . feature = 'transcript' #NEWLINE# y . start = all_start #NEWLINE# y . end = all_end #NEWLINE# y . strand = strand #NEWLINE# y . transcript_id = transid #NEWLINE# y . gene_id = geneid #NEWLINE# options . stdout . write ( ( '%s\n' % str ( y ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( 'merge-genes' == options . method ) : #NEWLINE# #INDENT# gffs = GTF . iterator_sorted_chunks ( GTF . flat_gene_iterator ( GTF . iterator ( options . stdin ) ) , sort_by = 'contig-strand-start' ) #NEWLINE# def iterate_chunks( gff_chunks ) : #NEWLINE# #INDENT# last = gff_chunks . next ( ) #NEWLINE# to_join = [ last ] #NEWLINE# for gffs in gff_chunks : #NEWLINE# #INDENT# d = ( gffs [ 0 ] . start - last [ ( - 1 ) ] . end ) #NEWLINE# if ( ( gffs [ 0 ] . contig == last [ 0 ] . contig ) and ( gffs [ 0 ] . strand == last [ 0 ] . strand ) ) : #NEWLINE# #INDENT# assert ( gffs [ 0 ] . start >= last [ 0 ] . start ) , ( 'input file should be sorted by contig, strand and position: d=%i:\nlast=\n%s\nthis=\n%s\n' % ( d , '\n' . join ( [ str ( x ) for x in last ] ) , '\n' . join ( [ str ( x ) for x in gffs ] ) ) ) #NEWLINE# #UNINDENT# if ( ( gffs [ 0 ] . contig != last [ 0 ] . contig ) or ( gffs [ 0 ] . strand != last [ 0 ] . strand ) or ( d > 0 ) ) : #NEWLINE# #INDENT# ( yield to_join ) #NEWLINE# to_join = [ ] #NEWLINE# #UNINDENT# last = gffs #NEWLINE# to_join . append ( gffs ) #NEWLINE# #UNINDENT# ( yield to_join ) #NEWLINE# raise StopIteration #NEWLINE# #UNINDENT# for chunks in iterate_chunks ( gffs ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# if ( len ( chunks ) > 1 ) : #NEWLINE# #INDENT# gene_id = ( 'merged_%s' % chunks [ 0 ] [ 0 ] . gene_id ) #NEWLINE# transcript_id = ( 'merged_%s' % chunks [ 0 ] [ 0 ] . transcript_id ) #NEWLINE# info = ',' . join ( [ x [ 0 ] . gene_id for x in chunks ] ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# gene_id = chunks [ 0 ] [ 0 ] . gene_id #NEWLINE# transcript_id = chunks [ 0 ] [ 0 ] . transcript_id #NEWLINE# info = None #NEWLINE# #UNINDENT# intervals = [ ] #NEWLINE# for c in chunks : #NEWLINE# #INDENT# intervals += [ ( x . start , x . end ) for x in c ] #NEWLINE# #UNINDENT# intervals = Intervals . combine ( intervals ) #NEWLINE# strand = chunks [ 0 ] [ 0 ] . strand #NEWLINE# for ( start , end ) in intervals : #NEWLINE# #INDENT# y = GTF . Entry ( ) #NEWLINE# y . fromGTF ( chunks [ 0 ] [ 0 ] , gene_id , transcript_id ) #NEWLINE# y . start = start #NEWLINE# y . end = end #NEWLINE# y . strand = strand #NEWLINE# if info : #NEWLINE# #INDENT# y . addAttribute ( 'merged' , info ) #NEWLINE# #UNINDENT# options . stdout . write ( ( '%s\n' % str ( y ) ) ) #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'renumber-genes' ) : #NEWLINE# #INDENT# map_old2new = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# for gtf in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# if ( gtf . gene_id not in map_old2new ) : #NEWLINE# #INDENT# map_old2new [ gtf . gene_id ] = ( options . pattern % ( len ( map_old2new ) + 1 ) ) #NEWLINE# #UNINDENT# gtf . setAttribute ( 'gene_id' , map_old2new [ gtf . gene_id ] ) #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gtf ) ) ) #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'unset-genes' ) : #NEWLINE# #INDENT# map_old2new = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# for gtf in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# key = gtf . transcript_id #NEWLINE# if ( key not in map_old2new ) : #NEWLINE# #INDENT# map_old2new [ key ] = ( options . pattern % ( len ( map_old2new ) + 1 ) ) #NEWLINE# #UNINDENT# gtf . setAttribute ( 'gene_id' , map_old2new [ key ] ) #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gtf ) ) ) #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'renumber-transcripts' ) : #NEWLINE# #INDENT# map_old2new = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# for gtf in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# key = ( gtf . gene_id , gtf . transcript_id ) #NEWLINE# if ( key not in map_old2new ) : #NEWLINE# #INDENT# map_old2new [ key ] = ( options . pattern % ( len ( map_old2new ) + 1 ) ) #NEWLINE# #UNINDENT# gtf . setAttribute ( 'transcript_id' , map_old2new [ key ] ) #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gtf ) ) ) #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'transcripts2genes' ) : #NEWLINE# #INDENT# transcripts = set ( ) #NEWLINE# genes = set ( ) #NEWLINE# ignore_strand = options . ignore_strand #NEWLINE# for gtfs in GTF . iterator_transcripts2genes ( GTF . iterator ( options . stdin ) ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# for gtf in gtfs : #NEWLINE# #INDENT# if ignore_strand : #NEWLINE# #INDENT# gtf . strand = '.' #NEWLINE# #UNINDENT# options . stdout . write ( ( '%s\n' % str ( gtf ) ) ) #NEWLINE# transcripts . add ( gtf . transcript_id ) #NEWLINE# genes . add ( gtf . gene_id ) #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# noutput += 1 #NEWLINE# #UNINDENT# E . info ( ( 'transcripts2genes: transcripts=%i, genes=%i' % ( len ( transcripts ) , len ( genes ) ) ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method in ( 'rename-genes' , 'rename-transcripts' ) ) : #NEWLINE# #INDENT# map_old2new = IOTools . readMap ( IOTools . openFile ( options . filename_filter , 'r' ) ) #NEWLINE# if ( options . method == 'rename-transcripts' ) : #NEWLINE# #INDENT# is_gene_id = False #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'rename-genes' ) : #NEWLINE# #INDENT# is_gene_id = True #NEWLINE# #UNINDENT# #UNINDENT# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# if is_gene_id : #NEWLINE# #INDENT# if ( gff . gene_id in map_old2new ) : #NEWLINE# #INDENT# gff . setAttribute ( 'gene_id' , map_old2new [ gff . gene_id ] ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# E . debug ( ( 'removing missing gene_id %s' % gff . gene_id ) ) #NEWLINE# ndiscarded += 1 #NEWLINE# continue #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( gff . transcript_id in map_old2new ) : #NEWLINE# #INDENT# gff . setAttribute ( 'transcript_id' , map_old2new [ gff . transcript_id ] ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# E . debug ( ( 'removing missing transcript_id %s' % gff . transcript_id ) ) #NEWLINE# ndiscarded += 1 #NEWLINE# continue #NEWLINE# #UNINDENT# #UNINDENT# noutput += 1 #NEWLINE# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'filter' ) : #NEWLINE# #INDENT# keep_genes = set ( ) #NEWLINE# if ( options . filter_method == 'longest-gene' ) : #NEWLINE# #INDENT# iterator = GTF . flat_gene_iterator ( GTF . iterator ( options . stdin ) ) #NEWLINE# coords = [ ] #NEWLINE# gffs = [ ] #NEWLINE# for gff in iterator : #NEWLINE# #INDENT# gff . sort ( key = ( lambda x : x . start ) ) #NEWLINE# coords . append ( ( gff [ 0 ] . contig , min ( [ x . start for x in gff ] ) , max ( [ x . end for x in gff ] ) , gff [ 0 ] . gene_id ) ) #NEWLINE# gffs . append ( gff ) #NEWLINE# #UNINDENT# coords . sort ( ) #NEWLINE# last_contig = None #NEWLINE# max_end = 0 #NEWLINE# longest_gene_id = None #NEWLINE# longest_length = None #NEWLINE# for ( contig , start , end , gene_id ) in coords : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# if ( ( contig != last_contig ) or ( start >= max_end ) ) : #NEWLINE# #INDENT# if longest_gene_id : #NEWLINE# #INDENT# keep_genes . add ( longest_gene_id ) #NEWLINE# #UNINDENT# longest_gene_id = gene_id #NEWLINE# longest_length = ( end - start ) #NEWLINE# max_end = end #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( end - start ) > longest_length ) : #NEWLINE# #INDENT# ( longest_length , longest_gene_id ) = ( ( end - start ) , gene_id ) #NEWLINE# #UNINDENT# #UNINDENT# last_contig = contig #NEWLINE# max_end = max ( max_end , end ) #NEWLINE# #UNINDENT# keep_genes . add ( longest_gene_id ) #NEWLINE# invert = options . invert_filter #NEWLINE# for gff in gffs : #NEWLINE# #INDENT# keep = ( gff [ 0 ] . gene_id in keep_genes ) #NEWLINE# if ( ( keep and ( not invert ) ) or ( ( not keep ) and invert ) ) : #NEWLINE# #INDENT# noutput += 1 #NEWLINE# for g in gff : #NEWLINE# #INDENT# nfeatures += 1 #NEWLINE# options . stdout . write ( ( '%s\n' % g ) ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . filter_method in ( 'longest-transcript' , 'representative-transcript' ) ) : #NEWLINE# #INDENT# iterator = GTF . gene_iterator ( GTF . iterator ( options . stdin ) ) #NEWLINE# def selectLongestTranscript( gene ) : #NEWLINE# #INDENT# r = [ ] #NEWLINE# for transcript in gene : #NEWLINE# #INDENT# transcript . sort ( key = ( lambda x : x . start ) ) #NEWLINE# length = ( transcript [ ( - 1 ) ] . end - transcript [ 0 ] . start ) #NEWLINE# r . append ( ( length , transcript ) ) #NEWLINE# #UNINDENT# r . sort ( ) #NEWLINE# return r [ ( - 1 ) ] [ 1 ] #NEWLINE# #UNINDENT# def selectRepresentativeTranscript( gene ) : #NEWLINE# #INDENT# 'select a representative transcript.\n\n                The representative transcript represent the largest number\n                of exons over all transcripts.\n                ' #NEWLINE# all_exons = [ ] #NEWLINE# for transcript in gene : #NEWLINE# #INDENT# all_exons . extend ( [ ( x . start , x . end ) for x in transcript if ( x . feature == 'exon' ) ] ) #NEWLINE# #UNINDENT# exon_counts = { #NEWLINE# #INDENT# #NEWLINE# #UNINDENT# } #NEWLINE# for ( key , exons ) in itertools . groupby ( all_exons ) : #NEWLINE# #INDENT# exon_counts [ key ] = len ( list ( exons ) ) #NEWLINE# #UNINDENT# transcript_counts = [ ] #NEWLINE# for transcript in gene : #NEWLINE# #INDENT# count = sum ( [ exon_counts [ ( x . start , x . end ) ] for x in transcript if ( x . feature == 'exon' ) ] ) #NEWLINE# transcript_counts . append ( ( count , transcript [ 0 ] . transcript_id , transcript ) ) #NEWLINE# #UNINDENT# transcript_counts . sort ( ) #NEWLINE# return transcript_counts [ ( - 1 ) ] [ ( - 1 ) ] #NEWLINE# #UNINDENT# if ( options . filter_method == 'longest-transcript' ) : #NEWLINE# #INDENT# _select = selectLongestTranscript #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . filter_method == 'representative-transcript' ) : #NEWLINE# #INDENT# _select = selectRepresentativeTranscript #NEWLINE# #UNINDENT# #UNINDENT# for gene in iterator : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# transcript = _select ( sorted ( gene ) ) #NEWLINE# noutput += 1 #NEWLINE# for g in transcript : #NEWLINE# #INDENT# nfeatures += 1 #NEWLINE# options . stdout . write ( ( '%s\n' % g ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . filter_method in ( 'gene' , 'transcript' ) ) : #NEWLINE# #INDENT# if options . filename_filter : #NEWLINE# #INDENT# ids = IOTools . readList ( IOTools . openFile ( options . filename_filter , 'r' ) ) #NEWLINE# E . info ( ( 'read %i ids' % len ( ids ) ) ) #NEWLINE# ids = set ( ids ) #NEWLINE# by_gene = ( options . filter_method == 'gene' ) #NEWLINE# by_transcript = ( options . filter_method == 'transcript' ) #NEWLINE# invert = options . invert_filter #NEWLINE# ignore_strand = options . ignore_strand #NEWLINE# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# keep = False #NEWLINE# if by_gene : #NEWLINE# #INDENT# keep = ( gff . gene_id in ids ) #NEWLINE# #UNINDENT# if by_transcript : #NEWLINE# #INDENT# keep = ( gff . transcript_id in ids ) #NEWLINE# #UNINDENT# if ( ( invert and keep ) or ( ( not invert ) and ( not keep ) ) ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# if ignore_strand : #NEWLINE# #INDENT# gff . strand = '.' #NEWLINE# #UNINDENT# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# nfeatures += 1 #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if options . sample_size : #NEWLINE# #INDENT# if ( options . filter_method == 'gene' ) : #NEWLINE# #INDENT# iterator = GTF . flat_gene_iterator ( GTF . iterator ( options . stdin ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . filter_method == 'transcript' ) : #NEWLINE# #INDENT# iterator = GTF . transcript_iterator ( GTF . iterator ( options . stdin ) ) #NEWLINE# #UNINDENT# #UNINDENT# if options . min_exons_length : #NEWLINE# #INDENT# iterator = GTF . iterator_min_feature_length ( iterator , min_length = options . min_exons_length , feature = 'exon' ) #NEWLINE# #UNINDENT# data = [ x for x in iterator ] #NEWLINE# ninput = len ( data ) #NEWLINE# if ( len ( data ) > options . sample_size ) : #NEWLINE# #INDENT# data = random . sample ( data , options . sample_size ) #NEWLINE# #UNINDENT# for d in data : #NEWLINE# #INDENT# noutput += 1 #NEWLINE# for dd in d : #NEWLINE# #INDENT# nfeatures += 1 #NEWLINE# options . stdout . write ( ( str ( dd ) + '\n' ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# assert False , 'please supply either a filename ' #NEWLINE# 'with ids to filter with (--map-tsv-file) or a sample-size.' #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . filter_method in ( 'proteincoding' , 'lincrna' , 'processed-pseudogene' ) ) : #NEWLINE# #INDENT# tag = { #NEWLINE# #INDENT# 'proteincoding' : 'protein_coding' , #NEWLINE# 'processed-pseudogene' : 'processed_pseudogene' , #NEWLINE# 'lincrna' : 'lincRNA' , #NEWLINE# #UNINDENT# } [ options . filter_method ] #NEWLINE# rx = re . compile ( ( '"%s"' % tag ) ) #NEWLINE# if ( not options . invert_filter ) : #NEWLINE# #INDENT# f = ( lambda x : ( ( x . source == tag ) or rx . search ( x . attributes ) ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# f = ( lambda x : ( ( x . source != tag ) and ( not rx . search ( x . attributes ) ) ) ) #NEWLINE# #UNINDENT# for gff in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# if f ( gff ) : #NEWLINE# #INDENT# options . stdout . write ( ( str ( gff ) + '\n' ) ) #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'exons2introns' ) : #NEWLINE# #INDENT# for gffs in GTF . flat_gene_iterator ( GTF . iterator ( options . stdin ) ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# cds_ranges = GTF . asRanges ( gffs , 'CDS' ) #NEWLINE# exon_ranges = GTF . asRanges ( gffs , 'exon' ) #NEWLINE# input_ranges = Intervals . combine ( ( cds_ranges + exon_ranges ) ) #NEWLINE# if ( len ( input_ranges ) > 1 ) : #NEWLINE# #INDENT# last = input_ranges [ 0 ] [ 1 ] #NEWLINE# output_ranges = [ ] #NEWLINE# for ( start , end ) in input_ranges [ 1 : ] : #NEWLINE# #INDENT# output_ranges . append ( ( last , start ) ) #NEWLINE# last = end #NEWLINE# #UNINDENT# if options . intron_border : #NEWLINE# #INDENT# b = options . intron_border #NEWLINE# output_ranges = [ ( ( x [ 0 ] + b ) , ( x [ 1 ] - b ) ) for x in output_ranges ] #NEWLINE# #UNINDENT# if options . intron_min_length : #NEWLINE# #INDENT# l = options . intron_min_length #NEWLINE# output_ranges = [ x for x in output_ranges if ( ( x [ 1 ] - x [ 0 ] ) > l ) ] #NEWLINE# #UNINDENT# for ( start , end ) in output_ranges : #NEWLINE# #INDENT# entry = GTF . Entry ( ) #NEWLINE# entry . copy ( gffs [ 0 ] ) #NEWLINE# entry . clearAttributes ( ) #NEWLINE# entry . transcript_id = 'merged' #NEWLINE# entry . feature = 'intron' #NEWLINE# entry . start = start #NEWLINE# entry . end = end #NEWLINE# options . stdout . write ( ( '%s\n' % str ( entry ) ) ) #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# noutput += 1 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'set-score-to-distance' ) : #NEWLINE# #INDENT# for gffs in GTF . transcript_iterator ( GTF . iterator ( options . stdin ) ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# strand = Genomics . convertStrand ( gffs [ 0 ] . strand ) #NEWLINE# ( all_start , all_end ) = ( min ( [ x . start for x in gffs ] ) , max ( [ x . end for x in gffs ] ) ) #NEWLINE# if ( strand != '.' ) : #NEWLINE# #INDENT# t = 0 #NEWLINE# if ( strand == '-' ) : #NEWLINE# #INDENT# gffs . reverse ( ) #NEWLINE# #UNINDENT# for gff in gffs : #NEWLINE# #INDENT# gff . score = t #NEWLINE# t += ( gff . end - gff . start ) #NEWLINE# #UNINDENT# if ( strand == '-' ) : #NEWLINE# #INDENT# gffs . reverse ( ) #NEWLINE# #UNINDENT# #UNINDENT# for gff in gffs : #NEWLINE# #INDENT# options . stdout . write ( ( '%s\n' % str ( gff ) ) ) #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'remove-overlapping' ) : #NEWLINE# #INDENT# index = GTF . readAndIndex ( GTF . iterator ( IOTools . openFile ( options . filename_gff , 'r' ) ) ) #NEWLINE# for gffs in GTF . transcript_iterator ( GTF . iterator ( options . stdin ) ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# found = False #NEWLINE# for e in gffs : #NEWLINE# #INDENT# if index . contains ( e . contig , e . start , e . end ) : #NEWLINE# #INDENT# found = True #NEWLINE# break #NEWLINE# #UNINDENT# #UNINDENT# if found : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# noutput += 1 #NEWLINE# for e in gffs : #NEWLINE# #INDENT# nfeatures += 1 #NEWLINE# options . stdout . write ( ( '%s\n' % str ( e ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'intersect-transcripts' ) : #NEWLINE# #INDENT# for gffs in GTF . gene_iterator ( GTF . iterator ( options . stdin ) , strict = options . strict ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# r = [ ] #NEWLINE# for g in gffs : #NEWLINE# #INDENT# if options . with_utr : #NEWLINE# #INDENT# ranges = GTF . asRanges ( g , 'exon' ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# ranges = GTF . asRanges ( g , 'CDS' ) #NEWLINE# #UNINDENT# r . append ( ranges ) #NEWLINE# #UNINDENT# result = r [ 0 ] #NEWLINE# for x in r [ 1 : ] : #NEWLINE# #INDENT# result = Intervals . intersect ( result , x ) #NEWLINE# #UNINDENT# entry = GTF . Entry ( ) #NEWLINE# entry . copy ( gffs [ 0 ] [ 0 ] ) #NEWLINE# entry . clearAttributes ( ) #NEWLINE# entry . transcript_id = 'merged' #NEWLINE# entry . feature = 'exon' #NEWLINE# for ( start , end ) in result : #NEWLINE# #INDENT# entry . start = start #NEWLINE# entry . end = end #NEWLINE# options . stdout . write ( ( '%s\n' % str ( entry ) ) ) #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( 'rename-duplicates' == options . method ) : #NEWLINE# #INDENT# assert ( options . duplicate_feature in [ 'gene' , 'transcript' , 'both' ] ) , "for renaming duplicates, --duplicate-feature must be set to one of 'gene', transcript' or 'both'" #NEWLINE# gene_ids = list ( ) #NEWLINE# transcript_ids = list ( ) #NEWLINE# gtfs = list ( ) #NEWLINE# for gtf in GTF . iterator ( options . stdin ) : #NEWLINE# #INDENT# gtfs . append ( gtf ) #NEWLINE# if ( gtf . feature == 'CDS' ) : #NEWLINE# #INDENT# gene_ids . append ( gtf . gene_id ) #NEWLINE# transcript_ids . append ( gtf . transcript_id ) #NEWLINE# #UNINDENT# #UNINDENT# dup_gene = [ item for item in set ( gene_ids ) if ( gene_ids . count ( item ) > 1 ) ] #NEWLINE# dup_transcript = [ item for item in set ( transcript_ids ) if ( transcript_ids . count ( item ) > 1 ) ] #NEWLINE# E . info ( ( 'Number of duplicated gene_ids: %i' % len ( dup_gene ) ) ) #NEWLINE# E . info ( ( 'Number of duplicated transcript_ids: %i' % len ( dup_transcript ) ) ) #NEWLINE# gene_dict = dict ( zip ( dup_gene , ( [ 0 ] * len ( dup_gene ) ) ) ) #NEWLINE# transcript_dict = dict ( zip ( dup_transcript , ( [ 0 ] * len ( dup_transcript ) ) ) ) #NEWLINE# for gtf in gtfs : #NEWLINE# #INDENT# if ( gtf . feature == 'CDS' ) : #NEWLINE# #INDENT# if ( options . duplicate_feature in [ 'both' , 'gene' ] ) : #NEWLINE# #INDENT# if ( gtf . gene_id in dup_gene ) : #NEWLINE# #INDENT# gene_dict [ gtf . gene_id ] = ( gene_dict [ gtf . gene_id ] + 1 ) #NEWLINE# gtf . attributes = gtf . attributes . strip ( ) #NEWLINE# gtf . setAttribute ( 'gene_id' , ( ( gtf . gene_id + '.' ) + str ( gene_dict [ gtf . gene_id ] ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# if ( options . duplicate_feature in [ 'both' , 'transcript' ] ) : #NEWLINE# #INDENT# if ( gtf . transcript_id in dup_transcript ) : #NEWLINE# #INDENT# transcript_dict [ gtf . transcript_id ] = ( transcript_dict [ gtf . transcript_id ] + 1 ) #NEWLINE# gtf . attributes = gtf . attributes . strip ( ) #NEWLINE# gtf . setAttribute ( 'transcript_id' , ( ( gtf . transcript_id + '.' ) + str ( transcript_dict [ gtf . transcript_id ] ) ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# options . stdout . write ( ( '%s\n' % gtf ) ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method in ( 'merge-exons' , 'merge-introns' , 'merge-transcripts' ) ) : #NEWLINE# #INDENT# for gffs in GTF . flat_gene_iterator ( GTF . iterator ( options . stdin ) , strict = options . strict ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# cds_ranges = GTF . asRanges ( gffs , 'CDS' ) #NEWLINE# exon_ranges = GTF . asRanges ( gffs , 'exon' ) #NEWLINE# strands = set ( [ x . strand for x in gffs ] ) #NEWLINE# contigs = set ( [ x . contig for x in gffs ] ) #NEWLINE# if ( len ( strands ) > 1 ) : #NEWLINE# #INDENT# raise ValueError ( ( "can not merge gene '%s' on multiple strands: %s" % ( gffs [ 0 ] . gene_id , str ( strands ) ) ) ) #NEWLINE# #UNINDENT# if ( len ( contigs ) > 1 ) : #NEWLINE# #INDENT# raise ValueError ( ( "can not merge gene '%s' on multiple contigs: %s" % ( gffs [ 0 ] . gene_id , str ( contigs ) ) ) ) #NEWLINE# #UNINDENT# strand = Genomics . convertStrand ( gffs [ 0 ] . strand ) #NEWLINE# utr_ranges = [ ] #NEWLINE# if ( cds_ranges and options . mark_utr ) : #NEWLINE# #INDENT# ( cds_start , cds_end ) = ( cds_ranges [ 0 ] [ 0 ] , cds_ranges [ ( - 1 ) ] [ 1 ] ) #NEWLINE# midpoint = ( ( ( cds_end - cds_start ) / 2 ) + cds_start ) #NEWLINE# utr_ranges = [ ] #NEWLINE# for ( start , end ) in Intervals . truncate ( exon_ranges , cds_ranges ) : #NEWLINE# #INDENT# if ( ( end - start ) > 3 ) : #NEWLINE# #INDENT# if ( strand == '.' ) : #NEWLINE# #INDENT# feature = 'UTR' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( strand == '+' ) : #NEWLINE# #INDENT# if ( start < midpoint ) : #NEWLINE# #INDENT# feature = 'UTR5' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# feature = 'UTR3' #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( strand == '-' ) : #NEWLINE# #INDENT# if ( start < midpoint ) : #NEWLINE# #INDENT# feature = 'UTR3' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# feature = 'UTR5' #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# utr_ranges . append ( ( feature , start , end ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# try : #NEWLINE# #INDENT# biotypes = [ x [ 'gene_biotype' ] for x in gffs ] #NEWLINE# biotype = ':' . join ( set ( biotypes ) ) #NEWLINE# #UNINDENT# except ( KeyError , AttributeError ) : #NEWLINE# #INDENT# biotype = None #NEWLINE# #UNINDENT# def output_ranges( ranges , gffs , biotype = None , use_geneid = False ) : #NEWLINE# #INDENT# result = [ ] #NEWLINE# for ( feature , start , end ) in ranges : #NEWLINE# #INDENT# entry = GTF . Entry ( ) #NEWLINE# entry . copy ( gffs [ 0 ] ) #NEWLINE# entry . clearAttributes ( ) #NEWLINE# entry . feature = feature #NEWLINE# if use_geneid : #NEWLINE# #INDENT# entry . transcript_id = entry . gene_id #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# entry . transcript_id = 'merged' #NEWLINE# #UNINDENT# if biotype : #NEWLINE# #INDENT# entry . addAttribute ( 'gene_biotype' , biotype ) #NEWLINE# #UNINDENT# entry . start = start #NEWLINE# entry . end = end #NEWLINE# result . append ( entry ) #NEWLINE# #UNINDENT# return result #NEWLINE# #UNINDENT# result = [ ] #NEWLINE# if ( options . method == 'merge-exons' ) : #NEWLINE# #INDENT# if options . with_utr : #NEWLINE# #INDENT# if options . mark_utr : #NEWLINE# #INDENT# result . extend ( output_ranges ( utr_ranges , gffs , biotype , options . use_geneid ) ) #NEWLINE# r = [ ( 'CDS' , x , y ) for ( x , y ) in Intervals . combineAtDistance ( cds_ranges , options . merge_exons_distance ) ] #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# r = [ ( 'exon' , x , y ) for ( x , y ) in Intervals . combineAtDistance ( exon_ranges , options . merge_exons_distance ) ] #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# r = [ ( 'CDS' , x , y ) for ( x , y ) in Intervals . combineAtDistance ( cds_ranges , options . merge_exons_distance ) ] #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'merge-transcripts' ) : #NEWLINE# #INDENT# if options . with_utr : #NEWLINE# #INDENT# r = [ ( 'exon' , exon_ranges [ 0 ] [ 0 ] , exon_ranges [ ( - 1 ) ] [ 1 ] ) ] #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if cds_ranges : #NEWLINE# #INDENT# r = [ ( 'exon' , cds_ranges [ 0 ] [ 0 ] , cds_ranges [ ( - 1 ) ] [ 1 ] ) ] #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# continue #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'merge-introns' ) : #NEWLINE# #INDENT# if ( len ( exon_ranges ) >= 2 ) : #NEWLINE# #INDENT# r = [ ( 'exon' , exon_ranges [ 0 ] [ 1 ] , exon_ranges [ ( - 1 ) ] [ 0 ] ) ] #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# ndiscarded += 1 #NEWLINE# continue #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# result . extend ( output_ranges ( r , gffs , biotype , options . use_geneid ) ) #NEWLINE# result . sort ( key = ( lambda x : x . start ) ) #NEWLINE# for x in result : #NEWLINE# #INDENT# options . stdout . write ( ( '%s\n' % str ( x ) ) ) #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'find-retained-introns' ) : #NEWLINE# #INDENT# for gene in GTF . gene_iterator ( GTF . iterator ( options . stdin ) ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# found_any = False #NEWLINE# for intron in find_retained_introns ( gene ) : #NEWLINE# #INDENT# found_any = True #NEWLINE# options . stdout . write ( ( '%s\n' % str ( intron ) ) ) #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# if found_any : #NEWLINE# #INDENT# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( options . method == 'genes-to-unique-chunks' ) : #NEWLINE# #INDENT# for gene in GTF . flat_gene_iterator ( GTF . iterator ( options . stdin ) ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# for exon in gene_to_blocks ( gene ) : #NEWLINE# #INDENT# options . stdout . write ( ( '%s\n' % str ( exon ) ) ) #NEWLINE# nfeatures += 1 #NEWLINE# #UNINDENT# noutput += 1 #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# raise ValueError ( ( "unknown method '%s'" % options . method ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# E . info ( ( 'ninput=%i, noutput=%i, nfeatures=%i, ndiscarded=%i' % ( ninput , noutput , nfeatures , ndiscarded ) ) ) #NEWLINE# E . Stop ( )