Алгоритмы коррекции неравномерности освещения в изображениях системы параллельного секвенирования ДНК

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

Существуют различные публикации, в которых обсуждаются причины и отдельные составляющие ошибок, вызванных неравномерностью освещения при измерениях сигналов флуоресценции [1, 2, 3]. В работе [1] показано, что флуоресцентное излучение обычно пропорционально интенсивности освещающего света (за исключением случаев, когда происходит истощение основного состояния флуорофора). Поэтому, если однородный флуоресцентный образец освещается неравномерно, результирующая флюоресценция обычно также будет неравномерной. Неравномерное освещение может быть крайне вредным для количественных измерений, поскольку оно может привести к тому, что интенсивность объекта в одной области поля зрения будет измеряться иначе, чем интенсивность объекта с равной концентрацией флуорофора в другой области поля зрения. Из-за присущих процессу формирования изображения несовершенств микроскопические изображения часто искажаются изменениями интенсивности, проявляющимися в виде больших площадей градиентов интенсивности, отсутствующих в неискаженных изображениях [5]. Это явление еще называют затенением или неоднородностью интенсивности [6], а также называют виньетированием, т. е. ослаблением яркости вдали от центра изображения, что часто приводит к тому, что внешние края изображения оказываются значительно темнее, чем центр [7, 8]. Неравномерное освещение может быть вызвано неточной установкой реакционной ячейки, неравномерность фонового освещения, неравномерностью пространственной чувствительности видеокамеры или ее нелинейности [6]. Согласно работе [8] несколько механизмов могут быть ответственны за эффекты виньетирования. Некоторые из них возникают из-за оптических свойств объективов камеры, наиболее известным из которых является закон падения освещенности по оси. Этот вклад в виньетирование усиливается при съемке под большими углами отклонения от оптической оси [1]. Другие источники виньетирования имеют геометрическую природу. Например, свет, падающий под косыми углами к оптической оси, может быть частично заблокирован диафрагмой или оправой линзы [8]. Виньетирование может быть обусловлено и другими факторами, от освещения, конструкции светового пути между камерой и микроскопом или поведения устройства формирования изображений [1, 8]. Компенсация затенения изображения осуществляется как аппаратными, так и программными методами [10].

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

Данная публикация является развитием работ по разработкам алгоритмов первичной обработки изображений сигналов флуоресценции для программного обеспечения системы параллельного секвенирования «Нанофор СПС» [11 – 12]. В этих работах и в ссылках в них дается описание принципа действия системы параллельного секвенирования. Для уточнения терминологии, применяемой в следующих подпараграфах данной статьи, приведем определение нескольких терминов и дадим краткое описание эксперимента, проводимого в секвенаторе. Реакционная ячейка, входящая в состав секвенатора «Нанофор СПС», представляет собой микрофлюидное одноканальное устройство, в которое вводится исследуемое вещество и химические реагенты, участвующие в химических реакциях, происходящих в процессе эксперимента. В секвенаторе «Нанофор СПС» имеются четыре видеоканала: A, C, G, T, по числу нуклеотидов, входящих в состав ДНК. Под действием лазерного излучения возбуждается флуорофор, длина волны которого соответствует определенной букве. В процессе секвенирования генома, содержащего несколько десятков миллионов пар нуклеотидов, снимается серия из нескольких десятков изображений по всей длине реакционной ячейки. Каждое такое изображение называется полем зрения. Съемка изображений на всей длине реакционной ячейки называется циклом или сканом. В нижеследующем тексте номер цикла обозначен буквой с. После завершения текущего цикла производится переход к следующей стадии. В этой стадии через реакционную ячейку пропускают реагенты, отщепляющие флуорофор и останавливающие процесс синтеза. Затем в реакционную ячейку подаются другие реагенты, с помощью которых начинается новый процесс синтеза – новый цикл ( с = с+ 1 ) , где к каждому иммобилизированному кластеру нуклеиновых кислот добавляются новые нуклеотиды. Кластер – группа идентичных ДНК-молекул, прикрепленных к поверхности реакционной ячейки и амплифицированных в одном месте. В секвенаторе «Нанофор СПС» проводится от нескольких десятков до 250 циклов в зависимости от анализируемой пробы и поставленных задач эксперимента.

