Сборка генома из ридов на пангеномном графе с помощью открытого ПО

Предыдущие материалы цикла о полногеномном секвенировании: 1, 2, 3, 4.

Один из авторов статьи, перечисленных в конце, предложил исследовать геном Вадима Зорина — врача-биофизика, кандидата биологических наук, которому в 67 лет удаётся сохранять поразительную выносливость и добиваться впечатляющих спортивных результатов. Другие его регалии также приведены в конце статьи. Мы собирались рассчитать все доступные полигенные оценки рисков и особенностей и найти те, которые заметно отличают Вадима от среднего европейца: резко выраженные либо подтверждённые сразу несколькими исследованиями. Геном уже был отсеквенирован, а компания передала Вадиму список вариантов и исходные риды. Для анализа мы решили загрузить варианты в открытый проект Just-DNA-Lite, подробно описанный Сергеем в третьей и четвёртой частях.

Однако ничего не вышло. Геном собрали относительно устаревшего референса hg19, стандарта 2009 года, который Just-DNA-Lite не поддерживает. Поскольку исходные риды сохранились, оставалось заново пересобрать геном. И уж если браться за эту задачу, то основательно. В третьей части Сергей использовал для своего генома связку bwa-mem2 + DeepVariant и линейный референс. На этот раз мы выбрали граф-пангеном и наиболее мощные открытые инструменты из тех, что удалось найти.

Попутно мы пересобрали ещё два генома — Сергея Каменева и его жены, опробовали несколько способов вызова вариантов на одних и тех же выравниваниях и сопоставили результаты с опубликованными эталонными сравнениями и исходными ридами. Главный вывод: оптимальный рецепт зависит от секвенатора. Для Illumina лучшим оказался pangenome-aware DeepVariant, для MGI — обычный DeepVariant.

Над задачей работали не только люди. Конвейер, сравнения и проверки выполнял Opus 5.5, а выводы дополнительно проверяла Astra. Мы формулировали задачи и задавали неудобные вопросы.

Сложность материала возрастает по мере чтения. В первой половине доступно объясняются устройство генома и метрики качества сборки, во второй разбираются технические подробности: оборудование, команды, время выполнения и возникшие проблемы.

Как из ридов собирают геном

Геном человека содержит около трёх миллиардов «букв» — нуклеотидов A, C, G и T. Каждая хромосома представлена двумя копиями: материнской и отцовской. Одну такую копию генома называют гаплотипом.

Обычный секвенатор не способен прочитать хромосому целиком. ДНК разрезают на фрагменты, после чего прибор считывает отрезки длиной примерно 150 букв. Каждый такой отрезок называется ридом. Для полного генома нужны сотни миллионов ридов: у Вадима их было 431,5 млн пар. Парные риды считываются с двух концов одного фрагмента.

В результате получается головоломка из сотен миллионов деталей. Собирают её по аналогии с картинкой на коробке — референсу. Это эталонная последовательность, с которой сравнивают данные. Она не принадлежит какому-то «идеальному человеку»: референс составлен из ДНК нескольких доноров, главным образом одного. Его периодически обновляют. GRCh37, известный также как hg19, появился в 2009 году, GRCh38 — в 2013-м.

Выравнивание. Первый этап сборки. Для каждого рида программа-выравниватель, например bwa-mem2 или giraffe, ищет подходящее место в референсе, размещает его там и оценивает уверенность в выборе. Эта оценка называется MAPQ. Высокое значение указывает на однозначное положение. Ноль означает, что существует как минимум другое столь же подходящее место, откуда рид тоже мог происходить.

Вызов вариантов. На втором этапе программа-вызователь вариантов, или caller, анализирует распределённые по геному риды и определяет отличия от референса. Каждое такое отличие называют вариантом. SNV — замена одной буквы; индел — вставка или удаление нескольких букв. Например, при синдроме Жильбера в гене UGT1A1 встречается дополнительный повтор TA: индел C>CAT.

Генотип. Это сочетание значений в конкретной позиции двух копий хромосомы. В файле его обозначают цифрами: 0 — буква совпадает с референсом, 1 — первый альтернативный вариант, 2 — второй:

• 0/0 — обе копии совпадают с референсом;

• 0/1 — гетерозигота: вариант есть в одной копии, но отсутствует в другой;

• 1/1 — гомозигота по варианту: он присутствует в обеих копиях;

• 1/2 — в двух копиях находятся разные варианты.

Выходные файлы. VCF содержит список вариантов и генотипов; именно его загружают в сервисы анализа. gVCF дополнительно включает блоки, где последовательность подтверждённо совпадает с референсом. Это позволяет отличить отсутствие варианта от непрочитанного участка. CRAM хранит сжатое выравнивание всех ридов.

dbSNP и rsID. dbSNP — международный каталог известных вариантов. Каждому присвоен идентификатор rsID: например, rs7412 обозначает один из вариантов APOE. Такие номера добавляют к вариантам в сборке, чтобы их было удобно искать в публикациях и базах данных.

Недостатки исходной сборки

Лаборатория передала Вадиму VCF с 5,05 млн записей и два файла ридов, полученных на Illumina NovaSeq 6000, 2×151. У VCF обнаружились три проблемы.

Референс hg19. Он же GRCh37. В следующей версии, GRCh38, исправлены тысячи ошибок сборки и заполнена часть пробелов. Современные базы и инструменты ориентированы именно на неё. Одна и та же буква генома может иметь разные координаты в hg19 и GRCh38, поэтому файл на основе hg19 нельзя напрямую сопоставить с базой для GRCh38. Just-DNA-Lite, рассчитывающий полигенные оценки, не поддерживает GRCh37 и завис при загрузке генома.

