今年の終わりに

今年の秋から始めたブログですが、私の思っていた以上に見ている方がいらっしゃるようで。
ありがたいことです。

第二世代シーケンサーがようやく日本でも本格的に使われるようになって、来年はきっとどこかで第三世代を使った結果がでてくるんでしょうね。
結果といっても、「やってみた」程度かもしれませんが。
かく言う私も、大いに興味があります。
来年は忙しくなりそう。

さて、今年を締める軽い話題として、日本で次世代シーケンサーの３強のノベルティグッズについて。
セミナーとかに行くと、たいてい、ボールペンやノートがもらえる。　



上から、ロッシュ、ライフテック、イルミナ　のボールペン。

ロッシュのは2色プラスシャーペンで便利そう。　ライフテックのは金属性で高級感あり。　イルミナのは軽くて持ち運びやすい。　

ライフテックのは実は今年のものでは無い。2年くらい前にもらったやつだ。

イルミナのは実はもう一種類あって、そっちはボディが透明で、これもプラスチック製のやつ。　いずれもインクの色は黒のみ。

ちなみに、アフィメトリクス社のボールペンは黒・赤・青の3色。　僕はこれが一番好き。

できればみんな、3色か4色ボールペンにして欲しい。

 

さて、ノートは最近は少なくなった。ロッシュだけが今年配っていた。

 

それ以外のグッズとして、ライフテックはアプライドバイオの時代からぬいぐるみ系が多い。

 

右のくまさんは、2年前にSOLiDセミナーで全員に配っていたもの。
左の犬くんは、今年配っていたもの。
犬以外にもいろいろな動物があった気がする。
いずれにしても、ぬいぐるみ系はおじさんがたがもらってもどうかな～って気がするけど、みんな喜んでもらっていた。　子供や孫にあげるのかな。

さてさて、来年はどんなグッズがもらえるか？　
こんなことも楽しみのひとつ。

ではみなさん、良いお年を！

BWA マッピングツール

BWA，Burrows-Wheeler Alignment Tool　を、ひとから勧められたので試してみた。
このアライメント（マッピング）ツールは、ギャップを許してアライメントすることができる。
また、リード配列は全部が同じ長さである必要が無い。
ということで、SNPなどのバリアントを見つけるにはいいと、そのひとは言っていたのだが、ここは議論の余地が残るところだ。

さて、そんなことはさておき、さっそく実行ファイルをダウンロード
http://sourceforge.net/projects/bio-bwa/
http://sourceforge.net/projects/bio-bwa/files/

インストールするマシンはLinuxの64ビット。
ダウンロードしたら、以下のコマンドを投げていこう！
$ bunzip2 bwa-0.5.9rc1.tar.bz2　これは解凍ね
$ tar xvf bwa-0.5.9rc1.tar
$ cd bwa-0.5.9rc1/　ディレクトリに移動して、
$ make　メークして、
$ sudo apt-get install bwa　スーパーユーザー権限でインストール！

だ、だ、だ、、、、ってインストールが完了したら、準備OK

失敗したら、近くのLinuxマニアに助けを呼ぼう。

ここからが本題。

BWAはBurrows-Wheeler Transform (BWT)を基本としているアライメントツールで、ショートリード向けのものと、ロングリード向けのものの2種類がある。
ショートリード向けは200塩基未満のリード長を3%未満のエラー率で、ギャップを許すグローバルアライメントを行い、ペアエンドにも対応している。
ロングリード向けのものはBWA-SWといい、Smith-Watermanアルゴリズムに似たヒューリステッィクなアルゴリズムで高いスコアを見つけながらローカルアライメントを行う。　

ショートリードでもこのBWA-SWを実行できなくはないが、エラーは多くなる。

・・・なんて、いかにも前から知っていたような書き方をしたが、ここに書いてあるのを訳しただけ。

もっと詳しく知りたいひとは文献までさかのぼってみよう。

さて、ランの流れだ。　こんな感じで進む。

