In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from scipy.stats import norm

Задание 1¶

Задание 1.1 Генерация данных¶

Ниже имеются гистограммы, отражающие распределение сгенерированных данных.

In [ ]:
np.random.seed(42)
n1 = int(input())
n2 = int(input())
mu1, sd1, mu2, sd2 = 0, 1, 1, 1
x1 = np.random.normal(mu1, sd1, n1)
x2 = np.random.normal(mu2, sd2, n2)
x = np.concatenate([x1, x2])
answers = [0 for i in range(n1)]
answers.extend([1 for i in range(n2)])
100
5

Задание 1.2 EM-оценка¶

Формула логарифма правдоподобия, с помощью которого подбирается оптимальное k на каждом шаге оптимизации

image.png

Шаг 0: начинаем со случайных разбиений
Шаг 1: оценка параметров разбиений (шаг оценки)
Шаг 2: оптимизация разбиений с помощью логарифма правдоподобия (шаг максимизации) --> вернуться к шагу 1

In [ ]:
# @title EM алгоритм
def maximize_likelihood(arr, exp1, exp2, disp1, disp2):
  likelihoods = []
  sorted_arr = sorted(arr)
  for i in range(2, len(arr) - 2):
    arr1, arr2 = sorted_arr[:i], sorted_arr[i + 1:]
    likelihood_log = - np.sum((arr1 - exp1) ** 2) / (2 * disp1 ** 2) - np.sum((arr2 - exp2) ** 2) / (2 * disp2 ** 2)
    likelihoods.append((likelihood_log, i))
  return max(likelihoods, key=lambda x: x[0])


def expectation_maximization(arr, convergence_threshold=0.1, n_convergence=15):

  k = np.random.randint(2, len(arr) - 2)
  x1_p, x2_p = sorted(arr)[:k], sorted(arr)[k + 1:]
  prev_curr_likelihood, cnt = 0, 0
  while cnt < n_convergence:
    # Expectation step
    e1, e2, d1, d2 = np.mean(x1_p), np.mean(x2_p), np.std(x1_p, ddof=1), np.std(x2_p, ddof=1)
    # Maximization step
    current_likelihood, k = maximize_likelihood(arr, e1, e2, d1, d2)
    change = current_likelihood - prev_curr_likelihood
    prev_curr_likelihood = current_likelihood
    # Convergence check
    if change < convergence_threshold:
      cnt += 1
    else:
      cnt = 0
    x1_p, x2_p = sorted(arr)[:k], sorted(arr)[k + 1:]
  if e1 > e2:
    e1, d1, e2, d2 = e2, d2, e1, d1
  return k, e1, d1, e2, d2

Подобрав оптимальное k для разбиений, найдем необходимые для ответа величины. Во-первых, мат.ожидания μ1, μ2.

In [ ]:
res_1 = expectation_maximization(x)
print(f'μ1 is {res_1[1]:.2f}, correct is {mu1}\nμ2 is {res_1[3]:.2f}, correct is {mu2}')
μ1 is -0.66, correct is 0
μ2 is 1.13, correct is 1

Рассмотрим, как выглядят разбиения, построенные с помощью expectation maximization и реальные данные. Заметнее всего то, что не получилось точно предсказать стандартное отклонение.

In [ ]:
# @title Рисование графиков
def draw_EM_results(arr, x1_orig, x2_orig, EM_results):

  k, e1, d1, e2, d2 = EM_results
  part1, part2 = sorted(arr)[:k], sorted(arr)[k + 1:]
  fig, ax = plt.subplots(1, 2, figsize=(12, 4))
  ax[0].set_title('Original distributions')
  ax[1].set_title('Predicted distributions')
  for a in [x1_orig, x2_orig]:
      sns.histplot(a, ax=ax[0], stat='density', bins=np.arange(min(arr), max(arr), 0.3), alpha=0.7, ls='None')

  x_axis = np.arange(min(arr), max(arr), 0.3)
  y_axis_1 = norm.pdf(x_axis, e1, d1)
  y_axis_2 = norm.pdf(x_axis, e2, d2)

  ax[1].plot(x_axis, y_axis_1)
  ax[1].plot(x_axis, y_axis_2)
  for a in [x1_orig, x2_orig]:
    sns.histplot(a, ax=ax[1], stat='density', color='teal', alpha=0.3, ls='None')
