Моделирование первичных данных одномолекулярного секвенатора. Ч. 1. Симуляция детектирования излучения волновода нулевой моды

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

Разработки компаний Pacific Biosciences (PacBio) и Oxford Nanopore Technologies (ONT)

представляют собой две основные технологии секвенирования третьего поколения, по которым сегодня развивается одномолекулярное секвенирование. Первая технология, названная SMRT (Single Molecule Real Time sequencing), вначале обеспечивала получение непрерывных длинных чтений (continuous longreads, или CLR) с уровнем ошибок 15%, а затем получила развитие в методе Hi-Fi ридов, снижающем уровень ошибок до величины < 1%. Вторая технология дает возможность получать ультрадлинные риды до 4 Mbp, но с высоким уровнем ошибок, которые можно снизить почти до уровня метода Hi-Fi, используя технологии "дуплексного" секвенирования [1]. Обе технологии послужили для завершающего этапа расшифровки генома человека — секвенирования

Y-хромосомы, имеющей сложную повторяющуюся структуру [2].

В открытом доступе существует несколько симуляторов длинных последовательностей, имитирующих данные секвенаторов третьего поколения. Это, например, программы SimLoRD, NPBSS, PBSIM, PaSS [3–6]. Они позволяют получить последовательности длинных ридов (longreads) в формате FASTA-файлов, моделируя различные ошибки чтения, что дает возможность использовать эти данные в различных задачах биоинформатики для отладки алгоритмов вторичной обработки. Такие симуляторы дают высокую степень приближения к результатам, полученным при реальном секвенировании. Однако они не позволяют получить первичные, или сырые (raw), данные, т.е. изображения (movies) и трейсы (traces) для SMRT-секвенаторов. От них отличается симулятор Nanopore SimulatION, который специально предложен как симулятор первичных данных для технологии ONT-сигналов ионного тока [7].

Симуляторы первичных данных для секвенаторов, регистрирующих сигналы флуоресценции, отсутствуют в публичном доступе. PacBio в своих внутренних тестах, очевидно, использует собственные генераторы симуляций, которые не публикуются в открытом виде. Это представляет проблему в случае необходимости имитации первичных данных секвенаторов PacBio, например, для разработки и оптимизации собственных алгоритмов первичной обработки данных типа movie и trace для оптимизации процедуры поиска баз (base-calling).

Таким образом, для решения задач распознавания флуоресцирующих объектов на исходных изображениях и дальнейшей интерпретации сигналов в методе одномолекулярного секвенирования ввиду недостатка и дороговизны получения экспериментальных данных и отсутствия симуляторов становятся актуальными задачи имитационного моделирования и программирования.

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

СИСТЕМА ДЕТЕКТИРОВАНИЯ ОДНОМОЛЕКУЛЯРНОГО СЕКВЕНАТОРА

Систему детектирования при одномолекулярном секвенировании по технологии, близкой к PacBio [8–12], можно представить как специализированный четырехканальный, по числу спектральных диапазонов, флюоресцентный микроскоп с цифровой регистрацией изображений флуоресцирующих объектов ВНМ1 (волноводов нулевой моды). Каждый спектральный диапазон соответствует одному из четырех флуоресцентно меченных нуклеотидов, традиционно обозначаемых A, T, G и C. ВНМ представляет собой наноструктуру, в которой протекает реакция секвенирования. Это цилиндрический наноколодец с глубиной и шириной около 100 нм, сформированный в тонкой пленке алюминия, объем которого не превышает 10–19 л. В каждом ВНМ в процессе встраивания полимеразой очередного меченого нуклеотида из реакционного буферного раствора происходит эмиссия флуоресценции соответствующей метки под действием возбуждающего лазерного излучения, падающего на прозрачную поверхность измерительной ячейки снизу. Ячейка содержит большое число ВНМ и находится в поле зрения микрообъектива высокого разрешения с целью одновременного наблюдения флуоресценции во всех наноструктурах.

Выходящие из объектива пучки света несут информацию о пространственном распределении ВНМ и спектральном составе излучения. Для разделения спектральных диапазонов регистрируемой флуоресценции световые пучки проходят через дисперсионную призму [13]. При построении изображения излучение от каждого ВНМ будет разделяться в направлении дисперсии призмы в соответствии со спектральным составом флуоресценции четырех излучающих меток (рис. 1).

