Как понять, что две последовательности похожи?

Я однажды сидел на каком-то не очень интересном докладе, что-то про организацию какой-то базы, как вдруг после слов докладчика "Если две последовательности похожи друг на друга более, чем 95%" весь зал вдруг проснулся и сразу же прозвучал вопрос "Как вы сравниваете последовательности?". Было видно, что весь зал внезапно проснулся, и, замерев, уставился на бедного докладчика. Он смутился и ответил "Ну, стандартным способом". - "Каким это еще таким стандартным способом?" - сразу же спросил кто-то в зале. Докладчик задумался и после долгой паузы выдал "Через blast". В зале поднялся шум, все стали обсуждать это.

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

Это только в теории все просто. Вот у нас есть две последовательности
CGGACACACAAA 
и 
CGACAACACAAT

Мы просто считаем минимальное количество вставок-замен и удалений, суммируем их и выясняем, что это 3. А три - это и есть наше искомое расстояние между ними.

Но ДНК - это формат хранения данных. Данные триплетами транслируются в аминокислоты, из которых уже собираются белки, которые и являются теми микромеханизмами, из которых все построено.

Что означает "триплетами транслируются в аминокислоты"? Нуклеотидная последовательность разбивается на тройки, которые потом по таблице конвертируются в аминокислоты

CGG ACA CAC AAA

 |      |      |     |

R    T    H    K

Или так: CGG ACA CAC AAA --->  RTHK

Эти тройки называют кодонами. Все время забываю, сколько этих аминокислот, толи 20, толи 21. Существуют еще управляющие сигналы начала и конца трансляции кода и сама таблица трансляции слегка отличается у некоторых групп организмов.

А теперь представьте себе, что в коде произошла вставка, на третью позицию вклинился нуклеотид А.

CGGACACACAAA
CGAGACACACAAA

И оттранслируем его заново: CGAGACACACAAA ---> RDTQ

Все поломалось. Мы получили вообще другой набор аминокислот, которые, собравшись в белок, точно будут делать что-то совсем другое, чем от него ожидалось. Вставка и удаление всего одного нуклеотида это почти всегда означает полную поломку механизмов работы этой части генома, и, вполне возможно, даже смерть всего организма, если это произошло в критической и недублированной его части. Одна маленькая вставка или удаление ДНК вызывает огромную разницу в аминокислотной последовательности, если оно произошло в кодирующей части. То есть той, которая будет транслирована в аминокислоты. А если не в кодирующей, то, может, ничего и страшного. Буквой больше, буквой меньше, там все равно может быть забит давно не работающий ген, типа хвоста у человека или просто что-то для "форматирования" кода.

А вот замена нуклеотид вызывают куда меньше эффектов. Если посмотрите на таблицу, то видно, что меньше всего на выбор аминокислоты влияет третий нуклеотид. Он часто можно быть вообще любым. Можно даже подобрать совсем другую последовательность ДНК, которая будет транслирована в точно такой же набор аминокислот, то есть физически белок, который будет что-то делать в организме, будет точно такой же.

Например,

CGGACACACAAA --->  RTHK
AGAACGCATAAG --->  RTHK

одинаковые для организма, одинаковые для биолога и разные только для программиста. Поэтому биологи любят сравнивать оттранслированные последовательности. Сначала перевести ДНК в аминокислоты и сравнивать уже аминокислоты. Они называют это "сравнением в пространстве протеинов" - protein space. Что немного терминологически неверно и сбивает с толку. Белки - это уже следующий уровень, а тут мы пока смотрим только на аминокислоты.

Кстати, даже если замена нуклеотида приводит к замене аминокислоты, то реальный белок, состоящий из их цепочки, все равно может получиться вполне рабочим. В некоторых случаях замещение одной аминокислоты другой делает белок нерабочим, в некоторых - проходит вообще незамеченным. Еще в начале 90ых некоторые люди собрали статистику на тему какие аминокислоты могут замещать какие, а какие нет, и захардкодили это в некоторые утилиты сравнения аминокислотых последовательностей, но с тех пор, насколько мне известно, никто так глубоко не копал, чтобы перепроверить эти таблицы согласно современной статистике. Оно в общем-то как-то работает.

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

CGG ACA CAC AAA --->  RTHK
GGA CAC ACA AA --->  GHT

GAC ACA CAA A --->  DTQ

Их называют соответственно, первый, второй и третий фрейм. А ведь есть еще обратно комплементарная сторона. Итого 6 разных вариантов - 6 фреймов. Какой выбрать? Проще всего не думать и сравнивать все со всеми. Итого 36 сравнений вместо одного. Так часто и делают.

