![](https://habrastorage.org/files/932/cfb/b58/932cfbb58df74fa1ae2ecad6e2f313d1.png)
Перевод поста Пола-Жана Летурно (Paul-Jean Letourneau) "Searching Genomes with Mathematica and HadoopLink".
Код, приведенный в статье, можно скачать здесь.
Примечание: этот пост написан как продолжение поста Большие массивы данных в Mathematica с HadoopLink.
Примечание переводчика: автор данной статьи под термином геном понимает всю совокупность генов некоторого структурного элемента живой материи. Это несколько отличается от стандартных определений, близких по смыслу, в которых подразумевается либо вся совокупность генов конкретного вида (Ridley, M. (2006). Genome. New York, NY: Harper Perennial), либо полный набор генетических инструкций, которые можно найти в клетке (http://www.genome.gov/Glossary/index.cfm?id=90). В данном посте будем пользоваться представлением автора.
В моём предыдущем посте я описал, как писать алгоритмы MapReduce (вики) в Mathematica с помощью пакета HadoopLink. Теперь давайте копнём немного глубже и напишем более серьёзный алгоритм MapReduce.
Я уже писал раньше о некоторых занятных возможностях в сфере геномики в Wolfram|Alpha. Если вам это интересно, вы даже можете осуществлять поиск по человеческому геному определённых последовательностей ДНК. Биологам часто требуется найти расположение фрагмента ДНК, которые они нашли в лаборатории, для определения того, какому животному принадлежит этого фрагмент, или из какой он хромосомы. Давайте используем HadoopLink для создания геномной поисковой системы!
Как и прежде, мы загружаем пакет HadoopLink:
![](https://habrastorage.org/getpro/habr/post_images/9ce/335/9aa/9ce3359aa4c33eed76b2c017f679731d.png)
И устанавливаем ссылку на главном узле Hadoop:
![](https://habrastorage.org/getpro/habr/post_images/31f/416/65a/31f41665a74f40ee9c73a2b931a46fea.png)
Чтобы проиллюстрировать идею, возьмём небольшой по размерам человеческий митохондриальный геном из GenomeData:
![](https://habrastorage.org/getpro/habr/post_images/aca/5ba/8c3/aca5ba8c3dac04ba08a9f4c33ea64280.png)
Сперва разделим геном на отдельные основания (A, Т, С, G):
![](https://habrastorage.org/getpro/habr/post_images/ea6/f0a/a30/ea6f0aa30d25932c720ff3de99374b83.png)
Они будут нашими парами ключ-значение (k1, v1). Значения — начальные позиции каждого основания генома:
{k1, v1} = {основание, позиция}
То есть мы только что создали индекс для генома.
Экспортируем этот индекс в Hadoop Distributed File System (HDFS):
![](https://habrastorage.org/getpro/habr/post_images/9da/244/f7f/9da244f7fe9c3d9febd2b742b0bb8212.png)
Для запросов мы будем использовать последовательность из 11 оснований, о которых мы знаем, что они содержатся в геноме:
![](https://habrastorage.org/getpro/habr/post_images/449/ba7/032/449ba7032269c0773df8e660af3dcd99.png)
Мы используем последовательность с большим количеством повторений в ней для того, чтобы усложнить задачу нашим алгоритмам.
Наш геномный поисковик на этот запрос должен вернуть позицию 515:
![](https://habrastorage.org/getpro/habr/post_images/afe/80c/821/afe80c82196f1d248702e8d1ad28e4ab.png)
Теперь нам нужны mapper и reducer.
Как было сказано в первой части — mapper принимает пару ключ-значение (шаг 1), и выдаёт другую пару (шаг 2):
![](https://habrastorage.org/getpro/habr/post_images/ef7/392/506/ef7392506b9e71d630c1e92e93b42553.png)
Mapper получает основание из индекса генома в качестве входного значения и выдаёт пары ключ-значение для каждого расположения в запросе, где встречаются индексы оснований (скоро вы узнаете почему):
(1) Ввод: {k1, v1} = {индекс основания, позиция в геноме}
![GenomeSearchMapper GenomeSearchMapper](https://habrastorage.org/getpro/habr/post_images/5ab/d34/e4c/5abd34e4ce4591834c67c0abd9cab3dd.png)
Ключ на выходе — положение в геноме, значение на выходе — положение в запросе:
(2) Вывод: {k2, v2} = {позиция в геноме, позиция в запросе}
Какая разница между позицией в геноме и позицией в запросе? Положение в запросе есть положение основания в запросе, в то время как положение в геноме есть положение в целом геноме.
К примеру, скажем, mapper получает пару ключ-значение с основанием A на позиции 517:
![](https://habrastorage.org/getpro/habr/post_images/fed/e6c/92f/fede6c92fbf4da140a95d5d73294599e.png)
Позиции в запросе для основания А в запрашиваемой последовательности GCACACACACA есть 3, 5, 7, 9, и 11:
![](https://habrastorage.org/getpro/habr/post_images/03f/2c8/692/03f2c8692b5ba9690b81bd52f4138ed3.png)
Вот эта последовательность с подсвеченными позициями:
![](https://habrastorage.org/getpro/habr/post_images/306/a54/8da/306a548da6e60dd79af7cd4d7a95379f.png)
Mapper имеет только одну пару ключ-значение с одной базой индексов вместе с запрашиваемой последовательностью. У него не имеется остальной части генома для сравнения, так что ему требуется найти все возможные способы выстраивания запроса с основанием A на позиции 517:
![](https://habrastorage.org/getpro/habr/post_images/30e/e74/58b/30ee7458b835e2266bdcf881efee9713.png)
Здесь цвета соответствуют каждой A в запросе (по горизонтали) с их получающимися положениями в геноме (по вертикали). Возьмем, к примеру A в третьем основании в запросе (показано зеленым). Если его встроить с A на позиции 517, то запрашиваемая последовательность будет начинаться в геноме с позиции 515 (517 — 3 + 1 = 515) (также показано зеленым).
Точно так же, выделенное красным основание (пятое положение в запросе) заставит запрашиваемую последовательность выстраиваться, начиная с позиции 513 в геноме (также показано в красным). И то же самое для седьмой позиции в запросе с 511 позицией в геноме, (фиолетовый) девятой позиции в запросе с 509 позицией в геноме (оранжевый), одиннадцатой позиции в запросе с 507 позицией в геноме (коричневый).
Лишь одна из этих расстановок является правильной. В данном случае — для третьей позиции в запросе (зеленый цвет), которая позволяет запрашиваемой последовательности встроиться в геном. Но mapper этого не знает, и потому рассматривает все возможные совпадения.
Так так reducer собирает ключи, он соберёт и все основания, которые соответствуют одним и тем же позициям в геноме:
Ввод: {k2, {v2 …}} = {позиция в геноме, {позиции в запросе…}}
Теперь для данного положения в геноме reducer будет искать те случаи, когда значения формируют полную последовательность из позиций запросов:
![](https://habrastorage.org/getpro/habr/post_images/ccc/db1/e00/cccdb1e005a0c48870bf1a498d82d0dd.png)
Если reducer находит полное соответствие, он выдаёт положение генома:
Вывод: {k3, v3} = {запрашиваемая последовательность, положение в геноме}
Давайте теперь запросим у нашего геномного поисковика последовательность GCACACACACA:
![](https://habrastorage.org/getpro/habr/post_images/c08/20e/508/c0820e5084b8a0cfeff26d96949ebfaa.png)
(См. часть 1, где приводилось описание функции HadoopMapReduceJob).
И импортируем соответствия в геноме из HDFS:
![](https://habrastorage.org/getpro/habr/post_images/cf5/4a9/add/cf54a9adde6e76fa74840a401b3b51f0.png)
Соответствующее положение в геноме — 515, и это правильный ответ! Наша геномная поисковая система работает!
Давайте теперь осуществим поиск с запросом, которому должны соответствовать два разных положения в геноме:
![](https://habrastorage.org/getpro/habr/post_images/c17/442/d71/c17442d71a3e8c0da5b388a479c93ffd.png)
Этому запросу должны соответствовать позиции 10 и 2277:
![](https://habrastorage.org/getpro/habr/post_images/6dd/1c9/58b/6dd1c958b52162193e5cfe2846443915.png)
Да, он нашел оба совпадения!
Теперь давайте отмасштабируем всё это для всего человеческого генома. Первым шагом является создание индексации, на этот раз для всего генома, а не только митохондриального. Чтобы это сделать, я скачал весь человеческий геном в виде текстовых файлов с государственного сервера и импортировал их в HDFS:
![](https://habrastorage.org/getpro/habr/post_images/140/686/103/140686103c33b56c5fda7c9464de5bb1.png)
![](https://habrastorage.org/getpro/habr/post_images/89b/73b/cb6/89b73bcb627913fd34da47405ec0b2b8.png)
![](https://habrastorage.org/getpro/habr/post_images/8bb/d3b/43e/8bbd3b43e6cdf7ba565ed3c269c9b67b.png)
Получается по одному текстовому файлу на хромосому, содержащие необработанные последовательности хромосом:
![](https://habrastorage.org/getpro/habr/post_images/fa2/b4d/b2e/fa2b4db2ef32f4df9e9d069654e2102d.png)
Затем я применил MapReduce для создания пар ключ-значение в индексе на HDFS, который выглядит так:
[hs_ref_GRCh37.p13_alts.fa, 121] G
[hs_ref_GRCh37.p13_alts.fa, 122] A
[hs_ref_GRCh37.p13_alts.fa, 123] A
[hs_ref_GRCh37.p13_alts.fa, 124] T
[hs_ref_GRCh37.p13_alts.fa, 125] T
[hs_ref_GRCh37.p13_alts.fa, 126] C
[hs_ref_GRCh37.p13_alts.fa, 127] A
[hs_ref_GRCh37.p13_alts.fa, 128] G
[hs_ref_GRCh37.p13_alts.fa, 129] C
[hs_ref_GRCh37.p13_alts.fa, 130] T
Одно небольшое отличие от предыдущего примера в том, что ключ сейчас представляет собой {хромосома, позиция в геноме}, а значение — основание на этой позиции. Так что теперь хромосома будет в ключе. Я немного изменил mapper, чтобы он мог работать с новым представлением ключа:
![](https://habrastorage.org/getpro/habr/post_images/ce2/622/8f8/ce26228f86ead80614157e468421b585.png)
Reducer остался без изменений.
Давайте запустим поиск для той же последовательности:
![](https://habrastorage.org/getpro/habr/post_images/aa3/170/bd3/aa3170bd366fdc41198c66531954e4da.png)
На этот раз мы получаем совпадения для всего генома:
![](https://habrastorage.org/getpro/habr/post_images/401/aea/3b2/401aea3b20626194ea65575dcee0c948.png)
И мы можем ещё улучшить алгоритм. Как насчёт поиска неточных соответствий вместо точных? Слегка поменяем reducer, указав то, насколько большая часть запроса должна совпадать:
![](https://habrastorage.org/getpro/habr/post_images/1c4/174/c2b/1c4174c2b8688aa330cbecbe16ff880e.png)
Это не самый эффективный способ поиска в геноме, но он показывает, как легко писать прототипы алгоритмов на основе MapReduce и запускать их в Mathematica. Если вам интересны подробности, рекомендую посмотреть мои недавние сообщения. А на интересующие вопросы вам с радостью ответят в HadoopLink GitHub Repo!
Комментарии (3)
kluwert
28.07.2015 11:25Норамальные люди юзают Матлаб + Bioinformatics toolbox. А всё, что в статье написано — для любителей экзотики.
Mistx
Очень интересно, спасибо за статью! А на домашнем ПК можно развернуть инфраструктуру MapReduce, и в обозримый срок осуществить поиск по геному, или нужен кластер?
ninch
На практике люди пользуются либо общедоступными облачными сервисами вроде BLAST и Genome Browser (когда вещей для поиска мало, но надо внимательно посмотреть на результат), либо готовыми высокооптимизированными программами (когда искать надо очень много, но и над результатами не трясутся), которые можно запускать где угодно и почти как угодно (хоть в консоли, хоть из личного вебгуя, распараллеливание тривиально разделением входящих данных).
Как и все остальные статьи про вольфрам, здесь представлено крайне специфическое, заточенное строго под вольфрам, решение в лоб для весьма далёкой от реальности задачи. Это скорее пособие для студентов-биоинформатиков для понимания работы базовых алгоритмов. Приведённый тут алгоритм, например, неоптимален по памяти и сложности, для поиска чего-то длиннее 10 нуклеотидов в чём-то длиннее митохондриального генома и правда потребуется разворачивать MapReduce.