Как посчитать биологические данные и не уронить сервер и ноутбук?

Наверняка вы слышали о биоинформатике. Звучит перспективно, приятно и полезно. Часто, ввиду всеобщих рассказов о перспективности и возможностях направления, некоторые люди из IT или из «мокрой» биологии (так называют область биологии, где работают в лаборатории с бактериями и прочими возможными объектами живой и не очень природы и реагентами) хотят перейти в биоинформатику. Однако далеко не все понимают, что же это за область такая и почему с ней сложно работать.

Почему е̶л̶к̶и̶ серверы падают и как решить такие проблемы? Разберем сегодня. 

Какой кажется работа биоинформатика

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

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

То есть тут важно усвоить, что секвенатор выдает только куски последовательностей цельного ДНК (методика такая, о ней несколько позже). Кстати, почитать о механизме подробнее можно здесь. Задача биоинформатика же — собрать все воедино. Или найти какую-то часть во всей горе, или объединить некоторые кусочки для дальнейшего анализа (и проводить этот самый дальнейший анализ). 

Кажется сложным? Тогда смотрите. Чтобы найти один ген (специфический) в геноме риса Oryza Sativa могут потребоваться недели, если не месяцы. Причем длина гена может составлять как несколько сотен нуклеотидов, так и несколько тысяч. Почему так? По различным базам данных разбросаны только куски отсеквенированного генома риса. При этом базы данных не интегрированы между собой. И вот ищи да собирай пазл воедино. А если нужно найти ген в неизвестном еще человечеству виде, то придется с нуля разбирать последовательность и пытаться найти нужный ген.

Зачем искать отдельный ген? Просто потому что он может нести ценную информацию в самых разных смыслах. Это и устойчивость к антибиотикам, и устойчивость к различным заболеваниям. Мы уже отметили, что анализ данных может быть очень долгим. И тем не менее его востребованность со временем только растет. С каждым годом разрабатываются новые модели секвенаторов, более производительные. С каждым годом востребованность и популярность молекулярно-генетических манипуляций становится выше. Необходимость обработки и интерпретации полученных данных растет быстрее, чем производительность современных вычислительных комплексов. 

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

Секвенируют геном. И получают десятки тысяч генов (человеческий геном всего 20-25 тысяч генов). Сравнивают гены здоровых и больных. И пытаются найти закономерности. Вроде алгоритм попарного сравнения не так сложен и страшен. Однако гены могут работать в группах, образовывать кластеры, регулировать другие группы или один ген и так далее. Да даже мутации могут появиться не ввиду заболевания. А просто другой стрессовый фактор воздействовал. Необходимо не только провести множественные сравнения (с поправками, конечно), но и выявить, где отличия — не шум (об этой проблеме ниже). Реально удается разобрать в лучшем случае несколько процентов всей информации. На большее не хватает ни сил людей, ни киловатт энергии (гранты мелкие, электричество дорогое, извините), ни серверов.  

Сам же сигнал бывает очень слабый и его поймать очень трудно. Помните эпоху аналогового ТВ? Иногда поймать канал и настроить его было непросто. Биологические сигналы в этом плане очень похожи. Они очень слабы, их легко испортить неосторожным движением. Многие сигналы с различного диагностического оборудования получаются на грани их пороговой чувствительности. И в этом отношении очень легко получить ложные результаты и запороть исследование.

Часто это загадка для многих биологов, несмотря на то, что они могут быть отличными специалистами. Стоит добавить еще и то, что часто измерить какой-то биологический параметр иногда бывает просто невозможно, из-за его связи с другими. Т.е. нужно производить анализ сразу буквально всего, а современные аналитические методы в своем большинстве позволяют в один момент времени проанализировать только что-то одно. Что делать с остальным? Опять же, сложные математические модели, проведение симуляций, тесты в лабе, и т. п. Сложно все.

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

Так как оптимизировать вычисления на сервере (или на рабочем ноутбуке/компьютере) для биоинформатика? Существует несколько подходов.  Мы можем разделить задачи на несколько пакетов и поставить пакетную обработку по очереди. Дальше возможно множество не очень понятных терминов. Поэтому поясняем.

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