Но интересно, что правильный фрейм можно угадать. Точно так же, как и в естественном языке, частоты троек нуклеотид сильно разные. Какие-то тройки используются сильно чаще других.

К примеру, слово "правильный". Имеет тройки пра, рав, ави, вил, иль, льн, ьны, ный. А слово "щшкуоатвцаь" тройки щшк, шку, куо, оат, атв, твц, вца, цаь.
Взглянув только на "цаь" можно сразу сказать, что такого слова в русском языке нет. Но можно более формально подойти. Взять большой текст на русском, посчитать статистику троект по нему и получить, допустим, такие цифры частот для данных троек

(где 10%% это 0.010)

пра 10%%
рав 20%%
ави 10%%
вил 5%%
иль 2%%
льн 1%%
ьны 5%%
ный 50%%

щшк 0%%
шку 4%%
куо 1%%
оат 1%%
атв 3%%
твц 1%%
вца 3%%
цаь 0%%

и ститаем (где логарифм двоичный)

complexity ("правильный") =
log(1/0.010) + log(1/0.020) + log(1/0.010) + log(1/0.005) + log(1/0.002) + log(1/0.001) + log(1/0.005) + log(1/0.050) =
log(100) + log(50) + log(100) + log(200) + log(500) + log(1000) + log(200) + log(20) ~= 7 + 6 + 7 + 8 + 9 + 10 + 8 + 4 = 59

А теперь complexity("щшкуоатвцаь") =
log(1/0) + log(250) + log(1000) + log(1000) + log(333) + log(1000) + log(333) + log(1/0) ~= 8 + 10 + 10 + 8 + 10 + 8 + 2log(1/0) = 54 + две бесконечности.

Даже если вдруг щшк и цаь будут иметь частоту в одну промилле, то мы получим минимум 74.

Если есть вопросы откуда здесь вообще берутся частоты и логарифмы, то советую погуглить на тему кодирования Хаффмана и его более современное и развитое применение - intra кодирования в кодеках.

Можно пойти даже еще на уровень выше и заметить, что тройки аминокислот тоже распределены неравномерно. Скорее всего из трех вариантов RTH, GHT или DTQ встречается в реальных белках в разы чаще остальных двух. Я здесь когда-нибудь вставлю реально посчитанную таблицу частот, а пока поверьте мне на слово.

В итоге можно посчитать сложность (complexity) кодирования в битах как двоичный логарифм единицы на частоту и можно просто выбрать фрейм с минимальной complexity. Это реально работает. Как правило, один из фреймов выглядит гораздо менее сложным, то есть состоит из гораздо более часто встречающихся троек, чем остальные. Но это не везде так. Есть области генома, которые не транслируются в аминокислоты. Там такого может и не быть. Но, может, тогда и не надо такие области сравнивать в пространстве аминокислот?

Биологи поступают обычно немного по другому. Они знают, что трансляция может начаться только с определенных мест. Их называют ORF - Open Reading Frame. И заканчивается трансляция тоже в определенном месте. Если с окончанием трансляции более-менее понятно - там просто идет стоп-триплет (стоп-кодон), то начало трансляции не гарантируется старт-кодоном и существует некий набор эвристик, который оценивает, действительно ли трансляция начинается в этом месте и заканчивается в другом. Очень часто биологи из всего генома выбирают только эти небольшие регионы и исследуют только их, потому что с остальным просто непонятно что делать. Отсюда и идет поверье, что, мол, в ДНК много мусора, который ничего не значит.

Пока что вся история с превращением ДНК в некие реально работающие молекулы-запчасти, которые, плавая, сами собираются, например, в клеточные стенки, выглядит достаточно стройной - есть определенные места в ДНК, откуда начинается трансляция в аминокислоты, цепочки из которых сворачиваются в белки, из которых все и состоит. Но тут есть один очень неприятный нюанс, который настолько все усложняет, что его обычно просто игнорируют. Называется он сплайсинг. Для программиста это выглядит как будто трансляция началась с какого-то места, идет-идет, потом, бац и перепрыгнула куда-то в другое место, идет, потом опять перепрыгнула куда-то и продолжилась с него. В основном такое творится только у эукариот, у бактерий этого обычно нет. Причина этого в том, что трансляция идет не напрямую с ДНК, потому что она двойная и просто не пролазит в машину для трансляции, а с её одинарной временной копии - РНК, которая нестабильная из-за отсутствия комплементарной последовательности, спутывается сама с собой, в результате чего из нее выпадают некоторые части.
Часто вы не знаете, читаете ли кусок ДНК или РНК, из которой некоторые части уже отвалились. Иногда в описании к данным написано что там, хотя не стоит этому сильно доверять, а, бывает, что просто смесь - и куски ДНК и РНК. Обе вещи полезны. С одной стороны ДНК - ровно то, что записано. С другой стороны рабочие белки будут оттранислированы именно с РНК. У бактерий проблем с вываливающимися кусками из РНК обычно нет, и это одна из многих причин, по которой приходится писать специализированный софт для бактерий и эукариот и почему бактерии исследовать гораздо проще.