Вызователь вариантов. Судя по полям файла, использовался freebayes — классический статистический вызователь. Современные нейросетевые инструменты ошибаются существенно реже; далее мы приведём конкретные цифры.

Фильтрация. В заголовке указано одно условие, QUAL > 5, но столбец FILTER пуст. По такому файлу сложно отделить надёжные варианты от сомнительных.

Координаты можно программно перевести с одного референса на другой. Этот процесс называется лифтовером. Однако для нашей задачи такой подход неудовлетворителен: технологии выравнивания развиваются, GRCh37 устарел, а лифтовер лишь перерабатывает прежний результат. Он теряет варианты в местах различий между сборками и не устраняет ошибки выравнивания и вызова. Пересборка по исходным ридам относительно GRCh38 позволяет решить все три проблемы одновременно.

Какие участки генома сложнее всего анализировать

Около 90% генома составляют относительно «лёгкие» участки с уникальной последовательностью. Риды размещаются там однозначно, а современные конвейеры дают почти одинаковые результаты. Различия сосредоточены в оставшихся 10%.

Один из главных источников сложностей — паралоги, почти одинаковые копии участка ДНК в разных местах генома. Они появляются вследствие дупликации — редкой мутации, при которой участок копируется и дополнительная копия остаётся в геноме. Обычно это происходит из-за сбоя при образовании половых клеток: парные хромосомы обмениваются участками, и при их смещённом совмещении одна хромосома получает участок дважды. Подобные события происходили у наших далёких предков, после чего копии наследовались потомками. Со временем в них накопились небольшие различия, но совпадение может по-прежнему составлять 99% и более. Длинные копии, от тысячи букв, называют сегментными дупликациями. Для коротких ридов они особенно проблематичны: рид из одной копии может выровняться на другую, и установить его истинное происхождение очень трудно.

Сложные области можно разделить на четыре группы:

• Сегментные дупликации и паралоги. Это гены-близнецы: трипсиногены, включая PRSS1 с ложным «панкреатитом» из третьей части, CYP2D6 и его псевдогены — неработающие копии гена, а также SMN1/SMN2. Риды соседней копии могут попасть не туда и создать целые кластеры ложных «мутаций».

• Центромеры и сателлиты. Центромера — перетяжка хромосомы, за которую её растягивают при делении клетки. Она состоит из сателлитной ДНК: коротких последовательностей, тысячи раз повторённых подряд. Основной человеческий повтор называется альфа-сателлитом и содержит около 171 буквы. Представьте книгу, в которой один абзац напечатан на сотнях страниц: по обрывку в 150 букв невозможно определить страницу. Долгое время истинную последовательность центромер не удавалось собрать, поэтому в GRCh38 на её месте находится модель. Здесь концентрируется множество ложных вызовов.

• Тандемные повторы и гомополимеры. Тандемный повтор — короткий мотив, повторяющийся подряд, например CACACACA…. Гомополимер — многократное повторение одной буквы: AAAAAAA…. В таких местах секвенатор может ошибаться в числе повторов, а один и тот же индел допустимо записать в разных позициях.

• MHC/HLA на 6-й хромосоме. Это наиболее изменчивый участок генома с генами тканевой совместимости. По ним иммунная система отличает собственные клетки от чужих; их проверяют при трансплантации органов. Последовательность каждого человека здесь значительно отличается от референса, поэтому выравниватель неохотно размещает риды с многочисленными отличиями. Возникает смещение к референсу: программа обнаруживает референсную последовательность там, где её на самом деле нет.

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

Как оценивают качество сборки

При наличии эталона: точность и полнота

Для объективной оценки конвейера нужен геном, правильный ответ для которого уже известен. Консорциум GIAB (Genome in a Bottle, NIST) на протяжении многих лет исследует несколько геномов, HG001–HG007, всеми доступными методами и публикует списки достоверных вариантов. T2T-Q100 для HG002 — геном, полностью собранный от теломеры до теломеры. Теломеры находятся на концах хромосом, укорачиваются при каждом делении клетки и служат своеобразным счётчиком оставшихся делений. Точность этого эталона составляет примерно одну ошибку на миллион букв. Он охватывает и сложные области, недоступные прежним эталонам.

Конвейер запускают на ридах эталонного генома и вычисляют три показателя. Точность показывает, какая доля обнаруженных вариантов действительно существует. Полнота — какую долю настоящих вариантов удалось найти. F1 — их среднее гармоническое. F1, равный 0,9979, соответствует приблизительно двум ошибкам на тысячу.

Именно так получены все приведённые ниже таблицы с числом ошибок: авторы методов обрабатывали эталонный геном и подсчитывали пропущенные и лишние вызовы.

Без эталона: транзиции и трансверсии

Для геномов Вадима и Сергея, как и для генома большинства людей, эталона нет. Создание эталона занимает годы: геном читают разными методами, включая длинные риды PacBio HiFi и Oxford Nanopore, а затем изучают вручную. Проверить конвейер можно на открытых ридах эталонного генома, но это не покажет, где именно он ошибся в вашем. Поэтому используют косвенные признаки. Первый из них — Ts/Tv, отношение числа транзиций к числу трансверсий среди однобуквенных замен. Чтобы понять, почему этот показатель помогает оценить шум, понадобится немного химии.