ㅤㅤㅤㅤㅤУсловная схема гена

Еще нужно отметить, что биоинформатики — настоящие форматные маньяки (точнее, становиться такими их заставляет жизнь). Вместо одного формата файла — тысячи. Вместо написания одного алгоритма — тысячи различных инструментов для обработки данных на каждом шагу. Это все необходимо ввиду сложности представления данных. 

В итоге получается следующее. Секвенатор выдает файлы формата fasta, архивированные при помощи gz. Затем биоинформатик д̶е̶л̶а̶е̶т̶ ̶в̶и̶д̶,̶ ̶ч̶т̶о̶ ̶ч̶т̶о̶-̶т̶о̶ ̶д̶е̶л̶а̶е̶т̶ ̶в̶ ̶т̶е̶р̶м̶и̶н̶а̶л̶е̶ ̶и̶ ̶м̶е̶н̶я̶е̶т̶ ̶ф̶о̶р̶м̶а̶т̶ ̶ф̶а̶й̶л̶о̶в̶ делает выравнивание. Это такое сравнение с эталоном (ведь геном человека уже давно отсеквенирован. И сделан «усредненный» вариант здорового человека). В итоге появляется уже новый файл в формате bam/sam.

После часто нужно выискать определенные мутации. Причем мутации, которые очень условно могут иметь значение. То есть это не случайное неправильное прочтение единичного рида. После такого сканирования появляется еще один файл — vcf.

ㅤㅤㅤㅤㅤКонец вставки. 

Например, у вас есть 10 образов экзома (отсеквенированы только экзоны) людей со злокачественной опухолью. Пусть даны исходно .fasta файлы. Фаста формат достаточно тяжеловесный и неструктурированный. Там все последовательности еще похожи на тот самый «взрыв журналов». Наша базовая задача — выровнять геном. То есть собрать страницы журнала и строчки в нем. Человеческий геном давным-давно собран. Потому выравниваем на референс. После выравнивания у нас есть относительно структурированный файл, с которым можно в дальнейшем делать различные манипуляции: обрезать нужную часть, индексировать, выискивать мутации и другое. Запускать каждый из 10 образцов по очереди на выравнивание, потом на индексирование, обрезку можно. Но зачем? Гораздо лучше поставить пакетную обработку и выпить чаю с коллегами. Или поспать.

Если доступны несколько процессоров, то запустить обработку одновременно на всех процессорах, или же распределить обработку между несколькими системами (распределенные вычисления). Допустим, у нас много образцов (уже не 10 пациентов, а нормальная выборка из сотен и хорошо, если это отдельные файлы). Их нужно обработать определенным образом, например, как было сказано выше. Можно ли каждый образец обрабатывать независимо? Определенно. 

Что если работа уже не разделяется на независимые задачи, а следующий шаг зависит от результата предыдущего?

Берем все те же 10 образцов генома людей. После выравнивания каждый файл нужно проиндексировать. Это значит добавить «координаты» на последовательности для быстрого поиска. Как будто у вас есть ленточный накопитель, на который записали «треки», которые нужно разметить перед последующей обработкой. Потом найти значимые точковые мутации (когда меняется только один нуклеотид (SNP, снип)) в образцах генома. И сравнить эти образцы между собой. 

Конечно, можно запустить несколько циклов пакетной обработки. Но в норме в лаборатории нет времени, и результаты должны быть вчера.

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

ㅤㅤㅤㅤㅤВычислительные стратегии

Очевидно, что далеко не все инструменты для распределенных вычислений достаточно гибкие, чтобы адаптироваться к индивидуальным требованиям без значительных временных и трудозатрат (ну и деняк). А часть таких инструментов в принципе низкоуровневые и совсем нечеловекочитаемые (OpenMP). Поэтому хотелось бы предоставить некоторые готовые гибкие решения, которые вполне можно адаптировать для биологических данных и снизить проблемы функциональной совместимости и трудозатраты исследователей (в том числе затраты на понимание работы с инструментом).  Это была долгая прелюдия, чтобы подвести вас к самим вариантам оптимизации. Давайте разберем их.

