Предыдущие части цикла про полногеномное секвенирование: 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 — замена одной буквы. Индел — вставка или удаление нескольких букв. Пример: синдром Жильбера — это лишний повтор TA в гене UGT1A1, индел 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 исправлены тысячи ошибок сборки и закрыта часть пробелов. Современные базы и инструменты ориентируются на 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. С учётом направления транзиций возможно 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 раз.
Референс | Ошибок |
|---|---|
персональный граф (выборка гаплотипов) | 19 480 |
частотный граф d9 | 21 660 |
линейный GRCh38 (BWA) | 31 697 |
полный граф без выборки | 36 387 |
GATK 4 (для сравнения) | 74 718 |
Нейросети против формул
Граф улучшает первый шаг сборки, выравнивание. Второй шаг, вызов вариантов, улучшили нейросети.
Классический вызыватель вариантов 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. Число ошибок (пропуски плюс лишние вызовы):
Конвейер | Весь геном | SNV | Инделы | Сегм. дупликации | Гены (CDS) |
|---|---|---|---|---|---|
GATK 4.6 | 236 359 | 138 328 | 98 031 | — | — |
BWA + DeepVariant | 125 309 | 60 075 | 65 234 | 30 145 | 593 |
giraffe + DeepVariant | 107 102 | 40 457 | 66 645 | 16 697 | 312 |
BWA + pangenome-aware DV | 119 050 | 54 117 | 64 933 | 26 755 | 508 |
giraffe + pangenome-aware DV | 101 646 | 35 210 | 66 436 | 13 297 | 248 |
DRAGEN 4.3 (Illumina, коммерческий) | 73 745 | 26 797 | 46 948 | 14 287 | 347 |
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.
make_examplesчитает риды, находит места-кандидаты и рисует для каждого «картинку». Работает только на CPU, версии для GPU у этого шага нет. Поэтому даже при включённой видеокарте воет вентилятор процессора.call_variants— нейросеть смотрит на картинки и решает, какой генотип. Этот шаг можно считать на GPU.postprocess_variantsсобирает ответы в VCF и gVCF, тоже на CPU.
Модели DeepVariant. Google поставляет отдельные модели под разные данные, выбор флагом --model_type: WGS (полный геном, короткие риды), WES (экзом), PACBIO и ONT_R104 (длинные риды) и другие. У модели WGS есть ещё «малая модель»: простые однозначные места она решает сама, без большой нейросети. Pangenome-aware — отдельная модель и отдельный образ, малой модели в ней нет.
Где запускать
Этому софту очень нужна оперативная память, да и ядер нужно много. Оптимально железо на 256 ГБ и 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. Трёх дней хватило на несколько пересборок со всеми экспериментами и сравнениями. При выборе сервера главное — большая оперативная память.
Память важнее ядер. На этом мы и споткнулись. Запустили 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. Они пишутся и читаются последовательно, а места занимают десятки гигабайт.
Сколько потоков ставить
Число потоков выбирается под программу, а главное ограничение — не уходить в своп.
Программа | Во что упирается | Потоков |
|---|---|---|
| память; второй поток на ядре заполняет простои в ожидании RAM | до 120 |
| CPU, лёгкий, идёт параллельно с giraffe | 24, |
DeepVariant | CPU + ~1,5–2 ГБ на процесс | до 60 |
pangenome-aware | CPU + ~4,2 ГБ на процесс + 12 ГБ граф | (RAM − 12 − 20) / 4,2 ≈ 52 |
| вычисления; TensorFlow сам занимает ядра, сверх физических ядер не ускоряется | по умолчанию |
Формула для многопроцессных задач: 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×).
Шаг | Что | Время |
|---|---|---|
1 | FASTQ на NVMe (64 ГБ) | ~5 мин |
2 | граф HPRC v1.1 | ~15 мин, однократно |
2b | сборка индекса гаплотипов | 41 мин, однократно |
3 | KMC: k-меры ридов (k=29) | 6,5 мин |
4 | список путей GRCh38 в графе | <1 мин |
5 | giraffe: выборка гаплотипов + персональные индексы | ~27 мин |
5 | giraffe → суржекция (проекция на линейный GRCh38) → сортировка → дубликаты → BAM | ~1 ч 27 мин |
6 | пол по покрытию X и Y | <1 мин |
7 | DeepVariant (обычный, CPU) / pangenome-aware (CPU, 52 процесса) | ~2,5 ч / 7 ч 16 мин |
8 | rsID из dbSNP | ~25 мин |
10 | BAM → CRAM + сверка контрольных сумм | ~15 мин |
Итого ~6–7 часов с обычным DeepVariant и ~11–12 часов с pangenome-aware. CRAM получился 16,7 ГБ против 41,6 ГБ у BAM.
Грабли, на которые мы наступили
Готовый
.haplот HPRC не читается современным vg. Файл 2023 года записан в формате версии 1, а vg 1.75 ждёт 5–6 («Expected version 5 to 6, got version 1»). Более новой версии на сервере HPRC нет. Обновление vg не помогло бы: устарел сам файл. Собрали сами из графа:vg index -j … --no-nested-distance,vg gbwt -r,vg haplotypes -H. 41 минута, 40 ГиБ памяти.Amazon S3 зависал на TLS-рукопожатии по адресу
s3-us-west-2.amazonaws.com/human-pangenomics/…, а адресhttps://human-pangenomics.s3.amazonaws.com/…работал. Качали в 8 потоков по диапазонам и сверяли S3 ETag (это MD5 частей по 8 МиБ). Большие файлы иногда заметно быстрее качать напрямую, без VPN.Ноль rsID. Файл соответствия имён хромосом NCBI записан с переводами строк Windows, и хромосомы назвались
chr1\r. Аннотация отработала «успешно» и не нашла ни одного rsID. Теперь этап падает, если rsID меньше чем у половины вариантов.Терялись 60% инделов при поиске rsID. dbSNP пишет варианты в повторах «целиком» (
rs112329108:TACACACACACACAC → T, TAC, …), а DeepVariant кратко (TACAC → T, TAC). Сравнение строк видит разные записи. Решение: нормализовать обе стороны (bcftools norm -m -any), искать rsID по нормализованным аллелям и возвращать их в исходные записи. Найдено 99,4% инделов вместо 38%.Полный dbSNP (29,5 ГБ, 1,1 млрд записей) bcftools разбирает часами из-за огромного поля INFO. Один раз делаем облегчённую копию: хромосома, позиция, rsID, аллели.
«Ошибка контрольной суммы» CRAM оказалась ложной.
samtools checksumпишет в первую строку имя файла, а оно у BAM и CRAM разное. Сравнивать надо без строк-комментариев.Образ pangenome-aware DeepVariant не поместился на 30-гигабайтный отдельный раздел с Docker. Запустили через apptainer (контейнерная среда, популярная в биоинформатике), распаковав образ в каталог на btrfs.
В запускалке pangenome-aware нет флагов гаплоидных X/Y. Взяли команды этапов через
--dry_runи дописали--haploid_contigs chrX,chrY --par_regions_bedв шагpostprocess_variants.Временные файлы сортировки. При FASTQ от ~100 ГБ
samtools sortсбрасывает на диск десятки гигабайт. На третьем геноме NVMe заполнился почти до конца. Каталог для временных файлов лучше вынести на большую партицию.
Сборка генома на основе ридов Illumina
Выравнивание дало 98,52% выровненных ридов, правильно спаренных 97,41%. Пол по покрытию мужской (chrX ×0,51, chrY ×0,32 от аутосом), поэтому X и Y вызываются гаплоидно. Одна странность: дубликатов 28,5%. Дубликаты — это несколько ридов, прочитанных с одного и того же исходного фрагмента ДНК. Они возникают, когда при подготовке образца ДНК размножают ПЦР, а на проточных ячейках NovaSeq ещё и как «оптические» дубликаты соседних кластеров. Новой информации они не несут и могут размножить ошибку ПЦР, поэтому их помечают (samtools markdup), и DeepVariant их не учитывает. Эффективное покрытие получилось ~30× вместо ~42×. Для генома этого достаточно.
На одних и тех же выравниваниях мы запустили два варианта:
обычный DeepVariant WGS с рецептом Google для giraffe;
pangenome-aware DeepVariant.
обычный DeepVariant | pangenome-aware | |
|---|---|---|
PASS-вариантов | 5 019 977 | 5 054 315 |
Новые SNV (нет в dbSNP) | 46 129 | 49 546 |
Ts/Tv всех | 1,90 | 1,88 |
После фильтра (QUAL ≥ 20, DP ≥ 10): известные SNV | 3 697 460 | 3 757 697 |
После фильтра: известные инделы | 855 508 | 888 908 |
После фильтра: Ts/Tv | 2,03 | 2,02 |
Время DeepVariant (CPU) | ~2,7 ч | 7,3 ч |
Разница во времени объясняется малой моделью. У обычного DeepVariant через большую нейросеть прошли ~3 млн примеров, остальное решила малая модель. У pangenome-aware через сеть прошли все ~9 млн, отсюда 5 часов на шаге call_variants вместо полутора.
По метрикам pangenome-aware находит больше. Но «есть в dbSNP» и Ts/Tv (отношение транзиций к трансверсиям, у настоящих вариантов около 2) в трудных зонах ненадёжны: dbSNP сам полон отличий между копиями паралогов, почти одинаковых участков в разных местах генома, записанных как SNP. Поэтому мы полезли в сами риды.
Где расходятся сборки. Разных генотипов на общих вариантах ~35 тыс., почти все в тандемных повторах, сегментных дупликациях и сателлитах, то есть в трудных зонах из первой половины статьи. В «лёгких» участках генома их 636. Чтобы понять, кто прав, мы взяли по 300 случайных мест с разными генотипами и посчитали прямо в CRAM, сколько ридов несут каждый аллель. Для сравнения: у гетерозиготы, в которой программа уверена, ~16 ридов с вариантом и ~14 с референсом, у такой же гомозиготы ~30 с вариантом.
Расхождение | Риды показывают | Риды показывают |
|---|---|---|
обычный | 34% — прав обычный | 18% — прав pangenome-aware |
обычный | 44% — прав pangenome-aware | 17% — прав обычный |
Остальное — неясные места, где мало однозначно выровненных ридов. В пересчёте на число мест 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, на геном Сергея не переносится.
После фильтра (QUAL ≥ 20, DP ≥ 10) | BWA (июль) | giraffe d9 (июль) | новая, обычный DV | новая, pangenome-aware |
|---|---|---|---|---|
известные SNV | 3 604 475 | 3 630 598 | 3 651 002 | 3 718 379 |
новые SNV | 4 324 | 7 467 | 7 424 | 7 208 |
гетерозиготы на chrX | 1 408 | 1 304 | 0 | 0 |
Для генома Сергея мы оставили персональный граф + обычный DeepVariant с порогом MAPQ 0, как в рецепте Google для giraffe. Про сам порог оговорюсь: мы не доказали, что 0 лучше 5. Там, где программа уверена в генотипе, оба порога дают одно и то же. Разница видна только в паралогах, где уверенность низкая. Как мы выбирали и как сначала ошиблись — в спойлере.
Порог MAPQ: 0 или 5
У DeepVariant есть параметр min_mapping_quality: риды с MAPQ ниже порога нейросеть не видит. По умолчанию он равен 5, а в рецепте Google для giraffe стоит 0 (в версии рецепта r1.6 был 1, с r1.8 стал 0). Рабочий процесс команды vg для обычной модели оставляет 5, для pangenome-aware ставит 0. Между vg 1.67 у Google и нашей 1.75 к тому же переписали оценку выравниваний, так что шкала MAPQ могла сдвинуться. Мы пересчитали DeepVariant с порогом 5 на тех же выравниваниях обоих геномов.
Сначала метрики сбили с толку. MQ 5 выглядел чище: вдвое меньше новых SNV, Ts/Tv 1,96 вместо 1,90. Нашлось ~21–23 тыс. мест, где MQ 0 даёт гетерозиготу, а MQ 5 гомозиготу, причём покрытие при MQ 0 там в 2,3 раза выше. Напрашивался вывод: это риды из паралога, MQ 5 чинит ложные гетерозиготы. Opus уже предложил заменить выданные файлы на MQ 5. Мы показали выводы Astra, и она возразила: повышенное покрытие ещё не доказывает, что MQ 5 прав, Ts/Tv нельзя переводить в долю ошибок, и надо смотреть на генотипы с высокой уверенностью.
Расхождения есть только там, где программа не уверена в генотипе. На позициях, где обе сборки уверены в генотипе (GQ ≥ 20), конфликтов почти нет: у Вадима 12 из 4,3 млн, у Сергея 13 из 4,3 млн. Вся разница приходится на вызовы с низкой уверенностью в сателлитах и сегментных дупликациях.
Риды показали обратное первому впечатлению. Проверяли на геноме Вадима. В 400 спорных местах (0/1 у MQ 0, 1/1 у MQ 5) медиана ридов с вариантом и высоким MAPQ равна 13, как у гетерозиготы (у гомозигот с высокой уверенностью 30). «Гомозиготная» картина нашлась лишь в 16% мест, «гетерозиготная» в 71%. Механизм: в паре почти одинаковых копий риды с вариантом отличимы и получают высокий MAPQ, а риды с референсной буквой совпадают с обеими копиями и получают MAPQ 0, включая собственные референсные риды человека. Порог 5 их выбрасывает и превращает настоящую гетерозиготу в «гомозиготу».
Итог: по ридам в спорных местах чаще прав порог 0, и Google измерял точность связки giraffe + DeepVariant тоже с ним. Строгим доказательством это не является: в таких местах возможно разное число копий, а уверенность обеих сборок низкая. На геноме Сергея, с MGI, проверку по ридам мы не повторяли, там есть только сравнение по GQ с той же картиной.
Мораль для тех, кто сравнивает сборки без эталона: «есть в dbSNP» и Ts/Tv в трудных зонах ведут в неверную сторону. Решают риды.
Третий геном: куда делись 12% ридов
Третий геном — жены Сергея, тоже с MGI, 398,6 млн пар. Собрали тем же рецептом, что и геном Сергея, за 6,5 часа. PASS-вариантов 5,02 млн, rsID нашлись у 99%. На X 73 тыс. гетерозигот, и здесь это норма: у женщины две X, они вызываются диплоидно.
Удивила доля выровненных ридов: 87,47%. Конвейер тут ни при чём: в файле выравнивания от самой лаборатории, сделанном через BWA, невыровненных столько же, около 12%.
Мы посмотрели, что это за риды. Качество у них высокое, адаптеров и мусора нет. Отличается состав: доля букв G и C у них 55,8% с узким пиком, а у человеческих ридов 42%. Это ДНК бактерий. Образцом была слюна, так что, скорее всего, это микрофлора полости рта.
На результат это не влияет: человеческих ридов осталось на покрытие ~35×. Но за секвенирование бактерий вы тоже платите, и низкий процент выравнивания ещё не означает плохую сборку.
Что в итоге
Геном | Секвенатор | Сборка | Почему |
|---|---|---|---|
Вадим | Illumina NovaSeq | giraffe (персональный граф) + pangenome-aware DeepVariant | больше вызовов, подтверждённых ридами; меньше |
Сергей | MGI DNBSEQ | giraffe (персональный граф) + обычный DeepVariant, MAPQ 0 | pangenome-aware на MGI шумит у центромер |
жена Сергея | MGI DNBSEQ | то же, что у Сергея | тот же секвенатор |
Задача, с которой всё началось, решена. Вместо файла по 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. Шаги отмечаются метками, и после сбоя скрипт продолжает с того же места.
Две главные команды конвейера
Выравнивание на персональном графе, суржекция на GRCh38, сортировка и пометка дубликатов:
# 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
Обычный DeepVariant с рецептом Google для giraffe. Для мужчины добавляются гаплоидные X и Y:
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
Литература
Asri M., Chang P.-C., Mier J. C. et al. Pangenome-aware DeepVariant. bioRxiv, 2025. doi:10.1101/2025.06.05.657102
Sirén J., Eskandar P., Ungaro M. T. et al. Personalized pangenome references. Nature Methods, 2024. Препринт: doi:10.1101/2023.12.13.571553
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
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
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 ООО «МТ ФИНАНС»