Буквы ДНК относятся к двум группам: пурины A и G имеют двухкольцевую структуру, а пиримидины C и T — однокольцевую. Транзиция — замена внутри одной группы: A↔G или C↔T. Трансверсия — замена пурина на пиримидин или наоборот: A↔C, A↔T, G↔C, G↔T. Если учитывать направление, возможны четыре транзиции и восемь трансверсий. Следовательно, для случайных замен, в том числе ошибок секвенатора и другого шума, Ts/Tv ≈ 4/8 = 0,5.

Среди настоящих мутаций транзиции встречаются значительно чаще. Основная причина связана с метилированием. У человека C в паре CG часто несёт химическую метку — метильную группу. Такие метки в начале гена обычно выключают его. В частности, благодаря этому в нейронах и клетках печени работают разные наборы генов, хотя их геном одинаков.

У метки есть побочный эффект. Молекула C содержит аминогруппу — атом азота с двумя атомами водорода. Иногда при взаимодействии с водой она самопроизвольно отщепляется. Обычный C при этом превращается в урацил, которого в ДНК быть не должно; клетка замечает повреждение и исправляет его. Метилированный C превращается в полноценный T, поэтому клетка нередко пропускает ошибку. Так транзиция C→T закрепляется в геноме.

Кроме того, ферменты восстановления ДНК легче путают буквы сходной формы. В результате для настоящих человеческих вариантов Ts/Tv составляет примерно 2,0–2,1 по всему геному и около 2,8–3,0 в белок-кодирующих областях.

Отсюда способ самопроверки: чем больше ложных вызовов, тем ближе Ts/Tv к 0,5. Если для 4 млн SNV, то есть однобуквенных замен, получено 1,90 вместо ожидаемых 2,0, грубая модель случайных ошибок указывает на несколько десятков тысяч ложных вариантов. Но модель приблизительная: систематические ошибки, например различия между копиями паралогов, имеют собственное значение Ts/Tv, поэтому напрямую перевести показатель в процент ошибок нельзя. Однажды мы уже ошиблись на этом основании; подробности — в спойлере о пороге MAPQ.

Без эталона: X-хромосома, dbSNP и риды

X и Y у мужчин. У мужчины одна X- и одна Y-хромосома, поэтому гетерозигота — вариант только в одной из двух копий — на X невозможна. Исключение составляют псевдоаутосомные участки PAR на концах X и Y, общие для обеих хромосом. Каждая гетерозигота на X вне PAR является ошибкой. Это бесплатный способ проверить качество.

Доля известных вариантов. Настоящий вариант, вероятно, уже встречался у других людей и внесён в dbSNP. Большое число отсутствующих в каталоге вариантов — повод насторожиться. В сложных областях этот признак ненадёжен; на примере генома Вадима мы покажем почему.

Подсчёт ридов. Если сборки расходятся в генотипе, можно открыть выравнивание и подсчитать риды с каждой буквой. Это наиболее трудоёмкий, но и наиболее убедительный способ проверки.

Пангеном: референс на основе множества людей

Линейный референс представляет собой геном почти одного человека. Если у вас есть вариант, отсутствующий в референсе, содержащий его рид выровняется хуже. А когда таких вариантов подряд много, как в MHC, рид может вообще не выровняться или попасть в паралог — почти идентичную копию участка в другом месте генома.

Консорциум HPRC (Human Pangenome Reference Consortium) собрал полные геномы десятков людей разного происхождения и объединил их в граф. Каждый узел графа соответствует фрагменту ДНК длиной от одной до тысяч букв. Ребро соединяет фрагменты, которые хотя бы у одного человека следуют друг за другом.

Граф похож на карту со всеми возможными дорогами. В одинаковых у всех людей участках проходит одна трасса. В местах различий она разделяется на несколько объездов, которые затем снова соединяются. Например, после общего фрагмента …GATT один путь проходит через A, другой — через G, а затем оба ведут к общему CCTA…. Такую развилку с последующим слиянием называют пузырём, или bubble: это общепринятый термин в работах о пангеноме. Геном конкретного человека представляет собой один маршрут по карте от начала хромосомы до конца. Хромосома — линейная молекула с двумя равноправными концами, но по соглашению координаты отсчитывают от конца короткого плеча.

Как по обычной дороге можно двигаться в обе стороны, так и узел графа можно проходить в двух направлениях. ДНК состоит из двух цепей, каждая из которых имеет химическое направление. Клеточные ферменты и секвенатор читают цепь только в одну сторону. Цепи направлены навстречу друг другу, поэтому рид с одной цепи идёт по хромосоме слева направо, а с другой — справа налево; при этом буквы заменяются комплементарными: A на T, C на G. Кроме того, у некоторых людей отдельные участки целиком перевёрнуты. Поэтому узлы можно проходить прямо и обратно, а граф называют двунаправленным.

Для выравнивателя это означает следующее: рид с вариантом может найти путь, где тот уже присутствует, и быть размещён с высокой уверенностью. Мы использовали граф HPRC v1.1, построенный относительно GRCh38. Он содержит 88 гаплотипов 44 человек; гаплотип — одна из двух копий генома, материнская или отцовская.

Графовый выравниватель называется giraffe. Обычные инструменты не умеют работать с графом, поэтому giraffe проецирует результаты обратно на линейный GRCh38. Такая линейная проекция называется суржекцией. В итоге VCF получает привычные координаты GRCh38, необходимые для Just-DNA-Lite.