In [ ]:
draw_EM_results(x, x1, x2, res_1)
No description has been provided for this image

Определим для каждого числа вероятность принадлежать распределению 1 и 2, основываясь на построенной модели, и сравним с исходными распределениями с помощью метрики accuracy score.

In [ ]:
from sklearn.metrics import accuracy_score

def probabilities(arr, m1, s1, m2, s2, n1, n2):
  prior_p1, prior_p2 = n1 / (n1 + n2), n2 / (n1 + n2)
  res = []
  for el in arr:
    p1, p2 = norm.pdf(el, m1, s1), norm.pdf(el, m2, s2)
    p = prior_p1 * p1 + prior_p2 * p2
    res.append((prior_p1 * p1 / p, prior_p2 * p2 / p))
  return res

probs = probabilities(x, *res_1[1:], len(x1), len(x2))
predicted = [0 if el[0] > el[1] else 1 for el in probs]
round(accuracy_score(answers, predicted), 4)
Out[ ]:
0.686

Сами вероятности выглядят так, первый элемент - вероятность принадлежать первому распределению, второй - второму.

In [ ]:
probs[:2]
Out[ ]:
[(0.33309915402406204, 0.6669008459759379),
 (0.8768575471952835, 0.12314245280471656)]

Также определим количество чисел, пришедших из распределений 1 и 2, согласно предсказаниям модели.

In [ ]:
print(f'Принадлежат 1-му распределению: {len(predicted) - sum(predicted)}, второму: {sum(predicted)}')
Принадлежат 1-му распределению: 751, второму: 749

Проведя вышеперечисленные действия для датасетов размеров 10 и 5, 100 и 50, 1000 и 500, я пришла к следующему выводу: точность модели растет при повышении размеров датасетов. Она меньше подвержена влиянию выбросов при увеличении размеров выборки.


Также хочется отметить, что точность возрастает и при снижении пересечения двух распределений (см.ниже -- увеличим параметры исходного распределения (исходно это (1, 1)), а именно матожидание).

In [ ]:
test_x2, test_x3 = np.random.normal(mu2 + 1, sd2, n2), np.random.normal(mu2 + 2, sd2, n2)
test1 = np.concatenate([x1, test_x2])
test2 = np.concatenate([x1, test_x3])
res_2 = expectation_maximization(test1)
res_3 = expectation_maximization(test2)
print(f'μ1 is {res_2[1]:.2f}, correct is {mu1}\nμ2 is {res_2[3]:.2f}, correct is {mu2 + 1}')
print(f'μ1 is {res_3[1]:.2f}, correct is {mu1}\nμ2 is {res_3[3]:.2f}, correct is {mu2 + 2}')
μ1 is -0.36, correct is 0
μ2 is 1.92, correct is 2
μ1 is -0.29, correct is 0
μ2 is 2.61, correct is 3
In [ ]:
draw_EM_results(test1, x1, test_x2, res_2)
No description has been provided for this image
In [ ]:
draw_EM_results(test2, x1, test_x3, res_3)
No description has been provided for this image
In [ ]:
probs_1 = probabilities(test1, *res_2[1:], len(x1), len(x2))
predicted_1 = [0 if el[0] > el[1] else 1 for el in probs_1]
round(accuracy_score(answers, predicted_1), 4)
Out[ ]:
0.846
In [ ]:
probs_2 = probabilities(test2, *res_3[1:], len(x1), len(x2))
predicted_2 = [0 if el[0] > el[1] else 1 for el in probs_2]
round(accuracy_score(answers, predicted_2), 4)
Out[ ]:
0.9127

Действительно, точность предсказаний возросла.

Задание 2. Последовательности и их модели¶

2.1 Создаем последовательности¶

Сгенерируем 1000 последовательностей по 3000 bp над алфавитом {a, b}. Вероятности переходов выберем случайно (независимых из них только две).

In [ ]:
a = ['a', 'b']
prob_1, prob_2, start_prob = np.random.random(), np.random.rand(), np.random.rand()
starting_prob = [start_prob, 1 - start_prob]
transition_matrix = [[prob_1, 1 - prob_1],
                     [prob_2, 1-prob_2]]
