Elixir в биоинформатике

В этой статье я расскажу о своей попытке использования библиотеки GenStage, а в частности модуля Flow, для реализации одного из алгоритмов биоинформатики. На протяжении последних двух лет я занимался разработкой комплексной системы хранения и поиска результатов метагеномного анализа (метагеномика) углеводородного сырья. Наверное, для многих это китайская грамота. Фактически такой анализ означает выявление всех типов микроорганизмов, обитающих, к примеру, в залежах нефти. Некоторые из этих микроорганизмов, преимущественно бактерии, способны разъедать стальные трубы и создавать множество других неблагоприятных эффектов.

Первый этап метагеномного анализа – секвенирование генома нескольких тысяч особей, найденных в углеводородном сырье. Так как в результате секвенирования «считывается» не весь геном, а только отдельные его участки, то определить, какой особи принадлежит тот или иной участок – вычислительно сложная задача.

Что касается компьютерной реализации, анализ заключается в передаче исходных данных, так называемых азотистых оснований (A, T, U, C, и G), в цепочку процессов. Одной из известных программ для решения этой задачи является Qiime  (читается «чайм»). Она состоит из множества связанных друг с другом приложений, написанных на Python. Сначала я разработал свой фреймворк потоковой обработки данных для Elixir, но, как только появился GenStage, мне стало интересно оценить его возможности в проведении исследований подобного рода.

Переход на GenStage

Вместо того, чтобы пытаться переписать Qiime на Elixir (задание не из самых простых) я решил начать с малого, а именно выяснить, каким образом можно реализовать на Elixir простейшие алгоритмы биоинформатики и как пользоваться преимуществами параллельного исполнения GenStage (возможно, он сработает быстрее). С этой целью я прошёл замечательный онлайн-курс на сайте Coursera под названием «Находим скрытые в ДНК сообщения (Биоинформатика I)», в котором описываются задачи биоинформатики и пути их решения. Для их реализации можно выбрать любой язык программирования.

Одна из таких задач – поиск повторяющихся участков генома. Определим условие задачи. k-мер – это последовательность из k нуклеотидов, например, CTGGAT – это 6-мер. Дана последовательность, которая может состоять из миллионов нуклеотидов. Требуется найти k-меры, появляющиеся в ней неоднократно. Вместо прочёсывания всей последовательности, обратим внимание на k-меры, содержащиеся в отдельных её частях. В частности, в последовательности из 1500 нуклеотидов, k-меры будем искать только на участке из 500 нуклеотидов.

Пример из онлайн-курса

Дано: CGGACTCGACAGATGTGAAGAACGACAATGTGAAGACTCGACACGACAGAGTGAAGAGAAGAGGAAACATTGTAA
Задание: Найти 5-мер, присутствующий хотя бы 4 раза на участке из 50 нуклеотидов.
Решение: CGACA GAAGA

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

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

Реализация

А вот и моя реализация приведённого алгоритма с использованием модуля Flow из GenStage:

defmodule Bio do
alias Experimental.Flow

def clump(seq, k_mer_len, subseq_len, times) do
  |> seq
  |> String.to_charlist
  |> Stream.chunk(subseq_len, 1)
  |> Flow.from_enumerable
  |> Flow.partition
  |> Flow.map(fn e -> Stream.chunk(e, k_mer_len, 1) end)
  |> Flow.map(
        fn e ->
          Enum.reduce(e, %{},
            fn w, acc ->
              Map.update(acc, w, 1, & &1 + 1)
            end)
        end)
  |> Flow.flat_map(
        fn e ->
          Enum.reject(e, fn({_, n}) -> n < times end)
        end)
  |> Flow.map(fn({seq, _}) -> seq end)
  |> Enum.uniq
end
end

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

1. Последовательности нуклеотидов записываются в поток данных Flow.
2. Далее туда же помещаются k-меры, подсчитывается количество их повторений, и это значение записывается в поле Map.
3. Удаляются те элементы (k-меры), которые появляются меньше необходимого количества раз.
4. И, напоследок, вызывается функция Enum.uniq, которая отсеивает повторяющиеся элементы (важно не то, сколько раз появилась последовательность, а то, что она встретилась определённое количество раз). 

Обратите внимание на то, что я использую функцию Stream.chunk/4. Эта функция реализована в модулях Enum и Stream, но в Flow её нет. Будучи в замешательстве по поводу того, нужна ли отдельная реализация функции chunk для модуля Flow, я отправил этот вопрос в список рассылки Elixir. Создатель языка, Жозе Валим, любезно на него ответил, и более того, предоставил улучшенную реализацию функции clump (см. ниже)!

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

Реализация функции от Жозе:

def clump(seq, k_mer_len, subseq_len, times) do
  seq
  |> String.to_charlist
  |> Stream.chunk(subseq_len, 1)
  |> Flow.from_enumerable
  |> Flow.flat_map(&find_sequences(&1, k_mer_len, times))
  |> Enum.uniq
end
def find_sequences(subseq, k_mer_len, times) do
  subseq
  |> Stream.chunk(k_mer_len, 1)
  |> Enum.reduce(%{}, fn w, acc ->
       Map.update(acc, w, 1, & &1 + 1)
     end)
  |> Enum.reject(fn({_, n}) -> n < times end)
  |> Enum.map(fn({seq, _}) -> seq end)
end

Заключение от переводчика

Вы можете найти оригинал статьи по ссылке: GenStage and Bioinformatics
К сожалению, в сети ещё слишком мало информации об Эликсире на русском, поэтому, если вам понравилась статья, поддержите её плюсами и репостами, а если вы хотели бы периодически читать об Эликсире что-то новое, то присоединяйтесь к Вунш – русскоязычному сообществу Эликсира – и подписывайтесь на рассылку, чтобы получать переводы самых интересных статей!

Автор: tresstensel

Источник