Есть важная тонкость: полный граф содержит слишком много вариантов, из-за чего риды могут выравниваться даже хуже, чем на линейный референс. Дополнительные пути сбивают выравниватель. Раньше граф прореживали по частоте: в версии d9 оставлены варианты, встречающиеся минимум у девяти гаплотипов. В июле Сергей пробовал пересобрать свой геном на d9. Теперь существует более удачный способ — выборка гаплотипов (Sirén et al., Nature Methods, 2024). Программа анализирует ваши риды и для каждого небольшого участка отдельно выбирает наиболее похожие из 88 гаплотипов HPRC. Из этих фрагментов она собирает два маршрута, соответствующие двум копиям человеческого генома, и добавляет GRCh38. Получается персональный диплоидный граф. Среди 44 доноров может не быть человека, похожего на вас целиком, но для отдельных участков подходящие «родственники» находятся.

Значение этого подхода видно в дополнительных таблицах статьи Sirén et al.: HG002, NovaSeq 30×, эталон GIAB v4.2.1, число ошибок. Обозначение 30× означает, что каждая буква генома прочитана в среднем 30 раз.

Нейросети вместо формул

Граф повышает качество первого этапа — выравнивания. Второй этап, вызов вариантов, удалось улучшить с помощью нейросетей.

Классический вызователь GATK HaplotypeCaller опирается на статистику и вычисляет вероятности по формулам. Он локально пересобирает гаплотипы, оценивает правдоподобие с учётом модели ошибок секвенатора, а затем варианты фильтруются по вручную подобранным порогам или с помощью дополнительного машинного обучения. Модель ошибок описывают люди, поэтому каждый новый тип артефактов приходится сначала изучать, а затем программировать. К тому же семейству относится freebayes, использованный для файла Вадима.

DeepVariant от Google работает иначе. Каждую подозрительную позицию он преобразует в «картинку» — стопку окружающих ридов. Буква, её качество, MAPQ и направление рида кодируются цветами, как в обычном изображении. Затем картинку классифицирует свёрточная нейросеть, подобная тем, что распознают кошек на фотографиях, но обученная на эталонных геномах. Она самостоятельно усваивает признаки артефактов секвенатора, ридов из паралогов и ошибок в гомополимерах. Для нового типа данных сеть переобучают, не переписывая формулы.

Pangenome-aware DeepVariant (Asri et al., 2025, UCSC и Google) идёт дальше: к изображению ридов добавляются гаплотипы пангенома HPRC для той же области. Сеть видит, например, что такой же индел присутствует у 30 из 88 гаплотипов, и ей проще отличить настоящий вариант от артефакта. При этом инструмент не заменяет giraffe: выравниватель, BWA или giraffe, размещает риды, а обычный или pangenome-aware DeepVariant определяет варианты.

Следующие цифры взяты из дополнительной таблицы 4 препринта Asri et al.: HG002, Illumina NovaSeq 35×, эталон T2T-Q100. Указана сумма пропусков и лишних вызовов:

CDS — участки генов, по которым клетка синтезирует белки. Они занимают около 1,5% генома, но изучены лучше остальных: здесь сосредоточена большая часть известных мутаций, вызывающих наследственные заболевания. Остальная часть включает регуляторные участки, включающие и выключающие гены, промежутки внутри генов, повторы и области с пока неизвестными функциями. Варианты для полигенных оценок чаще находятся именно в регуляторных участках, поэтому нам важна точность по всему геному.

Нейросетевой DeepVariant допускает вдвое меньше ошибок, чем GATK. Граф дополнительно снижает их число на 15%: для SNV — на треть, для сегментных дупликаций и генов — почти вдвое. Pangenome-aware даёт ещё около 5% улучшения, а в дупликациях и генах — 20%. Коммерческий DRAGEN показывает лучший общий результат почти исключительно за счёт инделов в тандемных повторах; в дупликациях и генах связка giraffe + pangenome-aware его догоняет и даже превосходит.

Так сложился наш план: использовать giraffe с выборкой гаплотипов и DeepVariant, а pangenome-aware проверить отдельно.

Научно-популярная часть на этом заканчивается. Далее — подробности для тех, кто захочет повторить эксперимент.

Термины для чтения VCF

В технической части понадобятся ещё несколько определений.

Шкала Phred. Способ выразить вероятность ошибки: Q 10 соответствует 10%, Q 20 — 1%, Q 30 — 0,1%. Чем больше значение, тем выше уверенность.

MAPQ. Выраженная по шкале Phred уверенность в правильном размещении рида. Например, MAPQ 20 соответствует оценочной вероятности ошибки 1%. MAPQ 0 означает, что есть другое место, подходящее не хуже.

GQ и QUAL. GQ отражает уверенность программы в генотипе, например 0/1. QUAL показывает уверенность в том, что позиция вообще отличается от референса. DeepVariant ограничивает GQ значением 99.

DP, AD, VAF. Характеристики покрытия. DP — число ридов, покрывающих позицию. AD — число ридов с референсной буквой и с вариантом. VAF — доля ридов с вариантом. Для настоящей гетерозиготы VAF близок к 0,5, для гомозиготы — к 1. Пример для повтора TATATA…: 0/1, GQ 20, DP 40, AD 29,8, VAF 0,2. Здесь гетерозигота определена с уверенностью всего около 99%, а доля варианта подозрительно мала: в повторах риды часто ошибаются на единицу повтора. Такую находку следует перепроверить.

