<s> '''<STR_LIT>''' <EOL> import sys <EOL> import re <EOL> import CGAT . Experiment as E <EOL> import numpy <EOL> def Old ( ) : <EOL> ref_nucleotides = numpy . zeros ( max_end - min_start , numpy . int0 ) <EOL> target_nucleotides = numpy . zeros ( max_end - min_start , numpy . int0 ) <EOL> for x in targets : <EOL> target_nucleotides [ x . start - min_start : x . end - min_start ] = <NUM_LIT:1> <EOL> for x in references : <EOL> ref_nucleotides [ x . start - min_start : x . end - min_start ] = <NUM_LIT:1> <EOL> not_ref = numpy . logical_not ( ref_nucleotides ) <EOL> not_target = numpy . logical_not ( target_nucleotides ) <EOL> v = numpy . logical_and ( ref_nucleotides , target_nucleotides ) <EOL> tp = numpy . dot ( v , v ) <EOL> v = numpy . logical_and ( ref_nucleotides , not_target ) <EOL> fn = numpy . dot ( v , v ) <EOL> v = numpy . logical_and ( not_ref , target_nucleotides ) <EOL> fp = numpy . dot ( v , v ) <EOL> v = numpy . logical_and ( not_ref , not_target ) <EOL> tn = numpy . dot ( v , v ) <EOL> def GetCounts ( nucleotides , val ) : <EOL> """<STR_LIT>""" <EOL> t = <NUM_LIT:0> <EOL> for x in nucleotides : <EOL> if x == val : <EOL> t += <NUM_LIT:1> <EOL> return t <EOL> def AnalyseOverlaps ( references , targets ) : <EOL> """<STR_LIT>""" <EOL> if len ( references ) == <NUM_LIT:0> or len ( targets ) == <NUM_LIT:0> : <EOL> return ( ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) ) <EOL> targets_start = min ( map ( lambda x : x . start , targets ) ) <EOL> targets_end = max ( map ( lambda x : x . end , targets ) ) <EOL> refs_start = min ( map ( lambda x : x . start , references ) ) <EOL> refs_end = max ( map ( lambda x : x . end , references ) ) <EOL> min_start = min ( targets_start , refs_start ) <EOL> max_end = max ( targets_end , refs_end ) <EOL> nucleotides = numpy . zeros ( max_end - min_start , numpy . int0 ) <EOL> for x in targets [ : <NUM_LIT:100> ] : <EOL> nucleotides [ x . start - min_start : x . end - min_start ] = <NUM_LIT:1> <EOL> for x in references [ : <NUM_LIT:100> ] : <EOL> nucleotides [ x . start - min_start : x . end - min_start ] = numpy . array ( [ <NUM_LIT:2> ] * ( x . end - x . start ) , numpy . int0 ) <EOL> tp = GetCounts ( nucleotides , <NUM_LIT:3> ) <EOL> fp = GetCounts ( nucleotides , <NUM_LIT:1> ) <EOL> tn = GetCounts ( nucleotides , <NUM_LIT:0> ) <EOL> fn = GetCounts ( nucleotides , <NUM_LIT:2> ) <EOL> return ( tp , fp , tn , fn ) <EOL> def CalculateSpecificitySensitivity ( tp , fp , tn , fn ) : <EOL> if fp + tp == <NUM_LIT:0> : <EOL> spec = <NUM_LIT:1.0> <EOL> else : <EOL> spec = float ( tp ) / ( fp + tp ) <EOL> if tp + fn == <NUM_LIT:0> : <EOL> sens = <NUM_LIT:0.0> <EOL> else : <EOL> sens = float ( tp ) / ( tp + fn ) <EOL> return spec , sens <EOL> def CalculateCorrelationCoefficient ( tp , fp , tn , fn ) : <EOL> x = ( ( tp + fn ) * ( tn + fp ) * ( tp + fp ) * ( tn + fn ) ) <EOL> if x == <NUM_LIT:0> : <EOL> return <NUM_LIT:0.0> <EOL> else : <EOL> return ( tp * tn - fn * fp ) / x <EOL> def GetFirstOverlaps ( gffs , index ) : <EOL> """<STR_LIT>""" <EOL> this_name = gffs [ index ] . mName <EOL> this_strand = gffs [ index ] . strand <EOL> l = len ( gffs ) <EOL> min_start = gffs [ index ] . start <EOL> max_end = gffs [ index ] . end <EOL> overlaps = [ gffs [ index ] ] <EOL> index += <NUM_LIT:1> <EOL> while index < l : <EOL> if this_name != gffs [ index ] . mName or this_strand != gffs [ index ] . strand or max_end < gffs [ index ] . start : <EOL> return overlaps , index , min_start , max_end <EOL> max_end = max ( max_end , gffs [ index ] . end ) <EOL> overlaps . append ( gffs [ index ] ) <EOL> index += <NUM_LIT:1> <EOL> return overlaps , index , min_start , max_end <EOL> def CountMatchesPerGene ( gffs , <EOL> rx_gene , <EOL> rx_other , <EOL> write = ( ) , <EOL> outfile = sys . stdout ) : <EOL> """<STR_LIT>""" <EOL> map_gene2matches = { } <EOL> for x in gffs : <EOL> try : <EOL> gene = rx_gene . search ( x . mInfo ) . groups ( ) [ <NUM_LIT:0> ] <EOL> except AttributeError : <EOL> print "<STR_LIT>" % x . mInfo <EOL> sys . stdout . flush ( ) <EOL> sys . exit ( <NUM_LIT:1> ) <EOL> if gene not in map_gene2matches : <EOL> map_gene2matches [ gene ] = [ <NUM_LIT:0> , <NUM_LIT:0> , { } , { } ] <EOL> info = map_gene2matches [ gene ] <EOL> info [ <NUM_LIT:0> ] += <NUM_LIT:1> <EOL> try : <EOL> if x . mStatus == "<STR_LIT>" : <EOL> info [ <NUM_LIT:1> ] += <NUM_LIT:1> <EOL> elif x . mStatus in ( "<STR_LIT>" , "<STR_LIT>" ) : <EOL> for y in x . mMatches : <EOL> other_gene = rx_other . search ( y . mInfo ) . groups ( ) [ <NUM_LIT:0> ] <EOL> if other_gene not in info [ <NUM_LIT:2> ] : <EOL> info [ <NUM_LIT:2> ] [ other_gene ] = <NUM_LIT:0> <EOL> if other_gene not in info [ <NUM_LIT:3> ] : <EOL> info [ <NUM_LIT:3> ] [ other_gene ] = <NUM_LIT:0> <EOL> if GTF . Identity ( x , y , max_slippage = options . max_exon_slippage ) : <EOL> info [ <NUM_LIT:2> ] [ other_gene ] += <NUM_LIT:1> <EOL> info [ <NUM_LIT:3> ] [ other_gene ] += <NUM_LIT:1> <EOL> elif GTF . HalfIdentity ( x , y , max_slippage = options . max_exon_slippage ) : <EOL> info [ <NUM_LIT:3> ] [ other_gene ] += <NUM_LIT:1> <EOL> except AttributeError : <EOL> print "<STR_LIT>" , str ( x ) <EOL> total_extra = <NUM_LIT:0> <EOL> total_match = <NUM_LIT:0> <EOL> total_partial_match = <NUM_LIT:0> <EOL> total = len ( map_gene2matches ) <EOL> for gene , info in map_gene2matches . items ( ) : <EOL> l , nextra , matches , partials = info <EOL> if nextra == l : <EOL> total_extra += <NUM_LIT:1> <EOL> if "<STR_LIT>" in write : <EOL> outfile . write ( "<STR_LIT>" % ( gene , ) ) <EOL> else : <EOL> for x , y in matches . items ( ) : <EOL> if y == l : <EOL> total_match += <NUM_LIT:1> <EOL> if "<STR_LIT>" in write : <EOL> outfile . write ( "<STR_LIT>" % ( gene , x ) ) <EOL> break <EOL> else : <EOL> for x , y in partials . items ( ) : <EOL> if y == l : <EOL> total_partial_match += <NUM_LIT:1> <EOL> if "<STR_LIT>" in write : <EOL> outfile . write ( "<STR_LIT>" % ( gene , x ) ) <EOL> break <EOL> return total , total_match , total_partial_match , total_extra <EOL> def main ( argv = None ) : <EOL> """<STR_LIT>""" <EOL> if argv is None : <EOL> argv = sys . argv <EOL> parser = E . OptionParser ( <EOL> version = "<STR_LIT>" , usage = globals ( ) [ "<STR_LIT>" ] ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , <EOL> help = "<STR_LIT>" , action = "<STR_LIT:store_true>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , <EOL> help = "<STR_LIT>" , action = "<STR_LIT:store_true>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , action = "<STR_LIT:store_true>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , action = "<STR_LIT:store_true>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , dest = "<STR_LIT>" , action = "<STR_LIT>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , dest = "<STR_LIT>" , action = "<STR_LIT>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , dest = "<STR_LIT>" , action = "<STR_LIT>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . add_option ( "<STR_LIT>" , dest = "<STR_LIT>" , type = "<STR_LIT:string>" , <EOL> help = "<STR_LIT>" ) <EOL> parser . set_defaults ( <EOL> remove_redundancy = False , <EOL> max_exon_slippage = <NUM_LIT:9> , <EOL> write_missed_exons = False , <EOL> write_matched_exons = False , <EOL> write_missed_genes = False , <EOL> write_wrong_exons = False , <EOL> write_wrong_genes = False , <EOL> do_nucleotides = True , <EOL> do_exons = True , <EOL> do_genes = True , <EOL> regex_reference = None , <EOL> regex_target = None , <EOL> outfile_pattern = "<STR_LIT>" , <EOL> ) <EOL> ( options , args ) = E . Start ( parser ) <EOL> if len ( args ) != <NUM_LIT:2> : <EOL> print USAGE <EOL> print "<STR_LIT>" <EOL> sys . exit ( <NUM_LIT:1> ) <EOL> input_filename_target , input_filename_reference = args <EOL> if options . loglevel >= <NUM_LIT:1> : <EOL> print "<STR_LIT>" % input_filename_target <EOL> print "<STR_LIT>" , <EOL> sys . stdout . flush ( ) <EOL> gff_targets = GTF . readFromFile ( open ( input_filename_target , "<STR_LIT:r>" ) ) <EOL> if options . loglevel >= <NUM_LIT:1> : <EOL> print "<STR_LIT>" % ( len ( gff_targets ) ) <EOL> sys . stdout . flush ( ) <EOL> if options . loglevel >= <NUM_LIT:1> : <EOL> print "<STR_LIT>" % input_filename_reference <EOL> print "<STR_LIT>" , <EOL> sys . stdout . flush ( ) <EOL> gff_references = GTF . readFromFile ( open ( input_filename_reference , "<STR_LIT:r>" ) ) <EOL> if options . loglevel >= <NUM_LIT:1> : <EOL> print "<STR_LIT>" % ( len ( gff_references ) ) <EOL> sys . stdout . flush ( ) <EOL> if options . remove_redundancy : <EOL> gff_targets = GTF . CombineOverlaps ( gff_targets ) <EOL> gff_references = GTF . CombineOverlaps ( gff_references ) <EOL> if options . loglevel >= <NUM_LIT:1> : <EOL> print "<STR_LIT>" % ( len ( gff_targets ) , len ( gff_references ) ) <EOL> if options . loglevel >= <NUM_LIT:1> : <EOL> print "<STR_LIT>" , <EOL> sys . stdout . flush ( ) <EOL> gff_targets . sort ( lambda x , y : cmp ( ( x . mName , x . strand , x . start , x . end ) , <EOL> ( y . mName , y . strand , y . start , y . end ) ) ) <EOL> gff_references . sort ( lambda x , y : cmp ( ( x . mName , x . strand , x . start , x . end ) , <EOL> ( y . mName , y . strand , y . start , y . end ) ) ) <EOL> ntargets = len ( gff_targets ) <EOL> nreferences = len ( gff_references ) <EOL> if options . loglevel >= <NUM_LIT:1> : <EOL> print "<STR_LIT>" <EOL> sys . stdout . flush ( ) <EOL> if options . do_nucleotides : <EOL> print """<STR_LIT>""" <EOL> headers = ( <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) <EOL> print "<STR_LIT:\t>" . join ( headers ) <EOL> first_r , first_t = <NUM_LIT:0> , <NUM_LIT:0> <EOL> r , t = <NUM_LIT:0> , <NUM_LIT:0> <EOL> ttp , tfp , ttn , tfn = <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> <EOL> while r < nreferences and t < ntargets : <EOL> this_name = gff_references [ r ] . mName <EOL> this_strand = gff_references [ r ] . strand <EOL> while r < nreferences and gff_references [ r ] . mName == this_name and gff_references [ r ] . strand == this_strand : <EOL> r += <NUM_LIT:1> <EOL> while t < ntargets and ( gff_targets [ t ] . mName != this_name or <EOL> gff_targets [ t ] . strand != this_strand ) : <EOL> t += <NUM_LIT:1> <EOL> first_t = t <EOL> while t < ntargets and gff_targets [ t ] . mName == this_name and gff_targets [ t ] . strand == this_strand : <EOL> t += <NUM_LIT:1> <EOL> tp , fp , tn , fn = AnalyseOverlaps ( gff_references [ first_r : r ] , <EOL> gff_targets [ first_t : t ] ) <EOL> spec , sens = CalculateSpecificitySensitivity ( tp , fp , tn , fn ) <EOL> cc = CalculateCorrelationCoefficient ( tp , fp , tn , fn ) <EOL> print "<STR_LIT>" % ( this_name , this_strand , tp , fp , tn , fn , spec , sens , cc ) <EOL> ttp += tp <EOL> tfp += fp <EOL> ttn += tn <EOL> tfn += fn <EOL> first_r , first_t = r , t <EOL> spec , sens = CalculateSpecificitySensitivity ( ttp , tfp , ttn , tfn ) <EOL> cc = CalculateCorrelationCoefficient ( ttp , tfp , ttn , tfn ) <EOL> print "<STR_LIT>" % ( "<STR_LIT:all>" , "<STR_LIT:all>" , ttp , tfp , ttn , tfn , spec , sens , cc ) <EOL> sys . stdout . flush ( ) <EOL> if options . do_exons or options . do_genes : <EOL> print """<STR_LIT>""" % ( options . max_exon_slippage ) <EOL> headers = ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) <EOL> print "<STR_LIT:\t>" . join ( headers ) <EOL> r , t = <NUM_LIT:0> , <NUM_LIT:0> <EOL> next_r , next_t = r , t <EOL> tp , fp , tn , fn = <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> <EOL> ttp , tfp , ttn , tfn = <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> <EOL> ptp , pfp , ptn , pfn = <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> <EOL> tptp , tpfp , tptn , tpfn = <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> <EOL> missed_exons , wrong_exons = <NUM_LIT:0> , <NUM_LIT:0> <EOL> tmissed_exons , twrong_exons = <NUM_LIT:0> , <NUM_LIT:0> <EOL> last_partial_overlap = False <EOL> last_increased_ref = False <EOL> while r < nreferences and t < ntargets : <EOL> this_name = gff_references [ r ] . mName <EOL> this_strand = gff_references [ r ] . strand <EOL> if next_r == r : <EOL> ref_overlaps , next_r , ref_start , ref_end = GetFirstOverlaps ( <EOL> gff_references , r ) <EOL> if next_t == t : <EOL> target_overlaps , next_t , target_start , target_end = GetFirstOverlaps ( <EOL> gff_targets , t ) <EOL> if options . loglevel >= <NUM_LIT:3> : <EOL> print "<STR_LIT>" <EOL> for x in ref_overlaps : <EOL> print "<STR_LIT:#>" , str ( x ) <EOL> for x in target_overlaps : <EOL> print "<STR_LIT:#>" , str ( x ) <EOL> do_summary = False <EOL> if next_r < nreferences and ( this_name != gff_references [ next_r ] . mName or <EOL> this_strand != gff_references [ next_r ] . strand ) : <EOL> if options . loglevel >= <NUM_LIT:3> : <EOL> print "<STR_LIT>" <EOL> do_summary = True <EOL> last_increased_ref = False <EOL> last_partial_overlap = False <EOL> next_name = gff_references [ next_r ] . mName <EOL> next_strand = gff_references [ next_r ] . strand <EOL> while next_t < ntargets and next_name != gff_targets [ next_t ] . mName or next_strand != gff_targets [ next_t ] . strand : <EOL> fp += <NUM_LIT:1> <EOL> pfp += <NUM_LIT:1> <EOL> target_overlaps , next_t , target_start , target_end = GetFirstOverlaps ( <EOL> gff_targets , next_t ) <EOL> for x in gff_targets [ t : next_t ] : <EOL> x . mStatus = "<STR_LIT>" <EOL> for x in gff_references [ r : next_r ] : <EOL> x . mStatus = "<STR_LIT>" <EOL> r , t = next_r , next_t <EOL> elif next_t < ntargets and ( this_name != gff_targets [ next_t ] . mName or <EOL> this_strand != gff_targets [ next_t ] . strand ) : <EOL> if options . loglevel >= <NUM_LIT:3> : <EOL> print "<STR_LIT>" <EOL> do_summary = True <EOL> last_increased_ref = False <EOL> last_partial_overlap = False <EOL> next_name = gff_targets [ next_t ] . mName <EOL> next_strand = gff_targets [ next_t ] . strand <EOL> while next_r < nreferences and next_name != gff_references [ next_r ] . mName or next_strand != gff_references [ next_r ] . strand : <EOL> fn += <NUM_LIT:1> <EOL> pfn += <NUM_LIT:1> <EOL> reference_overlaps , next_r , references_start , references_end = GetFirstOverlaps ( <EOL> gff_references , next_r ) <EOL> for x in gff_targets [ t : next_t ] : <EOL> x . mStatus = "<STR_LIT>" <EOL> for x in gff_references [ r : next_r ] : <EOL> x . mStatus = "<STR_LIT>" <EOL> r , t = next_r , next_t <EOL> else : <EOL> ref_status , target_status = None , None <EOL> if options . loglevel >= <NUM_LIT:3> : <EOL> print "<STR_LIT>" <EOL> if min ( ref_end , target_end ) - max ( ref_start , target_start ) > <NUM_LIT:0> : <EOL> last_increased_ref = False <EOL> last_partial_overlap = False <EOL> found = False <EOL> for rr in ref_overlaps : <EOL> xfound = False <EOL> for tt in target_overlaps : <EOL> if GTF . Identity ( rr , tt , max_slippage = options . max_exon_slippage ) : <EOL> xfound = True <EOL> break <EOL> if xfound : <EOL> found = True <EOL> break <EOL> if found : <EOL> ref_status = "<STR_LIT>" <EOL> target_status = "<STR_LIT>" <EOL> tp += <NUM_LIT:1> <EOL> ptp += <NUM_LIT:1> <EOL> if options . write_matched_exons : <EOL> print "<STR_LIT>" <EOL> for x in ref_overlaps : <EOL> print "<STR_LIT:#>" , str ( x ) <EOL> for x in target_overlaps : <EOL> print "<STR_LIT:#>" , str ( x ) <EOL> else : <EOL> fn += <NUM_LIT:1> <EOL> for rr in ref_overlaps : <EOL> xfound = False <EOL> for tt in target_overlaps : <EOL> if GTF . HalfIdentity ( rr , tt , max_slippage = options . max_exon_slippage ) : <EOL> xfound = True <EOL> break <EOL> if xfound : <EOL> found = True <EOL> break <EOL> if found : <EOL> ptp += <NUM_LIT:1> <EOL> code = "<STR_LIT>" <EOL> ref_status = "<STR_LIT>" <EOL> target_status = "<STR_LIT>" <EOL> else : <EOL> pfn += <NUM_LIT:1> <EOL> code = "<STR_LIT>" <EOL> ref_status = "<STR_LIT>" <EOL> target_status = "<STR_LIT>" <EOL> if options . write_missed_exons : <EOL> print "<STR_LIT>" % code <EOL> for x in ref_overlaps : <EOL> print "<STR_LIT:#>" , str ( x ) <EOL> for x in target_overlaps : <EOL> print "<STR_LIT:#>" , str ( x ) <EOL> if ref_end == target_end : <EOL> r , t = next_r , next_t <EOL> elif ref_end < target_end : <EOL> r = next_r <EOL> last_increased_ref = True <EOL> last_partial_overlap = True <EOL> else : <EOL> t = next_t <EOL> last_increased_ref = False <EOL> last_partial_overlap = True <EOL> else : <EOL> if ref_end < target_start : <EOL> if not ( last_partial_overlap and not last_increased_ref ) : <EOL> if options . write_missed_exons : <EOL> print "<STR_LIT>" <EOL> for x in ref_overlaps : <EOL> print "<STR_LIT:#>" , str ( x ) <EOL> missed_exons += <NUM_LIT:1> <EOL> fn += <NUM_LIT:1> <EOL> pfn += <NUM_LIT:1> <EOL> ref_status = "<STR_LIT>" <EOL> r = next_r <EOL> else : <EOL> if not ( last_partial_overlap and last_increased_ref ) : <EOL> if options . write_wrong_exons : <EOL> print "<STR_LIT>" <EOL> for x in target_overlaps : <EOL> print "<STR_LIT:#>" , str ( x ) <EOL> wrong_exons += <NUM_LIT:1> <EOL> fp += <NUM_LIT:1> <EOL> pfp += <NUM_LIT:1> <EOL> target_status = "<STR_LIT>" <EOL> t = next_t <EOL> last_partial_overlap = False <EOL> if options . loglevel >= <NUM_LIT:3> : <EOL> print "<STR_LIT>" % ( ref_status , target_status ) <EOL> if ref_status : <EOL> for rr in ref_overlaps : <EOL> rr . mStatus = ref_status <EOL> if ref_status in ( "<STR_LIT>" , "<STR_LIT>" ) and options . do_genes : <EOL> for rr in ref_overlaps : <EOL> rr . mMatches = target_overlaps <EOL> if target_status : <EOL> for tt in target_overlaps : <EOL> tt . mStatus = target_status <EOL> if target_status in ( "<STR_LIT>" , "<STR_LIT>" ) and options . do_genes : <EOL> for tt in target_overlaps : <EOL> tt . mMatches = ref_overlaps <EOL> if do_summary or r >= nreferences or t >= ntargets : <EOL> ttp += tp <EOL> tfp += fp <EOL> ttn += tn <EOL> tfn += fn <EOL> tptp += ptp <EOL> tpfp += pfp <EOL> tptn += ptn <EOL> tpfn += pfn <EOL> tmissed_exons += missed_exons <EOL> twrong_exons += wrong_exons <EOL> if tp + fn != <NUM_LIT:0> : <EOL> pmissed_exons = "<STR_LIT>" % ( float ( missed_exons ) / ( tp + fn ) ) <EOL> else : <EOL> pmissed_exons = "<STR_LIT:0>" <EOL> if tp + fp != <NUM_LIT:0> : <EOL> pwrong_exons = "<STR_LIT>" % ( float ( wrong_exons ) / ( tp + fp ) ) <EOL> else : <EOL> pwrong_exons = "<STR_LIT>" <EOL> spec , sens = CalculateSpecificitySensitivity ( tp , fp , tn , fn ) <EOL> cc = ( spec + sens ) / <NUM_LIT> <EOL> print "<STR_LIT>" % ( this_name , this_strand , <EOL> tp , fp , tn , fn , <EOL> spec , sens , cc , <EOL> missed_exons , wrong_exons , <EOL> pmissed_exons , pwrong_exons ) <EOL> spec , sens = CalculateSpecificitySensitivity ( <EOL> ptp , pfp , ptn , pfn ) <EOL> cc = ( spec + sens ) / <NUM_LIT> <EOL> print "<STR_LIT>" % ( this_name , this_strand , <EOL> ptp , pfp , ptn , pfn , <EOL> spec , sens , cc , <EOL> missed_exons , wrong_exons , <EOL> pmissed_exons , pwrong_exons ) <EOL> tp , fp , tn , fn = <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> <EOL> ptp , pfp , ptn , pfn = <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> <EOL> missed_exons , wrong_exons = <NUM_LIT:0> , <NUM_LIT:0> <EOL> if t < ntargets : <EOL> for x in gff_targets [ t : ntargets ] : <EOL> x . mStatus = "<STR_LIT>" <EOL> if r < nreferences : <EOL> for x in gff_references [ r : nreferences ] : <EOL> x . mStatus = "<STR_LIT>" <EOL> spec , sens = CalculateSpecificitySensitivity ( ttp , tfp , ttn , tfn ) <EOL> cc = ( spec + sens ) / <NUM_LIT> <EOL> print "<STR_LIT>" % ( "<STR_LIT:all>" , "<STR_LIT:all>" , ttp , tfp , ttn , tfn , <EOL> spec , sens , cc , <EOL> tmissed_exons , twrong_exons , <EOL> float ( tmissed_exons ) / ( ttp + tfn ) , <EOL> float ( twrong_exons ) / ( ttp + tfp ) ) <EOL> spec , sens = CalculateSpecificitySensitivity ( tptp , tpfp , tptn , tpfn ) <EOL> cc = ( spec + sens ) / <NUM_LIT> <EOL> print "<STR_LIT>" % ( "<STR_LIT:all>" , "<STR_LIT:all>" , tptp , tpfp , tptn , tpfn , <EOL> spec , sens , cc , <EOL> tmissed_exons , twrong_exons , <EOL> float ( tmissed_exons ) / ( ttp + tfn ) , <EOL> float ( twrong_exons ) / ( ttp + tfp ) ) <EOL> if options . do_genes and options . regex_reference and options . regex_target : <EOL> print """<STR_LIT>""" % ( options . max_exon_slippage ) <EOL> out_options = [ ] <EOL> if options . write_missed_genes : <EOL> out_options . append ( "<STR_LIT>" ) <EOL> if options . loglevel >= <NUM_LIT:2> : <EOL> print "<STR_LIT>" <EOL> sys . stdout . flush ( ) <EOL> ( ref_total , ref_match , ref_partial , ref_extra ) = CountMatchesPerGene ( gff_references , <EOL> re . compile ( options . regex_reference ) , <EOL> re . compile ( options . regex_target ) , <EOL> write = out_options , <EOL> outfile = open ( options . outfile_pattern % "<STR_LIT>" , "<STR_LIT:w>" ) ) <EOL> if options . loglevel >= <NUM_LIT:2> : <EOL> print "<STR_LIT>" <EOL> sys . stdout . flush ( ) <EOL> ( target_total , target_match , target_partial , target_extra ) = CountMatchesPerGene ( gff_targets , <EOL> re . compile ( options . regex_target ) , <EOL> re . compile ( <EOL> options . regex_reference ) , <EOL> write = out_options , <EOL> outfile = open ( options . outfile_pattern % "<STR_LIT:target>" , "<STR_LIT:w>" ) ) <EOL> if options . loglevel >= <NUM_LIT:1> : <EOL> print "<STR_LIT>" % ( ref_total , ref_match , ref_partial , ref_extra ) <EOL> print "<STR_LIT>" % ( target_total , target_match , target_partial , target_extra ) <EOL> headers = ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ) <EOL> print "<STR_LIT:\t>" . join ( headers ) <EOL> tp = ref_match <EOL> fp = target_extra <EOL> tn = <NUM_LIT:0> <EOL> fn = ref_total - ref_match <EOL> wrong_genes = target_extra <EOL> missed_genes = ref_extra <EOL> spec , sens = CalculateSpecificitySensitivity ( tp , fp , tn , fn ) <EOL> cc = ( spec + sens ) / <NUM_LIT> <EOL> if tp + fp == <NUM_LIT:0> : <EOL> fp = nreferences <EOL> print "<STR_LIT>" % ( tp , fp , tn , fn , <EOL> spec , sens , cc , <EOL> missed_genes , wrong_genes , <EOL> float ( missed_genes ) / ( tp + fn ) , <EOL> float ( wrong_genes ) / ( tp + fp ) ) <EOL> tp = ref_match + ref_partial <EOL> fp = target_extra <EOL> tn = <NUM_LIT:0> <EOL> fn = ref_total - ref_match - ref_partial <EOL> wrong_genes = target_extra <EOL> missed_genes = ref_extra <EOL> spec , sens = CalculateSpecificitySensitivity ( tp , fp , tn , fn ) <EOL> cc = ( spec + sens ) / <NUM_LIT> <EOL> print "<STR_LIT>" % ( tp , fp , tn , fn , <EOL> spec , sens , cc , <EOL> missed_genes , wrong_genes , <EOL> float ( missed_genes ) / ( tp + fn ) , <EOL> float ( wrong_genes ) / ( tp + fp ) ) <EOL> E . Stop ( ) <EOL> if __name__ == "<STR_LIT:__main__>" : <EOL> sys . exit ( main ( sys . argv ) ) </s>
