Собираем геном из ридов на графе-пангеноме открытым софтом

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

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

Не вышло. Геном оказался собран по устаревшему референсу 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 — замена одной буквы. Индел — вставка или удаление нескольких букв. Пример: синдром Жильбера — это лишний повтор TA в гене UGT1A1, индел C>CAT.

Генотип ^[6]. То, что стоит на двух копиях хромосомы в данной точке. В файле он записан цифрами: 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 исправлены тысячи ошибок сборки и закрыта часть пробелов. Современные базы и инструменты ориентируются на GRCh38. Одна и та же буква генома в hg19 и GRCh38 стоит на разных позициях, поэтому файл по hg19 нельзя сравнивать с базой по GRCh38. Проект для подсчёта полигенных оценок Just-DNA-Lite не умеет работать со сборкой GRCh37: при загрузке генома он завис.

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

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

Координаты из одного референса в другой можно перевести программно, это называется лифтовер. Но это идеологически неверно и даёт плохой результат: технологии выравнивания развиваются, референс GRCh37 устарел, и лифтовер — это пережёвывание жёваного. Он теряет варианты там, где сборки различаются, и не исправляет ошибки ^[7] выравнивания и вызова. Пересборка генома из ридов в 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 — геном, собранный полностью, от теломеры до теломеры. Теломеры ^[8] — структуры на обоих концах хромосомы; при каждом делении клетки они укорачиваются и работают как счётчик оставшихся делений. Точность эталона — около одной ошибки на миллион букв. Он покрывает и трудные зоны, куда старые эталоны не дотягивались.

Конвейер прогоняют на ридах эталонного генома и считают три числа. Точность — какая доля найденных вариантов настоящие. Полнота — какую долю настоящих вариантов нашли. 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. С учётом направления транзиций возможно 4, трансверсий 8. Поэтому у случайных замен, то есть у ошибок секвенатора и прочего мусора, 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 оставлены варианты, которые встречаются минимум у 9 гаплотипов. В июле Сергей пробовал пересобрать свой геном на 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 работает на статистических методах: вероятности он считает по формулам. Он локально пересобирает гаплотипы, считает правдоподобия по модели ошибок секвенатора, а потом варианты фильтруют вручную подобранными порогами или машинным обучением ^[9] поверх. Модель ошибок описана людьми, и каждый новый тип артефакта приходится понимать и кодировать. К тому же семейству относится freebayes, которым собран файл Вадима.

DeepVariant от Google устроен иначе. Каждое подозрительное место он превращает в «картинку»: стопку ридов вокруг позиции. Буква, её качество, MAPQ (уверенность в месте рида) и направление рида закодированы в ней как цвета в обычном изображении. Картинку классифицирует свёрточная нейросеть — того же типа, что распознают кошек на фотографиях, — обученная на эталонных геномах. Сеть сама выучивает, как выглядит артефакт секвенатора, рид из паралога или заикание ^[10] в гомополимере. Под новый тип данных сеть переобучают, формулы переписывать не приходится.

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 — отдельная модель и отдельный образ, малой модели в ней нет.

Где запускать

Этому софту очень нужна оперативная память ^[11], да и ядер нужно много. Оптимально железо на 256 ГБ и 32–64 ядра. Покупать такой сервер дорого, поэтому мы арендовали его на три дня у RUVDS.

Сколько нужно памяти.

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

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

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

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

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

У RUVDS есть выделенные серверы ^[13] с посуточной оплатой. Мощный сервер на Xeon Gold стоит от 5700 рублей в сутки, есть и линейка на AMD EPYC ^[14].

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

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

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

Если пересобирать на личном железе, где есть HDD и небольшой SSD, данные придётся делить:

• только SSD или NVMe: копии FASTQ на время выравнивания, персональный граф с индексами, BAM и временные файлы DeepVariant. Giraffe читает оба файла ридов одновременно, и на HDD это превращается в бесконечные скачки головки: в июле giraffe четыре часа почти ничего не выдавал. Копирование 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 взять видеокарту, если она видна в контейнере. На третьем геноме шаг нейросети сам ушёл на медленную карту и занял около трёх часов вместо полутора. Чтобы считать на процессоре, видеокарту в контейнер не надо передавать вовсе: без --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 млн, отсюда 5 часов на шаге 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 ч и ~3 тыс. лишних сомнительных 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) с новыми, обычным DeepVariant и pangenome-aware, на одних и тех же ридах.

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

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

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 часа. PASS-вариантов 5,02 млн, rsID нашлись у 99%. На X 73 тыс. гетерозигот, и здесь это норма: у женщины две X, они вызываются диплоидно.

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

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

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

Что в итоге

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

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

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

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

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

Скрипты

Полные скрипты лежат в репозитории sukamenev/pangenome-rebuild ^[16]:

• 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 ^[17]

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

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 ^[19]

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 ^[20]

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 ^[21]

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

Авторы:

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

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

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

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