def print_stats(genomes):
    """
    print substitution rate data to table
    genomes[genome][contig][sample] = \
            {'bp_stats':{}, 'sub_rates'[locus] = {ref PnPs, consensus PnPs}}
    """
    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]))