FILTER: PASS, RefCall, NoCall. Итоговое решение вызователя по записи. PASS — вариант принят; RefCall — нейросеть сочла последовательность референсной; NoCall — определить генотип не удалось. Интерпретаторы дополнительно применяют пороги: Just-DNA-Lite по умолчанию сохраняет записи с QUAL ≥ 20 и DP ≥ 10. Далее мы называем это результатом «после фильтра».

Нормализация вариантов. Один индел в повторе можно записать несколькими способами: удаление CA из CACACA допустимо поставить в любую позицию повтора. Нормализация с помощью bcftools norm приводит записи к единому виду: сдвигает вариант максимально влево, обрезает общие буквы по краям и разделяет записи с несколькими вариантами. Без неё сравнение двух VCF обнаруживает мнимые расхождения.

Внутреннее устройство конвейера

vg и giraffe. vg — набор программ для работы с графами-пангеномами, а giraffe — одна из его команд. Граф HPRC v1.1 построен с помощью Minigraph-Cactus.

Выборка гаплотипов. KMC подсчитывает во всех ридах подстроки длиной 29 — k-меры. Giraffe использует индекс гаплотипов в файле .hapl: хромосомы разделены там на участки примерно по 10 тысяч букв, и для каждого указаны k-меры, отличающие гаплотипы друг от друга. Выбор происходит в два этапа. Сначала для каждого участка отбирают 32 кандидата, чьи k-меры лучше всего совпадают с ридами. Затем перебирают пары кандидатов и сохраняют ту, которая совместно лучше остальных объясняет данные. Этот режим называют диплоидной выборкой. Именно так построен персональный граф в приведённой выше таблице Sirén et al.

Три этапа DeepVariant.

1. make_examples читает риды, находит позиции-кандидаты и создаёт для каждой «картинку». Этап работает только на CPU; версии для GPU нет. Поэтому процессор активно нагружен даже при использовании видеокарты.

2. call_variants передаёт картинки нейросети, которая определяет генотип. Этот этап можно выполнять на GPU.

3. postprocess_variants собирает результаты в VCF и gVCF и тоже работает на CPU.

Модели DeepVariant. Google предоставляет отдельные модели для разных данных. Их выбирают флагом --model_type: WGS для полного генома с короткими ридами, WES для экзома, PACBIO и ONT_R104 для длинных ридов и другие. У WGS есть дополнительная «малая модель», которая самостоятельно обрабатывает простые однозначные позиции без большой нейросети. Pangenome-aware поставляется как отдельная модель и отдельный образ; малой модели у неё нет.

Где запускать обработку

Этим программам необходимо много оперативной памяти и процессорных ядер. Оптимальная конфигурация — 256 ГБ RAM и 32–64 ядра. Покупка такого сервера обходится дорого, поэтому мы арендовали его у RUVDS на три дня.

Требования к памяти.

• подготовка индекса гаплотипов — 40 ГиБ;

• giraffe с персональным графом — около 65 ГБ, параллельная сортировка требует ещё 50–70 ГБ;

• pangenome-aware DeepVariant — около 4,2 ГБ на параллельный процесс и ещё 12 ГБ для графа.

На домашнем компьютере с 16–32 ГБ конвейер не запустится: первый ресурсоёмкий этап завершится из-за нехватки памяти. Если памяти достаточно, но доступно только восемь ядер, обработка одного генома займёт много дней.

Почему аренда. В 2026 году серверная DRAM подорожала примерно на 90 % за квартал. Для ограниченной по времени задачи, подобной нашей, аренда выгоднее покупки.

RUVDS предлагает выделенные серверы с оплатой по дням. Производительный сервер на Xeon Gold стоит от 5700 рублей в сутки; также доступна линейка на AMD EPYC.

Наш выбор. AMD EPYC 7713 (Zen 3, Milan): 64 физических и 128 логических ядер, 256 ГБ DDR4-3200. Трёх дней хватило на несколько пересборок, эксперименты и сравнения. При выборе сервера главное — достаточный объём RAM.

Память важнее числа ядер. Мы убедились в этом на собственном опыте. Запуск pangenome-aware с 80 процессами привёл к использованию 123 ГБ свопа при пике 244 из 251 ГБ RAM. С 52 процессами обработка прошла без свопа, а каждый процесс ускорился на 43%: около 1320 кандидатов в минуту вместо 920. Общая скорость почти не изменилась: 28 «лишних» процессов в основном ожидали подкачки с диска. При этом CPU был загружен лишь примерно на 41%.

Диск. На арендованном сервере беспокоиться о нём не пришлось: серверы RUVDS на AMD EPYC оснащены NVMe-накопителями суммарным объёмом от 9 ТБ.

На личном компьютере с HDD и небольшим SSD данные придётся распределить:

• только SSD или NVMe: временные копии FASTQ для выравнивания, персональный граф с индексами, BAM и временные файлы DeepVariant. Giraffe одновременно читает оба файла ридов, поэтому на HDD головка постоянно перемещается: в июле программа четыре часа почти не выдавала результатов. Копирование 64 ГБ на NVMe занимает минуты, а после выравнивания копии можно удалить;

• HDD подойдёт для хранения исходных FASTQ, итоговых CRAM и VCF, а также временных файлов samtools sort. Последние читаются и записываются последовательно, но занимают десятки гигабайт.

Как выбрать число потоков

Число потоков зависит от программы. Главное ограничение — не допускать активного использования свопа.

