Скорость

При проектировании архитектуры платформы 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  млрд. прочтений) выглядят так:

Время выполненияШаг пайплайнаVMCPURAM
2ч 41мУдаление адаптеров и прочтений низкого качества188 Гб
0ч 55мВыравнивание BWA-MEM12817921920 Гб
2ч 18мДедупликация и сортировка11632 Гб
0ч 36мРекалибровка BQSR4008001588 Гб
1ч 01мПоиск вариантов HaplotypeCaller4008002400 Гб
3ч 29мПоиск вариантов Strelka213264 Гб
1ч 14мПоиск ошибок свёрточной нейросетью4008002400 Гб
0ч 52мПрименение фильтров к VCF-файлу124 Гб

Единственный шаг, который сложно распараллелить на несколько виртуальных машин – это этап удаления адаптеров, время его выполнения прямо пропорционально размеру входной пары 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  млн. прочтений) выглядит так:

Время выполненияШаг пайплайнаVMCPURAM
0ч 11мУдаление адаптеров и прочтений низкого качества188 Гб
0ч 15мВыравнивание BWA-MEM8128120 Гб
0ч 14мДедупликация и сортировка11632 Гб
0ч 44мРекалибровка BQSR173468 Гб
0ч 47мПоиск вариантов HaplotypeCaller3939117 Гб
0ч 23мПоиск вариантов Strelka211616 Гб
0ч 49мПоиск ошибок свёрточной нейросетью3978156 Гб
0ч 12мПрименение фильтров к VCF-файлу114 Гб

Общее время до получения VCF-файла составило 3ч 16м (арифметическая сумма времени всех шагов больше, так как некоторые из них выполняются параллельно). Главным преимуществом облачного подхода к биоинформационному анализу является независимость обработки образцов друг от друга. То есть, общее время одновременного анализа 10 или 20 экзомов на платформе seq24 тоже будет составлять 3-4 часа.