2010年11月28日日曜日

こんなブログでも・・・

毎日50人くらいが訪れている。
そしてその数は右肩上がりで増えている。
話題はかなりマニアックなのに。 画像も少ないし。

これからももっと続けようと感じている。
ターゲットは 「次世代シーケンサーを学び始めた学生」 を意識して、わかりやすく書こうと思う。
だから、厳密には正しくない ところもあるかと思う。
ご了承願いたい。
コメント・意見も歓迎する。



 

2010年11月26日金曜日

マッピングツールの種類

ショートリードの解析は、つまるところ
1.リファレンスが無いデノボアセンブリ と
2.リファレンスがあるそれ以外
に分かれる。
1は、シーケンサー本来の使い方であり、王道なのだが、僕はちょっと退屈だ。 いまいちゲノムプロジェクトの魅力が感じられないのは、僕がまだ半人前の証拠か。

それ以外の解析、トランスクリプトームやエピゲノム、SNPやAllelic Imbalance などの方が興味を惹かれる。
これらの解析では通常リファレンスがあるので、マッピングが行われる。以前紹介したBowtieは、僕のお気に入りのマッピングツールだが、それ以外にも、BFASTやSOAPなどが良く論文に現れる。
けれどもやっぱり、ユーザーが最初に使ってみようと思うのは、シーケンサーについているソフトのマッピングツールではないだろうか。

SOLiDなら付属ソフト「BioScope」のGUIから選べる "Map Data" とか "bFast Map Data"
Illuminaなら同じく付属ソフト「CASAVA(キャッサバと発音)」の "Eland"

シーケンサーについているソフトなので、配列データをベースコールした後、そのままマッピングできるのが便利…だろう。 いずれにせよ、マッピングは精度と速さ、目的に応じた使い勝手、で選ぶべきなので、
1.シーケンサー付属のマッピングツール
2.Bowtieなどのフリーのマッピングツール
3.CLC-BioやNextGeneなどの有償パッケージソフトのマッピングツール

を試すのが最も良いと思う。
マッピングのアルゴリズムは、今やどれも遜色ない。というか、精度の差は優越付けがたい。
僕の場合、マッピングの目的は、短い配列をできるだけユニークにアラインさせることにある。
大抵の場合、ヒトゲノムに対して、なので、リファレンス情報に不自由はない。
ターゲットを決めたディープシーケンスのときがほとんどだが、たまにホールゲノムに対しての時もある。 
マッピングの後の解析は、SNPやInDelの検出、新規転写産物の検出、たまにChIP-Seqやメチレーション、だ。
転座解析は、まだやってないが、興味はある。

これらに必要なマッピングアルゴリズムは何だろう?
短い配列だから、ユニークにアラインさせることは難しいか?

短いと言っても、リード長、今や40、50は当たり前。 Pair Endで実験すればForwardとReverseの間のフラグメントの長さが決まっているから、その情報を頼りにユニークにアラインさせることは難くない。 以前、メーカーの人に聞いてみたが同じような答えだった。

読んだ配列に、ギャップがあった場合、リファレンス配列にちゃんとマッピングされるか?
大丈夫、短いInDelならほとんどのツールが判別してくれる。
短いInDelとは例えば、IlluminaのCASAVAのEland v2 では20塩基未満のギャップなら認識してくれる。v2 以前はギャップは見落としていたらしい。

短い配列のマッピングには、ローカルアライメントが使われ、例えば連続する12塩基が完全に一致する箇所をリードの中に先ず見つけ、その後両端を伸長してアライメントさせるアルゴリズムがある。 
この方法は昔からよくChIP-Seqの論文で見つけた。 
伸長していくとき、同時に、ミスマッチの箇所、1塩基ギャップの箇所を見つけ、あればこれらでアライメントの精度を評価する。
評価にはスコアリングが使われる。 わかりやすく例えると、最初、50ベースのリードが、50点持っていたとする。 ミスマッチがあればマイナス3点、ギャップスがあればマイナス1点、などと50から減点していく。 ある程度の点数以下になればこのリードは信用なし、としてアライメントに使われない。 
正しくはツールによってスコアの種類も異なり、最初の12塩基完全一致という前提も異なり、計算式もそれぞれだ。 でもわかりやすく言うとこんな感じで、スコアリングについては、どのマッピングツールも大体似たような考えを用いて精度を上げている。

