code: def variation(fastas, sams, s2b, annots, qual, min_cov, min_per, no_sub, bp_table, print_con):
      
      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
docs: quantify genomic variation from mapping
      genomes[genome][contig][sample] = \
              {'bp_stats':[], 'sub_rates'[locus] = {ref PnPs, consensus PnPs, ...}}