Метод обработки сигналов нестационарных процессов спектрометрии в задаче секвенирования ДНК с использованием обратных функций

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

Наиболее характерным примером таких спектральных сигналов является совокупность сигналов сиквенс-анализа фрагментов ДНК по Сэнгеру методом капиллярного электрофореза [1, 2]. Фореграмма — зарегистрированные в результате электрофореза данные, состоит из четырех спектров, или четырех цветовых каналов, в единой шкале времени. Основные параметры сигналов — ширина спектральных пиков и интервал времени между выходом пиков соседних (в цепи ДНК) нуклеотидов — одинаковы для всех каналов. Однако при анализе длинных фрагментов ДНК, порядка 1000 нуклеотидов, значения параметров существенно изменяются во времени. Так, в конце рабочего цикла анализа ширина пиков увеличивается, а расстояние между ними уменьшается, и спектральные пики практически сливаются друг с другом.

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

96 обостряющего фильтра [3–5], используемого для разделения слившихся пиков, включается образ искомого полезного пика, параметрами которого являются ширина пика и параметры или математическое выражение формы пика (например, гауссиан). Параметры образа напрямую определяют точность и корректность разделения слившихся пиков сигнала.

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

Для достоверной оценки параметров обрабатываемых спектральных сигналов нестационарных случайных процессов требуется найти и использовать такой метод коррекции спектров, который позволил бы оценивать параметры спектров с такой же точностью, как в условиях стационарных процессов. В работе предлагается один из возможных методов оценки параметров спектральных сигналов нестационарных процессов, который основан на использовании обратных функций [6], [7, pp. 103–120]. Метод апробирован в задачах обработки сигналов в программном обеспечении отечественного секвенатора ДНК.

ПРЕДЛАГАЕМОЕ РЕШЕНИЕ

Преобразование нестационарного спектра в стационарный

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

Преобразование шкалы времени

Задача решается при следующих условиях:

• - Параметр нестационарности , который в результате преобразования требуется сделать постоянным, представляет собой положительный, нену-

• левой интервал времени на всей шкале. Например, такой параметр, как ширина спектрального пика, не может быть отрицательным или нулевым — это не имеет физического смысла.

• - Шкалирующая функция , описывающая изменение параметра по шкале, — медленная и гладкая (см. далее). Например, для спектра в десятки тысяч отсчетов такая функция вполне приемлемо может быть представлена полиномом 5-го порядка, при том что ширина спектральных пиков (параметр нестационарности) обычно имеет порядок 5-30 отсчетов.

Цифровые оценки параметров приведены применительно к задаче капиллярного электрофореза для анализа ДНК.

Ниже использованы следующие обозначения:

t — "физическое" время, определяет шкалу ис

*) ходного, зарегистрированного спектра; )

τ — "виртуальное" время, определяет шкалу "виртуального", идеального спектра, где интервал параметра нестационарности постоянный (constant);

δ — параметр нестационарности на шкале t, например, ширина спектрального пика, интервал времени;

Δ — значение параметра нестационарности на "идеальной", виртуальной шкале τ , постоянный (constant) интервал времени.

Представим себе следующую мысленную модель системы, в которой необходимо решить обсуждаемую задачу.

Входной процесс представляет собой "идеальный" спектр с постоянной шириной спектральных пиков Δ в течение всего времени τ . Однако при его регистрации прибором возникает искажение шкалы времени t выходного зарегистрированного спектра, например, за счет "плавания" частоты задающего генератора, формирующего отсчеты оцифровки выходного спектра. Соответственно меняется в количестве отсчетов ширина спектральных пиков δ ( t ) на этой шкале.

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

Искажение шкалы описывается некоторой шкалирующей функцией t = f ( τ ) (см. рис. 1).

*) Здесь использовано такое имя, чтобы не путать, например, с таким устоявшимся термином, как "реальное" время (on-line). Следует отметить, что обработка спектрометрической информации, как правило, производится как раз в "нереальном" времени (off-line), когда весь спектр уже полностью накоплен и зарегистрирован.

Рис. 1. Пояснения к построению шкалы с помощью обратной функции.

а — шкалирующая и обратная функции, б — зависимость ширины пиков от времени

Ставится задача восстановления исходного идеального спектра по зарегистрированному спектру на шкале t .

Если функция f(τ) строго монотонная, то для нее существует обратная функция τ = f -1(t), которая также является строго монотонной. А из теории обратных функций [6], [7, pp. 103– 120] известно, что f-1 (fT)) = t f (fЛt)) = t. (1)

Таким образом, для того чтобы в нашей модели восстановить, а в реальности создать, "идеальную" шкалу τ , необходимо подвергнуть физическую шкалу t обратному преобразованию функцией f ^-1( t ).

В реальности в нашем распоряжении имеется зарегистрированный физический спектр на шкале t. Этот спектр можно обработать, найти спектральные пики и построить по ним зависимость ширины пиков от времени δ ( t ) (рис. 1).

На рис. 1 величины Δ и δ изображены условноиллюстративно, чтобы можно было что-то разгля- деть, реально же эти величины на 2–3 десятичных порядка меньше длины шкалы. Поэтому их вполне приемлемо использовать как дифференциалы. Тогда из рис. 1 производная τ по t будет dt(t) _ A dt “ 5(t).

Функция f –1(t) рассчитывается интегрированием ее производной f "(')_ T( t)_A^ d t. (3)

Функция 1/ δ ( t ) является ненулевой положительной функцией, и интеграл от нее будет строго монотонным. Поэтому выполняется приведенное выше условие, чтобы функция f ^- ( t ) была строго монотонной.

