Из-за универсальности ДНК крайне важно определить, за что отвечают определенные части гена. Глядя на транскрипцию РНК-полимеразы, становится ясно, что в генах есть три основных участка:

  • Промоторные области: области выше или рядом с местом, где будет транскрибироваться ген, с которым связывается РНК-полимераза.
  • Кодирующие области: часть ДНК или РНК, которая кодирует белки.
  • Область терминации: области ниже по течению, которые сигнализируют обрыве РНК-полимеразы.

Промоторные области являются жизненно важной частью экспрессии генов, поскольку они контролируют связывание РНК-полимеразы с ДНК. В результате для приложений в области генной инженерии и анализа генов крайне важно определить, какие гены кодируют какие белки и как изменить эти белки.

Скрытые марковские модели

Начнем с базовой скрытой марковской модели. Предполагается, что эти скрытые марковские модели являются марковскими процессами со скрытыми состояниями. Определим следующие термины:

  • Наблюдения/выбросы. Это наблюдения, которые пользователь передает модели. Для этой конкретной модели это будут «Прогулка», «Магазин» и «Чистка».
  • Состояния. Это «скрытые» значения, которые мы пытаемся вывести из наблюдений с помощью вероятностного анализа. Для этой конкретной модели это будут «Дождливый» и «Солнечный».

Каждая стрелка представляет другую вероятность. Например, вероятность перехода от дождливого к солнечному — это стрелка, указывающая от дождливого к солнечному с вероятностью 0,6. Это означает, что существует 60-процентная вероятность того, что из состояния «Дождь» он перейдет в состояние «Солнечно» на следующем шаге. Это различные вероятности:

  • Начальные вероятности: это вероятности того, что состояния будут иметь место на первом этапе. Одним из таких примеров является вероятность 30% начать игру в дождливом состоянии. Другими словами, P(дождь) = 0,3.
  • Вероятности перехода. Это вероятности в различных состояниях. Одним из таких примеров является упомянутый выше переход от Дождливого к Солнечному. Предполагая, что t_k — это текущее состояние, мы можем сказать, что P(t_k+1 = «Солнечно» | t_k = «Дождь») = 0,6.
  • Вероятности выбросов: это вероятности того, что состояние вызовет это конкретное излучение. Например, P(Магазин | Дождь) = 0,4.

Преобразование в код

Bio.HMM.MarkovModel предоставляет полезный инструмент для инициализации этих скрытых марковских моделей для анализа.

from Bio.HMM.MarkovModel import HiddenMarkovModel
transition_alphabet = ['Rainy', 'Sunny']
emission_alphabet = ['W', 'S', 'C'] # W - walk, S - shop, C - clean
initial_prob = {'Rainy': 0.3, 'Sunny': 0.7}
transition_prob = {('Rainy', 'Rainy'): 0.4,
                 ('Rainy', 'Sunny'): 0.6,
                 ('Sunny', 'Rainy'): 0.2,
                 ('Sunny', 'Sunny'): 0.8}
emission_prob = {('Rainy', 'W'): 0.1,
                 ('Rainy', 'S'): 0.4,
                 ('Rainy', 'C'): 0.5,
                 ('Sunny', 'W'): 0.6,
                 ('Sunny', 'S'): 0.3,
                 ('Sunny', 'C'): 0.1}
model = HiddenMarkovModel(transition_alphabet, emission_alphabet, initial_prob, transition_prob, emission_prob, 1, 1)

Алгоритм Витерби

Viterbi позволяет нам выводить наиболее оптимальные скрытые состояния из данного наблюдения. Более подробный анализ Viterbi будет приведен ниже, но с помощью BioHMM мы можем использовать мощь Viterbi с помощью одной строки кода.

model.viterbi(Seq('WSC'), transition_alphabet)

Это обеспечивает вывод:

(Seq('SunnySunnyRainy'), -4.5972020163389145)

Если задана последовательность «Прогулка», «Магазин» и «Уборка», вам могут быть предоставлены скрытые состояния «Солнечно», «Солнечно» и «Дождь» соответственно для каждого из этих шагов. Мы видим, что вероятность этого равна 10^(-4,597). Это может показаться небольшим числом, пока вы не заметите, что происходит многократное умножение десятичных чисел, а это означает, что меньшее число весьма вероятно.

