def replaceReads(targetbam, donorbam, outputbam, nameprefix=None, excludefile=None, allreads=False, keepqual=False, progress=False, keepsecondary=False, seed=None):
    ''' targetbam, donorbam, and outputbam are pysam.Samfile objects
        outputbam must be writeable and use targetbam as template
        read names in excludefile will not appear in final output
    '''

    if seed is not None: random.seed(int(seed))

    # check whether references are compatible
    if not compare_ref(targetbam, donorbam):
        sys.exit("Target and donor are aligned to incompatable reference genomes!")

    RG = getRGs(targetbam) # read groups

    exclude = {}
    if excludefile:
        exclude = getExcludedReads(excludefile)

    # load reads from donorbam into dict 
    sys.stdout.write("loading donor reads into dictionary...\n")
    nr = 0
    rdict = {}
    secondary = [] # track secondary alignments, if specified
    excount = 0 # number of excluded reads
    nullcount = 0 # number of null reads

    for read in donorbam.fetch(until_eof=True):
        if read.seq and not read.is_secondary: # sanity check - don't include null reads, secondary alignments, supplementary alignments.
            if read.qname not in exclude:
                pairname = 'F' # read is first in pair
                if read.is_read2:
                    pairname = 'S' # read is second in pair
                if not read.is_paired:
                    pairname = 'U' # read is unpaired
                if nameprefix:
                    qual = read.qual # temp
                    read.qname = nameprefix + read.qname # must set name _before_ setting quality (see pysam docs)
                    read.qual = qual
                extqname = ','.join((read.qname,pairname))
                rdict[extqname] = read
                nr += 1
            else: # excluded
                excount += 1
        elif (read.is_secondary and keepsecondary):
            secondary.append(read)
        else: # no seq!
            nullcount += 1

    sys.stdout.write("loaded " + str(nr) + " reads, (" + str(excount) + " excluded, " + str(nullcount) + " null or secondary --> ignored)\n")

    excount = 0
    recount = 0 # number of replaced reads
    used = {}
    prog = 0
    for read in targetbam.fetch(until_eof=True):

        prog += 1
        if progress and prog % 10000000 == 0:
            sys.stdout.write("processed " + str(prog) + " reads.\n")

        if not read.is_secondary and read.qname not in exclude and bin(read.flag & 2048) != bin(2048):
            pairname = 'F' # read is first in pair
            if read.is_read2:
                pairname = 'S' # read is second in pair
            if not read.is_paired:
                pairname = 'U' # read is unpaired
            if nameprefix:
                qual = read.qual # temp
                read.qname = nameprefix + read.qname
                read.qual = qual

            extqname = ','.join((read.qname,pairname))
            if extqname in rdict: # replace read
                if keepqual:
                    try:
                        rdict[extqname].qual = read.qual
                    except __HOLE__ as e:
                        sys.stdout.write("error replacing quality score for read: " + str(rdict[extqname].qname) + " : " + str(e) + "\n")
                        sys.stdout.write("donor:  " + str(rdict[extqname]) + "\n")
                        sys.stdout.write("target: " + str(read) + "\n")
                        sys.exit(1)

                rdict[extqname] = cleanup(rdict[extqname],read,RG)
                outputbam.write(rdict[extqname])  # write read from donor .bam
                used[extqname] = True
                recount += 1
            else:
                read = cleanup(read,None,RG)
                outputbam.write(read) # write read from target .bam
        else:
            excount += 1

    sys.stdout.write("replaced " + str(recount) + " reads (" + str(excount) + " excluded )\n")

    if keepsecondary:
        for secread in secondary:
            outputbam.write(secread)

    sys.stdout.write("kept " + str(len(secondary)) + " secondary reads.\n")

    nadded = 0
    # dump the unused reads from the donor if requested with --all
    if allreads:
        for extqname in rdict.keys():
            if extqname not in used and extqname not in exclude:
                rdict[extqname] = cleanup(rdict[extqname],None,RG)
                outputbam.write(rdict[extqname])
                nadded += 1
        sys.stdout.write("added " + str(nadded) + " reads due to --all\n")
dataset/ETHPy150Open adamewing/bamsurgeon/bamsurgeon/replacereads.py/replaceReads