Простой

13 мин

272

В прошлой статье я рассказал, как получил Starter Kit для сбора материалов для секвенирования и отправил его в США. А также о встроенном в Nebula функционале для его анализа.

В этой статье я приведу минимальную базу о том, как выполняется секвенирование, а потом разберём типы и содержание конкретных файлов, полученных мной от Nebula.

База базовая

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

Как машина вообще читает ДНК

Расскажу на примере доминирующей технологии — секвенаторов Illumina.

Сначала — копии на стекле

Фрагмент ДНК попадает на проточную ячейку — стеклянную пластину, густо засеянную короткими синтетическими зацепками, к которым ДНК прилипает концом. Дальше из одной молекулы делают тысячи копий, и делают их прямо на стекле.

Механизм называется мостиковой амплификацией. Амплификация — это размножение, наработка множества копий; делают её средствами ПЦР, только здесь весь процесс протекает прямо на поверхности стекла. Молекула цепляется за зацепку, ДНК-полимераза достраивает её копию, а исходную нить смывают. Оставшаяся копия изгибается дугой, дотягивается свободным концом до соседней зацепки, полимераза копирует её через этот «мостик», две нити разъединяют — и копий уже две. Цикл повторяется, число копий удваивается раз за разом. Та же логика, что в ПЦР-тесте на вирус: берут несколько фрагментов вирусного генома и экспоненциально размножают, пока из неуловимых единиц не выйдут миллиарды, которые уже видит прибор. Разница в том, что обычная ПЦР идёт в пробирке и копии плавают вперемешку, а здесь она приклеена к стеклу.

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

Зачем вообще тысяча копий? Одна молекула с одной светящейся меткой для камеры слишком тускла. А тысяча одинаковых копий, которые на каждом шаге зажигаются синхронно и одним цветом, складываются в яркую точку — её уже видно.

На этом этапе никакого света ещё нет: размножение идёт обычными, непомеченными буквами, а полимеразу и реактивы потом вымывают. Свет включится дальше.

Чтение по одной букве за такт

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

Дальше — цикл:

1. залили буквы, в каждом кластере пристроилась одна, реакция встала;

2. сфотографировали ячейку — каждый кластер светится цветом своей буквы;

3. смыли метку и сняли «стоп»;

4. повторили.

Сколько циклов — столько прочитанных букв. Стандарт сегодня — 150.

Порядок букв задаётся не пространством, а временем. Если бы все буквы фрагмента вспыхнули разом, получился бы винегрет, по которому видно только соотношение букв, но не их порядок. А так порядок вспышек кластера во времени даёт последовательность: синяя вспышка, потом красная, потом две зелёных — G, A, T, T.

«Цвет своей буквы» — упрощение из учебников. Современные машины (NovaSeq, NextSeq) не используют четыре разные краски. Там две краски и два снимка за цикл: C — красная, T — зелёная, A светится в обоих каналах, а G не светится вообще. Буква G — это темнота, отсутствие сигнала в обоих кадрах. Так вдвое меньше реактивов и съёмки.

А «сфотографировали ячейку» — это не один щелчок. Стекло огромно по сравнению с полем зрения объектива, поэтому внутри стоит флуоресцентный микроскоп с камерой, который проходит ячейку по крошечным участкам — их называют тайлами, — снимая за цикл тысячи кадров; на каждом — множество кластеров разом. Это больше похоже на работу микроскопа, методично обходящему предметное стекло, чем на планшетный сканет или цифровой фотоаппарат, который фотографирует всё поле сразу. За цикл ячейку обходят дважды — отдельно красный и зелёный кадр (на старых четырёхкрасочных машинах проходов было четыре).

Почему файлов два: парноконцевое чтение

Прочитав 150 букв с одного конца (это ушло в _1), машина делает химический разворот. Те же нити в том же кластере снова изгибаются мостиком, фермент-«ножницы» отрезает и вымывает уже прочитанные нити, и остаётся тот же кластер, но молекулы в нём развёрнуты другим концом вперёд. Машина читает ещё 150 букв — это уйдёт в _2.