Свет попадает на фотосенсоры видеокамеры, где в течение заданного времени экспозиции генерируется изображение. Для отображения всех ВНМ используется одна крупноформатная высокочувствительнаяs CMOS-видеокамера [14]. На фотосенсорах видеокамеры для каждого ВНМ можно выделить четыре смежных группы пикселей, в которых сигналы генерируются как световые пятна в соответствии с интенсивностью и спектральным составом попавшего на них света.

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

Рис. 1. Соответствие каналов детектирования дисперсии света на призме.

1 — индуцированное лазером излучение от ВНМ; 2 — дисперсионная призма;

3 — дисперсия света на сенсорах видеокамеры; 4 — сформированное изображение ВНМ в четырех спектральных каналах. Используется участок спектра 550–750 нм

Размеры изображения в пикселях зависят от дисперсии призмы, выбранного увеличения и физического размера сенсоров видеокамеры. Изображения световых пятен для разных нуклеотидов частично перекрываются в соответствии с перекрытием спектров флуоресцентных меток.

Полученное для каждого кадра изображение передается в компьютер. Серия таких изображений формирует набор сырых данных movies, которые подвергаются дальнейшей обработке. В процессе первичной обработки в компьютере проводится локализация изображений всех флуоресци- рующих ВНМ на общем изображении кадра и осуществляется трассировка сигналов флуоресценции для каждого спектрального канала, в результате чего формируются трейсы — наборы данных, где хранится информация об изменении интенсивности флуоресценции во времени. После предварительной фильтрации трейсов по выбранным критериям качества проводится basecalling-обнаружение и интерпретация импульсов в каждом спектральном канале, что в итоге обеспечивает формирование нуклеотидной последовательности (рис. 2).

Рис. 2. Обобщенная схема детектирования флуоресценции при одномолекулярном секвенировании с использованием дисперсионной призмы.

МОДЕЛИРОВАНИЕ БАЗОВЫХ ИЗОБРАЖЕНИЙ И СИГНАЛОВ ДЛЯ ОДНОГО ВНМ

Принцип формирования базового модельного изображения

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

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

Область симуляции светового пятна для одного спектрального диапазона представляет собой квадратное изображение, размер которого может быть установлен выборочно от 2 × 2 до 5 × 5 пикселей. Один пиксель такого базового изображения соответствует одному сенсору видеокамеры. Уровень яркости в одном пикселе соответствует числу фотоэлектронов, накопленных в соответствующем сенсоре за время однократной экспозиции. Сумма яркостей всех пикселей в области симуляции формирует интегральный базовый сигнал, который становится единичным отсчетом трейса в одном спектральном канале (рис. 3).

На рис. 3:

– Лазер 1 индуцирует флуоресценцию в ВНМ 2 с иммобилизованной на дне полимеразой, фрагментом цепи ДНК и меченными нуклеотидами в растворе.

– В процессе встраивания флуоресцентно-меченного нуклеотида из ВНМ под действием лазерного излучения испускается свет, который, попадая на сенсоры видеокамеры 3, в каждом формирует эквивалентный электрический сигнал в фотоэлектронах (ф. эл.).

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

– Суммарный сигнал всех пикселей базового изображения формирует интегральный базовый сигнал, соответствующий одному отсчету на трей-се 5 в моделируемом спектральном канале.

В качестве примера реального трейса на рис. 4 показан фрагмент длительностью 30 с. Данные получены при регистрации флуоресценции в канале G для одиночной молекулы стрептавидина, иммобилизованной на дне ВНМ.

Рис. 3. Формирование базового изображения и базового сигнала ВНМ в одном спектральном канале.

1 — лазер, индуцирующий флуоресценцию; 2 — волновод нулевой моды (ВНМ); 3 — сенсоры видеокамеры; 4 — пиксели возможных моделей базовых изображений; 5 — отсчеты трейса моделируемого спектрального канала

s

Рис. 4. Фрагмент трейса

Молекула стрептавидина помечена красителем Cy5. В ВНМ добавлены меченные этим же красителем нуклеотиды, свободно плавающие в буферном растворе. На отметке 6.5 с включили красный лазер и регистрировали выгорание красителя на стрептавидине, после 15 с добавили меченые нуклеотиды. Интегральный сигнал интенсивности свечения ВНМ определялся на участке изображения размером 3 × 3 пикселя. Интенсивность фоно- вого свечения ячейки после включения лазера составляет примерно 1200 фотоэлектронов, сигналы флуоресценции в ВНМ (базовые сигналы, или отсчеты трейса) превышают по амплитуде 200–400 фотоэлектронов.

