Метод оценки параметров спектральных пиков Текст научной статьи по специальности «Компьютерные и информационные науки»

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2017, том 27, № 3, c. 99-106

МАТЕМАТИЧЕСКИЕ МЕТОДЫ И МОДЕЛИРОВАНИЕ ^ =

В ПРИБОРОСТРОЕНИИ

УДК 543.426; 543.9 © Л. В. Новиков, В. В. Куркина

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

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

Кл. сл.: обработка сигналов, аналитические спектры, оценка параметров, идентификация пиков

ВВЕДЕНИЕ

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

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

Пакет XCMS [3], выполняет операции фильтрации, идентификации, определения площади и выравнивания времени удерживания пиков. В алгоритме обнаружения пиков используется метод согласованной фильтрации с применением

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

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

Программные пакеты MetAlign [5] и М2тте [6] выполняют коррекцию базовой линии, фильтрацию шума и артефактов в спектре масс, детектирование пиков, нормализацию, формирование таблиц пиков, визуализацию результатов.

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

В работе [8] авторы используют вторую производную гауссовой кривой, взятую с обратным знаком ("мексиканская шляпа"), в качестве ядра оператора для вычисления сглаженной второй производной хроматографического сигнала. Обнаружение пика производится по превышению второй производной порога, равному 10-кратной интенсивности шума. Интенсивность пиков определяется в точках максимума по интенсивности норми-

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

В работах [9, 10] используется функция "мексиканская шляпа" в качестве материнского вейв-лета для непрерывного вейвлет-преобразования (НВП) данных. Выполняя НВП при различных величинах масштаба материнского вейвлета, получают двумерный спектр на плоскости время/масштаб, содержащий на каждом масштабе (шкале) пики, порожденные исходными данными. Максимумы пиков, положения которых близки, формируют хребты, объединяющие эти максимумы по всем шкалам. Пик считается обнаруженным, если сумма его максимальных значений по шкалам (по хребту) превышает заданный порог. Положение пика определяется по отсчету при минимальном значении масштаба. Благодаря применению НВП в широком диапазоне масштабов метод позволяет разделить близкие пики без заметной седловины между ними. Точные значения параметров пиков определяют методом подгонки кривых. Применение "мексиканской шляпы" и вообще второй производной данных позволяет, кроме того, удалить базовую линию, линейную в пределах окна дифференцирующей функции.

Алгоритм детектирования пиков в тандеме (газовый хроматограф—спектрометр ионной мобильности) (GC-IMS) [11] основан на вычислении первой и второй частных производных двумерных данных. Величины интенсивности, гауссовой кривизны и матрицы Гессе в каждой точке двумерной поверхности сравниваются с соответствующими порогами, по результатам чего принимается решение о наличии пиков.

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

ТЕОРИЯ

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

|хк (7), к = 1,...,К|, например, гауссовой формы

1.5

0.5

и 8

-0.5

Номер отсчета i (х 10 )

Рис. 1. Модель аналитического сигнала, состоящий из ста пиков гауссовой формы при отношении сигнал/шум 2.5.

Ширина пиков линейно возрастает с увеличением i

и шума п (7) (предполагается, что базовая линия удалена):

у (7 )=Х Хк (7)+п (7),

к=1

где

хк (7) = А ехр|_(7 - 7к У/}; 7

непрерыв-

ное или дискретное время 7 = iSt (i = 1,2,...,N), N — число отсчетов данных, 57 — интервал дискретизации; Ак, 7к, цк — интенсивность, положение и среднеквадратическая ширина к-го пика. Пример такого аналитического спектра, состоящий из ста пиков гауссовой формы с добавлением белого шума, приведен на рис. 1.