Кластер за весь процесс не сдвигается с места, и камера читает второй конец в той же точке X, Y, где читала первый. Поэтому два конца одной молекулы связаны общим адресом на стекле, а программе-выравнивателю не нужно ничего сопоставлять — она берёт по строчке из _1 и _2 и сразу знает, что это пара. А зачем вообще читать с двух концов, а не сделать один рид подлиннее, станет видно на слепой зоне ниже.

Слепая зона посередине

Возникает логичный вопрос. Итак, фрагменты нарезают длиной около 400 букв, а читают по 150 с каждого конца. 150 + 150 = 300. Середину, около 100 букв, машина не читает вообще:

Почему просто не сделать 400 циклов и не прочитать всё? Мешает рассинхронизация кластера. Он работает, пока все тысяча копий шагают в ногу — пристраивают по одной букве за цикл. Но фермент неидеален: примерно в одном случае из тысячи он то отстаёт на шаг, то пристраивает сразу две буквы и убегает вперёд. С каждым циклом рассинхрон копится. На первом такте все копии вспыхивают одним цветом — сигнал чистый. К сто пятидесятому часть копий уже «поёт вразнобой», цвета подмешиваются, и оценка качества начинает падать. Сделай машина 300–400 циклов — на выходе был бы случайный мусор. Поэтому ставят жёсткий предел: дистанцию, на которой ещё гарантируется точность 99,9% (тот же показатель качества Phred, что в FASTQ).

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

Более того, пара концов работает якорем. Зная, что между концами одной молекулы должно быть около 400 букв, выравниватель проверяет, легли ли они на референс на нужном расстоянии. Если конец _1 сел нормально, а его пара _2 оказалась за 5000 букв или вообще на другой хромосоме — значит, в этом месте у человека выпал или переехал крупный кусок ДНК. Одиночными ридами такие крупные перестройки почти не поймать.

Сколько ДНК машина читает за раз

Цифры тут астрономические. Флагман Illumina NovaSeq X Plus за один запуск выдаёт до ~52 миллиардов ридов — это десятки миллиардов кластеров, одновременно светящихся на стекле.

Это доведённая до предела SIMD: один химический такт — одна «инструкция», которая применяется сразу к десяткам миллиардов независимых «потоков»-кластеров. Графический ускоритель параллелит задачу на десятки тысяч ядер; секвенатор оптически обрабатывает миллиарды дорожек за такт.

Прикинем, сколько надо на одного человека. Геном — около 3,1 миллиарда букв. Для тридцатикратного покрытия нужно прочитать ~90 миллиардов букв. При длине рида 150 это ~600 миллионов ридов, то есть ~300 миллионов пар. Вот они и весят два файла по ~55 ГБ. А раз один запуск машины — это десятки миллиардов ридов, за раз секвенируют не одного человека, а сразу несколько десятков.

Как же не перепутать, чья ДНК чья? А никак — машина и не разделяет. Кластеры всех людей лежат на одном стекле вперемешку, и секвенатор читает их подряд, не зная и не интересуясь, чей это фрагмент. Сортировка — потом, в софте.

Чтобы она была возможна, к каждому фрагменту ещё на этапе подготовки пришивают короткую метку-баркод (8–10 букв), причём в современных протоколах сразу с двух концов — своя слева и своя справа. Двойная метка страхует от путаницы: если из-за химии фрагмент случайно подхватит чужую метку с одного конца, со вторым концом она не совпадёт, и такой рид выбросят. Сами метки секвенатор читает отдельными короткими проходами, между чтением двух концов фрагмента.

А первый же шаг обработки сырых данных — сортировка: софт разбирает общую кучу ридов по баркодам и раскладывает их по индивидуальным файлам .fq.gz для каждого конкрентного человека.

Мой набор файлов после полногеномного секвенирования

55G NG101J7C24_1.fq.gz
56G NG101J7C24_2.fq.gz
57G NG101J7C24.mm2.sortdup.bqsr.cram
1,4M NG101J7C24.mm2.sortdup.bqsr.cram.crai
215M NG101J7C24.mm2.sortdup.bqsr.hc.vcf.gz
1,7M NG101J7C24.mm2.sortdup.bqsr.hc.vcf.gz.tbi

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

