<s> import sys <EOL> from itertools import ifilter , izip <EOL> from functools import partial <EOL> from operator import itemgetter <EOL> import numpy as np <EOL> from nucio import fileIterator , lineItemIterator <EOL> MAPQ_THRESH = <NUM_LIT:30> <EOL> name = itemgetter ( <NUM_LIT:0> ) <EOL> seq = itemgetter ( <NUM_LIT:2> ) <EOL> qual = itemgetter ( <NUM_LIT:3> ) <EOL> status = itemgetter ( <NUM_LIT:4> ) <EOL> ascore = itemgetter ( <NUM_LIT:5> ) <EOL> aqual = lambda r : int ( itemgetter ( <NUM_LIT:6> ) ( r ) ) <EOL> ref = itemgetter ( <NUM_LIT:7> ) <EOL> pos = lambda r : int ( itemgetter ( <NUM_LIT:8> ) ( r ) ) <EOL> strand = itemgetter ( <NUM_LIT:9> ) <EOL> if not len ( sys . argv ) == <NUM_LIT:3> : <EOL> sys . exit ( "<STR_LIT>" ) <EOL> names_eq = lambda name1 , name2 : name1 . split ( "<STR_LIT:/>" ) [ <NUM_LIT:0> ] == name2 . split ( "<STR_LIT:/>" ) [ <NUM_LIT:0> ] <EOL> filenames = sys . argv [ <NUM_LIT:1> : <NUM_LIT:3> ] <EOL> headfilt = lambda x : not x . startswith ( "<STR_LIT:#>" ) <EOL> filt_lii = partial ( lineItemIterator , filter_func = headfilt ) <EOL> filt_fits = map ( lambda fn : fileIterator ( fn , filt_lii ) , filenames ) <EOL> failrepeat = <NUM_LIT:0> <EOL> failmapq = <NUM_LIT:0> <EOL> failsameref = <NUM_LIT:0> <EOL> total = <NUM_LIT:0> <EOL> passed = <NUM_LIT:0> <EOL> insertNotRF = [ ] <EOL> insertRF = [ ] <EOL> for read1 , read2 in izip ( * filt_fits ) : <EOL> total += <NUM_LIT:1> <EOL> if not names_eq ( name ( read1 ) , name ( read2 ) ) : <EOL> sys . exit ( "<STR_LIT>" % ( name ( read1 ) , name ( read2 ) ) ) <EOL> if not status ( read1 ) == "<STR_LIT>" or not status ( read2 ) == "<STR_LIT>" : <EOL> failrepeat += <NUM_LIT:1> <EOL> continue <EOL> if not aqual ( read1 ) >= MAPQ_THRESH or not aqual ( read2 ) >= MAPQ_THRESH : <EOL> failmapq += <NUM_LIT:1> <EOL> continue <EOL> if ref ( read1 ) == ref ( read2 ) : <EOL> failsameref += <NUM_LIT:1> <EOL> if not ( strand ( read1 ) == "<STR_LIT:R>" and strand ( read2 ) == "<STR_LIT:F>" ) : <EOL> insertNotRF . append ( pos ( read2 ) - pos ( read1 ) ) <EOL> else : <EOL> insertRF . append ( pos ( read2 ) - pos ( read1 ) ) <EOL> continue <EOL> print "<STR_LIT:\n>" . join ( map ( "<STR_LIT:\t>" . join , [ read1 , read2 ] ) ) <EOL> passed += <NUM_LIT:1> <EOL> with open ( "<STR_LIT>" , "<STR_LIT:w>" ) as ofh : <EOL> ofh . write ( "<STR_LIT>" ) <EOL> ofh . write ( "<STR_LIT>" % failrepeat ) <EOL> ofh . write ( "<STR_LIT>" % ( MAPQ_THRESH , failmapq ) ) <EOL> ofh . write ( "<STR_LIT>" % failsameref ) <EOL> ofh . write ( "<STR_LIT>" % ( passed , total , float ( passed ) / total ) ) <EOL> ofh . write ( "<STR_LIT>" ) <EOL> h = np . histogram ( insertNotRF , bins = <NUM_LIT:1000> ) <EOL> for v , b in zip ( * h ) : <EOL> ofh . write ( "<STR_LIT>" % ( b , v ) ) <EOL> ofh . write ( "<STR_LIT>" ) <EOL> ofh . write ( "<STR_LIT>" ) <EOL> h = np . histogram ( insertRF , bins = <NUM_LIT:1000> ) <EOL> for v , b in zip ( * h ) : <EOL> ofh . write ( "<STR_LIT>" % ( b , v ) ) </s>