С появлением MapReduce от Google жить стало гораздо легче, жить стало веселее (с). MapReduce простыми словами — разделение общей задачи на элементарные, которые после выполнения на распределенных компьютерах (узлах кластера) сводятся в общий результат. Это и хорошая масштабируемость и великолепная производительность.

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

Hadoop находит свое применение в сборке геномов de novo. Сборка De novo — это альтернативный обычной сборке генома на референс подход. Иногда бывают неизвестные виды бактерий, кто-то принес нового жука с экспедиции и т. д. Это уже не просто сравнили эталон и подставили наиболее совпадающие риды (хотя это тоже не всегда просто). Это вычислительно затратный процесс и алгоритмически более сложный, нежели выравнивание на референс. 

Алгоритм следующий. При подготовке образца на секвенирование его фрагментируют. И это случайный процесс. В силу этого и при глубоком покрытии (это мера того, сколько раз кусок был секвенирован, т.е. сколько раз определенный рид был получен с фрагмента ДНК) обязательно найдутся перекрывающиеся риды. И таким образом можно совмещать последовательно перекрывающиеся риды и получать полную последовательность.

Последовательность, которая получилась из перекрывающихся ридов, называется контигом. Контиги затем используется для поиска таких ридов, которые максимально его покроют (т. е. это должен быть относительно длинный рид). Ну и в итоге объединяются контиги — получается скаффолд. Долго и упорно, пока не найдутся длинные повторы. Это, как правило, кончики хромосомки. Конечно, могут быть дыры в результате сборки. Но их потом «перепрочитывают». И становится жить хорошо.

MapReduce в контексте сборки де ново очень удобен. Его используют для построения графа де Брейна (та самая сборка без референса), которые очень удобны для работы с короткими ридами. Например, Contrail, инструмент для сборки последовательностей de novo на основе Hadoop, способен создавать списки смежности для всех k-мер в считываемых геномных последовательностях, а затем использовать «reduce», для сжатия перекрытий до единой последовательности.

k-меры в данном случае короткие геномные последовательности длиной в k нуклеотидов. Например, k=3. В таком случае для коротких последовательностей 
ACGT CGTA GTAC TACG возможные k-меры будут такими:

Списки смежности хранят k-мер и его ближайших соседей (другие k-меры, которые могут быть связаны с конкретным). На деле получается словарь с парами ключ-значение. Reduce «перерабатывает» данные, представляя на выходе в идеальном случае полную последовательность генома. 

Дифференциальную экспрессию (с использованием RNA-Seq) можно измерить с помощью программного конвейера Myrna (который запускается на hadoop) — конвейеры представляют собой потоки данных, состоящие из последовательных этапов, на которых к данным применяется программное обеспечение для биоинформатики. (Тут вопрос, поскольку дифференциальную экспрессию использовать можно много где).

У Hadoop есть, конечно, нюансы. Из-за использования в нем дисковой системы, промежуточные результаты не кешируются. А это может быть важно… Поэтому при использовании Хадупа возможна только пакетная обработка, при итеративной производительность будет падать. В таком случае гораздо оптимальнее использовать Apache Spark.

Наверное, самым неплохим примером будет сокращение обработки множества образцов секвенирования с 1000 часов до пары. Почитать об этом можно здесь.

В отличие от дисковых вычислений Hadoop, Spark выполняет вычисления с использованием памяти за счет внедрения абстракции resilient distributed dataset (RDD). Промежуточные результаты можно сохранять в памяти, и это более эффективно для итеративных операций. Spark может работать в 100 раз быстрее с точки зрения доступа к памяти, чем Hadoop. С точки зрения гибкости, Spark предоставляет SDK для Java, Scala, Python и R (последние два являются любимчиками биоинформатиков), а также интерактивную оболочку и оптимизированный движок, поддерживающий общие графики выполнения.

