<s> '''<STR_LIT>''' <EOL> import sys <EOL> import os <EOL> import pysam <EOL> from ngsutils . bam import bam_iter <EOL> def usage ( msg = None ) : <EOL> if msg : <EOL> sys . stdout . write ( '<STR_LIT>' % msg ) <EOL> sys . stdout . write ( __doc__ ) <EOL> sys . stdout . write ( '''<STR_LIT>''' ) <EOL> sys . exit ( <NUM_LIT:1> ) <EOL> def __flush_cur_reads ( cur_reads , outbam , inbam , countfile = None ) : <EOL> if cur_reads : <EOL> for k in cur_reads : <EOL> count = <NUM_LIT:0> <EOL> for i , ( mapq , idx , r ) in enumerate ( sorted ( cur_reads [ k ] ) [ : : - <NUM_LIT:1> ] ) : <EOL> count += <NUM_LIT:1> <EOL> if i > <NUM_LIT:0> : <EOL> r . is_duplicate = True <EOL> if outbam : <EOL> outbam . write ( r ) <EOL> if countfile : <EOL> countfile . write ( '<STR_LIT>' % ( inbam . references [ k [ <NUM_LIT:0> ] ] , k [ <NUM_LIT:1> ] , k [ <NUM_LIT:2> ] , count ) ) <EOL> def pcrdup_mark ( inbam , outbam , fragment = False , countfile = None ) : <EOL> cur_pos = None <EOL> cur_reads = { } <EOL> total = <NUM_LIT:0> <EOL> unique = <NUM_LIT:0> <EOL> duplicates = <NUM_LIT:0> <EOL> dup_list = set ( ) <EOL> def callback ( read ) : <EOL> return '<STR_LIT>' % ( total , unique , duplicates , read . qname ) <EOL> for read in bam_iter ( bamfile , callback = callback ) : <EOL> if not read . is_paired or read . is_read1 : <EOL> total += <NUM_LIT:1> <EOL> if read . is_unmapped : <EOL> __flush_cur_reads ( cur_reads , outbam , inbam , countfile ) <EOL> if outbam : <EOL> outbam . write ( read ) <EOL> continue <EOL> start_pos = ( read . tid , read . pos ) <EOL> if fragment : <EOL> dup_pos = ( read . tid , read . pos , '<STR_LIT>' ) <EOL> else : <EOL> dup_pos = ( read . tid , read . pos , read . isize ) <EOL> if not cur_pos or start_pos != cur_pos : <EOL> __flush_cur_reads ( cur_reads , outbam , inbam , countfile ) <EOL> cur_pos = start_pos <EOL> cur_reads = { } <EOL> idx = <NUM_LIT:0> <EOL> if not fragment and ( read . mate_is_unmapped or not read . is_paired or not read . is_proper_pair or read . isize < <NUM_LIT:0> ) : <EOL> if read . qname in dup_list : <EOL> read . is_duplicate = True <EOL> dup_list . remove ( read . qname ) <EOL> if outbam : <EOL> outbam . write ( read ) <EOL> elif dup_pos in cur_reads : <EOL> duplicates += <NUM_LIT:1> <EOL> if not fragment : <EOL> dup_list . add ( read . qname ) <EOL> cur_reads [ dup_pos ] . append ( ( read . mapq , - idx , read ) ) <EOL> else : <EOL> unique += <NUM_LIT:1> <EOL> cur_reads [ dup_pos ] = [ ( read . mapq , - idx , read ) , ] <EOL> idx += <NUM_LIT:1> <EOL> __flush_cur_reads ( cur_reads , outbam , inbam , countfile ) <EOL> sys . stdout . write ( '<STR_LIT>' % total ) <EOL> sys . stdout . write ( '<STR_LIT>' % unique ) <EOL> sys . stdout . write ( '<STR_LIT>' % duplicates ) <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> infile = None <EOL> outfile = None <EOL> countfname = None <EOL> fragment = False <EOL> last = None <EOL> for arg in sys . argv [ <NUM_LIT:1> : ] : <EOL> if arg == '<STR_LIT>' : <EOL> usage ( ) <EOL> elif last == '<STR_LIT>' : <EOL> countfname = arg <EOL> last = None <EOL> elif last == '<STR_LIT>' : <EOL> outfile = arg <EOL> last = None <EOL> elif arg in [ '<STR_LIT>' , '<STR_LIT>' ] : <EOL> last = arg <EOL> elif arg == '<STR_LIT>' : <EOL> fragment = True <EOL> elif not infile : <EOL> if os . path . exists ( arg ) : <EOL> infile = arg <EOL> else : <EOL> usage ( "<STR_LIT>" % arg ) <EOL> elif not outfile : <EOL> if not os . path . exists ( arg ) : <EOL> outfile = arg <EOL> else : <EOL> usage ( "<STR_LIT>" % arg ) <EOL> if not infile or not ( outfile or countfname ) : <EOL> usage ( ) <EOL> bamfile = pysam . Samfile ( infile , "<STR_LIT:rb>" ) <EOL> bamout = None <EOL> if outfile : <EOL> bamout = pysam . Samfile ( outfile , "<STR_LIT:wb>" , template = bamfile ) <EOL> if countfname : <EOL> countfile = open ( countfname , '<STR_LIT:w>' ) <EOL> else : <EOL> countfile = None <EOL> pcrdup_mark ( bamfile , bamout , fragment , countfile ) <EOL> bamfile . close ( ) <EOL> if bamout : <EOL> bamout . close ( ) <EOL> if countfile : <EOL> countfile . close ( ) </s>
