Поиск по геному с помощью Wolfram Language (Mathematica) и HadoopLink

http://blog.wolfram.com/2014/02/21/searching-genomes-with-mathematica-and-hadooplink/
  • Перевод

Перевод поста Пола-Жана Летурно (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:



И устанавливаем ссылку на главном узле Hadoop:



Чтобы проиллюстрировать идею, возьмём небольшой по размерам человеческий митохондриальный геном из GenomeData:



Сперва разделим геном на отдельные основания (A, Т, С, G):



Они будут нашими парами ключ-значение (k1, v1). Значения — начальные позиции каждого основания генома:

{k1, v1} = {основание, позиция}

То есть мы только что создали индекс для генома.

Экспортируем этот индекс в Hadoop Distributed File System (HDFS):



Для запросов мы будем использовать последовательность из 11 оснований, о которых мы знаем, что они содержатся в геноме:



Мы используем последовательность с большим количеством повторений в ней для того, чтобы усложнить задачу нашим алгоритмам.

Наш геномный поисковик на этот запрос должен вернуть позицию 515:



Теперь нам нужны mapper и reducer.

Как было сказано в первой частиmapper принимает пару ключ-значение (шаг 1), и выдаёт другую пару (шаг 2):



Mapper получает основание из индекса генома в качестве входного значения и выдаёт пары ключ-значение для каждого расположения в запросе, где встречаются индексы оснований (скоро вы узнаете почему):

(1) Ввод: {k1, v1} = {индекс основания, позиция в геноме}

GenomeSearchMapper

Ключ на выходе — положение в геноме, значение на выходе — положение в запросе:

(2) Вывод: {k2, v2} = {позиция в геноме, позиция в запросе}

Какая разница между позицией в геноме и позицией в запросе? Положение в запросе есть положение основания в запросе, в то время как положение в геноме есть положение в целом геноме.

К примеру, скажем, mapper получает пару ключ-значение с основанием A на позиции 517:



Позиции в запросе для основания А в запрашиваемой последовательности GCACACACACA есть 3, 5, 7, 9, и 11:



Вот эта последовательность с подсвеченными позициями:



Mapper имеет только одну пару ключ-значение с одной базой индексов вместе с запрашиваемой последовательностью. У него не имеется остальной части генома для сравнения, так что ему требуется найти все возможные способы выстраивания запроса с основанием A на позиции 517:



Здесь цвета соответствуют каждой A в запросе (по горизонтали) с их получающимися положениями в геноме (по вертикали). Возьмем, к примеру A в третьем основании в запросе (показано зеленым). Если его встроить с A на позиции 517, то запрашиваемая последовательность будет начинаться в геноме с позиции 515 (517 — 3 + 1 = 515) (также показано зеленым).

Точно так же, выделенное красным основание (пятое положение в запросе) заставит запрашиваемую последовательность выстраиваться, начиная с позиции 513 в геноме (также показано в красным). И то же самое для седьмой позиции в запросе с 511 позицией в геноме, (фиолетовый) девятой позиции в запросе с 509 позицией в геноме (оранжевый), одиннадцатой позиции в запросе с 507 позицией в геноме (коричневый).

Лишь одна из этих расстановок является правильной. В данном случае — для третьей позиции в запросе (зеленый цвет), которая позволяет запрашиваемой последовательности встроиться в геном. Но mapper этого не знает, и потому рассматривает все возможные совпадения.

Так так reducer собирает ключи, он соберёт и все основания, которые соответствуют одним и тем же позициям в геноме:

Ввод: {k2, {v2 …}} = {позиция в геноме, {позиции в запросе…}}

Теперь для данного положения в геноме reducer будет искать те случаи, когда значения формируют полную последовательность из позиций запросов:



Если reducer находит полное соответствие, он выдаёт положение генома:

Вывод: {k3, v3} = {запрашиваемая последовательность, положение в геноме}

Давайте теперь запросим у нашего геномного поисковика последовательность GCACACACACA:



(См. часть 1, где приводилось описание функции HadoopMapReduceJob).

И импортируем соответствия в геноме из HDFS:



Соответствующее положение в геноме — 515, и это правильный ответ! Наша геномная поисковая система работает!

Давайте теперь осуществим поиск с запросом, которому должны соответствовать два разных положения в геноме:



Этому запросу должны соответствовать позиции 10 и 2277:



Да, он нашел оба совпадения!

Теперь давайте отмасштабируем всё это для всего человеческого генома. Первым шагом является создание индексации, на этот раз для всего генома, а не только митохондриального. Чтобы это сделать, я скачал весь человеческий геном в виде текстовых файлов с государственного сервера и импортировал их в HDFS:





Получается по одному текстовому файлу на хромосому, содержащие необработанные последовательности хромосом:



Затем я применил 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, чтобы он мог работать с новым представлением ключа:



Reducer остался без изменений.

Давайте запустим поиск для той же последовательности:



На этот раз мы получаем совпадения для всего генома:



И мы можем ещё улучшить алгоритм. Как насчёт поиска неточных соответствий вместо точных? Слегка поменяем reducer, указав то, насколько большая часть запроса должна совпадать:



Это не самый эффективный способ поиска в геноме, но он показывает, как легко писать прототипы алгоритмов на основе MapReduce и запускать их в Mathematica. Если вам интересны подробности, рекомендую посмотреть мои недавние сообщения. А на интересующие вопросы вам с радостью ответят в HadoopLink GitHub Repo!
  • +18
  • 6,5k
  • 3
Wolfram Research 45,89
Wolfram Language, Mathematica, Wolfram Alpha и др.
Поделиться публикацией
Комментарии 3
  • 0
    Очень интересно, спасибо за статью! А на домашнем ПК можно развернуть инфраструктуру MapReduce, и в обозримый срок осуществить поиск по геному, или нужен кластер?
    • 0
      На практике люди пользуются либо общедоступными облачными сервисами вроде BLAST и Genome Browser (когда вещей для поиска мало, но надо внимательно посмотреть на результат), либо готовыми высокооптимизированными программами (когда искать надо очень много, но и над результатами не трясутся), которые можно запускать где угодно и почти как угодно (хоть в консоли, хоть из личного вебгуя, распараллеливание тривиально разделением входящих данных).

      Как и все остальные статьи про вольфрам, здесь представлено крайне специфическое, заточенное строго под вольфрам, решение в лоб для весьма далёкой от реальности задачи. Это скорее пособие для студентов-биоинформатиков для понимания работы базовых алгоритмов. Приведённый тут алгоритм, например, неоптимален по памяти и сложности, для поиска чего-то длиннее 10 нуклеотидов в чём-то длиннее митохондриального генома и правда потребуется разворачивать MapReduce.
    • 0
      Норамальные люди юзают Матлаб + Bioinformatics toolbox. А всё, что в статье написано — для любителей экзотики.

      Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

      Самое читаемое