[cite]10.1093/bioinformatics/btl355[/cite]

1 Introduction

피크를 찾는 일은 질량 분석 기반의 단백질 분석(Mass Spectrometry/MS-based Proteomic data analysis)에서 매우 중요한 전처리 과정입니다. 피크 검출의 결과과 다음에 이어지는 작업에 직접적으로 영향을 끼치므로 그 결과가 중요한데 반해 MS spectrum의 노이즈나 높은 오탐지율은 어려운 과제로 남아있습니다.

이 논문에선 특별히 surface enhanced laser desorption ionization-time of flight(SELDI-TOF) spectroscopy에 초점을 맞추고 있습니다. 이 분야에서는 관심 영역을 결정하기 위한 가장 첫번째 단계에서 피크 검출 방법이 사용됩니다. 현재 거의 대부분의 피크 검출 방법은 어떠한 임계값 이상의 신호대잡음비(SNR)를 이용하여 지역 최고점은 방법을 사용하고 있습니다. SNR은 피크의 주변 노이즈 수준과 피크의 비에 의해 결정되지만, 높은 값을 갖는 것이 꼭 실제 피크라고 할 수는 없습니다. 어떤 노이즈들은 피크와 비슷한 스파이크 형태로 나타나고, 이러한 경우 그보다 낮은 피크가 우리가 찾고자 하는 실제 피크일런지도 모릅니다.

기준선(baseline)의 제거와 스무딩(smoothing)은 전처리 단계로 보통 피크 검출 이전에 행해집니다. 하지만 이렇게 별도의 단계로 행해진 것들은 다음 단계에 영향을 끼치며 원래 값으로 복구하는 것 또한 불가능합니다. 한편 여러가지 방법들 중, Continuous wavelet transform (CWT) 기반의 패턴 매칭 알고리즘이 있습니다. 이를 이용하면 기준선은 자연스럽게 제거가 되며 스무딩도 필요치 않게 됩니다. 즉 이 방법을 이용하면 전처리 과정 없이 원 데이터(raw data)에 알고리즘을 적용할 수 있을 뿐 아니라 그 결과 또한 일관성있게 나올 것이라 예상할 수 있습니다.

MS data에서 실제 피크는 어떤 모양이나 패턴을 가지고 있습니다. 하지만 실제로는, 피크가 얼마만큼의 폭을 가지고 있는지, 그리고 샘플마다 이것의 수준이 크게 변화한다는 것이 어려운 점입니다. 뿐만 아니라 같은 샘플 내에서도 피크의 너비나 높이가 변화하기까지 합니다. CWT 기반의 방법은 이러한 문제점들을 해결하면서도 낮은 오탐지율을 보여줍니다. Lange가 2006년에 보여준 연구는 CWT를 이용하여 스펙트럼을 여러 조각들로 쪼갠 뒤에 CWT를 사용하여 특정 스케일에서의 피크와그 파라미터들을 탐지하였습니다. 하지만 데이터를 분석하기 이전에 여러개의 조각들에 대해서 적당한 스케일을 정하는 것은 쉽지 않은 문제입니다.

패천 매칭 방법을 좀더 경고하게 만들기 위해서 MS 피크 검출 방법에 CWT를 적용하였습니다. Lange와는 다르게 여기서는 전체 스펙터름에 CWT를 직접 적용한 뒤에 그 결과인 2D CWT 계수 행렬들을 이용하고자 합니다. 계수 행렬에서 나타나는 능선(ridge)들은 MS 스펙트럼에서의 피크와도 관련이 있기 때문에 계수 행렬에서의 능선들을 차아 연결하여 피크들을 찾고, 또한 계수 들로부터 효과적으로 SNR을 계산해낼 것입니다. 전처리 문제 없이 Wavelet 공간에서 계산을 직접 하기 때문에 기준선이나 스무딩 관련한 문제들을 생각하지 않아도 됩니다.

2 Methods

2.1 Continuous wavelet transform

Wavelet transformation은 크게 discrete wavelet transform (DWT)와 CWT로 나누게 됩니다. DWT는 스케일과 적용 지점을 2의 승수들로 사용하여 중복된 정보가 없이 최소한의 데이터로 원래의 신호를 복원할 수 있습니다. 하지만 CWT의 경우 연속된 공간에서의 모든 스케일에 대해서 계산하므로 중복이 발생할 수밖에 없지만 이러한 부분에서 피크의 모양이나 요소들을 쉽게 해석할 수 있는 정보를 제공합니다. 또한 CWT에서는 orthogonal wavelet을 사용할 필요가 없습니다. 이런 이유로 CWT는 패턴 매칭에서 주로 사용됩니다. CWT는 다음과 같이 표현할 수 있습니다.

