Исследование погрешностей моделирования сигмоидальной функцией сигналов полимеразной цепной реакции в реальном времени Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2011, том 21, № 4, c. 28-34

= ПРИБОРЫ И ИЗМЕРИТЕЛЬНЫЕ МЕТОДЫ ДЛЯ БИОХИМИИ

УДК 543.426; 543.9

© Ю. В. Белов, А. И. Петров, В. Е. Курочкин

ИССЛЕДОВАНИЕ ПОГРЕШНОСТЕЙ МОДЕЛИРОВАНИЯ СИГМОИДАЛЬНОЙ ФУНКЦИЕЙ СИГНАЛОВ ПОЛИМЕРАЗНОЙ ЦЕПНОЙ РЕАКЦИИ В РЕАЛЬНОМ ВРЕМЕНИ

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

Кл. сл.: ПЦР в реальном времени, ДНК, пороговый цикл, сигмоидальная функция

ВВЕДЕНИЕ

В статьях [1, 2] приведена классификация методов определения значения порогового цикла — величины, позволяющей судить об исходном количестве копий ДНК и сравнивать образцы между собой по результатам ПЦР. Показаны преимущества вычисления пороговых циклов методом моделирования на базе сигмоидальной функции (МСФ) [3, 4]. При рассмотрении погрешностей измерения порогового цикла [2] шумы моделировались последовательностью 50 случайных чисел, соответствующих 50 температурным циклам. Было показано, что такая шумовая последовательность одновременно действует на все участки сигнала ПЦР и вызывает случайное изменение порогового цикла. В настоящей статье для детального изучения особенностей МСФ предложено использовать нормированное изменение интенсивности флуоресценции в течение одного или нескольких циклов (импульсное возмущение).

ИССЛЕДОВАНИЕ ВЛИЯНИЯ НА ПАРАМЕТРЫ МСФ НОРМИРОВАННОГО ИЗМЕНЕНИЯ ИНТЕНСИВНОСТИ ФЛУОРЕСЦЕНЦИИ В ТЕЧЕНИЕ ОДНОГО ЦИКЛА

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

Для реализации нормированного изменения интенсивности использованы аналитическая фор-

ма и методика построения модернизированной МСФ, приведенные в статье [5]:

F = Fm / (1+exp ((Со - C) / k)) +

+fo + f • C + f • C2, (1)

где F — интенсивность флуоресценции, измеряется в относительных единицах флуоресценции (о.е.); Fm — максимальное изменение интенсивности флуоресценции в конце реакции; С — порядковый номер цикла амплификации; С0 — дробный номер цикла, при котором интенсивность флуоресценции достигает половины значения Fm; k — коэффициент, который определяется наклоном кривой; fO — смещение нулевой линии (фон детектора); f — коэффициент, который определяется линейным дрейфом нулевой линии; f — коэффициент, учитывающий квадратичную составляющую дрейфа нулевой линии.

Вычисления параметров МСФ выполнялись с помощью электронной таблицы Excel (табл. 1) в следующей последовательности.

1. В ячейках В11:В60 вычислялась исходная МСФ с параметрами, приведенными в ячейках В4 :В7, как зависимость от номера цикла в ячейках А11: А60 (рис. 1, а).

2. В ячейках С11:С60 вычислялась исходная МСФ с аналогичными параметрами, первоначально расположенными в ячейках С4: С7.

3. В ячейку D11 вносилось нормированное изменение интенсивности флуоресценции в течение первого цикла Fi = 100 о.е. (относительных единиц флуоресценции). В ячейках E11:E60 вычислялись суммы (B11+D11) : (B60+D60) . В ячейках F11:F60 вычислялись квадраты разностей (Е11-С11)2: (Е60-С60)2, а в ячейке С10 — сумма квадратов разностей Sk.

ИССЛЕДОВАНИЕ ПОГРЕШНОСТЕЙ МОДЕЛИРОВАНИЯ. Табл. 1. Электронная таблица для вычисления параметров МСФ

A B C D E F

1 i= 1 15

2 Уровень= 500 500 500

3 Cti = 23,500 23,490 23,515

4 C0i = 23,5 23,465 23,548