Мы все — практически 100% клоны!

Геномы двух случайных, неродственных людей совпадают примерно на 99,9%. Швед и кореец, вы и сосед по лестничной клетке — на уровне ДНК почти неотличимые копии одного и того же текста. Различается около одной тысячной, и в эту тысячную упаковано всё: цвет глаз, рост, переносимость лактозы, предрасположенность к болезням — всё, что делает человека уникальным.

От генома нам нужен только diff

Раз люди настолько одинаковы, хранить геном каждого целиком, бессмысленно — это терабайты ради того, что у всех совпадает. Поэтому договорились об одном общем образце: есть референсный геном (сборка GRCh38). Это не геном какого-то конкретного человека, а «типовой» геном, собранный из ДНК нескольких анонимных доноров, — усреднённый представитель вида Homo Sapience. А геном конкретного человека описывают как отличия от этого образца.

Поэтому ваш геном хранят не целиком, а как diff против эталона.

И эталон, и то, что снимают с вас, — это геном целиком, а не только гены. За кодирование белков — химических наномашин, на которых держится работа нашего тела, — отвечает всего около 1,5% ДНК; остальное долго называли «мусором», хотя там сидит регуляция и масса пока непонятого. Полногеномное секвенирование читает всё подряд, поэтому, если завтра какой-нибудь «мусорный» участок признают важным, пересдавать анализ не нужно — он уже снят. И даже когда сменится сам эталон (GRCh38 собран с пробелами примерно на 8%, а в 2022-м вышла бесшовная сборка T2T-CHM13), пересеквенировать ничего не придётся: те же данные просто заново сопоставят с новым эталоном. Сырьё с вас снимают один раз. Оговорка одна: совсем уж повторяющиеся куски короткие чтения берут плохо даже при 30× — но это предел длины чтений, а не деление на гены и «мусор».

Так же работает система контроля версий и инкрементальный бакап: есть базовый образ, дальше пишутся только изменения. Мой набор файлов — три уровня этой же идеи.

FASTQ  — сырые риды, как они вышли из машины       111 ГБ (_1 + _2)
CRAM   — те же риды, но записаны как отличия         57 ГБ
         от референса
VCF    — только сами отличия, готовый список        215 МБ

Дальше — по каждому уровню.

Разбираем файлы по типам

FASTQ — сырые чтения (риды)

NG101J7C24_1.fq.gz и NG101J7C24_2.fq.gz — первичные данные прямо из секвенатора, упакованные gzip (стандартный линуксовый архиватор). Внутри — сотни миллионов считанных коротких кусочков ДНК, их называют ридами, и для каждой прочитанной буквы оценка качества: насколько машина в ней уверена.

Файлов два, потому что каждый фрагмент читают с двух концов навстречу друг другу — это парноконцевое чтение (о нём мы говорили ранее). Левые концы всех фрагментов идут в _1, правые — в _2. Файлы синхронны построчно: рид номер N в первом и рид номер N во втором — это два конца одной молекулы.

CRAM — риды, выровненные по референсному геному

...mm2.sortdup.bqsr.cram — те же риды из FASTQ, но уже привязанные к конретным местам референса. «Выравнивание» — это стандартный термин (его же называют маппингом): для каждого короткого рида программа находит то место в референсе-эталоне, где он должен быть. Представьте, что книгу с референсным геномом пропустили через шредер, который порезал текст на кусочки по 150 букв подряд и перемешали; выравниватель/маппер для каждого обрывка находит страницу и строку, где он мог находиться.

А когда место найдено, рид можно сильно сжать. Представим, что наш рид на 100% совпал с эталоном. Зачем заново сохранять на диск все эти 150 букв A, T, G, C?

Вместо этого формат CRAM сохраняет метаданные: «Рид номер такой-то лежит на хромосоме 7 в позиции 1234, совпадение полное». Физически сами буквы из файла удаляются — программа прочитает их из эталонного генома, когда вы решите посмотреть этот участок. Если же в риде есть мутация, CRAM запишет её как патч к коду: «хромосома 7, позиция 1234, дальше всё совпадает, только буква 50 — A вместо G». На этом выносе букв в эталон и держится экономия формата. В отличие от старого формата BAM, который «в лоб» хранит каждую букву каждого рида, CRAM оставляет только координаты и отличия. Отсюда и арифметика: 111 ГБ исходных FASTQ ужались в 57 ГБ.