速度はどうか? 僕の体感ということを断わっておくと、
Bowtieはダントツに速い。
次にCLC-Bioが速い。
もちろん、他のフリーソフトも、条件を最適にすればそれなりなのだろうが、それを怠っているせいで、あまり速くは感じない。 
シーケンサー付属ツールは、実はガッツリ使っているわけではないのでわからない。

使い勝手はどうか?
これはその人それぞれ。
シーケンサーが外部のラボにあって、実際に触れない人なら、シーケンサー付属のソフトという選択肢は先ず消える。ということなら、フリーツールか有償パッケージソフトになる。
有償パッケージソフトが安定しているのは当然だ。
コマンドラインに格闘する楽しさは味わえないかもしれないが、SAMファイルを出してくる所までは簡単にできる。 パラメータ設定も十分種類がそろっている。 アルゴリズムもほとんど世界共通なものを使用しているので、心配はいらないと思う。
フリーツールも、コマンドラインを使って見せて、頭よさそーに自己陶酔できる、という利点もさることながら、SAMファイルまではちゃんと出してくれるだろう。 Linuxに詳しければシェルを組んでバッチで流すということもできるツールが多い(つまりマッピングの自動化)。

SAMファイルまたはBAMファイルまで出てくれば、マッピングツールとしての役割は終わりだ。
このファイルはその後のトランスクリプトーム解析、SNP解析、エピゲノム解析のスタートポイントになる。 まさに、マイクロアレイで言えば、ノーマライズ後のデータ、に似てはいないか。

2010年11月19日金曜日

情報革命

今日、夜テレビを付けたらソフトバンクの孫さんが出ずっぱりだった。
なんでも、NTの光回線部門の別会社化案を総務省が見送る判断を示したらしく、その反論と自身の経営者としての考えを述べていた。

さっきもWBSで、「これから10年20年後、携帯端末の性能は今の100万倍になる。 想像もつかないほどの情報革命、通信革命が起こり、それは体感しないとわからない。」と言っていた。
孫さんの言うことは正しいと思う。
携帯端末の性能は、その心臓部、エンジンであったりメモリであったりオペレーションソフトであったりするのだが、これが20年後には今の100万倍になるという。

DNAシーケンサーの性能は、シーケンスの読み取り技術と、計算機の性能と、ベースコール解析ソフトによって出される、配列の長さと数とコストで計算できると思う。
ほんの10年前、国際ゲノムプロジェクトV.S.セレラジェノミクスとの10年にも及ぶ熾烈な競争の結果、ヒトゲノムがなんとか読まれた。 この時使われたテクノロジーはサンガーキャピラリーシーケンサー。 使われた資源はどれだけのものだったのか。

キャピラリーで読める長さは大体500塩基としよう。16本キャピラリーだと、1ランで8000塩基。
SOLiD 5500xlのアップグレード版は、ショートリードとはいえ、1ランで300ギガ=3000億塩基
その差なんと 3750万倍 !!
そしてスループットはまだまだ上昇の余地があるという。

孫さんが情報革命に興奮しているのと同様、ライフサイエンスをかじっている僕も、現在のシーケンス革命に興奮せざるを得・・・ っていうのは言いすぎか。
でも、将来、個人ゲノムが診断目的に読まれるようになって、そのシーケンス情報をケータイできるようになって、それらが世界的につながって、自分と同じジェノタイプを持つひとがどれくらいいるかとかがわかるようになって、どんな病気に気を付けたらいいかとか、この薬は合わないとか、そんなことが瞬時にわかるような時代がすぐそこに来ている気がする。

たまにはこんな話も良いでしょう。

2010年11月15日月曜日

シーケンサーの種類 SOLiD (2)

今月の1日、ライフテクノロジーズ社はプレスリリースを出した。
http://www.lifetechnologies.co.jp/pr101102.html
それによると、SOLiDシリーズは、近々また一新するらしい。

