『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 上でもできてしまいますね。