Имитация фона, полезного сигнала и темнового шума

Если рассматривать базовый сигнал с точки зрения того, произошло или не произошло событие встраивания нуклеотида, базовое изображение может быть сформировано как:

• 1.    Симуляция наличия темнового тока в сенсорах видеокамеры в случае, если свет от ВНМ не излучается вообще (ситуация "a").

• 2.    Симуляция сигнала фонового излучения для случая, когда встраивания нуклеотида в данном спектральном диапазоне в течение времени текущей экспозиции не произошло (ситуация "b").

• 3.    Симуляция светового пятна от флуоресцентной вспышки, в том числе и сигнала фона, для имитации процесса встраивания нуклеотида (ситуация "c").

Эти ситуации представлены на рис. 5.

Рис. 5. Базовые изображения области симуляции размером 5×5 пикселей (а) и профили яркости для них (б).

На б вверху — отсчеты базовых сигналов при имитации темнового шума камеры (ситуация "a"). На б внизу — отсчеты базовых сигналов при имитации фона ячейки (ситуация "b") и событий встраивания нуклеотидов (ситуация "c").

За время симуляции, равное 310 мс, произошло два события встраивания. Каждое событие образовано пятью отсчетами с интенсивностью, превышающей 100 фотоэлектронов (ф.эл.). Отсчеты фоновых сигналов сформированы базовыми сигналами меньшей интенсивности

Рис. 6. Спектры возбуждения и эмиссии флуоресцентных красителей, используемых для мечения нуклеотидов.

Заштрихованы области взаимного спектрального перекрытия между каналами А и Т, T и G, G и C

Перекрестные спектральные помехи

При синтезе комплементарной цепи ДНК два разных нуклеотида не могут одновременно встроиться в первичную цепь, поэтому в одном ВНМ не могут одновременно светиться два флуорофора. Однако в момент свечения какого-либо одного нуклеотида в соседние с ним каналы детекции могут попадать перекрестные спектральные помехи (crosstalk). В модельном сигнале величина такой помехи может быть задана в процентах от интенсивности сигнала флуоресценции в канале, который является источником помехи. Величина помехи добавляется к интегральному сигналу яркости в моделируемом канале и зависит от формы спектров возбуждения и испускания используемых красителей (рис. 6). Помеху из канала A в канал T обозначим как CrTk ^AT и, наоборот, помеху из канала T в канал A обозначим как CrTk ^TA. Логика обозначений для других спектральных каналов аналогична.

Для моделирования разделения лучей от четырех флуоресцентных красителей при помощи дисперсионной призмы следует расположить четыре базовых изображения друг под другом по вертикали в порядке увеличения длины волны спектров эмиссии A (зеленый), T (желтый), G (красный) и C

(темно-красный). Каждое из четырех базовых изображений генерируется отдельно в соответствии с заданными для данного спектрального канала параметрами.

Алгоритм моделирования базового изображения Sbasic

Алгоритм моделирования (вариант 1) базового изображения в области симуляции для канала А (SbasicA) в случае встраивания меченного нуклеотида А (рис. 7):

• 1.    Задать размер области симуляции в пикселях и смещение моделируемого светового пятна внутри области.

• 2.    Задать амплитуду сигнала флуоресценции A в фотоэлектронах — Fl ^A.

• 3.    Задать амплитуду фонового сигнала BkGrd ^A, т .е. величину сигнала собственного свечения измерительной ячейки в отсутствие сигнала флуоресценции от встраиваемого нуклеотида в моделируемом канале.

• 4.    Сформировать двумерные изображения гауссова распределения яркости в области симуляции для фона и флуоресценции согласно заданным параметрам. Суммировать оба изображения.

  

  Рис. 7. Формирование базового изображения Sbasic для каждого из четырех спектральных каналов.

  В канале А моделируется полезный сигнал (наличие флуоресценции). Цифрами показаны соответствующие этапы алгоритма моделирования

  
  

• 5.    Моделировать шум, зависящий от величины сигнала, накопленного в сенсорах за время экспозиции, как сумму PRNU (Photo Response NonUniformity) и дробового шума (Shotgunnoise) — N dp . Для этого рассчитать величину шумов в каж-