Идентификация промоутеров

Возьмем более сложный пример, где два состояния — фоновые или промоторные области, а вероятности эмиссии — различные нуклеотиды (A, T, C, G).

Обратите внимание, что здесь не заданы начальные вероятности, поэтому начальная вероятность запуска в некотором состоянии s равна P(s) = 1/n для n состояний.

from Bio.HMM.MarkovModel import HiddenMarkovModel
transition_alphabet = ['B', 'P'] # B - background, P - promoter
emission_alphabet = ['A', 'T', 'C', 'G']
initial_prob = {'B': 0.5, 'P': 0.5}
transition_prob = {('B', 'B'): 0.85,
                  ('B', 'P'): 0.15,
                  ('P', 'P'): 0.75,
                  ('P', 'B'): 0.25}
emission_prob = {('B', 'A'): 0.25,
                ('B', 'T'): 0.25,
                ('B', 'C'): 0.25,
                ('B', 'G'): 0.25,
                ('P', 'A'): 0.15,
                ('P', 'T'): 0.13,
                ('P', 'C'): 0.42,
                ('P', 'G'): 0.3}
model = HiddenMarkovModel(transition_alphabet, emission_alphabet, initial_prob, transition_prob, emission_prob, 1, 1)

Как и выше, мы можем использовать Viterbi:

model.viterbi(Seq('AGTACACTGGT'), transition_alphabet)
model.viterbi(Seq('GCGCGCGCAA'), transition_alphabet)

Они обеспечивают соответствующие результаты:

(Seq('BBBBBBBBBBB'), -17.567574447856494)
(Seq('PPPPPPPPBB'), -15.314217188702496)

Как мы видим в первом примере, последовательность «AGTACACTGGT» имеет вероятное состояние всех фоновых областей. Это имеет смысл, поскольку при рассмотрении вероятностей эмиссии маловероятно, что промоторные области взаимодействуют с А или Т с вероятностью эмиссии 0,15 и 0,13 соответственно.

Однако во втором примере из-за большого количества C и G последовательность «GCGCGCGCAA» имеет вероятное состояние, состоящее из всех промоторных областей. Это имеет смысл, поскольку островки CpG являются обычными местами для областей промотора гена.

HMM с «Память»

Наконец, мы можем повысить качество нашего анализа, рассматривая последовательности динуклеотидов. Что, если бы мы могли повысить точность нашей оценки состояния, учитывая, каким был последний нуклеотид? Это позволило бы делать более точные прогнозы, поскольку это вероятный сценарий в нашей ДНК.

Наблюдается следующий шаблон состояний с 8 состояниями: A+, A-, C+, C-, G+, G-, T+, T-. + представляет промоторные области, а - представляет фоновые области. Буква перед символом представляет собой последний нуклеотид.

Вероятности перехода были взяты из Durbin, Eddy, et al. 1998, а вероятности эмиссии тривиальны (если буква в штате совпадает с буквой в эмиссии, то вероятность равна 1). Начальные вероятности были определены как 1/n, так что все состояния имели случайные начальные вероятности 0,125.

Вероятности перехода были инициализированы следующим образом:

array_promoter = [0.18, 0.274, 0.426, 0.12, 0.171, 0.368, 0.274, 0.188, 0.161, 0.339, 0.375, 0.125, 0.079, 0.355, 0.384, 0.182]
vals_pp = [0.75*val for val in array_promoter]
vals_pb = [0.25*val for val in array_promoter]
array_background = [0.3, 0.205, 0.285, 0.21, 0.322, 0.298, 0.078, 0.302, 0.248, 0.246, 0.298, 0.208, 0.177, 0.239, 0.292, 0.292]
vals_bp = [0.15*val for val in array_background]
vals_bb = [0.85*val for val in array_background]

В первом массиве рассматриваются вероятности, когда начальным состоянием является промоутер. Это выглядит следующим образом: АА, АС, АГ, АТ, …. Переход от промоторной области к промоторной области имеет вероятность перехода 75%, а переход к фоновой области составляет 25%. Таким образом, эти числа были соответственно умножены для представления истинных оценок. То же самое было сделано для случая, когда начальным состоянием был фон.