Ну хорошо, мы можем сравнивать последовательность в пространстве нуклеотид, можем - в пространстве аминокислот, но, особенно при сравнении вирусов приходится идти в пространство более высокого уровня. Я называю его пространством функциональных элементов. Просто каким бы вирус ни был, ему все равно приходится решать ряд проблем, как проникновение в клетку и построение оболочки. Как бы вирусы ни мутировали, в них все равно можно обнаружить все эти части и построить меру схожести даже с оооочень отдаленными вирусами, например, вирус группа связан с вирусами, заражающими растения. Сравнение в пространстве функциональных элементов это тема отдельного поста или даже постов.

Кто-то наверняка может сказать - можно использовать Deep Learning для идентификации функциональных частей в вирусах. Да, можно, никто не запрещает. Но точно так же можно забороть эту проблему обычным матаном, который будет во многом очень похож на Deep Learning, просто более явный его вид без разного рода скрытых и непонятных вещей и подгоночных коэффициентов.

Но вернемся к более простым вещам. Те, кто занимается data science скажут, что один из вариантов сравнения - это сравнение фич. Мы можем придумать и посчитать некоторые характеристики последовательностей и сравнивать их. Самая часто используемая характеристика последовательности называется GC content. Это просто процент букв C и G в последовательности. Экспериментально установлено, что это мера стабильности ДНК. Разные организмы имеют разные механизмы поддержания стабильности ДНК (да и другие причины наверняка есть) из-за чего только по одному этому числу часто можно сразу же предположить к какому примерно семейству организмов может относиться данная последовательность. Если две достаточно длинные последовательности имеют сильно разный GC Content, то организмы, к которым они принадлежат, скорее всего сильно разные. Это число даже обычно публикуют рядом с любой последовательностью, чтобы ученый мог на глаз оценить, что перед ним.

Можно все известные последовательности ДНК просто отсортировать по GC Content и получить очень неплохую группировку по степени похожести этих организмов в реальности. Когда мы пытаемся оценить на что похожа какая-то неизвестная последовательность, можно просто посчитать её GC Content и посмотреть на каком месте она окажется и в какую группу попадет. У меня где-то есть эта самая таблица и надо бы наверное её тут опубликовать как-нибудь будет.

Но GC Content это все технологии прошлого, есть расширение этой идеи - оценить распределение букв A C T и G или даже лучше всех троек этих букв и получить не одно число, а вектор и сравнивать вектора соотвествующие последовательностям. Data scientist'ы узнают тут операцию свертки на некоторые простые базовые элементы из троек нуклеотид. Вычислительно это легкая задача и работает сравнение и группировка последовательностей только на основании похожести распределениё троек просто замечательно хорошо даже для сильно мутирующих вирусов. Если вы получили какой-то новый неизвестный вирус, который ну просто никак не похож ни на что другое, самое простое, что можно сделать - посчитать статистику распределений троек нуклеотид и найти группу известных вирусов с похожими характеристиками.

Можно идти выше и собирать статистику не троек, а, например 16ек или 32ок. В анализе текста их называют n-граммы, а в бионформатике больше устоялся термин k-mer. Хотя 32mer очень красиво умещается в int64, но все равно уже возникает проблема нехватки памяти и скорости обработки, как и общей избыточности таких k-mer и как это красиво обходится minhash я напишу уже где-то в последующей статье.

Какой-то мегапост уже получился. Я закончу его на одной простой вещи. Если вы получили два вектора, у вас существует, конечно, много способов их сравнения. Можно посчитать сумму модулей разности, можно корень из суммы квадратов разностей, но самым удобным и точным по многим причинам способом является сравнение через косинус.

a * b = |a| * |b| * cos( угла между a и b)

Косинус угла между векторами как отношение скалярного произведения на произведение норм  (Евклидовых - корень из суммы квадратов) это и есть очень удобная мера схожести двух векторов. А задача сравнения двух последовательностей в итоге сводится к выборе пространства сравнения и вычисления координат последовательностей в этом пространстве.