CRAM не удаляет дубли. Геном прочитан с покрытием 30×, значит, на каждое место эталона приходится около тридцати ридов. CRAM честно сохранит их все в базе. Он просто тридцать раз запишет короткую фразу: «здесь лежит рид, совпадение полное». Это тот же полный набор из 600 миллионов ридов, что и в FASTQ. Именно поэтому файл CRAM невозможно раскодировать и посмотреть, если у вас на диске не лежит тот самый референсный геном — программе просто неоткуда будет взять удалённые буквы. При необходимости CRAM можно развернуть обратно в FASTQ.

Суффиксы в имени — журнал обработки: по нему виден весь пайплайн. Пайплайн (он же конвейер) на жаргоне биоинформатиков — это цепочка программ, где выход одной подаётся на вход следующей. Имя читается слева направо, шаг за шагом:

• mm2 — выравнивали minimap2, быстрым выравнивателем Хэна Ли (Heng Li).

• sortdup — риды отсортированы по координатам, а дубликаты, оставшиеся от копирования молекул в лаборатории, помечены, чтобы одну и ту же молекулу не считать несколько раз.

• bqsr — перекалибровка оценок качества. Секвенатор врёт не случайно, а систематически (например, после определённого сочетания букв). По известным точкам различий GATK строит модель этих систематических ошибок и поправляет оценки — сам Broad Institute называет это применением машинного обучения.

Сами суффиксы — не данные, а следы работы конвеера. Извлекать из них нечего: та же перекалибровка bqsr уже вшита в оценки качества внутри CRAM. Они могут помочь, когда нужно понять как был получен файл или когда нужно воспроизвести результат, сравнить с чужими данными или понять, почему два анализа разошлись.

Рядом лежит ...cram.crai на 1,4 МБ — индекс. Это оглавление к 57-гигабайтной книге: чтобы вытащить один ген на одной хромосоме, программе не нужно прочёсывать весь файл, она по индексу прыгает сразу к нужному месту. Без .crai большинство утилит CRAM просто не откроет — но это не данные: потеряется индекс — его можно пересоздать из CRAM одной командой за минуту.

VCF — собственно diff

...hc.vcf.gz — наконец, сама разница. Записаны только места, где моя ДНК отличается от референса: однонуклеотидные замены (одна буква поменялась на другую) и короткие вставки (инсерции) и выпадения букв (делеции) — так называемые инделы. Несколько миллионов строк, каждая в духе «в координате такой-то у тебя не как у эталона, а вот так».

Почему тогда CRAM (57 ГБ) в сотни раз больше VCF (215 МБ), раз оба — «отличия от референса»? Потому что отличия там разного масштаба. CRAM держит все 600 миллионов ридов: каждый рид со своими отличиями, с оценкой качества на каждую букву, и при тридцатикратном покрытии любая позиция представлена примерно тридцатью ридами — это сырые показания, со всеми ошибками чтения вперемешку. VCF — это уже чистовой результат: тридцать ридов в каждой точке свёрнуты в итоговый diff без дублирования, в файл попадают лишь те несколько миллионов позиций, где он разошёлся с эталоном, по строке на каждую.

Суффикс hc — это HaplotypeCaller из пакета GATK, программа, которая искала отличия. Как и mm2/sortdup/bqsr в имени CRAM, это пометка о том, чем сделан файл. Спорные места HaplotypeCaller разбирает дотошно: он не доверяет готовому выравниванию, а заново собирает участок из перекрывающихся ридов — строит граф и смотрит, какие варианты последовательности там вообще возможны. Поэтому он хорошо ловит даже сложные случаи, где рядом стоят несколько разных изменений.

