A statistical method for chromatographic alignment of LC-MS data (PETAL)

[Pei Wanga] [Department of Statistics, University of Chicago] [Biostatistics(2007), 8, 2, pp. 357–367]

・論文はこちら:

http://med.stanford.edu/tanglab/publications/PDFs/AStatisticalMethodForChromatographicAlignmentOfLC-MSData.pdf)

Peptide Element Alignment(PETAL)についての論文。

1. INTRODUCTION 略

2. LC-MS experiment

2-1. Generic LC-MS experiment 略

2-2. Existing alignment methods:アラインメントのゴールは、RTと実験のノイズの表出において異なる実験のmz-scan plotにおいて、対応するペプチドを一致させること。Bylund(2002)はtime warping method based on raw spectrum for alignment of LC-MS data, which is modification of the original correlated optimized warping(COW) algorithmを提案した。1つのファイルをテンプレートとして選び、RT axisを2つのイメージの間の類似性を最も高めるように2つのファイルを合わせる(warp(たゆませる))。これはdynamic time warping algorithm allowing every RT point to be movedでも使用された。しかし、古典的な1次元のchromatography profilesと比較して、LC-MSデータは質量情報も含まれているため、複雑である。言い換えると、1つの実験において同時に溶出した2つのペプチドは他の実験においてはそうとは限らない。したがって、2つのファイルのRTをmappingするだけでは各々のペプチドをアラインメントするには十分ではない。

raw spectrum dataを使う代わりに、Radulovic(2004) performed alignment based on the (mz, RT) values of detected features. mz domainをいくつかのintervalに分けて、fitted different piece-wise linear time warping functions for mz interval. その後、異なるファイルの最も隣接したピークを一致させるためにa "wobble" functionをトータルscan rangeの±1~2%内で移動することを許可する。ここに、mzの層は柔軟さと正確さを改善する。この方法論は探索されたペプチドの特性値(mz,RT)だけに依存するので、同位体分布のような元画像を利用することは出来ない。加えて、ヒトの血清のような複雑な形状を処理する際(区間の±1~2%以内に複数のペプチドが存在する)、"wobble"functionは多義の解釈をすることになるかもしれない。

msInspect(Bellow 2006), SpecArray(Li 2005), MZmine(Katajamaa&Oreic 2005)らは検知した形上リストの範囲で個々のペプチドのRTにバリエーションを許可するようなアラインメントソリューションを提供している。しかし、これらのmethodはpeak detectionの結果に頼っており、生のスペクトル情報を利用していないため、アラインメントの決定はpeak detectionのステップにおける不正確な見積もりに対して脆弱である。大半の方法はお互いにとても良く似たデータでは最もうまくはたらく。これらの方法は癌/非癌組織のような異なるdisease classのサンプルを解析するとき、バイアスを発生させがちである。

他にも関連するアルゴリズムとしてMALDI/SELDIスペクトルに対する階層的クラスター法を含むListgarten and Emili(2005)、a multi-scale wavelet decomposition approach data along the mz axis(Randolph&Yasui)、隠れマルコフモデル for multiple alignments of time series(Listgarten 2005)がある。Prakash(2006)はsignal mapping algorithm to perform comparison directly on the signal level of mass spectrometry experimentsを発表した。

現在の方法論の欠点に打ち勝つために、我々はPETALを提案する。これはペプチドのアラインメントのdetected peak features、raw spectrum dataのどちらの情報にも使える。

3. Peptide Element Alignment for LC-MS:LC-MS profilesにおいては、各々のペプチドは次の2点によって特徴づけられる。質量スペクトルとRTのレンジである。1スキャンの質量スペクトルはmz軸に沿ったIntensityの測定を記録したベクトルである。どのようなペプチドであっても、そのペプチドの1組の存在度だけを含む質量スペクトルを実験ノイズを伴わず定義される(Σl=1~L bl =1, Lはmz軸の合計ピクセル数、blはl番目のピクセルのintensity valueである)。あるスペクトル要素ベクトルbは対応するペプチドの同位体の位置とパターンによって一意的に決まりうる。加えて、我々はあるペプチドの理論的なRTレンジをRT = (rtbegin, rtend)のように表す。k番目のペプチドはPEPk=(bk, RTk)のように表すことが出来る。

