Past exam of the mathematical and Computational Biology course of the University of Cambridge / 2024 / Paper 2 / C8

from collections import Counter
A = np.asarray(msa_matrix)
gap_fraction = np.mean(A == ".", axis=0)
match_mask = gap_fraction < .30
match_columns = np.flatnonzero(match_mask)
print(match_columns + 1)  # human-readable, one-based alignment columns

def emission_counts(column):
    counts = Counter(str(x) for x in column if x != ".")
    n = sum(counts.values())
    if n == 0:
        raise ValueError("an all-gap column has no empirical emission distribution")
    return counts, n

match_emissions = {}
for state, j in enumerate(match_columns, 1):
    counts, n = emission_counts(A[:, j])
    match_emissions[state] = {aa: count/n for aa, count in counts.items()}
    print(state, j+1, match_emissions[state])

information = []
for j in range(A.shape[1]):
    counts, n = emission_counts(A[:, j])
    p = np.array(list(counts.values()), dtype=float) / n
    H = -np.sum(p * np.log2(p))
    information.append(float(np.log2(20) - H))
print(np.round(information, 4))

01–10: 1.9025, 1.5964, 1.9025, 3.8187, 1.9025, 1.5964, 1.9025, 3.5577, 1.9025, 1.9025
11–20: 2.4420, 3.3355, 1.5964, 4.3219, 3.8187, 1.6803, 3.0975, 3.3355, 3.0975, 2.4420
21–30: 2.5692, 2.8219, 3.5106, 3.3219, 4.3219, 4.3219, 2.8000, 2.0703, 2.6574, 1.8003
31–40: 2.1938, 1.9163, 3.7784, 3.2606, 1.6803, 1.6803, 2.6642, 2.5692, 2.8753, 2.2631
41–50: 2.8753, 3.3355, 1.6803, 2.9703, 3.3355, 2.7914, 2.2086, 2.5692, 4.3219, 1.3742
51–60: 2.2631, 4.3219, 4.3219, 3.8187, 4.3219, 3.5577, 2.8753, 1.9025, 2.7914, 2.0409
61–70: 3.3355, 2.7370, 2.4309, 2.4309, 2.5692, 2.0409, 2.5692, 3.3355, 2.9298, 3.5577
71–79: 1.9025, 1.8187, 4.3219, 2.2631, 3.8187, 2.4309, 1.9025, 1.9025, 1.9025

insert_columns = np.flatnonzero(~match_mask)
runs = np.split(insert_columns, np.flatnonzero(np.diff(insert_columns) > 1) + 1)
insert_emissions = []
for run in runs:
    if len(run):
        counts, n = emission_counts(A[:, run].ravel())
        insert_emissions.append({aa: count/n for aa, count in counts.items()})
        print(run+1, n, insert_emissions[-1])