#NEWLINE# def buildNormalizedBAM( infiles , outfile , normalize = True ) : #NEWLINE# #INDENT# 'build a normalized BAM file.\n\n    Infiles are merged and duplicated reads are removed. \n    If *normalize* is set, reads are removed such that all \n    files will have approximately the same number of reads.\n    ' #NEWLINE# min_reads = getMinimumMappedReads ( glob . glob ( '*.readstats' ) ) #NEWLINE# samfiles = [ ] #NEWLINE# num_reads = 0 #NEWLINE# for ( infile , statsfile ) in infiles : #NEWLINE# #INDENT# samfiles . append ( pysam . Samfile ( infile , 'rb' ) ) #NEWLINE# num_reads += getMappedReads ( statsfile ) #NEWLINE# #UNINDENT# threshold = ( float ( min_reads ) / num_reads ) #NEWLINE# pysam_out = pysam . Samfile ( outfile , 'wb' , template = samfiles [ 0 ] ) #NEWLINE# ( ninput , noutput , nduplicates ) = ( 0 , 0 , 0 ) #NEWLINE# ( last_contig , last_pos ) = ( None , None ) #NEWLINE# for pysam_in in samfiles : #NEWLINE# #INDENT# for read in pysam_in . fetch ( ) : #NEWLINE# #INDENT# ninput += 1 #NEWLINE# if ( ( read . rname == last_contig ) and ( read . pos == last_pos ) ) : #NEWLINE# #INDENT# nduplicates += 1 #NEWLINE# continue #NEWLINE# #UNINDENT# if ( normalize and ( random . random ( ) <= threshold ) ) : #NEWLINE# #INDENT# pysam_out . write ( read ) #NEWLINE# noutput += 1 #NEWLINE# #UNINDENT# ( last_contig , last_pos ) = ( read . rname , read . pos ) #NEWLINE# #UNINDENT# pysam_in . close ( ) #NEWLINE# #UNINDENT# pysam_out . close ( ) #NEWLINE# logs = open ( ( outfile + '.log' ) , 'w' ) #NEWLINE# logs . write ( ( '# min_reads=%i, threshold= %5.2f\n' % ( min_reads , threshold ) ) ) #NEWLINE# logs . write ( 'set\tcounts\tpercent\n' ) #NEWLINE# logs . write ( ( 'ninput\t%i\t%5.2f%%\n' % ( ninput , 100.0 ) ) ) #NEWLINE# nwithout_dups = ( ninput - nduplicates ) #NEWLINE# logs . write ( ( 'duplicates\t%i\t%5.2f%%\n' % ( nduplicates , ( ( 100.0 * nduplicates ) / ninput ) ) ) ) #NEWLINE# logs . write ( ( 'without duplicates\t%i\t%5.2f%%\n' % ( nwithout_dups , ( ( 100.0 * nwithout_dups ) / ninput ) ) ) ) #NEWLINE# logs . write ( ( 'target\t%i\t%5.2f%%\n' % ( min_reads , ( ( 100.0 * min_reads ) / nwithout_dups ) ) ) ) #NEWLINE# logs . write ( ( 'noutput\t%i\t%5.2f%%\n' % ( noutput , ( ( 100.0 * noutput ) / nwithout_dups ) ) ) ) #NEWLINE# logs . close ( ) #NEWLINE# if ( len ( samfiles ) > 1 ) : #NEWLINE# #INDENT# tmpfilename = P . getTempFilename ( ) #NEWLINE# pysam . sort ( outfile , tmpfilename ) #NEWLINE# shutil . move ( ( tmpfilename + '.bam' ) , outfile ) #NEWLINE# os . unlink ( tmpfilename ) #NEWLINE# #UNINDENT# pysam . index ( outfile ) #NEWLINE# P . info ( ( 'buildNormalizedBam: %i input, %i output (%5.2f%%), should be %i' % ( ninput , noutput , ( ( 100.0 * noutput ) / ninput ) , min_reads ) ) )