ピーク検出


R でバイナリを読む』で, R でピーク検出すればベースコールができるという話を書きました。 R には Peaks というパッケージがあり,ピーク検出が容易にできます。これを用いてベースコールを行ってみます。

いま chromat にはシークエンサーが読み取った蛍光強度の時系列データが入力されているとします。ピーク検出は SpectrumSearch で行えます。

# install.packages("Peaks")
library(Peaks)

peak.A <- SpectrumSearch(as.numeric(chromat$A))

ここで, peak.A$pos には,ベクトル chromat$A の検出されたピークの位置 (ベクトルのインデックス) が含まれます。したがって, AGCT のすべての塩基でこれを繰り返し, インデックス順にソートすれば,そのまま塩基配列になります。

詳しい実装はソースコードを参照してください。実質 6 行程度でクロマトグラムから塩基配列決定ができてしまいました。

実際の結果は以下の様になります。

ノイズがあるところではあまりきれいに読めませんが,ノイズが少ない場所ではきちんと塩基が読めています。クオリティは評価できませんが,少ない配列を目で確認しながら修正する,くらいの作業だったら R 上でもできてしまいますね。