def get_bp_stats(sam, genomes, s2b, qual):
    """
    save base pair frequency statistics for each base from sam file
    genomes = {} # genomes[genome][contig][sample] = {'bp_stats':{}}
    """
    for read in parse_sam(sam, qual):
        ref = read[2] # scaffold that read mapped to
        if ref not in s2b:
            continue
        genome = s2b[ref] # genome that scaffold belongs to
        refs = genomes[genome][ref][sam.name]['bp_stats']
        ref_start = int(read[3]) - 1 # position of start of alignment on reference
        ref_pos = int(read[3]) - 1 # for keeping track of reference region
        sequence = list(read[9]) # read sequence
        cigar = parse_cigar(read[5]) # parsed cigar string
        cigar_start = 0 # for keeping track of start of cigar regions
        bases = [] # bases to compare with reference
        for cigar_pos, status in cigar:
            if status == 'D': # deletion compared to reference
                refs[ref_pos - 1]['Del'].append(cigar_pos)
                for b in range(0, cigar_pos):
                    bases.append(False)
                ref_pos += cigar_pos
            else:
                cigar_stop = cigar_start + cigar_pos
                if status == 'M': # aligned to reference
                    for b in sequence[cigar_start:cigar_stop]: # bases for cigar region
                        bases.append(b)
                    ref_pos += cigar_pos
                elif status == 'I': # insertion compared to reference
                    refs[ref_pos - 1]['In'].append(cigar_pos)
                else:
                    print('# unrecognized cigar character: %s' % (status), file=sys.stderr)
                    exit()
                cigar_start += cigar_pos
        # add base to frequency at each position
        for base, position in zip(bases, list(range(ref_start, ref_start + len(bases)))):
            if base is False:
                continue
            try:
                refs[position][base.upper()] += 1
            except IndexError:
                continue
    return genomes
save base pair frequency statistics for each base from sam file
    genomes = {} # genomes[genome][contig][sample] = {'bp_stats':{}}