...vcf.gz.tbi — снова индекс, теперь к VCF. .gz у VCF — это bgzip, gzip-совместимый, но поблочный вариант: архив можно листать кусками, потому к нему и возможен такой индекс. И .crai, и .tbi — не данные, а служебные оглавления: каждый пересоздаётся из своего файла за секунды и нужен лишь чтобы быстро находить запись по координате, не читая весь файл целиком. Нужны они профессиональному софту, который прыгает по координатам: геномным браузерам вроде IGV, утилитам типа samtools и bcftools, онлайн-интерпретаторам, проверяющим конкретную позицию. Обычному пользователю они не нужны — если вы грузите данные в онлайн-сервис, он создаст индексы сам.

Итоговые 215 МБ — то, что энтузиасты загружают на сторонние сайты-интерпретаторы: на YFull для разбора Y-хромосомы и митохондриальной ДНК, на медицинские — для поиска вариантов, связанных со здоровьем и метаболизмом.

Рейтинг ценности файлов

Файлы неравноценны. Если из одного восстанавливается другой — беречь надо источник, а не производное. Главный источник — риды: только их сняли непосредственно с вашей ДНК, всё остальное из них выводится (были бы те же программы и тот же референс). По убыванию значимости:

1. Риды. Единственное, что нельзя получить заново. Лежат либо в FASTQ, либо внутри CRAM — там они не выброшены, а пересжаты, и при необходимости вытаскиваются из CRAM обратно в FASTQ. То есть FASTQ и CRAM во многом дублируют друг друга; если жмёт диск, логично держать что-то одно, и скорее CRAM — он меньше, а из него восстановимы и FASTQ, и VCF.

2. CRAM — риды плюс их места на референсе. Получается из FASTQ повторным выравниванием, но это часы счёта на все 600 миллионов ридов.

3. VCF — итоговый список отличий, ради которого всё и затевалось. Его грузят в сервисы; можно воспроизвести из CRAM можно, весит всего 215 МБ.

4. Индексы .crai и .tbi — маленькие файлы, чистые производные, пересоздаются за секунды.

Контрольные суммы и хранение на годы

Полногеномное секвенирование делается в идеале 1 раз на всю жизнь. И ваш геном вам может срочно понадобиться черее 20 лет. И тут нужно быть уверенным, что данные не испортились.

Проверить файлы на повреждения можно — контрольные суммы есть внутри. У .fq.gz это обычная gzip-сумма (CRC32) в конце потока; у .vcf.gz — CRC32 в каждом bgzip-блоке; CRAM считает CRC32 на свои блоки и контейнеры и вдобавок хранит MD5 референса, чтобы файл не раскодировали не тем эталоном. Тестируется в одну команду:

gzip -t  NG101J7C24_1.fq.gz                            # gzip-целостность FASTQ
bgzip -t NG101J7C24.mm2.sortdup.bqsr.hc.vcf.gz         # bgzip-блоки VCF
samtools quickcheck NG101J7C24.mm2.sortdup.bqsr.cram   # обрыв или порча CRAM

Если ваша секвенирующая лаба приложила к выгрузке список md5 — сверьтесь и по нему; сам список хешей тоже стоит положить в бакап, иначе потом будет не с чем сравнивать. Но контрольная сумма лишь ловит порчу, а не чинит её. Геном может понадобиться и через десять лет, а на таком сроке один диск — это не хранение. Что помогает:

• Минимум три копии в разных местах. У меня одна живёт на серверах Nebula, вторая — на ноуте, третья — на ПК.

• Желательна любая файловая система с контрольными суммами (Btrfs или ZFS).

• Холодная копия — желательно на магнитном диске, а не на SSD (это рекомендация, не требование). В SSD биты держатся зарядом в ячейках, и без питания заряд утекает: по стандарту JEDEC потребительский диск обязан хранить данные без питания всего около года (и тем меньше, чем теплее в шкафу и сильнее изношен диск). HDD держит магнитную запись куда дольше — порядка 10–15 лет, в хороших условиях и до двадцати (магнитное поле слабеет примерно на 1% в год). Про такой холодный архив легко забыть на годы, а достать его может понадобиться спустя десятилетие — под этот сценарий магнитная пластина надёжнее флеш-памяти.

Что со всем этим делать дальше

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

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