code: def print_stats(genomes):
      
      header = ['#genome', 'contig', 'locus', 'position', 'strand', 'length', \
                  'sample', 'coverage', \
                  'ref. Pn/Ps', 'ref SNP density', \
                  'consensus Pn/Ps', 'consensus SNP density', \
                  'product']
      print('\t'.join(header))
      for genome, contigs in list(genomes.items()):
          for contig, samples in list(contigs.items()):
              for sample, stats in list(samples.items()):
                  for locus, rates in list(stats['sub_rates'].items()):
                      length = rates['info']['length']
                      position, strand = rates['info']['position']
                      position = '%s-%s' % position
                      out = [genome, contig, locus, position, strand, length, \
                              sample, '%.2f' % (rates['cov']), \
                              rates['ref PnPs'], rates['ref SNP density'], \
                              rates['consensus PnPs'], rates['consensus SNP density'], \
                              rates['info']['product'][0]]
                      print('\t'.join([str(i) for i in out]))
docs: print substitution rate data to table
      genomes[genome][contig][sample] = \
              {'bp_stats':{}, 'sub_rates'[locus] = {ref PnPs, consensus PnPs}}