pd.DataFrame(transition_matrix, columns=a, index=a)
Out[ ]:
a b
a 0.518791 0.481209
b 0.703019 0.296981
In [ ]:
def create_seq_markov(alphabet, n, weights, starting_weights):
  res = np.random.choice(alphabet, p=starting_weights)
  for i in range(1, n):
    prev_symbol = res[i - 1]
    res += np.random.choice(alphabet, p=weights[alphabet.index(prev_symbol)])
  return res


sequences = []
for i in range(1000):
  sequences.append(create_seq_markov(a, 3000, transition_matrix, starting_prob))

2.2 (и заодно 2.4.1) оценка моделей¶

Теперь оценим несколько моделей, которые можно использовать для генерации таких последовательностей: бернуллиевскую и марковские 1-3 порядков. Для этого используем критерии AIC (Akaike information criterion) и BIC (Bayesian information criterion).
Начнем с биномиального распределения. Количество (независимых) параметров для нее = 1. Значения критериев рассчитываются по следующим формулам:
image.png, image.png
где k - число параметров модели, L - оптимальное правдоподобие, n - размер (обучающей) выборки

In [ ]:
# @title Информационные критерии для модели Бернулли
from collections import Counter

def calculate_criteria_binomial(seqs, alphabet, n, l):
  '''Модель -- биномиальное распределение'''
  frequencies = Counter(''.join(seqs))
  frequencies = dict(zip(alphabet, [frequencies[k] / (n * l) for k, v in frequencies.items()]))
  print('Frequencies:', frequencies)
  log_L = 0
  for seq in seqs:
    log_prob_seq = 0
    for el in seq:
      log_prob_seq += np.log(frequencies[el])
    log_L += log_prob_seq
  AIC_bernoulli = 2 * 1 - 2 * log_L
  BIC_bernoulli = 1 * np.log(l) - 2 * log_L
  print(f'AIC: {AIC_bernoulli:.04f}, BIC: {BIC_bernoulli:.04f}\nlog(L) = {log_L:.04f}')

calculate_criteria_binomial(sequences, a, 3000, 1000)
Frequencies: {'a': 0.4063873333333333, 'b': 0.5936126666666667}
AIC: 4478762.8393, BIC: 4478767.7471
log(L) = -2239380.4197

Далее применим критерии для марковских моделей 1, 2 и 3 порядка. Для построения матрицы перехода посчитаем количество (n + 1)-меров, где n - порядок и частоты самих символов. Оценим соответствующие вероятности, используя формулу условной вероятности: image.png

In [ ]:
# @title Получение матрицы переходов для модели Маркова 1 порядка
from itertools import product
from collections import OrderedDict


def calculate_kmers(seqs, alphabet, k):
  res = {''.join(el): 0 for el in product(''.join(alphabet), repeat=k)}
  res = OrderedDict(res)
  cnt = 0
  for s in seqs:
    for i in range(0, len(s) - k, k):
      kmer = s[i: i + k]
      res[kmer] += 1
      cnt += 1
  return {k: v / cnt for k, v in res.items()}


def create_transition_matrix(kmers, alphabet, frequencies):
  res = pd.DataFrame(columns=alphabet)
  for k, v in kmers.items():
    kmer, letter = k[:-1], k[-1]
    if len(kmer) < 1 or frequencies[kmer] == 0:
      res.loc[kmer, letter] = v
    else:
      res.loc[kmer, letter] = v / frequencies[kmer]
  return res


def markov_chain(seqs, alphabet, order):
    matrices = []
    prev_freq = None

    # оценка начальных условий из предположения об эргодичности
    # и также матрицы переходов

    for k in range(order + 1):
      kmers_freq = calculate_kmers(seqs, alphabet, k + 1)
      cur_matrix = create_transition_matrix(kmers_freq, alphabet, prev_freq)
      matrices.append(cur_matrix)
      prev_freq = kmers_freq

    return matrices