Формула для многопроцессных задач: N = min(60, (RAM − общая память − запас ~20 ГБ) / память на процесс). Если free или vmstat показывают рост свопа, N следует уменьшить.

История с видеокартой

Эту часть мы выполняли на домашнем компьютере Сергея с майнинговой NVIDIA CMP 90HX. В третьей статье он утверждал, что DeepVariant обработал его геном на этой карте. Позднее выяснилось, что это не так.

Во время работы нейросетевого этапа call_variants команда nvidia-smi показывала загрузку 0% и использование 1 МиБ памяти. Устройства /dev/nvidia* принадлежали группе video с правами 660. Контейнер DeepVariant запускался от обычного пользователя (-u uid:gid), но без этой группы. Внутри контейнера nvidia-smi работала, однако TensorFlow молча получал пустой список GPU и использовал процессор. Проблема решается флагом --group-add, но есть нюанс: --group-add video Docker интерпретирует по /etc/group контейнера, где video = 44, тогда как на хосте — 27. Поэтому нужно передавать номер группы устройства:

docker run --gpus all --group-add $(stat -c %g /dev/nvidia0) -u $(id -u):$(id -g) ... \
  run_deepvariant ... --call_variants_extra_args="execution_hardware=accelerator"

Флаг execution_hardware=accelerator заставляет DeepVariant завершиться с ошибкой, если GPU недоступен. Без него программа незаметно переключается на CPU.

Затем нас ждал сюрприз: на видеокарте обработка стала медленнее. У CMP 90HX шина ограничена PCIe 1.0 x4. Карта обрабатывала примерно 190–240 примеров в секунду против 520 у процессора.

Второй сюрприз оказался противоположным. Флаг execution_hardware=cpu не запрещает TensorFlow использовать видимую в контейнере видеокарту. При обработке третьего генома нейросетевой этап самостоятельно переключился на медленную карту и занял около трёх часов вместо полутора. Чтобы гарантированно использовать CPU, GPU не нужно передавать контейнеру вообще: запускайте его без --gpus all.

NVIDIA Parabricks

NVIDIA предлагает Parabricks — адаптированные для GPU инструменты геномики, включая BWA-MEM, HaplotypeCaller, DeepVariant, vg giraffe и pangenome-aware DeepVariant. На серверных GPU обработка полного генома от FASTQ до VCF занимает десятки минут вместо часов. Пакет бесплатен и распространяется как Docker-образ через NVIDIA NGC. Требуется видеокарта минимум с 16 ГБ памяти; некоторым инструментам по умолчанию нужно больше. У карты Сергея лишь 10 ГБ и ограниченная шина, поэтому мы использовали CPU. Если доступны A100/H100 или хотя бы карты с 24 ГБ памяти, начать стоит с Parabricks.

Этапы конвейера

vg и DeepVariant доступны в готовых контейнерных образах со всеми зависимостями: TensorFlow, CUDA и необходимыми версиями библиотек. Установка вручную занимает время и легко приводит к проблемам, поэтому оба инструмента запускаются в Docker. Контейнер видит только подключённые каталоги, а результаты сохраняет на диск хоста.

Время приведено для генома Вадима: 431,5 млн пар, покрытие около 42×.

Итого около 6–7 часов с обычным DeepVariant и 11–12 часов с pangenome-aware. Размер CRAM составил 16,7 ГБ против 41,6 ГБ у BAM.

Сборка генома по ридам Illumina

Удалось выровнять 98,52% ридов, правильно спаренных — 97,41%. По покрытию определён мужской пол: chrX — ×0,51, chrY — ×0,32 от покрытия аутосом. Поэтому X и Y обрабатываются как гаплоидные. Обнаружилась и особенность: 28,5% дубликатов. Это риды, прочитанные с одного исходного фрагмента ДНК. Они возникают при ПЦР-размножении материала и в проточных ячейках NovaSeq как «оптические» дубликаты соседних кластеров. Новых данных они не добавляют, но могут размножить ошибку ПЦР, поэтому их помечают с помощью samtools markdup, после чего DeepVariant их игнорирует. Эффективное покрытие составило около 30× вместо 42×, чего для генома достаточно.

На одних и тех же выравниваниях мы проверили два варианта:

1. обычный DeepVariant WGS с настройками Google для giraffe;

2. pangenome-aware DeepVariant.

Разницу во времени объясняет малая модель. В обычном DeepVariant большая нейросеть обработала около 3 млн примеров, остальные разрешила малая модель. В pangenome-aware через сеть прошли все примерно 9 млн примеров, поэтому этап call_variants занял пять часов вместо полутора.

По метрикам pangenome-aware обнаруживает больше вариантов. Однако наличие в dbSNP и показатель Ts/Tv — отношение транзиций к трансверсиям, составляющее около 2 для настоящих вариантов, — ненадёжны в сложных областях. Сам dbSNP содержит множество различий между паралогичными копиями, записанных как SNP. Поэтому мы проверили непосредственно риды.

Где сборки расходятся. Для общих вариантов найдено около 35 тысяч различающихся генотипов, почти все — в тандемных повторах, сегментных дупликациях и сателлитах. В «лёгких» областях их всего 636. Чтобы выяснить, кто прав, мы выбрали по 300 случайных позиций с разными генотипами и подсчитали в CRAM риды с каждым аллелем. Для ориентира: уверенно определённой гетерозиготе соответствуют около 16 ридов с вариантом и 14 с референсом, а уверенной гомозиготе — около 30 ридов с вариантом.

