#NEWLINE# def detectPolyA( query_id , matches , options , queries_fasta = None ) : #NEWLINE# #INDENT# 'detect PolyA tails and adjust coverage.\n\n    polyA detection. If there are ambiguous matches, the location\n    of the polyA tail is not straight-forward as suboptimal matches\n    might be declared to be a tail.\n\n    1. collect all matches with polyA tail - the remainder is unchanged\n    2. declare polyA tail by the shortest! unaligned segment and compute\n       coverage for each match appropriately. \n\n    The method checks whether the tails are consistent (always at the same end).\n    If not, an AssertionError is thrown\n    ' #NEWLINE# max_total = matches [ 0 ] . mQueryLength #NEWLINE# ( best_start , best_end , best_pA , best_pT , best_tail ) = ( 0 , 0 , 0 , 0 , '' ) #NEWLINE# tail_matches = [ ] #NEWLINE# new_matches = [ ] #NEWLINE# for match in matches : #NEWLINE# #INDENT# missing_start = match . mQueryFrom #NEWLINE# missing_end = ( match . mQueryLength - match . mQueryTo ) #NEWLINE# if ( missing_start < missing_end ) : #NEWLINE# #INDENT# smaller = missing_start #NEWLINE# larger = missing_end #NEWLINE# ( start , end ) = ( match . mQueryTo , match . mQueryLength ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# smaller = missing_end #NEWLINE# larger = missing_start #NEWLINE# ( start , end ) = ( 0 , match . mQueryFrom ) #NEWLINE# #UNINDENT# if ( not ( ( smaller < options . polyA_max_unaligned ) and ( larger > options . polyA_min_unaligned ) ) ) : #NEWLINE# #INDENT# new_matches . append ( match ) #NEWLINE# continue #NEWLINE# #UNINDENT# tail = queries_fasta . getSequence ( query_id ) [ start : end ] #NEWLINE# counts = { #NEWLINE# #INDENT# 'A' : 0 , #NEWLINE# 'T' : 0 , #NEWLINE# 'N' : 0 , #NEWLINE# #UNINDENT# } #NEWLINE# for c in tail . upper ( ) : #NEWLINE# #INDENT# counts [ c ] = ( counts . get ( c , 0 ) + 1 ) #NEWLINE# #UNINDENT# total = ( end - start ) #NEWLINE# pA = ( ( 100.0 * ( counts [ 'A' ] + counts [ 'N' ] ) ) / total ) #NEWLINE# pT = ( ( 100.0 * ( counts [ 'T' ] + counts [ 'N' ] ) ) / total ) #NEWLINE# if ( options . loglevel >= 5 ) : #NEWLINE# #INDENT# options . stdlog . write ( ( '# polyA detection: %s:%i-%i pA=%5.2f pT=%5.2f tail=%s\n' % ( query_id , start , end , best_start , pT , tail ) ) ) #NEWLINE# #UNINDENT# if ( max ( pA , pT ) < options . polyA_min_percent ) : #NEWLINE# #INDENT# new_matches . append ( match ) #NEWLINE# continue #NEWLINE# #UNINDENT# if ( total < max_total ) : #NEWLINE# #INDENT# max_total = total #NEWLINE# ( best_start , best_end , best_pA , best_pT , best_tail ) = ( start , end , pA , pT , tail ) #NEWLINE# #UNINDENT# if ( not ( ( best_start == start ) or ( best_end == end ) ) ) : #NEWLINE# #INDENT# if ( options . loglevel >= 1 ) : #NEWLINE# #INDENT# options . stdlog . write ( ( '# inconsistent polyA tails for %s: %i-%i, %i-%i - analysis skipped\n' % ( query_id , best_start , best_end , start , end ) ) ) #NEWLINE# #UNINDENT# return matches #NEWLINE# #UNINDENT# tail_matches . append ( match ) #NEWLINE# #UNINDENT# if tail_matches : #NEWLINE# #INDENT# for match in tail_matches : #NEWLINE# #INDENT# match . mQueryCoverage += ( ( 100.0 * float ( len ( best_tail ) ) ) / match . mQueryLength ) #NEWLINE# assert ( match . mQueryCoverage <= 100.0 ) , ( '%s: coverage %f > 100.0: incr=%f\n%s' % ( query_id , match . mQueryCoverage , ( float ( len ( best_tail ) ) / match . mQueryLength ) , str ( match ) ) ) #NEWLINE# new_matches . append ( match ) #NEWLINE# #UNINDENT# options . outfile_polyA . write ( ( '%s\t%i\t%i\t%i\t%5.2f\t%5.2f\t%s\n' % ( query_id , len ( tail_matches ) , best_start , best_end , best_pA , best_pT , best_tail ) ) ) #NEWLINE# #UNINDENT# assert ( len ( new_matches ) == len ( matches ) ) #NEWLINE# return new_matches