• дом пикселе для полученного в п. 4 суммарного изображения.2

• 6.    Моделировать шум, не зависящий от величины сигнала, накопленного в сенсорах, как сумму темнового шума DSNU (Dark Signal NonUniformity) и шума считывания камеры ReadOut — N undp . Для этого рассчитать величину шумов в каждом пикселе с помощью нормального распределения со средним, равным нулю, и СКО, указанным в паспорте видеокамеры для данных видов шумов.

• 7.    Задать смещение интенсивности сигнала относительно нуля для каждого пикселя — OffSet .³

• 8.    Добавить рассчитанные шумы и смещение к изображению фона и флуоресценции, полученному в п. 4.

Алгоритм моделирования (вариант 2) базового изображения в области симуляции для канала Т (SbasicT) (фон + спектральная помеха из А) в случае встраивания меченного нуклеотида А (рис. 7):

• 1.    Задать размер области симуляции в пикселях и смещение моделируемого светового пятна внутри области.

• 2.    Задать спектральную помеху (crosstalk) в канале T в процентах от амплитуды сигнала из канала A, создающего эту помеху — Fl ^A × CrTk ^AT.

• 3.    Задать амплитуду фонового сигнала BkGrd ^T — величину сигнала собственного свечения измерительной ячейки в отсутствие сигнала флуоресценции от встраиваемого нуклеотида в моделируемом канале.

• 4.    Сформировать двумерные изображения гауссова распределения яркости в области симуляции для фона и флуоресценции согласно заданным параметрам. В данном случае собственная флуоресценция в канале Т равна нулю, но формируется паразитный флуоресцентный сигнал от засветки из канала А. Суммировать оба изображения.

• 5.    Моделировать шум, зависящий от величины сигнала, накопленного в сенсорах за время экспозиции, как сумму PRNU (Photo Response NonUniformity) и дробового шума (Shotgunnoise) — N _dp. Для этого рассчитать величину шумов в каждом пикселе для полученного в п. 4 суммарного изображения.

• 6.    Моделировать шум, не зависящий от величины сигнала, накопленного в сенсорах, как сумму темнового шума DSNU (Dark Signal NonUniformity) и шума считывания камеры ReadOut — N undp . Для этого рассчитать величину шумов в каждом пикселе с помощью нормального распределения со средним, равным нулю, и СКО, указанным в паспорте видеокамеры для данных видов шумов.

• 7.    Задать смещение интенсивности сигнала относительно нуля для каждого пикселя — OffSet .

• 8.    Добавить рассчитанные шумы и смещение к изображению фона и флуоресценции, полученному в п. 4.

При необходимости моделировать темновой ток в сенсорах используются только пункты 6 и 7 алгоритма.

Пояснение к п. 4 алгоритма. Из ВНМ на видеокамеру попадает суммарный поток света, порождаемый как сигналом флуоресценции нуклеотида, если произошло встраивание, так и общим фоновым излучением ВНМ, обусловленным свечением поверхности ячейки и слабым сигналом от меченых нуклеотидов, плавающих в растворе. В обоих случаях изображение светового пятна в области симуляции может быть смоделировано двумерным гауссовым распределением яркости (1) с заданным σ ( σ x = σ y = σ ), которое соответствует радиусу светящегося объекта, и заданным случайным смещением центра пятна m x , m y . При этом сам оптический сигнал соответствует величине двумерного гауссового интеграла (рис. 8)

f ( x , y ) =

2 πσ _x σ _y

( x - mf ( y - m y ) 2 2 σ x ²       2 σ y ²

Смещение центра симметрии пятна задается опционально случайными величинами с равномерным распределением в пределах области симуляции и позволяет учесть проявление пространственной нестабильности изображения.

Общая яркость всех пикселей изображения в области симуляции определяется заданной величиной базового сигнала Fl (например, только для флуоресцентного сигнала). Программно базовый сигнал формируется как сумма яркостей всех пикселей в базовом изображении. Для этого каждый пиксель в области симуляции получает яркость A i,j с соответствующим весовым коэффициентом, как показано в (2):

A j =

777       (‘ - m x )2

Fl         2 2

e ^{2 σ}

2 πσ ²

⁽ j - m y ) 2

2 σ ²