$$C(a, b) = \int_R s(t) \psi_{a, b}(t) dt \\ \psi_{a,b}(t) = \frac{1}{\sqrt{a}}\psi(\frac{t-b}{a}), a \in R^+ - {0}, b \in R$$

여기서 $s(t)$는 원래의 신호이고 $a$는 스케일, $b$는 위치 이동(translation), $\psi(t)$는 mother wavelet, $\psi_{a,b}(t)$는 스케일과 위치 이동을 적용한 wavelet입니다. $C$는 2D 행렬로 wavelet 계수들을 요소로 갖습니다. 직관적으로 wavelet 계수는 신호와 wavelet간의 매칭 정보를 반영한다고 생각할 수 있는데, 높은 계수들이 높은 매칭도를 갖는다고 생각할 수 있습니다. 피크 검출에 이를 이용하려고 한다면 wavelet이 피크기 갖는 기본적인 특징을 가져야 할 것입니다. 때문에 이 연구에서는 Maxican Hat wavelet을 선택하였습니다. Maxican Hat wavelet의 유효 범위는 [-5, 5]이지만, $a=1$로 스케일 조절이 없었다고 하였을 때, 피크에 가장 잘 맞는 너비는 2개의 샘플 간격으로 나타났습니다. 즉, 스케일을 $a_1$로 조절한다면 피크의 너비는 $2 a_1$가 최적의 너비가 됩니다.

MS 스펙트럼을 CWT를 해보면, 피크의 중심 부분에서 CWT 계수가 지역 최고점이 나타납니다. 따라서 스케일 $a=1$로부터 시작하여 점차적으로 스케일을 늘려가면서 조사를 해보면, 점차적으로 지역 최고점의 계수의 값이 즈가하다가 피크의 너비와 스케일이 정확하게 들어맞았을 때 가장 최고점을 기록하게 되고, 다시 점점 낮아지게 될 것입니다. 이를 위치까지 고려하면 2D CWT 계수 행렬과 그 값을 갖는 3차원 공간에서 능선을 찾는 문제로 생각해볼 수 있습니다. 이것은 직접 CWT 계수들에서 최적값을 찾는 문제보다는 덜 민감하게 반응하는 장점이 있습니다. 결과적으로 능선에서 가장 최고 값을 갖는 CWT 스케일을 선택하면 피크의 너비를 구할 수 있습니다.

2.2 Removal of the baseline

CWT기반의 패턴 매칭 방법을 적용하면, 전처리 과정없이 원 데이터를 직접 계산하여 피크를 검출하게 됩니다. 더불어 기준선의 제거 또한 이루어집니다. 원 데이터의 피크를 $P_{raw}(t)$라고 한다면 이를 다음과 같이 표현할 수 있습니다.

$$P_{raw}(t) = P(t) + B(t) + C, t \in [t_1, t_2]$$

여기서 $P(t)$는 실제 피크를, $B(t)$는 평균값 0을 갖는 기준선 함수, $C$는 상수이고, $[t_1, t_2]$는 피크가 존재하는 영역입니다. 이를 앞의 CWT를 이용하여 나타내면 다음과 같아집니다.

$$C(a, b) = \int_R P(t) \psi_{a,b}(t)dt + \int_R B(t) \psi_{a,b}(t)dt + \int_R C \psi_{a,b}(t)dt$$

피크가 존재하는 영역 내에서 기준선은 천천히 변화하고 한 방향(monotonic)으로만 변한다고 가정한다면, 기준선은 상수 C에 odd function $B(t)$를 이용하여 추정할 수 있습니다. 먼저 Wavelet $\psi_{a,b}$는 평균 0을 가지므로 세번째 항은 0이 될 것입니다. 또한 Maxican Hat처럼 대칭인 wavelet을 사용하였으므로 두번째 항 또한 0에 가까워질 것입니다. 따라서 기준선 함수가 피크 영역 내에서 천천히 변화하는 한, CWT 계수를 구하는 과정에서 기준선은 자동적으로 제거가 될 것입니다.