我々はペプチド要素のライブラリー{PEPk}k=1~Kをターゲットサンプルにおいて起こりうる全てのペプチドの集合として定義する。ペプチド要素のライブラリによって、アラインメントのゴールは、各々のプロファイルのピーク特性を共通のライブラリと一致させることによって簡単に達成出来る。同じペプチド要素と一致した異なるプロファイルのピーク特性は同じペプチドの特性を表し、そのようにalignされなければならない。我々は損失関数を導入し、最適化問題を解くことによってalignmentの解を探す。

3-1. Loss Function:まず、1回のスキャンの質量スペクトルを考える。H個の異なるペプチド{PEPkh}h=1~Hが1回のスキャンで溶出したとする。そして、PEPkhの測定可能な存在度がβkhとする。完全な質量スペクトルは個々の溶出したペプチドの合計である。観測質量スペクトルYは、H個のペプチドのスペクトル要素ベクトルの線形結合で表すことが出来る。Y = Σh=1~H βkh・bkh + ε、εは実験機器のノイズ、bkhはPEPkhのスペクトル要素ベクトルである。

実際は、我々は観測されたYにおいてどのペプチドが溶出したのか知ることはない。しかし、ペプチド要素ライブラリ{PEPk}k=1~Kによって、我々はL1 penalized 最小2乗回帰モデル Y~ Σk=1~K βk bkにfitさせることによって、ペプチドの存在度を測定することが出来る。

{βk}k=1~K = argmin{βk}k,βk≧0 || Y - Σk=1~K βkbk||^2 + λ1 Σk=1~K w(t;RT^k)・|βk| (1)

λ1は非負パラメーター、tはscan YのRT、wは各々のペプチドのRTkの理論値だけでなく、scanのRTtに依存したweight functionである。weight function wは理論的なRT tがscanのRT tから大きく離れているペプチドに対して大きなペナルティを与える。そのとき、対応する予測値bkが頻度低く選ばれる。

(1)を解く際、βk(≠0)はペプチド要素ライブラリにおいてk番目に一致するYのシグナルの一部を表している。したがって、2つのスキャン(Y, Y' from two different profiles)に対して適した回帰モデルを得ていれば、これらのアラインメントは{βk}kと{β'k}kを比較することで行うことが出来る。

N個のプロファイルがあり、各々にMnのスキャンがあるとする。ペプチド要素ライブラリ{PEPk}k=1~Kが与えられたとすると、

{βmn,k}n,m,k = argmin{βmn,k:βmn,k>0} Σn=1~N Σm=1~Mn {||Ynm - Σk=1~K βmn,k bk||^2 +λ1 Σk=1~K βmn,k w(tmn;RTk)} (3)

を満たす{βmn,k}n,m,kを求めればよい。大半のケースではペプチド要素ライブラリ{PEPk}k=1~Kはavailableではない。我々は最も良く観測された全ての質量スペクトルを説明する{PEPk}k=1~Kを同定する必要がある。したがって、我々はすべてを含めた損失関数を定義する。

L~ = Σn=1~N Σm=1~Mn [ || Ynm - Σk=1~K βmn,k bk ||^2 +λ1 Σk=1~K βmn,k w(tmn; RTk) ] +λ2K (4)

and search for 

{{βmn,k}n,m, PEPk}k = argmin{βmn,k:βmn,k≧0;PEPk} L~ (5)

λ2Kは全体の複雑さをコントロールするペナルティ項。λ1とλ2の選択はsection5で議論。