Как геномные данные выглядят и откуда их берут

Ребенок смотрит на кошку и говорит "Как же повезло, что дырки на коже у кошки на том же месте, где глазки".

Нам в этом плане тоже очень повезло, что все живые организмы (и вирусы, которые ими формально не являются) имеют один и тот же вид генетического кода. Что говорит лишь о том, что когда-то у всех был общий предок, который потихоньку-помаленькой развился в вирус Эболы, плесень, птицу в небе, ну и лично в вас. У вас у всех одинаковый код. ACTG. Если сейчас немного забыть про вирусы, некоторые из которых немного по другому работают, то все остальные имеют молекулу ДНК, где друг за другом идет какая-то последовательность аденина, гуанина, цитозина и тимина, еще называемых нуклеотидами.

Это кусок вируса Эболы

CGGACACACAAAAAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGA

А это - человека

CATGTTTCCACTTACAGATCCTTCAAAAAGAGTGTTTCAAAACTGCTCTATGAAAAGGAAT

Молекулы состоящие из этих нуклеотид очень сложным и запутанным образом производят другие молекулы, часто очень сложной формы и обладающие множеством странных хических свойств. Эти молекулы уже собираются, например в оболочку клетки или же запускают какие-то процессы в организме.

Длина приведенной выше последовательности ДНК 61 bp. Сокращение bp - базовые пары. Парами их называют из-за того, что ДНК это двойная спираль. Параллельно этой последовательности идет такая же, но обратно комплементарная ей. 

CGGACACACAAAAAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGA

GCCTGTGTGTTTTTCTTTCTTCTTAAAAATCCTAGAAAACACACGCTTATTGATACTCCT

Под C идет G, под G - C, под A - T и наоборот. Но так как комплементарная последовательность однозначно восстанавливается из прямой, то ее обычно не пишут. Но у реальной ДНК нет ни верха, ни низа, и когда у вас есть некая необработанная сырая последовательность нуклеотид, вы обычно не знаете, прямая она или обратно комплементарная, и это нужно всегда держать в голове при написании кода. Оборудование, читающее ДНК, может прочитать или верхнюю строку справа налево или нижнюю задом наперед, но вы не можете это определить. По крайней мере оно не читает верхнюю задом наперед и нужно рассматривать 2 варианта расположения фрагмента, а не 4.

Но ДНК это "формат хранения" кода. Он идет такой двойной спиралью, потому что она куда более стабильна, чем одиночная последовательность. Но одиночные последовательности, без парной, тоже существуют. Некоторые вирусы даже существуют только в таком виде. Правда, они часто сворачиваются хотя бы кусками сами на себя. Когда последовательность вроде одинарная, но ее можно согнуть пополам, и ее одна половина будет комплементарна другой. Вирусы не смущаются постоянных поломок в коде, все равно они размножаются достаточно быстро. Такие одиночные последовательности называются РНК. Вот только РНК кодируется нуклеотидами ACUG вместо ACTG. Урацилом вместо тимина. Однако, чтобы не вносить путаницы, РНК тоже записывается как ACTG, чтобы не поддерживать два формата. 

Иногда встречаются какие-то странные символы в геномном коде типа R, которая означает A или G, но чаще всего код записывают именно верхним регистром и разных символов в нем, нет, к сожалению не 4, а 5. ACTG и N. Неизвестный символ. Однако N обычно обрабатывается особым образом и чаще всего хранить в памяти нужно именно 4 символа. Можно обойтись двумя битами на нуклеотид, что обычно и делают. А если вы посмотрите на битовое представление символов ACTG, то увидите интересную вещь:

A 65 1000001

C 67 1000011

T 84 1010100

G 71 1000111

Первый и второй биты образуют все различные комбинации. Можно быстро переводить нуклеотиды в двухбитный int подобными образом: (ch >> 1) & 3

Итак, откуда же их, эти последовательности ACTG берут? К сожалению нельзя просто так взять и прочитать последовательность ДНК какого-то организма. Вы, возможно, где-то уже слышали, что геном человека прочитан, причем достаточно давно. Ничего подобного. Существует уже 38ой "релиз" генома и целые комиссии сидят и постоянно его уточняют, выпускают "патчи" и полноценные "релизы". И даже в последнем релизе в нём все еще слишком много букв N, которые идут огромными регионами, не говоря уже про множество каких-то странных кусков генома кишечной палочки, которые вообще непонятно что там делают и скорее всего попали туда по ошибке.

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