def calculate_criteria_markov(seqs, matrices, alphabet):
  log_L = 0
  order = len(matrices) - 1
  for seq in seqs:
    log_L += np.log(float(matrices[0][seq[0]]))
    prev_el = seq[0]
    for k in range(1, order):
      kth_matrix = matrices[k]
      cur_el = seq[k]
      log_L += np.log(kth_matrix.loc[prev_el, cur_el])
      prev_el += cur_el
    t_matrix = matrices[-1]
    for i in range(order, len(seq)):
      cur_el, prev_el = seq[i], seq[i - order:i]
      log_L += np.log(t_matrix.loc[prev_el, cur_el])
  AIC = (len(alphabet) ** 2 - 1) * 2 - 2 * log_L
  BIC = (len(alphabet) ** 2 - 1) * len(seqs) - 2 * log_L
  return log_L, AIC, BIC


markov_matrices_1 = markov_chain(sequences, a, 1)
markov_matrices_1[-1]
Out[ ]:
a b
a 0.518835 0.480309
b 0.704073 0.297178

Теперь подсчитаем значения информационных критериев для этой модели. В результате получаем, что значения критериев для нее ниже, чем для биномиальной.

In [ ]:
criteria_markov_1 = calculate_criteria_markov(sequences, markov_matrices_1, a)
print(f'AIC: {criteria_markov_1[1]:.04f}, BIC: {criteria_markov_1[2]:.04f}\nlog(L) = {criteria_markov_1[0]:.04f}')
AIC: 3949634.6152, BIC: 3952628.6152
log(L) = -1974814.3076

Далее сделаем то же самое для моделей 2 и 3 порядков. Матрица переходов для моделей 2 порядка (матрицы с начальными распределениями не показаны).

In [ ]:
markov_matrices_2 = markov_chain(sequences, a, 2)
markov_matrices_2[-1]
Out[ ]:
a b
aa 0.518211 0.480115
ab 0.704337 0.295215
ba 0.521393 0.480098
bb 0.701788 0.300007

И значения информационных критериев:

In [ ]:
criteria_markov_2 = calculate_criteria_markov(sequences, markov_matrices_2, a)
print(f'AIC: {criteria_markov_2[1]:.04f}, BIC: {criteria_markov_2[2]:.04f}\nlog(L) = {criteria_markov_2[0]:.04f}')
AIC: 3949655.0948, BIC: 3952649.0948
log(L) = -1974824.5474

Как видим, значения критериев у модели 2 порядка возросли. И аналогично для модели 3 порядка, матрица переходов:

In [ ]:
markov_matrices_3 = markov_chain(sequences, a, 3)
markov_matrices_3[-1]
Out[ ]:
a b
aaa 0.518141 0.479585
aab 0.705886 0.295586
aba 0.519254 0.479809
abb 0.709289 0.300237
baa 0.519078 0.476121
bab 0.705236 0.296344
bba 0.522037 0.48012
bbb 0.697253 0.298549

И значения информационных критериев:

In [ ]:
criteria_markov_3 = calculate_criteria_markov(sequences, markov_matrices_3, a)
print(f'AIC: {criteria_markov_3[1]:.04f}, BIC: {criteria_markov_3[2]:.04f}\nlog(L) = {criteria_markov_3[0]:.04f}')
AIC: 3949625.2775, BIC: 3952619.2775
log(L) = -1974809.6387

Значение критериев наименьшее у модели 3 порядка. Какой вывод можно сделать из этого? Марковские модели на сгенерированных последовательностях показали себя лучше (значения информационных критериев меньше), чем бернуллиевская (что и ожидалось). При этом наименьшие значения критериев оказались у модели 3 порядка (3 < 1 < 2).

2.3 Применение различных моделей на реальных данных (геноме B.subtilis)¶

Для этого задания я взяла референсный геном B.subtilis с GenBank ID GCA_000009045.1. Его размер 4.2 Mb.

In [ ]:
!wget 'https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/009/045/GCA_000009045.1_ASM904v1/GCA_000009045.1_ASM904v1_genomic.fna.gz' -O genome.fna.gz
--2024-03-10 20:46:51--  https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/009/045/GCA_000009045.1_ASM904v1/GCA_000009045.1_ASM904v1_genomic.fna.gz
Resolving ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)... 130.14.250.11, 130.14.250.12, 2607:f220:41e:250::7, ...
Connecting to ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|130.14.250.11|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1248394 (1.2M) [application/x-gzip]
Saving to: ‘genome.fna.gz’