2.3 Peak identification process

먼저 MS 스펙트럼 데이터에 대해서 CWT를 1에서 64까지 2를 각녁으로 33개의 스케일에 대해서 수행합니다. 2D CWT 계수 행렬을 살펴보면 주된 피크는 길고 높은 능선을, 작은 피크들은 낮고 짧은 능선으로 나타납니다.

Identify the ridges by linking the local maxima

우리가 찾고자 하는 능선은 각 스케일에 대한 CWT 계수들의 지역 최고점들을 잇는 것으로 찾을 수 있습니다. 먼저 계수행렬의 각 스케일에서 지역 최고점들을 미리 찾아놓습니다. 최고점은 슬라이딩 윈도우를 이용해서 찾게 되는데, 슬라이딩 윈도우는 각 스케일에 비례하도록 설정합니다. 2D CWT 계수 행렬의 크기는 N * M인데, N은 스케일의 수, M은 MS 스펙트럼의 길이입니다. 능선을 찾는 알고리즘은 다음과 진행됩니다.

1. 가장 큰 스케일에서 능선의 시작점들을 여러개 찾아 초기화합니다. 행렬에서 가장 큰 스케일에 해당하는 n번째 행에서 최고점들을 찾고, 능선의 갭(gap)은 0으로 시작합니다.
2. 갭이 적절한 임계값보다 작은 각 능선에 대해서, 현재 스케일에서의 최고점 위치에 가장 가까운, 다음 스케일에서의 최고점을 찾습니다. 탐색 범위는 해당 스케일 수준의 슬라이딩 윈도우를 벗어나지 않도록 해야 합니다. 만약 가까운 지점이 없으면 갭을 1 증가시키고, 지점이 있다면 다시 갭을 0으로 초기화합니다.
3. 입계값보다 큰 갭 값을 갖는 능선이 있다면, 이를 저장하고 탐색 리스트에서 제거합니다.
4. 앞의 스케일에서 이어지지 않은 지점은 새로운 능선으로 다시 시작합니다.
5. 2-4과정을 가장 낮은 스케일(n=1)에 도달할 때 까지 반복합니다.

Definition of signal to noise ratio

가정에 따라 피크 부분의 신호 세기는 능선에서 가장 높은 계수를 갖는 부분에서 정의됩니다. 노이즈의 경우 기준선의 이동은 wavelet 공간으로 변환될 때 거의 사라졌을 것이므로, 가장 낮은 스케일의 CTW 계수에서 노이즈 수준을 추정할 수 있습니다. 가장 낮은 스케일에서의 CWT 계수 중, 피크 주변 윈도 내에서 계수의 절대값 95 백분위수를 이용하여 노이즈 수준을 정의합니다. 따라서 SNR은 피크의 신호 세기와 피크 부분에서의 노이즈 수준의 비로 계산할 수 있습니다.

Identify the peaks based on the ridge lines

주 피크(major peak)를 결정하기 위한 규칙은 아래와 같습니다.

1. 능선에서 가장 큰 값을 갖는 스케일은 피크의 너비에 비례하고, 너비는 적절한 범위 안이어야 합니다.
2. SNR은 적절한 임계값 이상이어야 합니다.
3. 능선의 길이는 임계값 이상이어야 합니다.

일반적으로 주 피크 주변에는 자잘한 다른 피크들이 존재하는데 이는 규칙 3으로 해결할 수 있습니다.

피크의 위치를 찾기 위한 방법으로 2가지 방법을 사용할 수 있습니다. 하나는 높은 스케일에서 낮은 스케일로 능선을 따라가다가, 높은 피크를 보이는 한 작은 스케일을 선택하는 방법이 있고 다른 방법으로는 능선의 적절한 스케일 범위 내에서 가장 높은 CWT 계수를 보고 피크의 중심점을 찾는 방법이 있습니다.

Refine the peak parameter estimation

계산의 편의를 위해서 미리 정해둔 스케일에 대해서 CWT를 계산하도록 되어있지만, 좀더 정확하게 피크의 폭과 중심 위치를 찾고자 한다면 기존의 찾아둔 파라미터를 기준으로 좀더 세세한 스케일에 대해서 계산을 반복함으로써 정확도를 높일 수 있습니다.

3 Results

MS 데이터셋은 CAMDA 2006 데이터 셋을 사용하였습니다.