1．データを用意する
リファレンス配列はFASTAフォーマット
リード配列はFASTQフォーマット　今回は例としてイルミナ・ペアエンドのショートリードを使う。

ほかのフォーマットに対応しているかは？

2．リファレンスデータベースの作成
リファレンス配列のデータベースを先に作る必要があり、FASTAファイルに、indexコマンドを使ってインデックスを貼る。

3．アライメントコマンド実行
アライメントはalnコマンドで実行する。Suffix Arrayという検索アルゴリズムの原理を使ってひとつひとつのリードのベストヒットを検索するらしい。

4．SAMファイルへ出力
alnコマンドで作成したsaiファイルを基に、マッピングファイルSAMを作成する。
samse/sampeコマンドはSuffix Arrayのインデックスを染色体の位置情報に変換する。
シングルリードならsamse、ペアエンドリードなら、sampeを実行。

sampeだよ。sampleと打ってしまいそうなので注意！

それでは実例を

まず、リファレンスに対し、インデックスを貼る（ヒト染色体1番のみを例に）
$ bwa index -p ./work/humanCh1 -a is ./Path_to_directory/hs_ref_GRCh37_chr1.fa

重要パラメータ
  -p: アウトプットするリファレンスデータベースの名前
  -a: インデックスの種類を指定　-a is または -a bwtsw
isは簡単で速い。ただしリファレンス塩基数の5.37倍のメモリを必要とし、2GB以上のサイズのデータベースは作れない。
bwtswは、BWT-SWのため。2GBの制限が無いので、ヒトの全ゲノムでデータベースを作ることができる。ただし10MB未満のデータベースは作れないし、速度もISに比べて遅い。

こんなインデックスファイルができたかな？



リファレンスデータセットの名前をhumanCh1にしたので、以下、alnの後のデータベース名はhumanCh1 になる。

リードファイルは、SRR027863_1.fastqとSRR027863_2.fastqで、これはペアエンドのファイルだ。

ひとつずつ、アラインさせる。

$ bwa aln ./work/humanCh1 ./work/SRR027863_1.fastq > ./work/SRR027863_1.sai
$ bwa aln ./work/humanCh1 ./work/SRR027863_2.fastq > ./work/SRR027863_2.sai

アラインして、拡張子saiのファイルができたら、いよいよマッピングファイルSAMの作成。
これはペアエンドなので、sampe

$ bwa sampe ./work/humanCh1 ./work/SRR027863_1.sai ./work/SRR027863_2.sai ./work/SRR027863_1.fastq ./work/SRR027863_2.fastq > ./work/SRR027863_Ch1.sam

ちゃんとパスを通しておけば、もっとすっきりするけどね。
基本は、
$ bwa sampe 「インデックスデータベース」 「アラインファイル1.sai」 「アラインファイル2.sai」 「リードファイル1.fastq」 「リードファイル2.fastq」 > 「マッピングファイル.sam」

これでSAMファイルは出来上がる。1.26Gくらいのでかいファイル。

今回は、シンプルにデフォルトで行った。

全部通しでも1時間くらいで終わった。

気をつけるとしたら、ショートリード用のリファレンスインデックスをつくるには、2GBのサイズ制限があるので、そのままヒト全ゲノムは使えないということ。

染色体ごとに分けるか、目的に応じて、いらない配列を削るかする必要があるかな。

サンプルデータの取得法 ２ NCBI

以前、「サンプルデータの取得法」というタイトルで書いた。
その後すぐに、NCBIのSRA（ここ）にて、データのフォーマットが変更されたので記しておく。

NCBIのSRAに行ってみて、何か検索してみると、例えばSRP000698を検索してみると、右にランの名前がリストされている。

以前はここから、FTPにリンクされて、そこからbzip2で圧縮されたリードファイルがダウンロードできた。
今はちょっと違う。
まず、FTPからは、sraまたはsra-liteという2種類のデータがダウンロードできる。
fastqファイルが欲しければsra-liteで良い。
Rocheの波形データを含むsffフォーマットが必要なら、sraから取得する。

