#NEWLINE# def trimUTR( args ) : #NEWLINE# #INDENT# "\n    %prog trimUTR gffile\n\n    Remove UTRs in the annotation set.\n\n    If reference GFF3 is provided, reinstate UTRs from reference\n    transcripts after trimming.\n\n    Note: After running trimUTR, it is advised to also run\n    `python -m jcvi.formats.gff fixboundaries` on the resultant GFF3\n    to adjust the boundaries of all parent 'gene' features\n    " #NEWLINE# import gffutils #NEWLINE# from jcvi.formats.base import SetFile #NEWLINE# p = OptionParser ( trimUTR . __doc__ ) #NEWLINE# p . add_option ( '--trim5' , default = None , type = 'str' , help = "File containing gene list for 5' UTR trimming" ) #NEWLINE# p . add_option ( '--trim3' , default = None , type = 'str' , help = "File containing gene list for 3' UTR trimming" ) #NEWLINE# p . add_option ( '--trimrange' , default = None , type = 'str' , help = ( 'File containing gene list for UTR trim back' + 'based on suggested (start, stop) coordinate range' ) ) #NEWLINE# p . add_option ( '--refgff' , default = None , type = 'str' , help = 'Reference GFF3 used as fallback to replace UTRs' ) #NEWLINE# p . set_outfile ( ) #NEWLINE# ( opts , args ) = p . parse_args ( args ) #NEWLINE# if ( len ( args ) != 1 ) : #NEWLINE# #INDENT# sys . exit ( ( not p . print_help ( ) ) ) #NEWLINE# #UNINDENT# ( gffile , ) = args #NEWLINE# gff = make_index ( gffile ) #NEWLINE# trim_both = ( False if ( opts . trim5 or opts . trim3 ) else True ) #NEWLINE# trim5 = ( SetFile ( opts . trim5 ) if opts . trim5 else set ( ) ) #NEWLINE# trim3 = ( SetFile ( opts . trim3 ) if opts . trim3 else set ( ) ) #NEWLINE# trimrange = dict ( ) #NEWLINE# if opts . trimrange : #NEWLINE# #INDENT# trf = must_open ( opts . trimrange ) #NEWLINE# for tr in trf : #NEWLINE# #INDENT# assert ( len ( tr . split ( '\t' ) ) == 3 ) , 'Must specify (start, stop) coordinate range' #NEWLINE# ( id , start , stop ) = tr . split ( '\t' ) #NEWLINE# trimrange [ id ] = ( int ( start ) , int ( stop ) ) #NEWLINE# #UNINDENT# trf . close ( ) #NEWLINE# #UNINDENT# refgff = ( make_index ( opts . refgff ) if opts . refgff else None ) #NEWLINE# fw = must_open ( opts . outfile , 'w' ) #NEWLINE# for feat in gff . iter_by_parent_childs ( featuretype = 'gene' , order_by = ( 'seqid' , 'start' ) , level = 1 ) : #NEWLINE# #INDENT# for c in feat : #NEWLINE# #INDENT# ( cid , ctype , cparent ) = ( c . id , c . featuretype , c . attributes . get ( 'Parent' , [ None ] ) [ 0 ] ) #NEWLINE# ( t5 , t3 ) = ( False , False ) #NEWLINE# if ( ctype == 'gene' ) : #NEWLINE# #INDENT# t5 = ( True if ( cid in trim5 ) else False ) #NEWLINE# t3 = ( True if ( cid in trim3 ) else False ) #NEWLINE# ( start , end ) = get_cds_minmax ( gff , cid ) #NEWLINE# trim ( c , start , end , trim5 = t5 , trim3 = t3 , both = trim_both ) #NEWLINE# fprint ( c , fw ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ctype == 'mRNA' ) : #NEWLINE# #INDENT# ( utr_types , extras ) = ( [ ] , set ( ) ) #NEWLINE# if any ( ( ( id in trim5 ) for id in ( cid , cparent ) ) ) : #NEWLINE# #INDENT# t5 = True #NEWLINE# trim5 . add ( cid ) #NEWLINE# #UNINDENT# if any ( ( ( id in trim3 ) for id in ( cid , cparent ) ) ) : #NEWLINE# #INDENT# t3 = True #NEWLINE# trim3 . add ( cid ) #NEWLINE# #UNINDENT# refc = None #NEWLINE# if refgff : #NEWLINE# #INDENT# try : #NEWLINE# #INDENT# refc = refgff [ cid ] #NEWLINE# refctype = refc . featuretype #NEWLINE# refptype = refgff [ refc . attributes [ 'Parent' ] [ 0 ] ] . featuretype #NEWLINE# if ( ( refctype == 'mRNA' ) and ( refptype == 'gene' ) ) : #NEWLINE# #INDENT# if cmp_children ( cid , gff , refgff , cftype = 'CDS' ) : #NEWLINE# #INDENT# reinstate ( c , refc , trim5 = t5 , trim3 = t3 , both = trim_both ) #NEWLINE# if t5 : #NEWLINE# #INDENT# utr_types . append ( 'five_prime_UTR' ) #NEWLINE# #UNINDENT# if t3 : #NEWLINE# #INDENT# utr_types . append ( 'three_prime_UTR' ) #NEWLINE# #UNINDENT# for utr_type in utr_types : #NEWLINE# #INDENT# for utr in refgff . children ( refc , featuretype = utr_type ) : #NEWLINE# #INDENT# extras . add ( utr ) #NEWLINE# for exon in refgff . region ( region = utr , featuretype = 'exon' ) : #NEWLINE# #INDENT# if ( exon . attributes [ 'Parent' ] [ 0 ] == cid ) : #NEWLINE# #INDENT# extras . add ( exon ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# refc = None #NEWLINE# #UNINDENT# #UNINDENT# except gffutils . exceptions . FeatureNotFoundError : #NEWLINE# #INDENT# pass #NEWLINE# #UNINDENT# #UNINDENT# ( start , end ) = get_cds_minmax ( gff , cid , level = 1 ) #NEWLINE# if ( cid in trimrange ) : #NEWLINE# #INDENT# ( start , end ) = range_minmax ( [ trimrange [ cid ] , ( start , end ) ] ) #NEWLINE# #UNINDENT# if ( not refc ) : #NEWLINE# #INDENT# trim ( c , start , end , trim5 = t5 , trim3 = t3 , both = trim_both ) #NEWLINE# #UNINDENT# fprint ( c , fw ) #NEWLINE# for cc in gff . children ( cid , order_by = 'start' ) : #NEWLINE# #INDENT# _ctype = cc . featuretype #NEWLINE# if ( _ctype not in utr_types ) : #NEWLINE# #INDENT# if ( _ctype != 'CDS' ) : #NEWLINE# #INDENT# if ( _ctype == 'exon' ) : #NEWLINE# #INDENT# eskip = [ range_overlap ( to_range ( cc ) , to_range ( x ) ) for x in extras if ( x . featuretype == 'exon' ) ] #NEWLINE# if any ( ( skip for skip in eskip ) ) : #NEWLINE# #INDENT# continue #NEWLINE# #UNINDENT# #UNINDENT# trim ( cc , start , end , trim5 = t5 , trim3 = t3 , both = trim_both ) #NEWLINE# fprint ( cc , fw ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# fprint ( cc , fw ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# for x in extras : #NEWLINE# #INDENT# fprint ( x , fw ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# fw . close ( )