В 2015 году ученые использовали Spark, чтобы развернуть набор GATK ( набор инструментов для генетического анализа, вроде картирования (определения положения генетической информации в геноме), поиска мутаций и др.) на основе кластеров. Чтобы сократить время выполнения, этот подход сохранял данные активными в памяти между фазами map и reduce (тут можно также передать привет Hadoop). На основе статистики времени выполнения активной рабочей нагрузки был разработан алгоритм динамической балансировки нагрузки. Этот алгоритм позволяет лучше использовать производительность системы. Результаты экспериментов показали, что этот метод позволяет увеличить скорость в 4,5 раза по сравнению с многопоточным конвейером GATK на одном узле. Например, делаем стандартный анализ GATK (это получение снипов, SNP, или полиморфизмов — точковых изменений в организме). 

Обычно происходит следующее. Сначала загружаем необходимые данные. Часто, к слову, на один образец приходится два файла — forward и reverse. У нас же ДНК двухцепочечная. Далее файлы индексируются с помощью bwa для подготовки к выравниванию. На выходе угадайте что? Новые форматы, конечно!

После — выравнивание и получение sam файла. Затем — сжатие sam в bam (перевод в двоичный формат). Подсчет количества ридов. Он нужен, чтобы предварительно понять, насколько качественное секвенирование получено.

Для дальнейших манипуляций с данными необходим сортинг. Это разбивка по хромосомам и индексам на них. После сортинга ищут дубликаты. Их можно соотнести с количеством ридов, полученных изначально, понять долю и еще больше убедиться в качестве секвенирования. Затем найти снипы (SNP), те самые точковые изменения нуклеотидов, которые отличают конкретный организм. 

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

Развернуть распределенные вычисления на компьютере или на неплохом сервере, конечно, здорово. Но не всегда это возможно в силу разных причин. Да и не всегда есть доступ к нужным компьютерам и серверам (звучит странно, но бывает). 

Она хороша по многим параметрам. Это и user-friendly интерфейс (это не запуск инструментов из терминала), и возможность добавлять инструменты в пайплайн. Это и набор самих инструментов. Биоинформатических. Он очень широк и подходит для работы на разных этапах и с самыми разными форматами файлов. Ведь Galaxy заточен под биоинформатику, под биологов без бэкграунда в программировании. Вариантов анализа — масса: выравнивание, поиск мутаций, филогенетический анализ. Да и в принципе любая обработка, о которой было рассказано выше.

Ознакомиться со всем спектром инструментов легче по ссылке. Также можно интегрировать NextFlow в Galaxy.

Nextflow — это языково-независимая платформа для разработки и выполнения распределенных конвейеров данных. Она достаточно популярна в биоинформатике ( и естественно-научной области в целом).

Поддерживаемые языки сценариев включают bash, csh, ksh, Python, Ruby и R. Любой язык сценариев, использующий стандартное объявление shebang Unix (#!/bin/bash), совместим с Nextflow. Это тоже удобство в освоении. Nextflow является гибким и масштабируемым. Однако необходимо будет хотя бы немного разбираться в DSL. 

Есть инструменты распределенных вычислений, подходящие под биоинформатические задачи. Однако далеко не всегда они подходят под конкретные исследования. Если позволяют навыки, лучше всего настроить Spark и на нем разворачивать самые разные тулы и конвейеры. Если не позволяют, или же если нет времени на обучение — на помощь приходят решения внутри Galaxy. Возможно, это будет менее производительно, но гораздо лучше ручной обработки. 

Конечно, биоинформатика — это всего лишь одна из научных дисциплин, получающих выгоду от распределенных вычислений. Список областей применения в исследованиях огромен. Распределенные вычисления — движущая сила развития криминалистики, агробизнеса, образования, страхования и анализа рисков, веб-архитектуры. Фактически, любая область, использующая анализ данных, сможет делать это намного лучше, намного быстрее и намного доступнее благодаря распределенным вычислениям. Может, именно внедрение распределенных вычислений и распараллеливания операций (в каждой лаборатории) приведет к открытию новых видов терапии неизлечимых на данный момент болезней (и это не только о людях). А потому в комментариях предлагаю дать советы по распределенным вычислениям для биологов и всех причастных!

Автор: Анастасия Лазукина.