さて、僕はsffはいらない。　fastqフォーマットが欲しい。
ということで、sra-lite　からファイルをダウンロードした。

ところが、このファイルを解凍するには、SRA Toolkit という特別なツールが必要だ。
こちらを参照する　http://www.ncbi.nlm.nih.gov/books/NBK49294/　と、Linux System上で動くとあるが最近Windows版も出たらしい。
ツール自体は、
http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?cmd=show&f=software&m=software&s=software
からダウンロードできる。

僕は64-bit Linux版を落としてきた。
次に、これを展開する。　
$ tar xvfz sratoolkit.2.0b4-3-centos_linux64.tar

ダーっと展開された。

実行コマンドはfastq-dump.2 らしいのでヘルプを見てみよう。

$ ./fastq-dump.2 -h

よくわからない。
さっきのNCBIのSRA　TOOLKITのサイトを見てみる。

早速試してみた。
基本コマンドは　$ ./fastq-dump.2 -A (OutPut fastqファイル名) -D (Downloadした.lite.sraファイル)　なので、
$ ./fastq-dump.2 -A SRR027865 -D /opt/(Downloadしてきた場所)/SRR027865.lite.sra

しばらくして、SRR027865_1.fastqとSRR027865_2.fastqというファイルができるはず。
このSRR027865はペアエンドのリードファイルなので、2つのfastqファイルができるのだ。
シングルエンドのファイルなら、1つのfastqが作成されるはず。

sra-liteからfastqへ変換完了！
ちょっと面倒くさくなった。　慣れればどうってことないが。
データ量が増えたから仕方ないか。
DDBJもいずれこうなるのかな？

BMB2010

分子生物学会が終わり、今年も終わりに近づいてきた。
暮れのこのイベントが、仲間の生存を確認するいい機会になっている。
今年も懐かしい顔に会えて、ほっとした。

昨年と比べ、次世代シーケンサーの結果を披露しているポスターが多かった。
と言っても全体に比べるとごくごく少数だが。　（全体数が多すぎるのか？）

種類別で言えば、まだまだ解析の方法論を述べているものが多い。
アンプリコンのシーケンスで特定配列のタイピングをしたもの。
2種類の生物のゲノムを読んでSNPを比べたもの。
エキソンスプライシングの差を特定していたもの。

個人的には、ゲノム構造（ヒストン構造などのエピゲノム）と遺伝子発現との関連をテーマにしている発表が興味深かった。　何十年も前から関連性は示唆されていたものの、これまでの技術では解明できなかった。　高速シーケンサー技術がこれを可能にするかもしれない。
僕は、1分子シーケンサーの登場こそが、真にシーケンサーによる解析のブレークスルーをもたらすと思う。
これについては、後日、ちゃんとまとめる。

さて、先日の書き込み、マッピングのところ、大切なことを忘れていたことに気が付いた。
僕はそれほど気にしていなかったのだが、ある方と話していて、気が付いた。

マッピングには、Gapを許すアルゴリズムと、Gapを許さないアルゴリズムがある。
Gapを許す方が、計算時間はかかるが、Insertion / Deletion の検出に向いている。
Gapを許さない方は、InDel検出には向いていないが、その分計算時間が短い。
フリーのツールでは、Bowtie、Soapが、Gapを許さない。
BWAがGapを許す。

また、リードをマッピングするときは、マッピングの前に、クオリティでトリミングすることがある。
塩基のクオリティは、リードの後ろほど低い。
50塩基のリードなら、ある程度、例えば40塩基まではクオリティが高くても、その後ガクンと落ちることがままある。
だから、マッピング前に、リードの後ろの方を削り落すことが大事なのだ。
大抵の有償ソフトには、デフォルトでその機能が付いているが、フリーのツールでは自分でパラメータを追加しないといけないことが多い。

