def get_bp_stats(sam, genomes, s2b, qual):
for read in parse_sam(sam, qual):<EOL><INDENT>ref = read[<NUM_LIT:2>] <EOL>if ref not in s2b:<EOL><INDENT>continue<EOL><DEDENT>genome = s2b[ref] <EOL>refs = genomes[genome][ref][sam.name]['<STR_LIT>']<EOL>ref_start = int(read[<NUM_LIT:3>]) - <NUM_LIT:1> <EOL>ref_pos = int(read[<NUM_LIT:3>]) - <NUM_LIT:1> <EOL>sequence = list(read[<NUM_LIT:9>]) <EOL>cigar = parse_cigar(read[<NUM_LIT:5>]) <EOL>cigar_start = <NUM_LIT:0> <EOL>bases = [] <EOL>for cigar_pos, status in cigar:<EOL><INDENT>if status == '<STR_LIT:D>': <EOL><INDENT>refs[ref_pos - <NUM_LIT:1>]['<STR_LIT>'].append(cigar_pos)<EOL>for b in range(<NUM_LIT:0>, cigar_pos):<EOL><INDENT>bases.append(False)<EOL><DEDENT>ref_pos += cigar_pos<EOL><DEDENT>else:<EOL><INDENT>cigar_stop = cigar_start + cigar_pos<EOL>if status == '<STR_LIT:M>': <EOL><INDENT>for b in sequence[cigar_start:cigar_stop]: <EOL><INDENT>bases.append(b)<EOL><DEDENT>ref_pos += cigar_pos<EOL><DEDENT>elif status == '<STR_LIT:I>': <EOL><INDENT>refs[ref_pos - <NUM_LIT:1>]['<STR_LIT>'].append(cigar_pos)<EOL><DEDENT>else:<EOL><INDENT>print('<STR_LIT>' % (status), file=sys.stderr)<EOL>exit()<EOL><DEDENT>cigar_start += cigar_pos<EOL><DEDENT><DEDENT>for base, position in zip(bases, list(range(ref_start, ref_start + len(bases)))):<EOL><INDENT>if base is False:<EOL><INDENT>continue<EOL><DEDENT>try:<EOL><INDENT>refs[position][base.upper()] += <NUM_LIT:1><EOL><DEDENT>except IndexError:<EOL><INDENT>continue<EOL><DEDENT><DEDENT><DEDENT>return genomes<EOL>
save base pair frequency statistics for each base from sam file
genomes = {} # genomes[genome][contig][sample] = {'bp_stats':{}}