def asm(chrom, start, end, bamfilename, reffile, kmersize, tmpdir, mutid='null', debug=False):
    bamfile  = pysam.Samfile(bamfilename,'rb')
    matefile = pysam.Samfile(bamfilename,'rb')

    readpairs = {}
    nreads = 0
    ndisc  = 0 # track discordant reads
    rquals = []
    mquals = []


    for read in bamfile.fetch(chrom,start,end):
        if not read.mate_is_unmapped and read.is_paired:
            try:
                mate = matefile.mate(read)
                readpairs[read.qname] = ReadPair(read,mate)
                nreads += 1
                if not read.is_proper_pair:
                    ndisc  += 1
                if nreads % 1000 == 0:
                    print "INFO\t" + now() + "\t" + mutid + "\tfound mates for", nreads, "reads,", float(ndisc)/float(nreads), "discordant."
                if read.is_read1:
                    if read.is_reverse:
                        rquals.append(read.qual[::-1])
                        mquals.append(mate.qual)
                    else:
                        rquals.append(read.qual)
                        mquals.append(mate.qual[::-1])
                else:
                    if read.is_reverse:
                        rquals.append(mate.qual)
                        mquals.append(read.qual[::-1])
                    else:
                        rquals.append(mate.qual[::-1])
                        mquals.append(read.qual)
            except __HOLE__:
                sys.stderr.write("WARN\t" + now() + "\t" + mutid + "\tcannot find mate for read marked paired: " + read.qname + "\n")

    sys.stdout.write("INFO\t" + now() + "\t" + mutid + "\tfound " + str(nreads) + " reads in region.\n")

    if nreads == 0:
        return []

    refseq = None
    if reffile:
        refseq = reffile.fetch(chrom,start,end)

    region = "%s:%d-%d" % (chrom, start, end)

    contigs = runVelvet(readpairs, region, refseq, kmersize, tmpdir, cov_cutoff=True, mutid=mutid, debug=debug)
    newcontigs = None

    for contig in contigs:
        contig.rquals = rquals
        contig.mquals = mquals
    return contigs
dataset/ETHPy150Open adamewing/bamsurgeon/bamsurgeon/asmregion.py/asm