После инициализации мы можем использовать Viterbi:

model.viterbi(Seq('AGTACACTGGT'), transition_alphabet)
model.viterbi(Seq('GCGCGCGCAA'), transition_alphabet)

Получаем следующие результаты:

(Seq('A-G-T-A-C-A-C-T-G-G-T-'), -17.77362274426406)
(Seq('G+C+G+C+G+C+G+C-A-A-'), -16.064944938458574)

Обратите внимание, что результаты остаются в основном такими же, с небольшой разницей во втором примере с 3 фоновыми нуклеотидами. Эти небольшие вариации показывают, что HMM с памятью имеют небольшое увеличение точности, но не так сильно из-за высокой точности мононуклеотидных HMM.

Алгоритм Витерби

Витерби был закодирован для определения его эффективности по сравнению с интерпретацией алгоритма Витерби в BioHMM.

Обратите внимание, что для понимания следующего раздела, скорее всего, потребуется предварительное знание алгоритма Витерби. За помощью обращайтесь сюда.

def viterbi_algorithm(observations, states, start_p, trans_p, emit_p, table):
    V = [{}]
    
    for st in states:
        V[0][st] = {"prob": start_p[st] * emit_p[st][observations[0]], "prev": None}

Сначала мы начнем с определения V, где все шаги нашего алгоритма Витерби будут храниться для печати в будущем. Затем мы перебираем все состояния и создаем словарь, инициализированный начальными вероятностями и предыдущими указателями (обратите внимание, что это первый шаг, поэтому предыдущих указателей нет).

for t in range(1, len(observations)):
        V.append({})
        for st in states:
            max_tr_prob = V[t - 1][states[0]]["prob"] * trans_p[states[0]][st]
            prev_st_selected = states[0] 
            for prev_st in states[1:]:
                tr_prob = V[t - 1][prev_st]["prob"] * trans_p[prev_st][st]
                if tr_prob > max_tr_prob:
                    max_tr_prob = tr_prob
                    prev_st_selected = prev_st
 
            max_prob = max_tr_prob * emit_p[st][observations[t]]
            V[t][st] = {"prob": max_prob, "prev": prev_st_selected}

Это создает указатели для программы Viterbi. Программа работает, определяя максимально возможную вероятность перехода путем умножения вероятности последнего шага на вероятность перехода.

Он обновляет это, перебирая каждое состояние, а затем умножает наилучшее переходное состояние на соответствующую вероятность испускания для окончательного сохранения в V. Он делает это для каждого отдельного состояния и для такого количества шагов, сколько определено длины последовательности.

opt = []
max_prob = 0.0
best_st = None
 
for st, data in V[-1].items():
    if data["prob"] > max_prob:
        max_prob = data["prob"]
        best_st = st
    opt.append(best_st)
    previous = best_st

Это определяет состояние максимальной вероятности для последнего шага запуска обратной трассировки.

for t in range(len(V) - 2, -1, -1):
        opt.insert(0, V[t + 1][previous]["prev"])
        previous = V[t + 1][previous]["prev"]
print (" ".join(opt) + ", probability = %s" % max_prob)

Затем он прослеживает матрицу указателей и смотрит на «prev», чтобы определить оптимальный путь. Наконец, он выводит окончательный результат компилятору с максимальной вероятностью.

Запустив эту модель на последовательности «GCGCGCGCAA», была создана следующая таблица DP со следующими результатами:

Сравнение результатов самокодированного алгоритма Витерби с реализацией HiddenMarkovModel дает следующее:

HiddenMarkovModel Implementation
G+C+G+C+G+C+G+C-A-A-, probability = -16.064944938458574

Own Implementation
G+C+G+C+G+C+G+C-A-A-, probability = 1.0545885728228732e-07

Как видно здесь, у них были одинаковые результирующие состояния, показывая, что модель действительно работала одинаково для обоих.

Просмотрите мой репозиторий Github для получения дополнительной информации и доступа к полному коду проекта!

TL;DR

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

Если хотите поговорить подробнее, запланируйте встречу: Календарно! Если вы заинтересованы в подключении, подпишитесь на меня в Linkedin, Github и Medium.