個々のペプチドによるRTのランダム性を除いては、LC-MS実験には構造的なRTのシフトがある。 したがって、我々は構造的傾向を調整するために全体の変換を行い、wを計算するための修正時間を使う。

3-2. Optimization Strategy:損失関数(5)から、もし{PEPk}kが与えられれば、最適解{βmn,k}をlasso/larsのようなL1-regression techniquesにより簡単に計算することが出来る。したがって、我々の主たる障害は適切なペプチド要素ライブラリ{PEPk}kを見つけることである。

要素ベクトルの全体空間を探すことは難しい。それゆえ、次の2stepをとる。まず、アラインメントに従った全てのプロファイルに基づいてペプチド要素のinitial collectionを構築する。これは実験における全てのペプチドを表すことが期待される。しかし、不要で不正確な要素を含むと思われる。その結果、我々は損失関数を最小化するようにinitial collectionの部分集合を探す。

3.2.1 Initial collection:initial setを構築するため、我々はプロファイル中で見出された各々のピーク特性のために1つのペプチド要素を含める。そのため、我々はまず理想的な同位体形上を見積もる必要がある。与えられた質量において、アミノ酸の配列の際から生じる同位体形上のバリエーションは実験ノイズによるバリエーションより小さいため、似た質量値と同じcharge statusにあるペプチドは同じ同位体パターンをもつと仮定する。経験的な同位体形上により、我々は測定されたmono-isotope positionと検出された特性のcharge valueに基づいてスペクトル要素bkを作成する。我々はこの集合をΩ0とする。

3.2.2 Subset selection:部分集合の選択には大きく2つの戦略がある。forward-stepwise and backward-stepwiseである。この選択においては、我々はbackward-stepwiseにfocusする。Backward-stepwiseはΩ0から始まり、不要な要素を繰り返し削除する。前述のように、各々のペプチドは異なるプロファイルに寄与されるinitial collectionにおいて1つ以上のペプチド要素に対応する。もしクラスターを適切な方法で形成出来れば、同じペプチドを表す要素をまとめ、無駄を省くことが出来る。

この目的のため、我々はElastic net(Zou&Hastie 2005)と呼ばれるスパース回帰を適用する。これは、損失関数 L(λ1,λ2,β) = |y - Xβ|^2 + λ2|β| + λ1|β|を最小化するものである。ridgペナルティはグループ化効果がある。backward-stepwiseは以下の通り。

1. 全てのM個のターゲットプロファイルの特性スキャンをとり、initial elementsΩ0とする。|Ω0| = N0

2. For j = 1 to Mに対して、Yj ~ {1/d(RTj, RTk) bk}k=1~N0 with a fixed number of maximam steps。つまり、各々の要素はM個の係数をM個のRegression modelsから得る。※1で"似た質量"で"同じcharge"であるペプチドは同じ同位体パターンを持つと仮定しており、N0<Mであるので、N0個のベクトルでM個のターゲットを表現することになる。

3. the coefficient vector (length of M)に基づいたクラスター要素(同じペプチドを表している要素)をひとまとめにする。N1個(<N0)の要素を得る。

4. 2-3を繰り返し、Nk=Nk-1となるまで続ける。

 

■PETALのアルゴリズムの理解(13/05/21現在)。

1. 全サンプルを重ねて、「似ている」ピークのみを抽出してイニシャルリストとする(N0個)

 2. 各測定毎に、ピークをイニシャルリストのピークに照らし合わせる。

照らし合わせは、各RT毎にmzの最小単位毎に行って、regressionモデルで近似。RTのずれが一定以上(RT±δの外)であるピークはペナルティを∞とし、アラインメントの基準から除外。

 3. 全測定データを用いてクラスターを作成。2で検出された「同じピーク」毎に、Intensityの値でクラスターを作成し、N1(<N0)個にまとめる

 4. 3が減少しなくなるまで、2,3を繰り返す。残ったピークでアラインメント実行。