genome.fna.gz       100%[===================>]   1.19M  --.-KB/s    in 0.09s   

2024-03-10 20:46:52 (13.4 MB/s) - ‘genome.fna.gz’ saved [1248394/1248394]

In [ ]:
!gunzip genome.fna.gz
In [ ]:
with open('genome.fna') as inp:
  # в первом строке название последовательности, его пропускаем
  inp.readline()
  genome = ''.join([line.strip() for line in inp])

Поделим всю последовательность на куски по 1000 нуклеотидов. Полный из них 4215, поэтому для простоты отбросим последний фрагмент (длина в районе 600 bp).

In [ ]:
genome_split = [genome[i:i + 1000] for i in range(0, len(genome), 1000)]
genome_split = genome_split[:-1]
len(genome_split)
Out[ ]:
4215

Оценим параметры бернуллииевской модели и посчитаем информационные критерии для нее.

In [ ]:
nucleotides = ['A', 'C', 'G', 'T']
calculate_criteria_binomial(genome_split, nucleotides, 1000, len(genome_split))
Frequencies: {'A': 0.2818268090154211, 'C': 0.283020640569395, 'G': 0.21806714116251483, 'T': 0.21708540925266903}
AIC: 11760609.4301, BIC: 11760615.7765
log(L) = -5880303.7151

Далее используем Марковские модели 1-5 порядков. Для удобства значения критериев организованы в таблицу. Матрицы переходов не показаны (если требуется, они записаны в переменную genome_matrices).

In [ ]:
model_results = pd.DataFrame(columns=[1, 2, 3, 4, 5], index=['LogL', 'AIC', 'BIC'])
genome_matrices = []
indices = ['LogL', 'AIC', 'BIC']
for m_order in range(1, 6):
  markov_matrices = markov_chain(genome_split, nucleotides, m_order)
  genome_matrices.append(markov_matrices[-1])
  criteria_markov = calculate_criteria_markov(genome_split, markov_matrices, nucleotides)
  for i in range(3):
    model_results.loc[indices[i], m_order] = criteria_markov[i]
In [ ]:
pd.options.display.float_format = '{:.4f}'.format
model_results['binomial'] = np.array([-5880303.7151, 11760609.4301, 11760615.7765])
model_results
Out[ ]:
1 2 3 4 5 binomial
LogL -5736721.2391 -5694604.5763 -5676147.3045 -5669468.4956 -5668169.7666 -5880303.7151
AIC 11473472.4781 11389239.1526 11352324.6090 11338966.9912 11336369.5331 11760609.4301
BIC 11536667.4781 11452434.1526 11415519.6090 11402161.9912 11399564.5331 11760615.7765

Визуализируем один из критериев, например, BIC. Заметим, что модели становятся лучше согласно обоим информационным критериям в порядке, указанном на графике (т.е. марковская модель 5 порядка -- самая лучшая). Интересно отметить, что после 3 порядка значения критериев изменяются не так сильно при переходе на более высокий порядок. В связи с этим можно было бы проверить модель на переобучение, создав тренировочную, тестовую и валидационную выборки.

In [ ]:
fig = sns.scatterplot(x=[str(x) for x in list(model_results.columns)], y=model_results.loc[['BIC']].transpose()['BIC'].tolist())
fig.set_title('BIC')
Out[ ]:
Text(0.5, 1.0, 'BIC')
No description has been provided for this image

2.4.2 Проверка эргодичности реальных данных¶

Выделим подпоследовательности из нарезки, начиная с индексов, равномерно распределенных на интервале [1, 3000]. Исключим последовательности, у которых начальный индекс не попал в промежуток [1, 1000] (т.е. пустые строки).

In [ ]:
substring_indices = np.random.randint(low=1, high=3000, size=len(genome_split))
substrings = [genome_split[i][substring_indices[i]:] for i in range(len(substring_indices))]
substrings = [x for x in substrings if len(x) > 0 ]
In [ ]:
fig1 = sns.histplot([len(x) for x in substrings])
fig1.set_title('Substring length distribution')
Out[ ]:
Text(0.5, 1.0, 'Substring length distribution')
No description has been provided for this image

На этом все, остальное не успела :)