calculate syn and non-syn sub. rates for each locus/gene

def sub_rates(stats, annots):
    """
    calculate syn and non-syn sub. rates for each locus/gene
    """
    nucs = ['A', 'T', 'G', 'C']
    subs = {} # subs[locus]{pn/ps reference, pn/ps consensus, snp density}
    for cds in annots:
        locus, location, info = cds
        subs[locus] = {}
        start, end, strand = location
        trans_table = info['transl_table'][0]
        wrt_reference = {'obs syn':0, 'pos syn':0, 'obs non-syn':0, 'pos non-syn':0}
        wrt_consensus = {'obs syn':0, 'pos syn':0, 'obs non-syn':0, 'pos non-syn':0} 
        for codon in parse_codons(stats, start, end, strand):
            # reference and consensus translations
            ref_codon = [i['ref'] for i in codon]
            con_codon = [i['consensus'][0] for i in codon]
            ref_aa = codon2aa(ref_codon, trans_table)
            con_aa = codon2aa(con_codon, trans_table)
            # count types of mutations with respect to reference
            wrt_reference = \
                    count_mutations(ref_codon, ref_aa, codon, wrt_reference, nucs, trans_table) 
            # count types of mutations with respect to consensus 
            wrt_consensus = \
                    count_mutations(con_codon, con_aa, codon, wrt_consensus, nucs, trans_table)
        # calculate pn/ps
        subs[locus]['ref PnPs']       = calc_PnPs(wrt_reference)
        subs[locus]['consensus PnPs'] = calc_PnPs(wrt_consensus)
        # SNP density
        locus_length = end - start + 1
        subs[locus]['ref SNP density'] \
                = float(wrt_reference['obs syn'] + wrt_reference['obs non-syn']) / locus_length
        subs[locus]['consensus SNP density'] \
                = float(wrt_consensus['obs syn'] + wrt_consensus['obs non-syn']) / locus_length
        # sequencing coverage
        subs[locus]['cov'] = calc_coverage(stats, start, end, locus_length, nucs)
        info['length'] = locus_length
        info['position'] = ((start, end), strand)
        # point to cds info
        subs[locus]['info'] = info
    return subs