一律に、後方10塩基を落とす、という方法なら、全部のリード長は同じになる。
でも、クオリティがXX未満の塩基があったらその後ろを落とす、という方法なら、リード長はまちまちだ。
そんなとき、リード長がまちまちでもマッピングできるのが、BWA。
BowtieとSoapは、リード長が全部同じでないとマッピングできない。
もちろん有償のソフトなら、リード長がまちまちでも問題無い。

さあ、ここまでが、マッピングで書き忘れていた大事なこと。

ところでDDBJ(http://www.ddbj.nig.ac.jp/)でも、シーケンスの登録データベースがあるのをご存じだろうか？
そこの担当者とも話す機会があった。
DDBJにシーケンスデータをアップして、そのまま解析ができるパイプライン(ここ)を作ったそうだ。
今は、マッピングとアセンブルの機能はあるらしい。
もっと広く認知されれば、徐々にサーバーのスペックも上げていくらしい。今はまだ、同時アクセスを制限している状態、ら　し　い。

マッピング、マッピング、マッピング

しつこいけれどまたマッピングの話。

ショートリードをリファレンス配列に対して貼り付けることを、マッピングと呼ぼう。
単に言葉の違いだと思うが、これをアライメントと呼ぶこともある。
はたまた、アセンブルと呼んでいる文献もある。

えっ!?　ショートリードを貼り付けることがアセンブルだって!?

と驚いた方もいるかと思うが、2つ以上の配列の一致する箇所を並べて長いコンセンサス配列を作ることがアセンブル。なら、ショートリードと、リファレンス配列をアラインさせることだって、広い意味でアセンブルだ。　というひともいるのだ。　

ま、それはさておいて、ここではマッピングと呼ぶ。　
正確には、「リファレンスへのショートリードのマッピング」

で、デノボアセンブルを除くと、マッピングが全ての出発点なのだが、このマッピングのアルゴリズム、一体どれだけ正確なのだろうか？
というか、ショートリードのマッピングって、どれだけ信頼できるのだろう。

Illumina社の付属ソフトCasavaで使われているELANDv2にしろ、CLC-BioやNextGENet等有償ソフトにしろ、良く使われているアルゴリズムは、クオリティの高い塩基を十分に含むリード i.e. ちゃんと読めているリード　を、リファレンス配列とのギャップがいくつ以上、ミスマッチがいくつ以上あるものを除いて、残りをマッピングしているに過ぎない。
つまり、ギャップがいくつあったらいけないよ。ミスマッチがいくつあったら同じと認めないよ。という決め事を先に設定し、それに沿ってショートリードの塩基配列がどこにヒットするかを求めているのである。

実に単純なアルゴリズムである！　
というか、これくらい単純にしないと数千万、数億のリード配列を今のコンピュータでは手軽にマッピングできない。
宿命ではあるが、その単純さゆえに、ミスアラインもあるだろう。　
本来マッピングされるべきでない場所にマッピングされる現象は、起こりうる。
それを防ぐには、Paired Endで読むのが良いとか言われるが、これは実験レベルの解決法だ。

ちなみにロッシュ454のロングリードは、僕は経験がないが、聞いた話ではBLATでもマッピングできるそうだ。　BLASTは時間がかかるという。　これは想像がつく。

ショートリードに話を戻すと、ミスアラインの問題は、実はそれほど気にしていない。
ショートリードのマッピングなんてそんなものだ、と思っている。
実験には、そういう覚悟も必要である。
パラメータをいくらいじったところで、ミスアラインが全く無くなることはありえなく、現在のショートリードのマッピングアルゴリズムが単純である限り、その信頼性はBLASTのe-valueほどはっきりと提示してくれない。
ショートリードのマッピングなんてそんなものだ。その代わり、大量のデータをParied Endで出力することで、ミスアラインの可能性を低くする。
しかしデータが大量になればなるほど、マッピングアルゴリズムは単純である必要があるのだ。
堂々巡りになってしまった。

このあたり、ロングリードのロッシュの担当者、ショートリードのライフテクノロジーズ担当者、イルミナ担当者と話をすると、それぞれ全く異なる見解が聞けて面白いと思う。