1.    Иллюстрации неравномерного освещения в секвенаторе «Нанофор СПС»

На рис. 1 показано изображение поля зрения из эксперимента на секвенаторе «Нанофор СПС». Центральная область изображения более яркая, а края более темные. Исходный размер изображения имел 1950 строк и 2000 столбцов.

Рис. 1. Типичное изображение поля зрения, демонстрирующее неравномерное освещение

На рис. 2 а показан результат обработки изображения, представленного на рис. 1. Изображение, показанное на рис. 1, было обработано с помощью морфологической операции эрозии со структурным элементом disk с радиусом 20 пикселей, а затем сглажено с помощью фильтра Гаусса со средним квадратичным отклонением (сигма), равным 55 пикселей. На рис. 2 б показано трехмерное изображение после обработки исходного изображения с помощью следующих процедур. Морфологическая обработка и сглаживание с помощью фильтра нижних частот позволили выделить из исходного изображения фоновую составляющую.

На рис. 3 а линиями показаны строки изображения, профили которых показаны на рис. 3 б . На этом рисунке показаны профили строк с номерами: 100, 1000 и 1850. Как видно из рис. 3 б , профили строк изображения имеют разную форму и амплитуду, что подтверждает результат неравномерного освещения. Профиль самой яркой линии соответствует строке 1000 и находится вверху рисунка. Профиль самой тусклой линии соответствует строке 1850 и находится внизу рисунка.

2.    Методы коррекции неравномерности освещения

Существует несколько различных подходов к коррекции неравномерного освещения. Методы коррекции могут быть перспективными, когда изображение фоновой составляющей, необходимое для калибровки, получено заранее, до проведения эксперимента, или ретроспективными, когда изображение фоновой составляющей, необходимое для калибровки, получается в ходе эксперимента, в режиме реального времени [13]. В настоящей работе прежде всего приводятся результаты исследования возможностей ретроспективных методов, так они могут быть использованы непосредственно в процессе выполнения эксперимента в режиме реального времени. О перспективных методах коррекции неравномерности освещения в изображениях, получаемых в ходе секвенирования на приборе «Нанофор СПС», будет сказано в конце настоящей работы.

Номер пикселя

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

а)

б)

Рис. 3. (а) Результат обработки изображения, представленного на рис. 1, (б) профили обработанных сигналов

Существуют различные ретроспективные методы. Большинство существующих методов коррекции предполагают, что поля смещения в яркостях пикселей изображений являются мультипликативными и медленно меняющимися [14].

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

Второй класс алгоритмов использует методы поверхностного выравнивания [15], требующие выбора ряда точек на фоне вручную или автоматически, а фон получается путем выравнивания параметрической поверхности [14]. Полиномиальный метод основан на предположении, что изменение интенсивности фонового изображения может быть получено путем подбора полиномиальной функции к значениям интенсивности ряда точек, выбранных на фоне изображения [1]. Он аппроксимирует изображение многочленом и использует отношение ортогональности многочленов Лежандра для расширения изображения как двойной суммы этих функций. Затем сумма оценивается для получения изображения, которое аппроксимирует проекцию на пространство многочленных изображений [1].

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

3.    Линейная модель формирования изображения

Распространенная линейная модель формирования изображения [1, 3, 4], описывающая связь между истинным изображением U(x; у) и полученным изображением N(x; у), выглядит следующим образом :

^N(x; у) = ^U(x; y)S_M^(x; у) + S_A^(x; у) ,                                                                  (1)

N(x; у) - полученное изображение, искаженное неравномерностью освещения;

U(x; у) - истинное изображение;

S_M(x; у) - мультипликативный компонент затенения;

S_A(x; У) — аддитивный компонент затенения .

Коррекция затенения связана с нахождением исправленного изображения U(x, у), которое восстанавливает истинное изображение U(x, у) из полученного изображения N(x, у) [3]

shading correction

N(x,у)----------------> U(x,у) ^ U(x,у) .

Изображение с коррекцией затенения U(x, у) можно рассчитать, использовав формулу (1) модели формирования изображения:

U(x> у) =

N(x,y)-S_A(x,y) $ м (*,У)    ,

где S_A(x, у) и S_M(x,у) - оценки аддитивной и мультипликативной компоненты затенения. Таким образом, проблему коррекции затенения можно рассматривать как проблему оценки аддитивной и мультипликативной компоненты затенения [3].

4.    Оценка аддитивной компоненты затенения

Для ретроспективной оценки соответствующего компонента аддитивного затенения S_A(x,у) исходное изображение обрабатывалось скользящим окном в виде квадрата сначала размером 11×11 и локальный фон был определен как среднее значение самых тусклых 5% пикселей этого пиксельного окна. Затем полученное среднее значение присваивалось всем пикселям окна. В результате для каждого поля зрения формировалось изображение фона, которое вычиталось из исходного изображения. В следующих анализах размер квадрата пиксельного окна тщательно изучался и его размеры составляли следующие значения: 13×13, 15×15, 17×17, 21×21 и 23×23. В качестве критерия сравнения различных размеров пиксельного окна для выявления наилучшей оценки компонента аддитивного затенения S_A(x, у) использовалось количество фрагментов нуклеиновых кислот (картированных ридов), которое было выбрано программой, осуществляющей сборку целого генома из отдельных его частей. Чем больше таких фрагментов, тем меньше ошибок в результатах построения генома. Перед началом такого сравнения были обработаны исходные файлы для 60 циклов в эксперименте по анализу последовательности нуклеотидов в стандартном геноме бактериофага Phix174. В исследуемые исходные файлы сначала не вносилась поправка на неравномерное освещение. Количество фрагментов нуклеиновых кислот, которое было выбрано программой, осуществляющей сборку целого генома из отдельных его частей, оказалось равным 147075. После этого осуществлялась оценка аддитивного компонента затенения и его вычитание из исходного, необработанного изображения.

На рис. 4 приведены результаты выявления количества фрагментов нуклеиновых кислот (картированных ридов), которое было выбрано программой, осуществляющей сборку целого генома из отдельных его частей исследования в зависимости от размеров пиксельного окна для оценки компонента аддитивного затенения S_A(x,у), а также от применяемого алгоритма. Кроме алгоритма со скользящим окном, применялись морфологические алгоритмы, реализующие операции эрозии и дилатации со структурными элементами rectangle и disk разных размеров: 11×11, 13×13, 15×15, 17×17, 19×19, 21×21 и 23×23.

На рисунке показано, что наибольшее количество картированных ридов при коррекции аддитивной составляющей затемнения имеет алгоритм, реализующий морфологическую обработку со структурным элементом rectangle размером 15×15. Для этого алгоритма количество картированных ридов оказалось равным 186249. Для алгоритма со скользящим окном наилучший результат оказался для окна 17×17, который составил 185477. Анализ чисел, стоящих на горизонтальной оси рисунка, показывает, что все рассмотренные алгоритмы дают близкие результаты, отличающиеся друг от друга примерно 1 %.

Таким образом, введение аддитивного компонента в коррекцию исходных файлов позволило увеличить количество картированных ридов до 186249, что больше чем на 25 % по сравнению с количеством картированных ридов для изображений, которые не подверглись коррекции (147075).

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

5.    Оценка мультипликативной составляющей затенения

Для ретроспективной оценки мультипликативного компонента затенения S _M(x, у) могут применяться различные методы, такие как фильтрация Гаусса [1, 16], гомоморфная фильтрация [17, 18] и другие.

Для обработки изображений секвенатора «Нанофор СПС» использовались три метода фильтрации нижних частот, основанных на предположении, что неравномерное освещение является сигналом низкой частоты. Поэтому фильтр нижних частот можно использовать для его извлечения из изображения. Интересующие объекты, т. е. сигналы флуоресценции кластеров ДНК, больше, чем вариация фона, а фон имеет другую интенсивность нежели флуоресценция кластеров. Размытие изображения с помощью фильтра нижних частот основано на том, что все измерения интенсивности являются смесью сигнала и фона. Полученное сглаженное изображение считается оценкой фона изображения [16]. Такая фильтрация может быть выполнена путем свертки исходного изображения N (x, у) с ядром фильтра F_m . В частности, для фильтра Гаусса такое ядро G(x, у) определяется как:

^G ° =^^exP(^-X-^,                                                          (4)

где параметр σ определяет степень размытия изображения. Кроме матрицы с ядром Гаусса, использовались матрицы average и disk разных размеров. С помощью матрицы типа average создается усредняющий фильтр нижних частот. Все компоненты этой матрицы – усредняющие коэффициенты являются одинаковыми числами, значения которых зависят от размера матрицы. С помощью матрицы типа disk создается круговой усредняющий фильтр нижних частот внутри квадратной матрицы размером равным удвоенному радиусу. В отличие матрицы типа average в матрице disk усредняющие коэффициенты находятся внутри круга диаметром равным размеру матрицы, а остальные компоненты этой матрицы равны нулю. Эффект свертки с фильтрующей матрицей заключается в сглаживании значений интенсивности. Цель состоит в том, чтобы сгладить изображение до тех пор, пока не будут сглажены исследуемые объекты флуоресценции. При этом полученное сглаженное изображение будет соответствовать фону.

6.    Влияние различных фильтрующих матриц на количество картированных ридов

Оценка S_M(x, у) мультипликативного компонента затенения осуществлялось путем выполнения операции свертки выражения N(x,у) — S_A(x,у) с ядром фильтра F_m:

S m (x, у) = (N(x, y)—S_A(x, y))*F m ,                                                               (5)

где знак * – оператор вычисления свертки.

Оценка откорректированного от затенения изображения I^s (x, у, с) для каждого цикла вычислялась по формуле (6), но с умножением на масштабирующий коэффициент W(c), зависящий от номера цикла с. Коэффициент W(c) вычислялся как среднее значение из максимальных интенсивностей в мультипликативных компонентах затемнения S_M(x, у, с) для каждого цикла. Таким образом, оценка изображения I^s (x, у, с), откорректированного от затенения для каждого цикла с , вычислялась по формуле:

Рис. 5. Количество картированных ридов в зависимости от типа фильтрующей матрицы, ее размера и значений а для фильтрации Гаусса при коррекции от затенения изображения с учетом адаптивного и мультипликативного компонентов затемнения §А(х, у) и SM(x, у)

7.    Результаты применения алгоритмов

Рассмотрим результаты применения алгоритмов, описанных в параграфах 4–6, для обработки данных эксперимента.

На рис. 6 показаны графики профилей 1000-й строки изображения поля зрения С эксперимента на секвенаторе «Нанофор СПС» для исходного файла (верхний график) и для файла после корректировки неравномерности освещения. На рисунке показано, что профиль строки в исходном файле находится на изогнутой вверх линии, тогда как нижний график находится очень близко к «нулевой» линии.

0            500          1000          1500          2000

Номер пикселя

Рис. 6. Графики профилей 1000 строки изображения для исходного файла (верхний график) и для файла после корректировки неравномерности освещения (нижний график) по горизонтальной оси номер пикселя. По вертикальной оси интенсивность флуоресценции в условных единицах

На рис. 7 показаны профили строк с номерами: 100, 1000 и 1850 для откорректированного изображения после того, как к нему были применены те же процедуры обработки, аналогичные показанным на рис. 3 б . Профили строк изображения близки к прямым линиям. Средние значения профилей строк 10, 1000 и 1850 близки к нулю и равны соответственно 0,005; 0,001 и 0,03, а их средние квадратичные отклонения равны 0,01; 0,004 и 0,04. Профиль самой тусклой линии соответствует строке 1850, находится внизу рисунка и немного отличается от прямой линии.

На рис. 8 показано трехмерное представление результата обработки изображения, откорректированного на влияние неравномерного освещения. Трехмерное представление имеет форму, близкую к плоскости, и eе среднее значение близко к нулю и равно 0,04. Это обстоятельство свидетельствует о том, что в откорректированном на влияние неравномерного освещения изображении удален фон.

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

Рис. 7. Профили строк с номерами 100, 1000 и 1850 для откорректированного изображения

Рис. 8. Трехмерное представление фона, выделенного из откорректированного на влияние неравномерного освещения изображения

Рис. 9. Трехмерное представление фона, выделенного из откорректированного на влияние неравномерного освещения изображения по программе, реализующей алгоритм BaSiC

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

По программе, реализующей алгоритм BaSiC, были обработаны те же данные эксперимента по построению генома, что и данные для проверки рассмотренных выше алгоритмов. В результате этой обработки было полу- чено 185083 картированных рида. При использовании рассмотренных выше алгоритмов количество картированных ридов оказывается примерно на 5% больше указанной величины, что важно для прецизионного генетического анализа. Одна из причин меньшего количества картированных ридов следующая: фон, выделенный из изображений, откорректированных по алгоритму BaSiC, близок к плоскости, но среднеквадратичное отклонение шума на его поверхности больше чем в 2 раза выше, чем среднеквадратичное отклонение шума на поверхности плоскости фона, выделенной по алгоритмам, рассмотренным выше. Большее значение шума на плоскости фона подтверждается сравнением его трехмерного изображения, представленного на рис. 9, с трехмерным изображением фона, представленным на рис. 8. Другая причина получения меньшего количества картированных ридов заключается в том, что некоторые объекты флуоресценции после коррекции неравномерности освещения оказались смещенными на 1–2 пикселя, что привело к появлению дополнительных ошибок в построении ридов.

8.    Перспективные методы коррекции неравномерного освещения в секвенаторе «Нанофор СПС»

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

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

Оценка откорректированного на влияние неравномерного освещения для перспективных методов вычислялась по формуле (7), аналогичной ретроспективным методам, только в качестве мультипликативной компоненты затемнения S_M (x, У) использовались значения интенсивностей из априорных изображений фоновых поверхностей для каждого канала, а в качестве весового коэффициента W использовались либо средние значения интенсивностей априорных изображений фоновых поверхностей для каждого канала, либо интенсивности, близкие к максимальным (96-процентный перцентиль). Такая операция использовалась для исключения аномальных результатов измерения флуоресцентных сигналов априорных изображение фоновой интенсивностей каждого из каналов (A, C, G, T).

Таким образом, формула для оценки интенсивностей, откорректированных на влияние неравномерного освещения, модифицированная для перспективных методов, корректировки, имеет вид:

U(x, у, с)

N(x,y,c)S_A(x,y,c) ^ S M (x,y)

В этой формуле, в отличие от формулы (6) мультипликативный компонент затемнения S_M(x, у) и весовой коэффициент W не зависят от номера цикла C .

Проверка возможностей использования априорных изображений фоновой составляющей для коррекции неравномерностей освещения показала, что количество фрагментов нуклеиновых кислот (картированных ридов) оказалось равным 178744. Что примерно на 10% меньше, чем при коррекции ретроспективными методами, описанными выше.

Выводы

• 1.    Использование перспективных и ретроспективных методов коррекции неравномерности освещения в изображениях системы параллельного секвенирования «Нанофор СПС» позволило увеличить количество фрагментов нуклеиновых кислот (картированных ридов), которое было выбрано программой, осуществляющей сборку целого генома из отдельных его частей более чем на 25 %, а в случае использования ретроспективных методов с наилучшими параметрами – до 33 %.

• 2.    При оценке аддитивного компонента в алгоритме коррекции неравномерности освещения в изображениях секвенатора выявлено, что наилучшие результаты дает алгоритм с использованием морфологической обработки со структурным элементом average размером 15×15 пикселей для построения изображения фона.

• 3.    При оценке мультипликативного компонента в алгоритме коррекции неравномерности освещения в изображениях секвенатора выявлено, что наилучшие результаты дает алгоритм фильтрации нижних частот на основе функции Гаусса со средним квадратичным отклонением 65 пикселей и размером фильтрующей матрицы 125×125 пикселей.

• 4.    Ретроспективные методы коррекции дают примерно на 10% лучший результат по коррекции изображений, искаженных неравномерностью освещения по сравнению с перспективными методами.

Работа выполнена при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию № 075–00444–25–00 (от 26.12.2024).