5 Fm = 1000 1024,976 976,255

6 k = 1,39 1, 439 1,338

7 fo = 0 25,951 -5,815

8 f1 = 0 -2,571 1,383

9 f2 = 0 0, 032 -0,016

10 Sk = 0 7663,986 8916,995

11 1 9,33E-05 23,414 100 100,0001 5865,656

25 15 2,204 2,57537 0 2,204367 22,845

33 23 411,029 414,080 0 411,029 9,306

34 24 588,971 589,243 0 588,971 0, 074

43 33 998,925 999,496 0 998,9251 0,326

60 50 1000 1002,188 0 1000 4,789

4. После установления условия минимального значения Sk и оптимизации параметров в ячейках С4:С9 по методу наименьших квадратов в ячейках С11:С60 получаются числовые значения оптимизированной модели при номере цикла в ячейке А11 i = 1.

5. Аналогичные действия выполнялись при установке величины Fi = 100 о.е. последовательно в ячейки D12 : D60 (i = 2-50). В ячейках С4 : С9 каждый раз получаются новые параметры МСФ.

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

В графическом виде результаты изменения величины С(ц в зависимости от номера цикла I, в который вносилось нормированное изменение флуоресценции, представлены на рис. 1, б (ряд 1). На этом графике нулевое значение на вертикальной оси соответствует величине С0 = 23.5 исходной модели. Наибольшие изменения величины С00, приблизительно равные 0.17, наблюдаются при положении нормированного изменения флуоресценции в 23-м и 24-м циклах, т. е. в средней части исходной МСФ, график которой приведен на рис. 1, а.

Если в ячейки D11:D60 табл. 1 последовательно вносить величину Fi = -100 о.е., то изображение графика, приведенного на рис. 1, становится зеркальным по отношению к горизонтальной оси.

Рис. 1. Графики зависимостей параметров МСФ от номера температурного цикла.

а — исходная МСФ при ^ = / = ^ = = 0 о.е., Fm = = 1000 о.е.; б — зависимости величин С0,- и Сй от номера цикла I внесения нормированного изменения флуоресценции

10 20 30 40 Номера циклов

50

60

О

О <и

0.2 0.15 0.1 0.05 0

■0.05 -0.1 -0.15 -0.2

■

(ЦТ ю 20 . 1 30 . 40 50

-Ряд1 Ряд2 РядЗ Ряд4

Номера циклов

а

0

б

Несколько точек этого графика, соответствующие 15-му, 23-му и 33-му циклам, приведены на рис. 1, б (ряд 2).

ОПРЕДЕЛЕНИЕ ПОРОГОВОГО ЦИКЛА МЕТОДОМ ПОРОГА

Для реализации метода порога в ячейку С3 табл. 1 добавлена величина порогового цикла С^. Величина С^ определяется на уровне 500 о.е., приведенном в ячейке C2 (оптимальный уровень порога [2]). При этих вычислениях использованы параметры МСФ в соответствии с формулой

С2 =

=$С$7+$С$5/ (1+ЕХР ( ($С$4-С3) /$С$6) ) + +$С$8*С3+$С$9*(С3)Л2. (2)

Применен инструмент "Поиск решения" из меню "Сервис"; при этом выполнялись следующие действия:

"Установить целевую ячейку С2", "Равной значению 5 00", "Выполнить", "Изменяя ячейки С3".

График изменения величины Сй в зависимости от номера цикла i, в который вносилось нормированное изменение флуоресценции ^ = 100 о.е., представлен на рис. 1, б (ряд 3). Несколько точек графика при величине ^ = -100 о.е. приведены на рис. 1, б (ряд 4). При сравнении графиков 3 и 1 рядов на рис. 1, б, можно отметить, что изменения величин Сй при внесении нормированного изменения флуоресценции меньше, чем изменения величин Со, при тех же условиях.

23.35 23.40 23.45 23.50 23.55 23.60 23.65

Дробный номер цикла

Рис. 2. Графики зависимости интенсивности модельных сигналов от дробного номера цикла. Пояснения в тексте

ЭФФЕКТ УМЕНЬШЕНИЯ ПОГРЕШНОСТЕЙ

ПРИ ИСПОЛЬЗОВАНИИ МЕТОДА ПОРОГА

С целью сравнения механизмов вычисления пороговых циклов методом порога Ct и методом "средней точки" С0 на рис. 2 приведены графики флуоресценции в диапазоне от 23.4 до 23.6 дробных единиц циклов (д.е.ц.); при этом интенсивность флуоресценции вычислялась с использованием формулы (2).

Средний график соответствует исходной модели с параметрами, приведенными в ячейках В4:В9 табл. 1. Точка 1 соответствует параметрам, приведенным в ячейках В2:В4, при этом номера пороговых циклов равны: С0 = Ct = 23.5. Нижний график построен при условии Fi5 = 100 о.е. (внесение нормированного изменения в величину флуоресценции 15-го цикла) при параметрах, приведенных в ячейках D4 :D9; точка 2 имеет координаты С015 = 23.548 д.е.ц. и F2 = 506.025 о.е., а точка 3 — СД5 = 23.515 д.е.ц. и F3 = 500 о.е. Верхний график построен при условии Fi5 = -100 о.е., точке 4 соответствует С0 = 23.450 д.е.ц., а точке 5 — Ct = 23.485 д.е.ц.

Сравнивая координаты точек 1 и 2, можно отметить, что нижний график смещен относительно среднего исходного графика, нормированного при f0 = 0 и Fm = 1000 о.е., по горизонтальной оси на величину 0.048 д.е.ц. и по вертикальной оси на величину 6.025 о.е.

Выгодно вместо величины СМ5 (точка 2) в качестве порогового цикла использовать величину СД5 (точка 3), поскольку в таком случае не учитывается перемещение графика по вертикальной оси, по-

этому влияние внесенного нормированного изменения флуоресценции на результат определения порогового цикла значительно уменьшается (в приведенном примере более чем в 3 раза). Аналогичный вывод можно сделать при рассмотрении координат точек 4 и 5 верхнего графика.

Усредненные вносимые изменения пороговых циклов можно сравнить путем вычисления значений стандартного отклонения величин Сй и С07 в пределах от 1 до 50 циклов: 5(СЙ) = 0.049 и 5(С0/) = = 0.061 соответственно для метода порога на уровне 500 о.е. и метода "средней точки".

Стандартное отклонение 5 величин Сй и С07 можно рассматривать как основную составляющую погрешностей измерения пороговых циклов, при этом метод порога имеет преимущество по сравнению с методом "средней точки" (5(СЙ) < < 50).

ИССЛЕДОВАНИЕ ВЛИЯНИЯ НА ПАРАМЕТРЫ МСФ НОРМИРОВАННОГО ИЗМЕНЕНИЯ ИНТЕНСИВНОСТИ ФЛУОРЕСЦЕНЦИИ В ТЕЧЕНИЕ НЕСКОЛЬКИХ ЦИКЛОВ

Если при каждых циклах / и (/+1) последовательно внести нормированные изменения флуоресценции Fi = Fi+l = 100 о.е., то с помощью табл. 1 можно вычислить новые параметры МСФ. Результаты вычисления зависимостей изменения величин С07 и Сй от номера цикла / приведены на рис. 3 (соответственно ряды 1 и 2). На этом графике нулевое значение на вертикальной оси соответствует величине С0 = С = 23.5 исходной модели.

О

О <и

<и г

со

-РЯД1 РЯД2 РядЗ

-0.4

Рис. 3. Графики зависимостей величин С0; и Сп от номера цикла / при воздействии нормированного изменения интенсивности флуоресценции в течение двух циклов

Номера циклов

0.08

О

О <и

<и г

(0

Если сравнить графики на рис. 3 и рис. 1, б (ряд1), то можно отметить, что максимальное изменение 23-го цикла с 2 импульсами, равное -0.338 д.е.ц., почти в 2 раза больше, чем максимальное изменение с 1 импульсом, которое равно -0.173 д.е.ц.

На рис. 3 (ряд3) приведен результат сложения со сдвигом на 1 цикл двух одинаковых зависимостей С0/, полученных при воздействии нормированного изменения интенсивности флуоресценции в течение одного цикла. Почти полное совпадение графиков рис. 3 (ряд1 и ряд3) позволяет сделать вывод об аддитивности результатов воздействия нормированных изменений интенсивности флуоресценции в течение нескольких циклов.

Этот вывод с учетом знака подтверждается зависимостями изменения величин С0, и Сй от номера цикла /, когда при каждом /-м цикле вносилось нормированное изменение флуоресценции ^ = = 100 о.е., а при следующем цикле вносилось нор-

мированное изменение флуоресценции = = -100 о.е. Графики этих зависимостей приведены на рис. 4.

Наибольшие изменения величин С0/, приблизительно равные ±0.06 д.е.ц., наблюдаются при / = = 21 и 25 циклам. Эти изменения почти в 3 раза меньше, чем максимальное изменение с одиночными импульсами, равное -0.173 д.е.ц. Изменения величин С0/ по амплитуде почти полностью совпадают с результатами вычитания двух числовых последовательностей изменений С0/ со сдвигом на 1 цикл.

СРАВНЕНИЕ ВЛИЯНИЯ НА ПАРАМЕТРЫ МСФ ШУМА И НОРМИРОВАННОГО ИЗМЕНЕНИЯ ИНТЕНСИВНОСТИ ФЛУОРЕСЦЕНЦИИ

Шумовые сигналы дискретно моделировались 10 последовательностями по 50 случайных чисел с максимальными амплитудами Г/м, равными ±1,

±10 и ±100 о.е. Эти последовательности поочередно вносились в ячейки D12 :D60, затем вычислялись величины С0 и Ct .

Значения стандартного отклонения S позволяют оценить усредненные погрешности измерения пороговых циклов, возникающие при воздействии на МСФ шума и нормированного изменения интенсивности флуоресценции в течение одного цикла. Значения S величин С0 и Ct для каждой серии из 10 последовательностей случайных чисел приведены в табл. 2 (столбцы 2-4). По результатам сравнения величин S, полученных при дополнительной обработке сигналов, можно оптимизировать измерения пороговых циклов.

В столбце 5 табл. 2 приведены значения S величин С0 и Ct для случая предварительной обработки шумовых сигналов с помощью фильтра 3С3 (трехкратное последовательное применение скользящего среднего значения — прямоугольное окно по 3 циклам [1]). В последнем столбце приведены величины S по 50 циклам при нормированном изменении интенсивности флуоресценции в течение одного цикла без фильтра.

Значение S при шуме ±100 о.е. можно сравнить со значением S при нормированном изменении интенсивности флуоресценции в течение одного цикла Fi = 100 о.е. в виде отношения

Е = S(CM) / S(Q>6) = 0.252 / 0.061 = 4.13.

Величину Е можно условно назвать "коэффициентом эффективности шума", поскольку эта величина характеризует увеличение воздействия нормированной шумовой последовательности в течение всех циклов МСФ по отношению к воздействию нормированного изменения интенсивности флуоресценции в течение одного цикла. После фильтрации эффективность шума уменьшается:

Еф = S(C05) / S(C06) = 0.213 / 0.061 = 3.49.

При сравнении полученных результатов видно, что величины приблизительно линейно зависят от уровня шума (столбцы 2-4 табл. 2). Значение 5 величины С0, равное 0.252, после фильтрации уменьшается на 15 % (0.213 / 0.252 = 0.85). Отношение значений 5 величин Сt и С0 меньше единицы на 12 % (0.222 / 0.252 = 0.88).

При одновременном использовании метода порога и фильтрации сигналов влияние шума уменьшается еще больше. Отношение значений 5 величин Сt (с фильтром) и С0 (без фильтра) меньше единицы на 22 %:

5(СЙ) / 5(0,4) = 0.197 / 0.252 = 0.78.

ЗАКЛЮЧЕНИЕ

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

Приведены результаты импульсного возмущения на МСФ в течение нескольких циклов. Показано, что эти результаты могут быть получены путем сложения результатов воздействия единичных импульсных воздействий с учетом знака. Аддитивность результатов позволяет рассматривать шум как эквивалент последовательности случайных импульсов.

Выполнено сравнение влияния на параметры МСФ шума и импульсного воздействия в течение одного цикла. Приведено определение "коэффициента эффективности шума".

Табл. 2. Результаты вычислений стандартного отклонения 5 величин С0 и С в зависимости от вида и величины воздействия

Стандартное отклонение S величин С0, Ct Вид и величина воздействия

Шум FiM = = ±1 о.е. Шум FiM = = ±10 о.е. Шум FiM = = ±100 о.е. Шум FiM = = ±100 о.е., фильтр 1 цикл F, = 100 о.е.

1 2 3 4 5 6

С0 С 0.00185 0.00176 0.019 0.018 0.252 0.222 0.213 0.197 0.061 0.049

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

Показано, что путем использования метода порога и фильтрации сигналов при моделировании сигналов ПЦР-РВ сигмоидальной функцией можно уменьшить влияние шума на погрешность измерения пороговых циклов (в приведенном примере величина стандартного отклонения понизилась на 22 %). Этот результат реализован в программном обеспечении анализатора нуклеиновых кислот АНК-32, который серийно выпускается в ИАП РАН [6].

Работа выполнена при поддержке Министерства образования и науки Российской Федерации в рамках федеральной целевой программы "Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 20072012 годы" при выполнении научно-исследовательской работы по лоту (шифр 2011-1.2-512-073) "Разработка макетов аппаратно-программных комплексов для лабораторной и функциональной диагностики" по теме "Разработка макета аппаратно-программного комплекса для лабораторной и функциональной диагностики широкого спектра признаков методом полиме-разной цепной реакции в реальном времени" (шифр заявки 2011-1.2-512-073-110).

2. Белов Ю.В., Петров А.И., Лавров В.В., Куроч-кин В.Е. Изучение влияния шумов детектора на погрешности количественных анализов нуклеиновых кислот на приборах ПЦР-РВ // Научное приборостроение. 2011. Т. 21, № 2, C. 27-33.

3. Pfaffl M.W. A new mathematical model for relative quantification in real-time RT-PCR // Nucleic Acids Research. Oxford University Press, 2001. V. 29, N 9. e45.

URL: (http://pathmicro.med.sc.edu/pcr/ realtime-home. htm).

4. Rutledge R.G., Stewart D. A kinetic-based sigmoidal model for the polymerase chain reaction and its application to high-capacity absolute quantitative real-time PCR.

URL: (http://www.biomedcentral.com/1472-6750/8/47).

5. Белов Ю.В., Петров А.И., Лавров В.В., Куроч-кин В.Е. Оптимизация параметров сигмоидальной функции при моделировании сигналов ПЦР в реальном времени // Научное приборостроение. 2011. Т. 21, № 3. C. 130-134.

6. Алексеев Я.И., Белов Ю.В., Варламов Д.А. и др. Приборы для диагностики биологических объектов на основе метода полимеразной цепной реакции в реальном времени (ПЦР-РВ) // Научное приборостроение. 2006. Т. 16, № 3. С. 132-136.

Институт аналитического приборостроения РАН, г. Санкт-Петербург

СПИСОК ЛИТЕРАТУРЫ

1. Белов Ю.В., Петров А.И., Лавров В.В., Куроч-кин В.Е. Особенности количественных измерений содержания нуклеиновых кислот методом полиме-разной цепной реакции в реальном времени // Научное приборостроение. 2011. Т. 21, № 1. С. 44-49.

Контакты: Белов Юрий Васильевич, bel3838@mail.ru

Материал поступил в редакцию 8.07.2011.

ERROR ANALYSIS OF REAL TIME PCR SIGNALS MODELED

BY SIGMOID FUNCTION

Yu. V. Belov, A. I. Petrov, V. E. Kurochkin

Institute for Analytical Instrumentation of RAS, Saint-Petersburg

A method of PCR signal model characteristic investigation is discussed. It is based on sigmoid function and uses fixed change of fluorescence intensity during one or several temperature cycles (pulse disturbance). A comparison of pulse and noise disturbance was performed. The conclusion was made that it is suitable to use additional signals filtration and threshold method to reduce noise impact on PCR signals modeling. It was also explained why errors of threshold cycles measuring are reduced.

Keywords: real time PCR, DNA, threshold cycle, sigmoid function