今までのSOLiD 4 は、そのスループットをさらに向上させ、SOLiD 5500xl に、さらに当初リリース予定だった低コスト版のSOLiD PI は、名前を変えて SOLiD 5500 になるらしい。
5500xl は、2プロ―チップまでのランが可能。 5500 は、1フローチップのランが可能。 って・・・。

で、そのスループットだが、ここにいい表がある。 ちょっと見て欲しい。

コスト/ランは、新しい5500xlシステムの方がSOLiD 4 より1.7倍以上高い!
今は180ギガベース/ランで、7500円/ギガベースだけど、将来的には(ソフトウェアの改善で)300ギガベース/ランになるらしいから、4700円/ギガベースか。 これはすごいかも。
でも、ソフトウェアの改善というところがちょっと気にかかる。 

5500/5500xl シリーズの機械は日立ハイテクノロジーズ社製。 
ということで、ハードについては信頼性があるかな? 

もうひとつ、これからの注目株は、今年の分子生物学会のランチョンセミナーで紹介される予定の、Ion Torrentシーケンサー。 蛍光色素を使わず、CCDカメラも無く、ヌクレオチドが取り込まれる時の水素イオンのリリースだけを検出してシーケンスする、この方法は画期的だ。
でもまだ製品の詳細は明らかにされていない。 
ここから先は想像だが、この新型シーケンサーは診断用シーケンサーの世界標準を狙っているのではないか。 低価格、低コスト/ラン、短いラン時間、レーザープリンターくらいの小さいサイズ、次世代の中では少ないデータスループット、長い読み取りリード長。 
どれをとっても大量スループットが必要無いラボにはちょうど良い。
あとは、使い勝手の良さ、データの精度、プロトコールの安定性、などが議論になるだろう。
いずれにせよ、まだサンプルデータなるものが無いので、これ以上のことは言わない。

それにしても、SOLiDを導入したラボは、しょっちゅう出されるアップグレードに苦労しているらしい。

2010年11月10日水曜日

SAMファイルをとりあえずビューワーで見る

Bowtieなどのマッピングツールでリード配列をリファレンス配列にマッピングすると、SAM (Sequence Alignment Map) というフォーマットのファイルができる。
SAMは実際、マッピングファイルの世界標準になりつつある。
サイズは、かなり大きい。 中身を見てどうこうということは無いけど、先ず、リファレンスに使った配列の情報がずらっと並ぶ。 @SQ SNから始まり、LN:の前までの、例えばgi|224589801|ref|NC_000010.10|が、リファレンス配列の名前だ。ちなみにこれはヒト10番染色体。
@HD VN:1.0 SO:unsorted
@SQ SN:gi|224589801|ref|NC_000010.10| LN:135534747
@SQ SN:gi|224589802|ref|NC_000011.9| LN:135006516
@SQ SN:gi|224589803|ref|NC_000012.11| LN:133851895
@SQ SN:gi|224589804|ref|NC_000013.10| LN:115169878
@SQ SN:gi|224589805|ref|NC_000014.8| LN:107349540
@SQ SN:gi|224589806|ref|NC_000015.9| LN:102531392
@SQ SN:gi|224589807|ref|NC_000016.9| LN:90354753
@SQ SN:gi|224589808|ref|NC_000017.10| LN:81195210
@SQ SN:gi|224589809|ref|NC_000018.9| LN:78077248
@SQ SN:gi|224589810|ref|NC_000019.9| LN:59128983
@SQ SN:gi|224589800|ref|NC_000001.10| LN:249250621
@SQ SN:gi|224589812|ref|NC_000020.10| LN:63025520
@SQ SN:gi|224589813|ref|NC_000021.8| LN:48129895
@SQ SN:gi|224589814|ref|NC_000022.10| LN:51304566
@SQ SN:gi|224589811|ref|NC_000002.11| LN:243199373
@SQ SN:gi|224589815|ref|NC_000003.11| LN:198022430
@SQ SN:gi|224589816|ref|NC_000004.11| LN:191154276
@SQ SN:gi|224589817|ref|NC_000005.9| LN:180915260
@SQ SN:gi|224589818|ref|NC_000006.11| LN:171115067
@SQ SN:gi|224589819|ref|NC_000007.13| LN:159138663
@SQ SN:gi|224589820|ref|NC_000008.10| LN:146364022
@SQ SN:gi|224589821|ref|NC_000009.11| LN:141213431
@SQ SN:gi|224589822|ref|NC_000023.10| LN:155270560
@SQ SN:gi|224589823|ref|NC_000024.9| LN:59373566
そしてBowtieのパスが書かれていて、そのあとから、各リード配列のマッピング情報がダーッと書かれるわけ。
1279_2_116_F3 4 * 0 0 * * 0 0 CCGTACGTCGTGGGTAGGGGCNGTGAGTCGCTTCGGGTCGAGGATCTGG =;?89,:4:=2,57-&*'&&/!,)&/.&),49,'1)-)'(,12+//-7) XM:i:0
これは人間が読むわけでは無いので、どうしても意味が知りたいひとはここを参照。
http://samtools.sourceforge.net/SAM1.pdf

