При проектировании архитектуры платформы seq24 целевым параметром являлось минимальное время обработки одного образца без потери качества анализа. Для различных типов панелей было приняты следующие временные параметры обработки от FASTQ до VCF:
- экзом – менее 4 часов (размер сырых данных 5 Гб, ~90 млн. прочтений, среднее покрытие 75x),
- геном – менее 9 часов (размер сырых данных 50 Гб, ~900 млн. прочтений, среднее покрытие 30x),
- геном – менее 1 часа (размер сырых данных 50 Гб, ~900 млн. прочтений, среднее покрытие 30x) при использовании пайплайна DRAGEN-GATK.
Главной проблемой анализа больших файлов является отсутствие линейного прироста производительности при масштабировании биоинформационных алгоритмов на многопроцессорных системах. Так, BWA-MEM эффективно распараллеливается не более чем на 12-16 ядер CPU, инструмент HaplotypeCaller из набора GATK4 всегда работает в однопоточном режиме (несколько потоков использует только алгоритм PairHMM), шаги дедупликации, сортировки и рекалибровки BQSR также имеют существенные ограничения в максимальном использовании доступных вычислительных мощностей – это приводит к тому, что значительно снизить время обработки образца невозможно даже при наличии мощного локального сервера без отказа от части вычислений или замены алгоритмов на более быстрые, но менее точные.
Например, по данным Chen et al. (Chen J., Li X., Zhong H., Meng Y., Du H. 2019. Systematic comparison of germline variant calling pipelines cross multiple next-generation sequencers. Sci Rep. 9 (1), 1–13) оптимизация пайплайна для 24 CPU обеспечивает для генома общее время выполнения двух шагов GATK4-BQSR + GATK4-HaplotypeCaller в диапазоне 35-40 часов и оценку полного времени выполнения всего пайплайна в диапазоне 70-90 часов с учетом сбора метрик качества, причем дальнейшее увеличение количества CPU не приводит к пропорциональному уменьшению времени обработки.
Для решения этой проблемы команда GATK (столкнувшаяся с необходимостью анализа десятков и сотен геномов в день) разработала ряд инструментов и подходов, реализованных в последней версии The Genome Analysis Toolkit, значительная часть которых была внедрена на платформе seq24. Все шаги, допускающие независимую обработку определенного набора прочтений (BWA-MEM, BQSR, HaplotypeCaller и т.д.) автоматически распределяются по необходимому количеству виртуальных машин с оптимальными для каждой подзадачи ресурсами и настройками, а шаги, требующие анализа полного набора данных (дедупликация, сортировка, Strelka2), выполняются в кластерном окружении с подстройкой размера кластера под объем данных.
В результате проведенной оптимизации тайминг и выделяемые на анализ одного генома ресурсы (размер 84 Гб, ~1.5 млрд. прочтений) выглядят так:
Время выполнения | Шаг пайплайна | VM | CPU | RAM |
---|---|---|---|---|
2ч 41м | Удаление адаптеров и прочтений низкого качества | 1 | 8 | 8 Гб |
0ч 55м | Выравнивание BWA-MEM | 128 | 1792 | 1920 Гб |
2ч 18м | Дедупликация и сортировка | 1 | 16 | 32 Гб |
0ч 36м | Рекалибровка BQSR | 400 | 800 | 1588 Гб |
1ч 01м | Поиск вариантов HaplotypeCaller | 400 | 800 | 2400 Гб |
3ч 29м | Поиск вариантов Strelka2 | 1 | 32 | 64 Гб |
1ч 14м | Поиск ошибок свёрточной нейросетью | 400 | 800 | 2400 Гб |
0ч 52м | Применение фильтров к VCF-файлу | 1 | 2 | 4 Гб |
Единственный шаг, который сложно распараллелить на несколько виртуальных машин – это этап удаления адаптеров, время его выполнения прямо пропорционально размеру входной пары FASTQ-файлов. Именно по этой причине мы рекомендуем не склеивать разные дорожки/прогоны одного образца, а загружать в систему несколько отдельных пар файлов (если они имеются). В таблице для примера приведено время обработки пары файлов по 42 Гб каждый. В случае, если бы этот набор данных был загружен в виде 4-8 пар файлов с разных дорожек или прогонов, время обработки на первом шаге сократилось бы в несколько раз.
Дальнейшие вычисления проводятся над примерно равными порциями исходных данных одновременно (например, 128 виртуальных машин для BWA-MEM или 400 виртуальных машин для HaplotypeCaller), это позволяет сократить затрачиваемое время ресурсоемких шагов с 30-40 до 1-2 часов, а общее время обработки генома уменьшить до 9-10 часов, при этом строго выполняются все рекомендации GATK Best Practices и дополнительно проводится поиск вариантов с помощью Strelka2, значительно увеличивающий характеристики точности выполняемого анализа. Следует отметить, что все алгоритмы применяются к полному набору данных, если это требуется для обеспечения максимальной точности результата, даже ценой потерь в скорости. Ряд сервисов оптимизирует некоторые ресурсоемкие шаги - например, публичные пайплайны платформы SevenBridges используют только фрагмент хромосомы 20 для рекалибровки BQSR, экстраполируя результат на весь геном, это значительно ускоряет анализ, но приводит к небольшому падению чувствительности и специфичности.
Эффект проведенной оптимизации вычислений представлен на графике:
Параллельно с обработкой без влияния на ее производительность производится сбор метрик качества с помощью инструментов CollectHsMetrics, CollectRawWgsMetrics, CollectAlignmentSummaryMetrics и др. в зависимости от типа анализа.
Аналогичная табица для экзома (размер 4.7 Гб, ~90 млн. прочтений) выглядит так:
Время выполнения | Шаг пайплайна | VM | CPU | RAM |
---|---|---|---|---|
0ч 11м | Удаление адаптеров и прочтений низкого качества | 1 | 8 | 8 Гб |
0ч 15м | Выравнивание BWA-MEM | 8 | 128 | 120 Гб |
0ч 14м | Дедупликация и сортировка | 1 | 16 | 32 Гб |
0ч 44м | Рекалибровка BQSR | 17 | 34 | 68 Гб |
0ч 47м | Поиск вариантов HaplotypeCaller | 39 | 39 | 117 Гб |
0ч 23м | Поиск вариантов Strelka2 | 1 | 16 | 16 Гб |
0ч 49м | Поиск ошибок свёрточной нейросетью | 39 | 78 | 156 Гб |
0ч 12м | Применение фильтров к VCF-файлу | 1 | 1 | 4 Гб |
Общее время до получения VCF-файла составило 3ч 16м (арифметическая сумма времени всех шагов больше, так как некоторые из них выполняются параллельно). Главным преимуществом облачного подхода к биоинформационному анализу является независимость обработки образцов друг от друга. То есть, общее время одновременного анализа 10 или 20 экзомов на платформе seq24 тоже будет составлять 3-4 часа.