Функцию f ( τ ) можно вычислить подобным же способом. Однако наиболее простым и надежным (в плане устранения возможного накопления ошибок вычислений при интегрировании) способом является формирование f как зеркальное отражение (или результат вращения на 180 градусов) функции f ^-1 относительно оси τ = t . Дело в том, что требование к точности вычисления этих функций, вообще говоря, вполне умеренное, но точность их соотношения, точность выполнения преобразований (1) должна быть абсолютной.

После создания шкалы τ в нее переносятся спектральные данные из шкалы t . Для этого на шкале τ задаются отсчеты τ i , как правило, с постоянным шагом. Для каждого из этих отсчетов при помощи функции f ^-1( t ) вычисляется соответствующее ему время t j и соответствующее этому времени значение спектральных данных D [ t ( t j )] переносится в отсчет i шкалы τ — D [ τ ( τ_i )]. Необходимо отметить, что время t j почти всегда не будет точно совпадать с временем какого-либо отсчета на шкале t . Поэтому значение D [ t ( t j )] необходимо вычислять при помощи интерполяции между соответствующими соседними отсчетами шкалы t .

Пример преобразования шкалы времени для модели спектра фореграммы, содержащей пики разной ширины

Рассмотренный в предыдущем подразделе способ преобразования временнóй шкалы был сначала проверен с помощью модели фореграммы из 25 000 отсчетов, в которой содержалось 160 пиков гауссовой формы. Пики находились на одинаковом расстоянии, но каждый из них отличался шириной. Изменение ширины пиков от времени происходило по линейному закону.

Следует обратить внимание на то, что в дискретной математике и в программировании элементы массивов данных идентифицируются номерами этих элементов в массиве. Поэтому для простоты на приведенных ниже рисунках шкалы t и τ оцифрованы в отсчетах и так же на рисунках названы. Шкалы ширины пиков оцифрованы в количестве отсчетов. Секундам эти оцифровки будут соответствовать при равенстве расстояния между отсчетами 1 с. Однако экспериментальные данные, использованные в данной работе, были получены на секвенаторе с заданием интервала между отсчетами, равным 0.25 с. То есть в секундах величина этих оцифровок будет в 4 раза меньше.

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

Как видно из рис. 2, г, зависимость ширины пиков от времени в диапазоне от 0 до 20000 представляет собой линию, близкую к прямой, параллельной горизонтальной оси.

а

б

о 5000    10000 15000   20000

номера отсчетов

Рис. 2. Фрагменты модели фореграммы и зависимости ширины пиков.

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

г

Среднее значение ширины пиков на этом участке равно 23.9, а среднее квадратичное отклонение — 0.1. На участке от 20 000 до 25 000 можно заметить небольшое увеличение ширины пиков, связанное с краевыми эффектами. Среднее значение ширины пиков во всем диапазоне модели спектра равно 24.0, а среднее квадратичное отклонение — 0.4.

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

Пример линеаризации шкалы времени для фореграммы, полученной на секвенаторе ДНК

На рис. 3, а, б, представлены фрагменты фореграммы, полученные на секвенаторе "Нано- фор-05": а — фрагмент в начале шкалы, б — фрагмент в конце шкалы. В конце шкалы ширина пиков увеличилась.

На рис. 4 представлена зависимость ширины пиков в количестве отсчетов на полувысоте пика. Сплошной линией показана аппроксимация этой зависимости полиномом 5-й степени. Аппроксимация этой зависимости была выполнена методом наименьших квадратов для построения линии регрессии. Параметры функции аппроксимации являются характеристикой нелинейности временнóй шкалы

На рис. 5 представлены прямая и обратная функции шкалы времени фореграммы. Прямая функция является результатом интегрирования функции аппроксимации ширины пиков.

Рис. 3. Фрагменты фореграммы: а — в начале шкалы, б — в конце шкалы

Номера отсчетов

Рис. 4. Зависимость ширины пиков от номеров от- счета.

Рис. 5 . Прямая t = f ( τ ) и обратная τ = f ^-1( t ) функции шкалы времени фореграммы

Сплошная линия — аппроксимация полиномом

Интеграл от функции аппроксимации времени между пиками характеризует шкалу времени и является t = f ( τ ). Следуя алгоритму, рассмотренному в предыдущем разделе, найдем обратную функцию τ = f ^-1( t ), заменяя зависимость переменной по вертикальной оси от переменной по горизонтальной оси зависимостью переменной по горизонтальной оси от переменной по вертикальной оси.

Вычисление значения D [ t ( t j )] производилось описанным выше способом с использованием для интерполяции D [ t ( t )] между соседними отсчетами полинома 2-го порядка.

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

до 1450 составляла 16 единиц, а в диапазоне от 16 120 до 16 280 — 23 единицы. После проведения преобразования фореграмм с использованием обратной функции ширина пиков в указанных, а также в других диапазонах оказалась одинаковой и равной 10 единицам.

На рис. 6 показаны фрагменты фореграммы, построенные в исходной (старой) и скорректированной (новой) временных шкалах.

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

Фрагмент фореграммы в новой шкале

Рис. 6. Фрагменты фореграммы.

а — в начале старой шкалы, б — в конце старой шкалы, в — в начале новой шкалы, г — в конце новой шкалы

МЕТОД ОБРАБОТКИ СИГНАЛОВ НЕСТАЦИОНАРНЫХ ПРОЦЕССОВ СПЕКТРОМЕТРИИ 101

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

РЕЗУЛЬТАТЫ И ВЫВОДЫ

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

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

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

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

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

Финансирование

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