К счастью, большая часть данных идет от машин компании Illumina. Счастье в том, что у них очень низкий уровень ошибок чтения. 1 ошибка примерно на 10 фрагментов. И каждый фрагмент длиной 100-200 нуклеотид. Очень часто они бывают парными. Когда аппаратура читает один фрагмент и на некотором расстоянии от него еще один. Это на самом деле очень помогает в некоторых исследованиях.

Так примерно выглядят парные фрагменты каждый длиной 100:

>SRR1919642.1
ATCTTCGTATAAAAACTAGACAGAATCATTCACAGAAACTAGTTTGTGATGTGTGTGTTCAACTCAAGGAGTTTAAACTTTCTTTTGATGGAGCAGTTTG
>SRR1919642.1
AAAGCGCTTGAAATCTCCAGCTGCAAATTCCACAAAAAGGGTGTTTAACATCTGCTCTTCTAAAGGAAAGTTCAACTCTATGAGTTGAATACACAAAGCA
>SRR1919642.2
NAAAATAAAAACACAAAGGTTCAATCTCTTCTGACCTTTGAAAGACACAGCACAGACAGTGGTCCTTAGGACGAAGAGCAGGAGACCCCTAATTCCGTCA
>SRR1919642.2
GGATTCAGTTCTTGAAACAAAACTCTGAGCCTTCAATGACCTTTCGGTCTATGTAAAAGCACTCCTGTCTTCCTGGCAGCAGTTGGACCTCACAATGTGG

Существует некоторая путаница в названии эти фрагментов. Исторически чтение происходило с некоторых точек на пластинах, поэтому саму пару называли спотами (spot - точка), а их половинки - ридами (read). Но в современных статьях спот называют ридом, а каждую из его половинок - фрагментом. Всегда есть небольшая путаница, что конкретный человек имеет под ридом - обе половинки или только одну. Мне лично привычней называть каждую половинку ридом, а в сумме они образуют спот. В приведенном выше примере 2 спота по 2 рида в каждом. Часто бывает, что есть всего 1 рид на спот, так что по сути это одно и то же.

Откуда я взял приведенный пример? Вот отсюда
https://www.ncbi.nlm.nih.gov/sra/?term=SRR1919642

Наряду с нуклеотидами в таких данных есть также данные об их "качестве" - уровне уверенности, что нуклеотид правильный. Интересно, что эти самые уровни занимают порядка 80% данных по обьему. Когда-то давно данных было очень мало и это никого не волновало. А сейчас правильность нуклеотид проверяют проводя многократное, десятки или сотни раз (бывает, что и тысячи. рекорд, который я видел - 1.2 млн) прочитывание всего генома кусочками, накладывая кусочков друг на друга и строя таким образом протяженные регионы - contiguous regions или contigs (контиги) на сленге. По историческим причинам данные уровня сигнала хранят для совместимости со старыми программами, но реально никого не волнует уровень сигнала при стократном перекрытии данных.

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

Прочитать целиком вирус - без проблем. Бактерия скорее всего прочтется большими такими кусками тысяч по 100 базовых пар. Не совсем будет понятно как они расположены и в одинственном ли они экземпляре, но этого более чем достаточно, для того, чтобы составить представление о полном геноме в несколько миллионов нуклеотид. Но если вы попытаетесь прочитать геном человека или вообще любого сложного организма (эукариота), то вы будете получать максимум непрерывные куски в несколько тысяч баз, вообще непонятно как распложенных относительно друг друга и безо всякой информации, весь ли вы геном можете прочитать или есть какие-то нечитаемые области.

Существуют митохондрии - бактерии, захваченные клетками предков человека много миллионов лет назад. Их геном прочитать легко. А вот получить непрерывный кусок человеческого ДНК хотя бы тысяч в 50 баз из 3 млрд имеющихся - мне пока что такой технологии неизвестно, и уже тем более, хотя у человека действительно есть хромосомы, представляющие непрерывную последовательность ДНК длинной в сотни миллионов, на самом деле людям известны только множество мелких непонятно как расположенных обрывков этих хромосом, которые могут для удобства храниться в одной последовательности и разделяться последовательностями NNNNNNNNNN.

Know your data!

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

Когда программист принимается за решение какой-то задачи, то ему обычно хочется начать с какой-то простой математической модели и к ней придумывать решение. Плохо то, что реальные бионформатические данные очень специфические и очень сильно отличаются от того, как это кажется "должно быть по идее". Я перелопатил множество таких данных и хочу поделиться своим личными наблюдениями о том, как они выглядят на самом деле, какие возникают задачи и как их можно решать.