quantify genomic variation from mapping
    genomes[genome][contig][sample] = \
            {'bp_stats':[], 'sub_rates'[locus] = {ref PnPs, consensus PnPs, ...}}

def variation(fastas, sams, s2b, annots, qual, min_cov, min_per, no_sub, bp_table, print_con):
    """
    quantify genomic variation from mapping
    genomes[genome][contig][sample] = \
            {'bp_stats':[], 'sub_rates'[locus] = {ref PnPs, consensus PnPs, ...}}
    """
    genomes = {}

    # save base-pair data structures for each genome-sample pair
    for sam in sams:
        genomes, s2b = save_refs(sam, fastas, genomes, s2b)

    # get base-pair frequencies for each genome-sample pair
    for sam in sams:
        genomes = get_bp_stats(sam, genomes, s2b, qual)

    # calculate base-pair frequencies
    genomes = calc_frequencies(genomes, bp_table, min_cov, min_per)

    # print consensus genome
    if print_con is True:
        print('# saving consensus sequences', file = sys.stderr)
        genomes = print_consensus(genomes)

    if no_sub is True:
        return genomes

    # calculate substitution rates
    for genome, contigs in list(genomes.items()):
        for contig, samples in list(contigs.items()):
            for sample in samples:
#                print(genomes.keys())
#                print(genomes[genome].keys())
#                print(genomes[genome][contig].keys())
#                print(annots.keys())
                genomes[genome][contig][sample]['sub_rates'] = \
                    sub_rates(genomes[genome][contig][sample]['bp_stats'], annots[contig])

    return genomes