Оставшиеся позиции неоднозначны: в них слишком мало однозначно выровненных ридов. В пересчёте на общее число позиций pangenome-aware оказывается прав примерно в 8,7 тысячи случаев, обычный — в 7,2 тысячи.

Варианты, найденные только одной моделью. После фильтра у pangenome-aware их 23,4 тысячи, из них 75% подтверждаются ридами. У обычной модели — 6 тысяч, подтверждаются 54%. Многие «дополнительные» вызовы pangenome-aware представляют собой ту же позицию, записанную иначе, например замену T>G вместо вставки T>TG, либо места, где обычный DeepVariant выдал NoCall.

Вывод: для генома Вадима pangenome-aware немного лучше. Это согласуется с эталонным сравнением Google, где также использовались данные Illumina NovaSeq. Цена улучшения — 7,3 часа работы DeepVariant вместо 2,7 и около трёх тысяч дополнительных сомнительных SNV. Вадиму мы передали сборку pangenome-aware.

Какие файлы передавать после пересборки.

• Обязательно: .vcf.gz и индекс .tbi. Это список вариантов для загрузки в сервисы анализа с указанием сборки GRCh38.

• Желательно: .g.vcf.gz с индексом, чтобы отличать отсутствие варианта от непрочитанной позиции, и .cram с индексом .crai. CRAM содержит все риды и позволяет без пересборки изучать отдельные участки, определять CYP2D6 и HLA специализированными программами и искать крупные делеции.

• Для справки: отчёты flagstat и stats, контрольные суммы md5.

CRAM содержит результат выравнивания, тогда как порог MAPQ и модель относятся к следующему этапу. Поэтому для всех вариантов вызова используется один и тот же CRAM.

Сборка геномов по ридам MGI

Геном Сергея секвенировали на MGI DNBSEQ: 403,5 млн пар. Химия и профиль ошибок этой платформы отличаются от Illumina. Мы сопоставили июльские сборки — BWA и giraffe на графе d9 — с новыми результатами обычного и pangenome-aware DeepVariant на тех же ридах.

Доля выровненных ридов почти не зависит от метода. На графе d9 она составила 94,61%, на персональном — 94,45%. Значит, 98,5% у Вадима связаны с Illumina, а не с графом. Около 5,5% невыровненных ридов Сергея — особенность образца.

Гаплоидный вызов X устранил ложные гетерозиготы. В июльской сборке на chrX вне PAR их было 2920, в сборке BWA — 3152. Все они ошибочны, поскольку у мужчины там только одна копия. Теперь их нет.

Pangenome-aware добавил шум на данных MGI. После фильтра модель выдала десятки тысяч дополнительных SNV с Ts/Tv 0,88, что близко к случайному шуму. Из них 39% сосредоточены примерно в 50 окнах размером в мегабазу возле центромер — повторяющихся областей в середине хромосом: chr1:122–124 Мб, chr13:16–17 Мб, chr11:52–53 Мб и других. Качество этих вызовов находится на границе фильтра. У Вадима на такие центромерные скопления приходится 15%. Вероятная причина — обучение модели на данных Illumina и Element, тогда как MGI для неё «чужая» платформа. Поэтому результаты эталонного сравнения Google на Illumina нельзя напрямую переносить на геном Сергея.

Для генома Сергея мы выбрали персональный граф и обычный DeepVariant с порогом MAPQ 0, как в рецепте Google для giraffe. При этом мы не доказали, что порог 0 лучше 5. В позициях с уверенно определённым генотипом результаты одинаковы; различия проявляются только в паралогах, где уверенность невысока. История выбора и первоначальной ошибки — в спойлере.

Третий геном: куда пропали 12% ридов

Третий геном принадлежит жене Сергея. Он тоже получен на MGI и содержит 398,6 млн пар. Мы обработали его по тому же рецепту за 6,5 часа. Получилось 5,02 млн PASS-вариантов, для 99% найдены rsID. На X обнаружено 73 тысячи гетерозигот, что нормально: у женщины две X-хромосомы, которые обрабатываются диплоидно.

Неожиданно низкой оказалась доля выровненных ридов — 87,47%. Однако дело не в конвейере: в лабораторном файле выравнивания, полученном с помощью BWA, невыровненными остались те же примерно 12%.

Мы изучили эти риды. Их качество высокое, адаптеров и мусора нет, но состав отличается: доля G и C составляет 55,8% с узким пиком против 42% у человеческих ридов. Это бактериальная ДНК. Поскольку образцом была слюна, вероятнее всего, речь о микрофлоре полости рта.

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

Итоги

Исходная задача решена. Вместо файла на hg19 Вадим получил VCF, gVCF и CRAM в координатах GRCh38, построенные на персональном графе с нейросетевым вызователем. Файл успешно загрузился в Just-DNA-Lite, а полигенные оценки были рассчитаны. Сами результаты мы не публикуем, поскольку это личные медицинские данные.

Одно предупреждение: мы хотели рассчитать все оценки из каталога. Через веб-интерфейс это происходит медленно, поэтому Сергей написал два скрипта для пакетного расчёта в нескольких процессах. Модели для полного набора занимают около 50 ГБ. Не рекомендуем отмечать все пункты подряд: выбирайте только необходимое. Установку Just-DNA-Lite, расчёт и интерпретацию полигенных оценок Сергей подробно разобрал в третьей и четвёртой частях.