Массив дискретных данных у (/'57 ) = у (i) образует N -мерный вектор

У = (у (1), у (2) у (N)) .

Сформируем текущую выборку из М +1 отсчетов с центром в точке / исходных данных:

У^^ =(у*, (1),ук(2),...,ук(М/2),...,ук (М +1)) =

= (у(/ -М/2),...,у(/),...,у(/ + М/2)),

где М выбирается всегда четным из условия 0.5 ^ 0.7 от ширины пика на половине высоты в отсчетах. При этом во избежание скачков на границах массива у , полагаем, что если / - М/2 < 0 ,

то у* (/ - М/2) = у (1), и если / + М/2 > N,

то

у* (/ + М/2) = у (N) .

0

-0,1 i .265

1.2655

1.266

1.2665 1.267

Номер отсчета i (х 104)

Рис. 2. Аппроксимация данных полиномом второго порядка в скользящем окне из 15 точек (М = 14) на сетке (¡= -7,..., j=0,..., j=7).

tl, t0, tr — точки аппроксимирующей параболы, в которых вычисляются первые производные

Аппроксимируем отсчеты в окне у^, полиномом второй степени (параболой) на сетке j = (-М/ 2,...,0,..., М/ 2):

У (j ) = Р^2 + Р2 j + Pз, где р1, р2, р3 — коэффициенты полинома.

М М~

Выберем в интервале — left, r

to = j = 0, tr = j =

t,, to, tr (l

' M

три точки

2 2

right) такие, что t, = j =

" M' " M'

4 , где 4 —

PS

ных данных в точке 1

У 0' ) = У ( 0) = Рз

ОПИСАНИЕ АЛГОРИТМА

(3)

4

целое число, большее или равное М/4 . На рис. 2 показан фрагмент спектра с окном данных из пятнадцати отсчетов (М = 14) и аппроксимирующая парабола с точками 11, t0 и tг.

Величина первой производной параболы в точке tl будет равна

dl = 2М + Р2 (1)

и в точке tr:

¿г = 2М + Р2. (2)

Величина второй производной параболы равна

Р1.

Величина параболы в точке j = 0 равна усредненному значению в скользящей выборке исход-

Используя формулы (1), (2), (3) рассмотрим алгоритм обнаружения пиков и их параметров. В начале обработки исходных данных по этим формулам в программе устанавливаются признаки обнаружения начала и вершины пика:

[0, пик отсутствует, [1, начало пика обнаружено.

10, вершина пика не обнаружена, [1, вершина пика обнаружена.

Условие начала пика имеет вид

^ > 0; ¿г > ^^1оре; Рз > Кш и PS = 0, (4)

где — порог нижнего уровня, выбираемого по стандартному отклонению шума и, ^тт = (1 + 2) и , К — порог по наклону фронта пика, для гауссовой формы пика выбирается из интервала К = 0.003 + 0.007 .

После обнаружения начала пика устанавливаются признаки PS = 1 и PSS = 0 .

- О-8

н

0

1 0.6

к

о

я

Й 0.4

! I

, * « * * * * ж Ж :

* tr * * k ** *** " it*

tSk t0 th " «min \ * * * * : * * *

i i 1

"°Йб5 1.266 1.267 1.268 1.269 1.27 1.271

Номер отсчета 1 (х 104)

Рис. 3. Фрагмент спектра с аппроксимирующей параболой (утолщенная линия) в начале пика.

— порог нижнего уровня; tsk — начало пика

- 0-8 н

0

1 0.6

к

о Я

г, , Ё 04 s

0.2 О

t0 . tr

tl Г t ь ** * s

*1 * th ir 'max # * **

* * * ** * * * * * * *

1.266 1.267 1.268 1.269 1.27 1.271

Номер отсчета 1 (х 104)

Рис. 4. Фрагмент спектра с аппроксимирующей параболой (утолщенная линия) на вершине пика.

^шах — порог максимального уровня; tm¡x к — положение вершины к-го пика

Начало k-го пика определяется по положению точки tj на шкале отсчетов i: tsk = arg (tlпри

условии PS = 1 (см. рис. 3).

Условие вершины пика имеет вид

PS = 1; dl > 0; dr < 0; p < 0 и max (y)> thmax , (5)

где thmax — порог максимального уровня, выбираемый из условия

^ = (4*6)^ .

Оценка амплитуды Лк и положение вершины tmax к к -го пика определяются по точке максимума 1тах параболы:

Лк = таХ (У) , tmax_к = (таХ (У )),. (6)

по шкале отсчетов 1 (см. рис. 4).

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

^ (Сах_к - 'С )< 3.

После обнаружения вершины пика устанавливаются признаки PS = 0 и PSS = 1.

Оценку ширины пика на половине высоты находим из условия равенства участка гауссовой модели в окрестности вершины пика и аппроксимирующей ее параболы (полагаем переменную t непрерывной):

(

Ак exp <{-

t - L

V

Mk

= Plt2 + P2t + P3 = У (t) .

В точке tj параболы имеем

Ak exP I-

ft -1 У

1 max k

Mk

= У (ti).

Прологарифмировав обе части равенства и выполнив преобразования, получим

Mk = ъ -1

1 max k

ln

y (ti)

Ширина пика на половине высоты будет равна

И1/2_к = >/21п2 Мк . (7)

Площадь к -го пика гауссовой формы определяется по формуле

s. =4п Ak

Mk

(8)

Условия конца пика имеют следующий вид: если

PSS = 1; dt <0; dr <0; ръ <thmm, (9)

то обнаружен одиночный (хорошо разделенный) пик.

При этом конец k-го пика определяется по положению точки tk на шкале отсчетов i:

tfk = arg (tk ).; устанавливается признак: пики разделены Rk = 1 (при условии Rk1 Ф 0) и признак PSS = 0.

Если условие (9) не выполняется, но выполняется условие

PSS = 1; dt < 0; dr > thsiope; Р3 > th^, (10)

то обнаружена седловина между соседними пиками.

При этом конец k-го пика tfk определяется по положению точки минимума параболы в седловине между пиками tsad = arg (min (y)) по шкале

отсчетов i: tfk = arg (tsad). (см. рис. 5).

1.268 1.269 1.27 1.271 Номер отсчета 1 (х 104)

Рис. 5. Фрагмент спектра с аппроксимирующей параболой (утолщенная линия) на седловине пика.

tsaй — положение седловины между пиками; tfk — конец пика

Если при этом выполняется условие, например,

У (^ )>= \ Лк,

т. е. если глубина седловины не превышает одной трети интенсивности пика, то можно считать, что соседние пики не разделены и устанавливается признак Rk = 0 . В противном случае устанавливается признак Як = 1.

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

- начало ts

к

- положение максимума tmax к,

- интенсивность Лк,

- площадь !Ик,

- конец tfk,

- признак разделения Rk .

МОДЕЛИРОВАНИЕ

пиков гауссовой формы с добавлением белого шума (рис. 1) путем сравнения результатов обработки с алгоритмом, приведенным в работе [12]. Относительная погрешность оценки амплитуды, ширины и площади пиков производилась по десяти реализациям шума путем вычисления СКО параметров каждого пика и деления его на истинное значение параметра. Погрешность положения пика определялась по его отклонению от истинного значения. Максимальные значения погрешностей для обоих методов обработки сведены в таблицу. Кроме того, вычислялась вероятность правильного обнаружения как отношение числа обнаруженных пиков по всем десяти реализациям, положение которых отличается от заданных в модели не более, чем на четыре отсчета, к общему числу модельных пиков. Вероятность ложного обнаружения вычислялась как отношение числа обнаруженных пиков по всем десяти реализациям, не заданных в модели, к общему числу модельных пиков. Сравнение столбцов 1 и 2 таблицы показывает явное преимущество предлагаемого алгоритма как по уровню погрешности оценки параметров пиков, так и по величине вероятности ложного обнаружения.

Проверка эффективности алгоритма производилась по модельному спектру, состоящему из 100

Абсолютные (для положения пиков) и относительные среднеквадратические ошибки (СКО) оценки параметров пиков при различных отношениях сигнал/шум для метода, изложенного в [12], — 1 и предлагаемого метода — 2

Параметр Отношение сигнал/шум

2.5 5 10 20

Метод расчета

1 2 1 2 1 2 1 2

Положение пика (в отсчетах) 3.98 3 1.68 2 0.75 1 0.38 0

Амплитуда 0.08 0.04 0.038 0.018 0.21 0.01 0.01 0.004

Ширина 0.90 0.34 0.46 0.23 0.19 0.13 0.08 0.047

Площадь 1.01 0.27 0.43 0.19 0.18 0.11 0.08 0.030

Вероятность правильного обнаружения 0.980 0.999 1 1 1 1 1 1

Вероятность ложного обнаружения 0.016 0.002 0.003 0.001 0.002 0 0.002 0

ЗАКЛЮЧЕНИЕ

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

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

1. Ewing B., Green Ph. Base-calling of automated sequencer traces using Phred. II. Error probabilities // Genome Res. 1998. Vol. 8, no. 3. P. 186-194. Doi: 10.1101/gr.8.3.186.

2. CookD.W., Rutan S.C. Chemometrics for the analysis of chromatographic data in metabolomics investigations // J. Chemometrics. 2014. Vol. 28, no. 9. P. 681-687.

3. Smith C.A., Want E.J., O'Maille G., Abagyan R., Siuz-dak G. XCMS: processing mass spectrometry data for metabolite profiling using nonlinear peak alignment, matching, and identification // Anal. Chem. 2006. Vol. 78, no. 3. P.779-787. Doi: 10.1021/ac051437y.

4. Wei X., Shi X., Kim S., Zhang L., Patrick J.S., Binkley J., McClain C., Zhang X. Data preprocessing method for liquid chromatography-mass spectrometry based metabo-lomics // Anal. Chem. 2012. Vol. 84, no. 18. P. 79637971. Doi: 10.1021/ac3016856.

5. Lommen A. MetAlign: interface-driven, versatile metabo-lomics tool for hyphenated full-scan mass spectrometry data preprocessing // Anal. Chem. 2009. Vol. 81, no. 8. P. 3079-3086. Doi: 10.1021/ac900036d.

6. Katajamaa M., Oresic M. Processing methods for differential analysis of LC/MS profile data // BMC Bioinfor-matics. 2005. Vol. 6. P. 179-190. Doi: 10.1186/14712105-6-179.

7. Vivo-Truyols G., Torres-Lapasio J.R., van Nederkas-sel A.M., Heyden Y.V., Massart D.L. Automatic program for peak detection and deconvolution of multi-overlapped chromatographic signals. Part I: Peak detection // Journal of Chromatography A. 2005. Vol. 1096, no. 1-2. P. 133145. Doi: 10.1016/j.chroma.2005.03.092.

8. FredrikssonM.J., Petersson P., Axelsson B.-O., BylundD. An automatic peak finding method for LC-MS data using Gaussian second derivative filtering // J. Sep. Sci. 2009.

Vol. 32, no. 22. P. 3906-3918. Doi: 10.1002/jssc.200900395. 9. Gregoire J.M., Dale D., van Dover B.. A wavelet transform algorithm for peak detection and application to powder x-ray diffraction data // Review of Scientific Instruments. 2011. Vol. 82, no. 1. 015105. Doi: 10.1063/1.3505103.

10. Du P., Kibbe W.A., Lin S.M. Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching // Bioinformatics. 2006. Vol. 22, no. 17. P. 2059-2065.

Doi: 10.1093/bioinformatics/btl355.

11. Slodzinski R., Hildebrand L., Vautz W. Peak detection algorithm based on second derivative properties for two dimensional ion mobility spectrometry signals // Integration of Practice-Oriented Knowledge Technology: Trends and Prospectives / Madjid Fathi (ed.). Springer-Verlag, Berlin, Heidelberg, 2013. P. 341-354.

12. O'Haver T. Interactive Signal Processing Tools. Peak Finding and Measurement. URL:

http://terpconnect.umd.edu/~toh/spectrum/SignalProcessi ngTools.html.

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

Санкт-Петербургский государственный технологический институт (Технический университет)

(Куркина В.В.)

Контакты: Новиков Лев Васильевич, novik3 8 @mail. ru

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

ISSN 0868-5886

NAUCHNOE PRIBOROSTROENIE, 2017, Vol. 27, No.3, pp. 99-106

THE METHOD FOR ESTIMATION OF SPECTRAL PEAK PARAMETERS

L. V. Novikov1, V. V. Kurkina2

1 Institute for Analytical Instrumentation of RAS, Saint-Petersburg, Russia 2St. Petersburg State Technological Institute (Technical University), Russia

A new economical algorithm for estimating the parameters of signals in mass spectrometry, chromatography, and other applications, representing a sequence of peaks against a noise background, is proposed. A traditional approach for the peak top detection by the intersection the zero line by the first derivative is developing. In order to increase the reliability of detection of the beginning, end, vertex and saddle between peaks, it is proposed to compare the values of the derivatives at three points of the sliding data window.

Keywords: signal processing, analytical spectra, estimation of parameters, identification of peaks

REFERENСES

1. Ewing B., Green Ph. Base-calling of automated sequencer traces using Phred. II. Error probabilities. Genome Res., 1998, vol. 8, no. 3, pp. 186-194. Doi: 10.1101/gr.8.3.186.

2. Cook D.W., Rutan S.C. Chemometrics for the analysis of chromatographic data in metabolomics investigations. J. Chemometrics, 2014, vol. 28, no. 9, pp. 681-687.

3. Smith C.A., Want E.J., O'Maille G., Abagyan R., Siuz-dak G. XCMS: processing mass spectrometry data for metabolite profiling using nonlinear peak alignment, matching, and identification. Anal. Chem., 2006, vol. 78, no. 3, pp.779-787. Doi: 10.1021/ac051437y.

4. Wei X., Shi X., Kim S., Zhang L., Patrick J.S., Binkley J., McClain C., Zhang X. Data preprocessing method for liquid chromatography-mass spectrometry based metabo-lomics. Anal. Chem., 2012, vol. 84, no. 18, pp. 79637971. Doi: 10.1021/ac3016856.

5. Lommen A. MetAlign: interface-driven, versatile metabo-lomics tool for hyphenated full-scan mass spectrometry data preprocessing. Anal. Chem., 2009, vol. 81, no. 8, pp. 3079-3086. Doi: 10.1021/ac900036d.

6. Katajamaa M., Oresic M. Processing methods for differential analysis of LC/MS profile data. BMC Bioinformat-ics, 2005, vol. 6. pp. 179-190. Doi: 10.1186/1471-2105-6179.

7. Vivo-Truyols G., Torres-Lapasio J.R., van Nederkas-sel A.M., Heyden Y.V., Massart D.L. Automatic program for peak detection and deconvolution of multi-overlapped chromatographic signals. Part I: Peak detection. Journal

Contacts: Novikov Lev Vasil'evich, novik3 8 @mail. ru

of Chromatography A, 2005. vol. 1096, no. 1-2, pp. 133145. Doi: 10.1016/j.chroma.2005.03.092.

8. Fredriksson M.J., Petersson P., Axelsson B.-O., By-lund D. An automatic peak finding method for LC-MS data using Gaussian second derivative filtering. J. Sep. Sci., 2009, vol. 32, no. 22, pp. 3906-3918.

Doi: 10.1002/jssc.200900395.

9. Gregoire J.M., Dale D., van Dover B. A wavelet transform algorithm for peak detection and application to powder x-ray diffraction data. Review of Scientific Instruments, 2011, vol. 82, no. 1, 015105. Doi: 10.1063/1.3505103.

10. Du P., Kibbe W.A., Lin S.M. Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching. Bioinformatics, 2006, vol. 22, no. 17, pp. 2059-2065.

Doi: 10.1093/bioinformatics/btl355.

11. Slodzinski R., Hildebrand L., Vautz W. Peak detection algorithm based on second derivative properties for two dimensional ion mobility spectrometry signals. Integration of Practice-Oriented Knowledge Technology: Trends and Prospectives, Madjid Fathi (ed.), Springer-Verlag, Berlin, Heidelberg, 2013, pp. 341-354.

12. O'Haver T. Interactive Signal Processing Tools. Peak Finding and Measurement. URL: http://terpconnect.umd.edu/~toh/spectrum/SignalProcessi ngTools.html.

Article received in edition: 20.06.2017