В итоге алгоритм формирования модельного двумерного базового изображения S bsc в области симуляции может быть так описан формулами.

• •    Моделируем наличие флуоресценции (пример для канала А):

s A_ ^Fl = ( Fl ^A * f ( x , y ) + BkGrd ^A * f ( x , y ) ) +

+ N dp + N undp + Offset .                         (3)

Рис. 8. Гауссова модель изображения сигнала флуоресценции ВНМ

• •    Моделируем только фон (пример для канала Т):

T_Bg

S bsc =

= ( Fl A * CrTk AT * f ( x , y ) + BkGrd T * f ( x , y ) ) +

⁺ N dp ⁺ N undp ⁺ Offset .                             (4)

• •    Моделируем темновой кадр:

s DC = N undp + Offset .                 (5)

Обоснование способа моделирования шумов

Пояснение к п. 5 алгоритма. Экспериментальное изучение зависимости шумов видеокамеры от величины регистрируемого сигнала позволило определить СКО сигнала яркости изображения измерительной ячейки, интегрированного в области 3 × 3 пикселя. На рис. 9 представлена зависимость относительного СКО (определенного в процентах по отношению к соответствующей величине сигнала) такого сигнала яркости для пустой ячейки, засвеченной зеленым и красным лазерами. Под "хорошими" пикселями здесь подразумевается основное множество пикселей камеры, для ко- торых измеренный разброс сигнала яркости близок к указанному в паспорте видеокамеры значению. Под "плохими" пикселями мы имеем в виду некоторые пиксели с уровнем разброса сигнала, превышающим паспортные значения в 5 и более раз. Попадание таких пикселей в область интегрирования сигналя яркости при построении трейса может увеличить СКО шума фона и затруднить последующий поиск сигналов-событий4.

Данные для построения зависимости получены в результате обработки серии кадров, зарегистрированных в режиме темного поля (при отсутствии внешнего источника света) и при включенных лазерах с последовательно увеличивающимся временем экспозиции от 10 мс до 250 мс.

Только "хорошие" пиксели, как в точке (210; 105)

Есть включения "плохих” пикселей, как в точке(210;111) Аппроксимация для "хороших" пикселей

Рис. 9. Зависимость относительного СКО сигнала яркости (интенсивности) изображения фона ячейки в области интегрирования 3×3 пикселя от величины этого сигнала с вычетом смещения уровня Offset

Это позволило измерить уровни яркости изображения собственного свечения ячейки в широком диапазоне значений. Для корректной обработки данных из величины интегрального сигнала был вычтен уровень смещения, равный 100 фотоэлектронов на пиксель, искусственно добавляемый в видеокамере для коррекции нулевого уровня. С помощью аппроксимации результатов измерений получена функциональная зависимость величины СКО интегрального сигнала N от уровня интенсивности этого интегрального сигнала S в области 3 × 3 пикселя:

N = 128.15 х S ^ ^;. (6)

Полученное выражение используется при моделировании базового изображения по формулам (3), (4) для определения величины шумов, зависящих от уровня сигнала N dp . Для этого вместо величины S в (6) следует подставить сумму сигналов фона и флуоресценции для моделируемого пикселя в области симуляции.

РЕЗУЛЬТАТЫ

Таким образом, на симулированном базовом изображении мы будем иметь дело с асимметричным пятном, заданным двумерным гауссовым распределением сигналов фона и флуоресценции, с добавлением спектральных помех и флуктуациями яркости в каждом пикселе. Интегральная яркость базового изображения — сумма яркости в пикселях в пределах области симуляции — также будет флуктуировать в разных реализациях даже при одинаковых параметрах симуляции, создавая таким образом стохастическую имитационную модель процесса детектирования излучения от ВНМ.

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

Рис. 10. Примеры базовых изображений нуклеотидов. Пояснения см. в тексте

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

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

На рис. 10:

• 1    — область 5 × 5 пикс.; а — флуоресценция в канале A; b — флуоресценция в канале A со спектральным перекрытием канала T; c — флуоресценция в канале G со спектральным перекрытием канала С.

• 2    — область 4 × 4 пикс.; a — флуоресценция в канале Т; b — флуоресценция в канале A.

• 3    — область 3 × 3 пикс.; a — флуоресценция в канале A; b — флуоресценция в канале G.

• 4    — область 2 × 2 пикс.; a — флуоресценция в канале G; b — флуоресценция в канале С.