さて、このファイルをとりあえずフリーのビューワーで見てみよう。
例として、IGV (Integrative Genomics Viewer)で見てみる。
登録(無料)してからダウンロードする。これはWindowsで動くので敷居が低い!

ところでさっきのSAMファイル、リファレンスの名前がgi|224589801|ref|NC_000010.10| のようになっているが、このままではIGVが認識しない。
ビューワーが認識できる染色体名に変換してあげなくてはいけない。
IGV なら、gi|224589801|ref|NC_000010.10| は、10 に変換する。以下24本の染色体も同じ。
やり方は色々あるだろうが、一番簡単なのは、Perlスクリプトを使って、
perl -pe "s/gi\|224589801\|ref\|NC_000010.10\|/10/g" SD_Agilent_Exome_F3.sam > temp1.sam
perl -pe "s/gi\|224589802\|ref\|NC_000011.9\|/11/g" temp1.sam > temp2.sam
perl -pe "s/gi\|224589803\|ref\|NC_000012.11\|/12/g" temp2.sam > temp.sam
perl -pe "s/gi\|224589804\|ref\|NC_000013.10\|/13/g" temp.sam > temp1.sam
のような置換コマンドを24本分繰り返すと、最終的に全ての染色体名をIGV用の染色体名に置き換えることができる。 こんな感じになった。

@SQ SN:10 LN:135534747
@SQ SN:11 LN:135006516
@SQ SN:12 LN:133851895
@SQ SN:13 LN:115169878
@SQ SN:14 LN:107349540
@SQ SN:15 LN:102531392
@SQ SN:16 LN:90354753
@SQ SN:17 LN:81195210
@SQ SN:18 LN:78077248
@SQ SN:19 LN:59128983
@SQ SN:1 LN:249250621
@SQ SN:20 LN:63025520
@SQ SN:21 LN:48129895
@SQ SN:22 LN:51304566
@SQ SN:2 LN:243199373
@SQ SN:3 LN:198022430
@SQ SN:4 LN:191154276
@SQ SN:5 LN:180915260
@SQ SN:6 LN:171115067
@SQ SN:7 LN:159138663
@SQ SN:8 LN:146364022
@SQ SN:9 LN:141213431
@SQ SN:23 LN:155270560
@SQ SN:24 LN:59373566
さて、そうしてSAMファイルをIGVでも認識できるようにしたら、インストールしたIGVを立ち上げてみよう。
そうしたら、次に、File > Run igvtools... で IGV tool を起動する。
そこで先ず、さっきのSAMファイルを染色体ごとにソートする。 CommandでSortを選び、リファレンスの名前をIGV用に置き換えたSAMファイルを指定して、Runする。

その次に、インデックスを作成する。同じくIGV tool上にて、CommandでIndexを選び、今ソートした結果ファイルを指定する。
そうすると拡張子 sai のファイルができる。
ここまできたら、さあ、IGVで取り込もう!

File > Load File で、ソート済みのSAMファイルを指定すると、自動的にインデックスファイルも認識されて取り込まれるはずだ。
しかし、最初は何も見えないだろう。 こんな風に
染色体を適当に選んでから、ズームインしてみよう!

こんな感じに見える。

まとめ: SAMファイルは、ちょっとリファレンスの名前を変更するめんどくささはあるけど、フリーのビューワーでちゃんと見ることができる。

ビューワーについてはまた今度。