CC BY
24ISSN 0868-5886 НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2011, том 21, № 3, c. 130-134
-- ОБРАБОТКА И ПРЕДСТАВЛЕНИЕ ДАННЫХ
УДК 543.426; 543.9
© Ю. В. Белов, А. И. Петров, В. В. Лавров, В. Е. Курочкин
ОПТИМИЗАЦИЯ ПАРАМЕТРОВ СИГМОИДАЛЬНОЙ ФУНКЦИИ ПРИ МОДЕЛИРОВАНИИ СИГНАЛОВ ПЦР В РЕАЛЬНОМ ВРЕМЕНИ
В статье предложена усовершенствованная модель сигналов ПЦР в реальном времени и методика ее реализации на основе сигмоидальной функции. С помощью численного эксперимента показано, что при соответствующем выборе уровня порога модели удается устранить влияние дрейфа нулевой линии на результаты моделирования. Приведены критерии приближения модели к фактическому росту интенсивности сигналов флуоресценции в ходе реакции. Показана возможность полной автоматизации вычислений пороговых циклов сигналов ПЦР в реальном времени.
Кл. сл.: ПЦР в реальном времени, пороговый цикл, количественные анализы, сигмоидальная функция
ВВЕДЕНИЕ
Первым шагом при обработке сигналов анализаторов нуклеиновых кислот методом полимераз-ной цепной реакции в реальном времени (ПЦР-РВ) при относительных количественных измерениях является определение пороговых циклов (далее — СО, при которых сигналы флуоресценции образцов достигают постоянного уровня, выбранного оператором. Высокая воспроизводимость определения пороговых циклов (среднеквадратичное отклонение менее чем 0.2 циклов) [1] позволяет сравнивать концентрацию исследуемых и калибровочных образцов [2]. При этом используется основное свойство ПЦР — рост сигналов на начальной стадии реакции происходит по экспоненциальному закону. В статье [3] был выполнен практический анализ погрешностей определения С путем расчета значений стандартного отклонения пороговых циклов. При малых объемах реактивов, которые экономически выгодны при выполнении массовых анализов, возникает необходимость изучения влияния шумов и дрейфа нулевой линии детектора на погрешности количественных измерений. В статье [4] было выполнено сравнение разных способов определения С В настоящей статье предложена усовершенствованная модель для обработки сигналов на основе сигмоидальной образной) функции (СФ) [5] и показана возможность практически полного устранения влияния дрейфа детектора.
АНАЛИТИЧЕСКАЯ ФОРМА СИГМОИДАЛЬНОЙ ФУНКЦИИ
Кинетическая модель сигналов ПЦР-РВ на базе нелинейной сигмоидальной функции (МСФ) мо-
жет быть выражена в следующем виде: F = Fm / (1+exp ((Со - C) / k)) +
+fo + f • C + f • C2, (1)
где F — интенсивность флуоресценции реакции при порядковом номере цикла амплификации С; Fm — максимальная интенсивность флуоресценции в конце реакции; С0 — дробный цикл, при котором интенсивность флуоресценции достигает половины значения Fm; k — коэффициент, который определяется наклоном кривой; f — смещение нулевой линии (фон детектора).
Традиционное выражение МСФ [3] дополнено нами еще двумя слагаемыми:
f • C + f • C2.
Коэффициент f определяется линейным дрейфом нулевой линии, а коэффициент f учитывает квадратичную составляющую дрейфа нулевой линии. Необходимость дополнительных слагаемых возникает в случаях, когда с целью экономии используются малые объемы реактивов, при этом интенсивность флуоресценции реакции становится соизмеримой с дрейфом нулевой линии.
Целью построения МСФ является уменьшение погрешности вычисления порогового цикла.
МЕТОДИКА ПОСТРОЕНИЯ И ОПТИМИЗАЦИИ ПАРАМЕТРОВ МСФ
С помощью электронной таблицы Excel можно построить и оптимизировать параметры МСФ (табл. 1).
В ячейках С4:С9 табл. 1 заданы исходные параметры модели. В терминах Excel интенсивность флуоресценции в ячейке С11 вычисляется следующим образом:
ОПТИМИЗАЦИЯ ПАРАМЕТРОВ СИГМОИДАЛЬНОЙ ФУНКЦИИ. Табл. 1. Электронная таблица для вычисления параметров МСФ
A B C D E F
1
2 Порог = 500
3 Ct = 23,500
4 Co = 23,5
5 Fm = 1000
6 k = 1,39
7 fo = 0
8 f1 = 0
9 f2 = 0
10 С= S = 0
11 1 9,33E-05 9,33E-05 0 9,33E-05 0
24 14 1, 07 1, 07 0 1, 07 0
33 23 411,029 411,029 0 411,029 0
34 24 588,971 588,971 0 588,971 0
43 33 998, 9 998, 9 0 998, 9 0
60 50 1000 1000 0 1000 0
С11 =$С$5/(1+EXP(($С$4-А11)/$С$6))+ $С$7+ $С$8*А11+$С$9*(А11)Л2,
где порядковый номера цикла С = 1 приведен в ячейке A11. Полностью МСФ строится в ячейках С11:С60 как зависимость от номера цикла, записанного в ячейках А11:А60 = 1^50.
В ячейках В11:В60 можно расположить полученную на анализаторе ПЦР зависимость флуоресценции от номера цикла в числовом выражении, желательно нормированную при Fm = 1000 и f0 = 0 о.е. (относительных единиц). Ячейки D11:D60 можно использовать для имитации дрейфа нулевой линии детектора и других видов помех. В ячейках E11:E60 вычисляются суммы (B11+D11):(B60+D60). В ячейках F11:F60 получаются квадраты разностей (Е11- С11)2:(Е60 - С60)2, а в ячейке С10 — остаточная сумма квадратов разностей S.
Оптимизация параметров модели выполняется путем использования инструмента Поиск решения из меню Сервис. Для этого необходимо выполнить следующий набор команд:
- "Установить целевую ячейку" -С10 - "Равной: минимальному значению" - "Выполнить" - "Изменяя ячейки" -С4 :С9.
В результате в ячейках С11 : С60 получаются чи-
(2)
словые значения зависимости оптимизированной модели от номера цикла, а в ячейках С4:С9 — новые параметры МСФ. Модель получается оптимизированной по методу наименьших квадратов ошибок.
В приведенном примере в ячейки В11:В60 внесена численная копия исходной нормированной МСФ из ячеек С11:С60, поэтому после выполнения вычислений параметры модели остались прежними.
Традиционно уровень порога выбирается оператором на участке экспоненциального роста флуоресценции [1]. Экспоненциальный участок может быть определен, как линейная часть графика сигнала флуоресценции на логарифмической шкале. На рис. 1 приведены почти полностью совпадающие графики реального сигнала и модели после нормировки = 1000 и ^ = 0 о.е.). Для определения пороговых циклов С в статье [4] был рекомендован уровень порога, равный 500 о.е. При таком уровне обеспечивается минимальное влияние шумов детектора на погрешности измерения, этот уровень условно можно считать концом участка экспоненциального роста флуоресценции (при линейном масштабе — это точка перегиба).
Дробный цикл С0 удобно принять в качестве величины порогового цикла С реального сигнала
Для имитации дрейфа нулевой линии детектора формировались числовые последовательности /0, / • С, /2 • С2 и /3 • С3 с величинами коэффициентов /0 = 100 о.е., £ = 2 о.е., £ = 0.04 о.е. и £з = 0.0008 о.е. При таких коэффициентах максимальное отклонение нулевой линии составляет 100 о.е., т. е. 10 % от максимальной нормированной интенсивности флуоресценции в конце реакции Fm = 1000 о.е. Дополнительно формировалась числовая последовательность суммы /0 + /1 • С + /2 • С2.
Числовые последовательности поочередно вносились в ячейки D11:D60. При этом в ячейках Е11:Е60 выполнялось суммирование исходной нормированной сигмоидальной функции с числовой последовательностью, имитирующей дрейф нулевой линии детектора, а в ячейке С10 получалась новая сумма квадратов разностей 5".
После выполнения действий по оптимизации параметров модели (с целью получения минимального значения 5) в ячейках С11:С60 получаются числовые значения зависимости новой оптимизированной модели от номера цикла. Новые параметры МСФ в ячейках С4:С9 приведены в табл. 2.
Сравнивая полученные в табл. 2 параметры (выделены жирным шрифтом) с заданными параметрами (/0 = 100, / = 2, /2 = 0.04 о.е.), можно сделать вывод о хорошем совпадении параметров моделей и исходных последовательностей с дрейфом нулевой линии. Особенно ценным качеством модели с предложенными дополнительными параметрами является малая зависимость от дрейфа величины С0, которую можно принять в качестве порогового цикла. Даже при дрейфе нулевой линии, который характеризуется кубической зависимостью, величина С0 отличается от исходной величины менее чем на 0.5 % амплификационного цикла.
Парамет- Исход- Имитация дрейфа нулевой линии детектора
ры МСФ ные значения параметров fc = = 100 о.е. (макс.) f • с = = 100 о.е. (макс.) f2 • C2 = = 100 о.е. (макс.) fo +f • C + +f2 • C2 = = 300 о.е. (макс.) fs • C3 = = 100 о.е. (макс.) fa • C3 = (от 14 до 33 циклов)
Co 23.5 23.499 23.499 23.500 23.500 23.496 23.500
Fm 1000 999.997 999.988 999.988 999.989 988.441 998.215
к 1.39 1.390 1.390 1.390 1.390 1.371 1.389
fo 0 99.999 -2.003 0.038 98.038 2.728 8.024
f1 0 0.000118 2.0004 -0.083 1.917 -0.900 -1.241
f2 0 -4.5E-07 -1.7E-06 0.042 0.042 0.060 0.057
S 0 2.07E-05 0.00016 0.00015 0.0002 49.041 0.010
R2 1 1 0.9994 0.9990 0.9973 0.9986 1
18 2 0 £ 2 2 4 г 6 г 8 3 о Н о м е г- ц [.-1 к л о
Рис. 1. Графики сигнала и модели после нормировки
ПЦР. Условие Ct = С0 выполняется, если выбрать уровень порога, при котором интенсивность флуоресценция достигает половины значения Fm (в этом примере Fm /2 = 500 о.е.).
ИССЛЕДОВАНИЕ ВЛИЯНИЯ ДРЕЙФА НУЛЕВОЙ ЛИНИИ ДЕТЕКТОРА
С помощью электронной таблицы Excel выполнен численный эксперимент, результаты которого доказывают, что даже при неблагоприятных условиях c помощью МСФ удается устранить влияние дрейфа нулевой линии.
Табл. 2. Результаты вычисления параметров МСФ
ОПТИМИЗАЦИЯ ПАРАМЕТРОВ СИГМОИДАЛЬНОЙ ФУНКЦИИ.
133
После завершения ПЦР, результатом которой является быстрый рост сигнала флуоресценции, может наблюдаться медленное увеличение или уменьшение сигнала флуоресценции. Этот эффект подобен дрейфу нулевой линии, хотя имеет другую причину. Уменьшить влияние этого эффекта можно путем использования укороченной модели сигмоидальной функции (УМСФ). Приведенная в ячейке В24 табл. 1 интенсивность модельного сигнала при 14-м цикле отличается от нуля примерно на 1 о.е., а интенсивность модельного сигнала при 33-м цикле (В43) отличается от максимального значения (1000 о.е.) также примерно на 1 о.е. Поэтому при допустимой погрешности интенсивности сигнала около 0.1 % можно применить УМСФ в интервале от 14 до 33 циклов. Результаты вычисления основных параметров (С0, k и Я) УМСФ при дрейфе нулевой линии с кубической зависимостью, приведенные в крайнем правом столбце табл. 2, мало отличаются от исходных параметров.
При использовании УМСФ можно сократить время, затраченное на выполнение ПЦР, если не регистрировать сигналы флуоресценции в течение нескольких первых циклов.
Для использования в качестве иллюстрационного материала можно построить упрощенные модели, используя данные табл. 1 с вновь вычисленными параметрами С0 и k при Fm =1000 о.е. и ве-
личинами f , fi, f2, равными нулю.
Критерием приближения прогнозируемой модели к исходной фактической последовательности может служить остаточная сумма квадратов разностей S. В безразмерной форме функциональную близость можно оценить с помощью квадрата коэффициента корреляции Я2 (функция Excel КВПИРСОН). Вычисленные значения этих критериев приближения приведены в табл. 2. Степень приближения можно считать очень высокой, если величина Я2 более 0.99.
Подобная методика построения МСФ реализована при модернизиации программного обеспечения анализатора нуклеиновых кислот АНК-32, который серийно выпускается в ИАП РАН [6]. Сигналы, имитирующие дрейф нулевой линии детектора (максимальное отклонение — до 300 о.е.), вводились в программу обработки информации в виде текстового файла. На рис. 2 одновременно изображены графики введенных сигналов и графики МСФ. Можно отметить почти полное попарное совпадение этих графиков. Различие величин С0 всех графиков не превышает 1 % порогового цикла.
В модернизированном программном обеспечении реализована нормировка (привязка графиков к нулю и к максимальному значению), а также выполняется коррекция наклона нулевой линии. Все
эти операции выполняются в ручном и автоматическом режимах.
ВЫВОДЫ
Предложена усовершенствованная модель сигналов ПЦР-РВ и методика оптимизации ее параметров на основе одного из наиболее эффективных способов обработки сигналов — способа аппроксимации с помощью сигмоидальной функции.
При определении порогового цикла предложено выбирать уровень порога, при котором интенсивность модельного сигнала достигает половины значения максимальной интенсивности в конце реакции. Выполнен численный эксперимент, результаты которого доказывают, что даже при неблагоприятных начальных условиях при таком выборе уровня порога удается устранить влияние дрейфа нулевой линии на результаты моделирования. Рассмотрены преимущества использования укороченной модели сигмоидальной функции.
Приведены критерии приближения прогнозируемой модели к фактическому росту интенсивности сигналов флуоресценции в ходе реакции.
Показана возможность полной автоматизации вычислений пороговых циклов сигналов ПЦР-РВ.
СПИСОК ЛИТЕРАТУРЫ
e45.
URL: (http://pathmicro.med.sc.edu/pcr/ realtime-home. htm).
3. Белов Ю.В., Петров А.И., Лавров В.В., Куроч-кин В.Е. Особенности количественных измерений содержания нуклеиновых кислот методом полиме-разной цепной реакции в реальном времени // Научное приборостроение. 2011. Т. 21, № 1. С. 44-49.
4. Белов Ю.В., Петров А.И., Лавров В.В., Куроч-кин В.Е. Изучение влияния шумов детектора на погрешности количественных анализов нуклеиновых кислот на приборах ПЦР-РВ // Научное приборостроение. 2011. Т. 21, № 2. C. 27-33.
5. 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).
6. Алексеев Я.И., Белов Ю.В., Варламов Д.А. и др. Приборы для диагностики биологических объектов на основе метода полимеразной цепной реакции в реальном времени (ПЦР-РВ) // Научное приборостроение. 2006. Т. 16, № 3. С. 132-136.
Институт аналитического приборостроения РАН, г. Санкт-Петербург
1. Stahlberg A., Zoric N., Aman P., Kubista M. Quantita- Контакты: Белов Юрий Васильевич, tive real-time PCR for cancer detection: the lymphoma Bel3838@mail ru
case // Expert Rev. Mol. Diagn. 2005. V. 5, N 2. P. 221-230.
URL:(http://tataa.com/iiles/documents/qpcr_review.pdi).
2. Pfaffl M.W. A new mathematical model for relative Материал поступил в редад^ 18042011. quantification in real-time RT-PCR // Nucleic Acids
Research. Oxford University Press, 2001. V. 29, N 9.
OPTIMIZATION OF SIGMOID FUNCTION PARAMETERS IN REAL-TIME PCR SIGNALS MODELING
Yu. V. Belov, A. I. Petrov, V. V. Lavrov, V. E. Kurochkin
Institute for Analytical Instrumentation of RAS, Saint-Petersburg
Advanced model of real-time PCR signals and methods of its realization on the sigmoid function base is discussed. Numerical experiment showed that the influence of the zero line drift on modeling results can be removed by choosing an appropriate level of threshold model. The criteria of model approximation to actual fluorescence signals growing during the reaction are proposed. Possibility of fully automatic computation of realtime PCR signals threshold cycles was shown.
Keywords: real-time PCR, threshold cycle, quantitative analysis, sigmoid function