Разница между лучшей и следующей за ней сборкой невелика и сосредоточена в сложных областях. При анализе любой сборки разумно оставлять лишь генотипы с высокой уверенностью, GQ ≥ 20, и осторожно интерпретировать находки в сегментных дупликациях. Надёжные результаты в геноме Сергея воспроизвелись во всех сборках: APOE, гомозигота Жильбера и отсутствие ложного кластера в PRSS1.

Ограничения остаются прежними. Конвейер обнаруживает только малые варианты — SNV и инделы. Для структурных перестроек, числа копий генов, CYP2D6 и HLA-типирования необходимы отдельные инструменты. Без эталона строго доказать превосходство одной сборки над другой нельзя. Окончательно сравнить рецепты позволит только обработка эталонного генома GIAB HG002 обоими способами.

Скоро чёрная пятница: в этом году она приходится на 27 ноября. Если вы давно планировали полногеномное секвенирование, можно вновь ожидать скидок. Сразу запросите у компании исходные риды FASTQ: они позволят пересобрать геном относительно нового референса и современным методом, независимо от того, насколько устареет выданный VCF.

Скрипты

Полные версии скриптов доступны в репозитории sukamenev/pangenome-rebuild:

• pipeline_v2_mq0.sh — giraffe с выборкой гаплотипов и DeepVariant WGS с MAPQ 0 для MGI;

• pipeline_v2_pangenome_aware.sh — аналогичный конвейер, но на шаге 7 используется pangenome-aware DeepVariant через apptainer, 52 процесса и гаплоидные X/Y для Illumina;

• annotate_rsid.sh — поиск rsID в dbSNP по нормализованным аллелям;

• compare_vcfs.sh, gq_discord.sh, pileup_mapq.py, pileup_sites.py — сравнение сборок без эталона, распределение расхождений по GQ и зонам GIAB, проверка спорных позиций по ридам;

• pdl.sh, s3etag.py — параллельное скачивание с S3 и проверка ETag.

Запуск: SAMPLE=… FQ1=… FQ2=… OUTDIR=… SEX=male|female|auto ./pipeline_v2_mq0.sh. Выполненные этапы помечаются, поэтому после сбоя скрипт продолжает работу с нужного места.

# k-меры ридов для выборки гаплотипов
kmc -k29 -m128 -okff -t120 -hp @fq.paths $SAMPLE kmctmp

docker run --rm -u $(id -u):$(id -g) -v $GRAPHDIR:/g:ro -v $STAGE:/s quay.io/vgteam/vg:v1.75.1 \
    vg giraffe -Z /g/$GBZ --haplotype-name /g/$HAPL --kff-name /s/$SAMPLE.kff \
      --index-basename /s/idx/$SAMPLE \
      -f /s/R1.fq.gz -f /s/R2.fq.gz -N $SAMPLE -R $SAMPLE \
      -o SAM --ref-paths /s/GRCh38.path_list.txt -P -L 3000 -t 120 \
  | sed -E 's/GRCh38#0#//g' \
  | samtools fixmate -m -@ 24 - - \
  | samtools sort -@ 24 -m 3G -T $SORTTMP/s - \
  | samtools markdup -@ 24 -O bam - $SAMPLE.bam

docker run --rm -u $(id -u):$(id -g) -v $WORK:/w -v $OUT:/pg -v $STAGE:/s google/deepvariant:1.10.0 \
    /opt/deepvariant/bin/run_deepvariant --model_type=WGS \
    --ref=/w/ref/GRCh38.fa --reads=/pg/$SAMPLE.bam \
    --make_examples_extra_args="min_mapping_quality=0,keep_legacy_allele_counter_behavior=true,normalize_reads=true" \
    --haploid_contigs=chrX,chrY --par_regions_bed=/w/ref/GRCh38_PAR.bed \
    --output_vcf=/pg/$SAMPLE.vcf.gz --output_gvcf=/pg/$SAMPLE.g.vcf.gz \
    --num_shards=60 --intermediate_results_dir=/s/dv_tmp

Литература

1. Asri M., Chang P.-C., Mier J. C. et al. Pangenome-aware DeepVariant. bioRxiv, 2025. doi:10.1101/2025.06.05.657102

2. Sirén J., Eskandar P., Ungaro M. T. et al. Personalized pangenome references. Nature Methods, 2024. Препринт: doi:10.1101/2023.12.13.571553

3. Liao W.-W., Asri M., Ebler J. et al. A draft human pangenome reference. Nature 617, 312–324 (2023). doi:10.1038/s41586-023-05896-x

4. Sirén J., Monlong J., Chang X. et al. Pangenomics enables genotyping of known structural variants in 5202 diverse genomes. Science, 2021 (выравниватель vg giraffe). doi:10.1126/science.abg8871

5. Poplin R. et al. A universal SNP and small-indel variant caller using deep neural networks. Nature Biotechnology, 2018 (DeepVariant). doi:10.1038/nbt.4235

Документация: рецепт Google «giraffe + DeepVariant», метрики pangenome-aware, рабочий процесс команды vg, выборка гаплотипов в vg, NVIDIA Parabricks.

Авторы:

• Сергей Каменев — алгоритмист, программист, популяризатор IT и AI;

• Вадим Леонидович Зорин — врач-биофизик, кандидат биологических наук, эксперт Сколково, лауреат премии РФ по науке и технике;

• Олег Буторин — специалист по консультированию в области внедрения ИИ в бизнес.

© 2026 ООО «МТ ФИНАНС»

Источник