Trinity: de novo Transcriptome の救世主？

前回の最後にチョイ出ししたTrinityというソフト、使ってみました。

http://trinityrnaseq.sourceforge.net/
Broad Instituteで開発されたこのソフトは、Perlで動かして使い、Linux上でテストされているそうです。
先のサイトからダウンロードしてインストールすると、

こんなディレクトリができあがり、/trinityrnaseq/sample_data/test_Trinity_Assembly の下に、サンプル用のペアエンドデータが作られているはずです。

とりあえずどんな結果が出てくるのかが待ち遠しいひとは、まずこのサンプルデータをランしてみましょう。　runMe.sh を打つようにマニュアルでは言っていますが、あえて（わざわざ）このシェルの中のコマンドを打ってみました。　（以下、一行です）

Trinity.pl --seqType fq --left reads.left.fq --right reads.right.fq --SS_lib_type RF --paired_fragment_length 280  --min_contig_length 305 --run_butterfly  --CPU 2 --bfly_opts "-V 10 --stderr"

• --seqType fq で、このリードがFastqだということを示し、
• --left / --right で、ペアエンドの配列ファイルを指定し、
• --SS_lib_type RF で、このリードがストランドスペシフィック（日本語で何と訳すんだろう）でリバース・フォワードの向きだと宣言し、
• --paired_fragment_length 280 はフラグメントの長さが280ベースだと言い、
• --min_contig_length 305　で、コンティグは最低305塩基必要とし、
• --run_butterfly --CPU 2 で、一緒にbutterfly もランして、その時はCPUを2個使えと言い、
• --bfly_opts "-V 10 --stderr"　これは何かよくわかりません

ちなみにシングルエンドリードのときは、--left / --right　を　--single にすればOKです。
さて、このコマンドを流すと、trinity_out_dir　というディレクトリが作られて結果が保存されますが、--output で出力先を決めることもできます。

結果で大事なのは、この中でTrinity.fastaです。　これが最終的なコンティグです。

これがそのファイルを開いたところです。　名前（ID）と、Contigの長さ、リードをマッピングした時のFPKM、がヘッダーに書かれているMulti-Fasta形式ですね。

 

Trinityは、Inchworm（シャクトリ虫）、Chrysalis（サナギ）、Butterfly（チョウ）という名前の3つに分かれています。

Inchwormでは、k-mer ベースでの de Bruijn graph でアセンブります。　Alternative Variantがある時でも、そのうちの代表的な転写産物をつくるようにアセンブるので、完全なTranscriptome、Isoformを反映したTranscriptomeができるわけではありません。　

デフォルトは k=25 です。　k-1 の重なりでContigを伸長していき、これ以上伸長しなくなったら止まるそうです。

Inchworm.K25.L48.fa という名前のファイルができます。　

メモリが足りないとこのファイルすらできません。

次のChrysalisは、名前の通り、動いているのか止まっているのかわからない、ユーザーを不安にさせるコマンドです。

冗談です。

Inchwormで作られたInchworm.K25.L48.fa のコンティグたちから似ている配列のContigを集め、Alternative Splicingや、Paralogっぽい配列をまとめる、極めて重要な作業をしている、そうです。

k-1 の重なり具合をもとに、Inchwormで出来上がったコンティグで、de Bruijn graph を作っていくつものパスを作るらしいです。

最後のButterflyは、いかにもチョウらしく、何百行ものコマンドラインが現れては流れていきます。　ChrysalisでまとめられたIsoform/Paralogのコンティグをもとに、リード配列を使って、あいまいさの残る部分を修正し、全体を再構築しているそうです。　

この過程でFPKMも算出されます。

文献によると、ベンチマークでは、ABySS、TransABySS、SOAP denovo、と比べると、完全長の遺伝子、Isoformの数ともにTrinityは優れているそうです。

しかーし、

このツールは恐らくショートリード向けでしょう。　k-mer が25と短いのがその理由です。

また、デモデータではなく、本番データをアセンブるには、結構な量のメモリが必要です。

76塩基程度のイルミナペアリード100万本あたり、1GB のメモリが必要だそうです。

とうことは、1,000万本で10GB、1億本で100GB !!!

彼らのお勧めは、256GB の、22コアCPUのサーバーとのこと。

とてもじゃないですが、常識的なお値段で買える代物ではありません。

と、いうことで、私は多分、リード数が2,000万本を超えるような de novo Transcriptome でこのツールを使うことはしないでしょう。　

分散化などをうまく施せば別ですが。

de novo Transcriptome; 454用のベストなアセンブラーはどれだ !? 【論文紹介】

Roche 454 と言えば、パイロシーケンスでロングリードを読める代表格です。
数百塩基も読めるロングリードは、未知の配列決定にも良く使われます。
de novo のシーケンスですね。
以前、de novo Transcriptome を話題にしましたが、リファレンスが未知の生物で、転写産物配列を決定したい場合、ロングリードの454が、真っ先に使うシーケンサーの候補に挙げられるでしょう。
もちろんショートリードのペアエンドで読む、ということも可能です。
しかし、長く読めるということは、未知の配列決定において非常に大きなアドバンテージですので、454を選ぶ方は多いのです。

ゲノムではなく、転写産物を読む場合、Isoformの存在が気になります。
スプライシングのバリアントを見たい場合、やっぱりロングリードが有利でしょうか。

実は、リファレンス未知の生物のTranscriptomeを行っている研究は、結構あります。
論文になっているものだけでも、昨年紹介した、ほかにも、

Coral larval (サンゴ): Meyer et al. BMC Genomics 10, 219 (2009).
Eucalyptus grandis (グランディスユーカリ): Novaes et al. BMC Genomics  9, 312 (2008).
Sarcophaga crassipalpis (ニクバエ) : Hahn et al. BMC Genomics 10, 234 (2009).
Populus trichocarpa　(ブラックコットンウツド): Geraldes et al. Mol.Ecol.Resour 11, 81 (2011).

なんかがあるようです。　（全部は読んでいませんが）

454のリードは一般的には、シーケンサーの機械に付属しているアセンブラー、Newblerでアセンブルすることが多いと思います。
でも、Newbler以外にも、ロングリードをアセンブルできるアセンブラーはいくつかあるんですよ。
そこで、454のリードをアセンブルできるアセンブラーを比較した論文を紹介します。

Kumar et al. Comparing de novo assemblers for 454 transcriptome data. BMC Genomics 11, 571 (2010).

この論文は2010年に出されていますから、今年は各ソフトもバージョンアップして、論文投稿時とは若干状況が違っていると思いますが、参考になると思います。

比較しているアセンブラーは、

1. Newbler 2.3
2. Newbler 2.5
3. CAP3
4. CLC Assembly Cell 3.0
5. MIRA 3.0
6. SeqMan NGen 2.1

の6つです。　CLCとSeqManは商用ソフト、Newblerはアカデミックはフリー、CAP3は非営利ならフリー、MIRAはどこでもフリー、だそうです。　
結論からいうと、Newbler 2.5が最も良いとのこと。
「454のシーケンスメーカーが作っているんだから当たり前だろ！」
私も思わず突っ込みましたが、ほかにも面白いことが書いてあるかと。

アルゴリズムの違いとして重要なところは、
CLCはde Bruijn graphを、その他の5つはOverlap-Layout-Consensus (OLC) を使っていること。

de Bruijn graphはVelvetやABySSといったアセンブラーでも使われています。
de Bruijn graphの特徴は、リードをk-merという決められた長さの塩基ブロックに切り、この部分の重なりをもとに、アセンブルをします。
k-mer = 31なら31塩基の重なりを見て、リード同士をつなげていきます。
k-merの範囲でのみ重なり具合を見るんですね。　本当はショートリードの、たくさんカバレージがあるデータに向いています。
一方OLCはクラシカルな方法？で、ペアワイズアライメントを元にしています。
より慎重に伸長していくのですね。
もちろん、6つのアセンブラーはそれぞれ、計算アルゴリズムを工夫しているので、結果はちがうのですが。

この論文では、彼らは線虫の一種をサンプルに、Roche 454 FLX を使って de novo Transcriptome をやっていました。
アダプタートリムした　741,387本のリード、約2億塩基のデータを、先の6つのアセンブラーでアセンブルしています。

最も速くアセンブルが終わったのは、CLCで4分、次がNewbler 2.5の45分。　一番遅いのがMIRAの3日。
1kb以上の長さのContigが一番多くできたのは、Newbler 2.5で7,661本、次がNewbler 2.3の6,320本。　一番少ないのはCLCの4,174本。
全Contigの合計塩基数は、一番多いのがSeqManの2,136万塩基、次がMIRAの2,134万塩基、3位はNewbler2.5の2,007万塩基。　最も少ないのはNewbler 2.5の1,446万塩基。
（詳細は同論文のTable4を参照）

Contigの合計塩基数でNewblerのバージョン2.5と2.3で成績が全然違うのにはびっくりします。
2,007万塩基と1,446万塩基ですからねえ。

6つのアセンブラーでできたContigをそれぞれBLATして、相似性・非相似性を調べ、新規配列がどれだけ作られたか、を彼らは次に見ています。
結果、MIRAとNewbler2.5は、より多くの（長い）ContigをCLCやCAP3より作ったが、その余分の配列は特に新規の配列ではなかった。　曰く、MIRAとNewbler2.5のContigには冗長性があるのではないか、と。
CLCのcontigは最も冗長性が低かったらしい。　ということは余分な配列は作らなかった。　
「待てよ、それって良いことなのか？　TranscriptomeではIsoformはつきもの。」

繰り返しですが、
de Burijn graph は高カバレージの大量リードアセンブルに向いています。　なぜならアセンブルがめっちゃ速いから。
OLCは低カバレージのアセンブルに向いています。　なぜなら正確だけど、ペアワイズアライメントなので計算量がかかるから。

私の経験でも、CLCの de Burijn graph アセンブラーは他と比べて速いです。
でも、速いがゆえに、犠牲にしていることがあります。

1. アセンブルに使われたリードの情報（どのリードがどのContig形成に用いられたか等の情報）は持たない。　他のアセンブラーではACEファイルというものを持つようです。
2. k-mer は最大でも31です。　これはCLCの設定です。　31塩基は、ショートリードには向いているかもしれませんが、400塩基を超すようなロングリードでは短いですし、結果、Contigが断片化してしまうことがあります。
3. Isoformの認識は失われる。　CLCのde Burijn graphの宿命です。　Contigが途中から枝分かれするようなとき、「枝」は別のContigとされます。　これを修正しているアルゴリズムは別途紹介します。

しかし、4分という速さは驚くべきことで、「とりあえずやってみよう」的な使い方には向いているでしょう。

現在、MiSeqやPGMに代表されるように、より多くの、ロングリードが出てくるようになってきています。　OLCアルゴリズムでは今後、計算速度が大きな壁になるでしょう。
ちなみにこの論文では、64bit Linuxの、普通のワークステーションレベルのPCを使っていました。


さて、同じ de Bruijn graph でも、Isoformを見つけよう！　パラログを認識しよう！というコンセプトのアセンブラーが最近出てきました。
454のようなロングリードではなく、ショートリード用ですが、Trinityというものです。

http://trinityrnaseq.sourceforge.net/

最近使い始めたばかりですので、このツールについてはまたの機会に。

MiSeqのデータセット

イルミナ社も、いよいよデスクトップ型シーケンサーをリリースするそうです。
それに先立ち、早くもデモデータセットが公開されています。　
http://www.illumina.com/systems/miseq/ecoli.ilmn

このページの、右上あたりにある、Download the data　から、FastqファイルとBAMファイルが落とせます。　E.coliのデータです。
サイズが大きく、ダウンロードに結構時間がかかるので、wget で取得すると良いと思います。

（wget http://.........ファイルのURL............）

ついでに、というかこっちの方が面白いかもしれませんが、BAMファイルの下の、プレゼンテーションファイルも落として見てみましょう。
びっくりしますが、スライドの後半部分に、MiSeqデータとIon Torrent PGMデータの比較、があるのです！

ここから先は実際にダウンロードしてからのお楽しみ？
当たり前ですが、MiSeqの方がいい！　と結論しています。

カバレージで言うと、1ランで1.7Gb読めるMiSeqは、393x
一方PGMは6~8ラン読んでも11~24Mbで、2~5x

クオリティはどうでしょう？
プレゼンにも書いてあるのですが、せっかくReadファイルを落とせるので、fastqファイルを取得し、FastQCにかけてみました。
ここから先は、イルミナさんのプレゼンではなく、私が行った内容です。
ちなみにPaired-Endです。以下、R1とR2の順に示しています。

Quality Score

GC Content per Sequence

その他、リード長は151、平均クオリティはQ35/Q34、GCはどちらも50%でした。

私が気になったのは、クオリティのばらつきです。
150塩基まで読めるとはいっても、リードによっては、後ろの方は結構Qが低い塩基も含まれるようですね。
とは言っても100塩基あたりまで、平均してQ30あったのは驚きました。


実際、このデータでアセンブルをすると、どれくらいつながるのでしょうか。
先のプレゼンの中に、じつは答えはあるのですが。

Ion Torrent データセット(2) クオリティチェック

先日、Ion Torrent PGM のデータセットが、Ion Torrent 社のホームページまたはBGIのFTPからダウンロードして誰でも使えるということを紹介しました。　Ion Torrent データセット
そのデータとは、E.coli O104-H4 の de novo sequence です。　
このことを書いた後、ある方と、クオリティの話をしました。


ということで、本日はクオリティチェックの結果を書きます。
使ったファイルは、Ion Torrent 社のサイトから取得した、LB226692株のfastq ファイルの、64.fastq ファイルです。　本当は64.fastq - 71.fastq まで、8ファイル全部やりましたが、ここに示すのは64.fastq ファイルのみということで。

使ったツールは、FastQC と、PRINSEQ　です。
両方とも同じようなツールです。
FastQCはPCの上で動くのに対し、PRINSEQはクラウド上で動くのが違いです。
ちなみにPRINSEQは、Internet Explore 9 ではまだうまく動かないみたいです。
FireFox、Google ChromeではOKでした。

FastQCの方が有名？

PRINSEQもいいですよ。
左下のUpload Fileをクリックして、ファイルを指定して、Upします。
でも難点は1つずつしかファイルをUPできないこと。　

ではここから、

1. リード長のばらつき
2. リードにおけるクオリティの変化
3. GCコンテンツ

 に絞って、結果を見ていきましょう。　2つのツールの結果画面を交互に示します。

リード長のばらつき　(Length Distribution)
FastQC

PRINSEQ

108塩基長のリードが最も多いようです。
ばらつき具合はグラフから一目了然ですね。

リードのクオリティ　(Read Quality Distribution)
FastQC

PRINSEQ

多くのシーケンサーデータ同様、Readの最初はクオリティが高いですね。　
後ろの方に行くに従ってだんだんと低くなり、90塩基付近では10を切ってしまいます。

GC含量　(GC Content Distribution)
FastQC

PRINSEQ

確かE.coliのゲノムはGCリッチで50%位だったと思います。まちがってたらすみません。

このデータもリードのCG含量が平均50%ですので、まあ、想定内でしょう。

FastQCとPRINSEQを例に、リードのクオリティチェックをしました。

本当はもっとメニューがあります。

FastQCはWindowsでも動きます。

PRINSEQはインターネットにデータを送り、結果はWebブラウザで閲覧します。

とっても簡単なので、是非一度、自身でお試しください。

Ion Torrent のデータセット

ヨーロッパで先月末に発生して1500人以上が感染、17人が死亡した、大腸菌O104-H4。
日本のニュースでも話題になりましたが、これは最初、新種の大腸菌ではないかと報道されました。
(読売新聞6月3日「欧州で広がるＯ１０４、新種の可能性…ＷＨＯ」)

やがて、ドイツと中国、その他のチームによってゲノムが読まれ、これが新種ではないものの、ハイブリッドな特徴を持つ厄介な菌であることがわかりました。　BGIによると、
このE.coliは、下痢原性大腸菌の一つのカテゴリーである腸管凝集性大腸菌 enteroaggregative E. coli (EAEC) の系統ですが、志賀毒素を産生するファージゲノムを自身のバクテリアゲノムに組み込み、さらに多薬剤耐性遺伝子をもゲノムに組み込んでいたことがわかりました。
環状ゲノムのサイズは　5，278 kbp　で、ほかに 88 kbp, 75 kbp, 1.5 kbp の3つのPlasmidから成るとのこと。
通常、志賀毒素を出す遺伝子、多剤耐性遺伝子は、ファージが関係するのですが、このE.coliは自身のゲノムにこれら遺伝子を組み込んでいることから、暫定的にShiga toxin-producing enteroaggregative Escherichia coli (STpEAEC)　と呼ばれているそうです。
腸管出血性大腸菌 enterohaemorrhagic E. coli (EHEC) に特徴は似ているけれど別の名前で呼ばれているんですね。知りませんでした。

さて、前置きが長くなりました。　
このドラフトシーケンシングに用いられた機器が、半導体シーケンサー　Ion Torrent PGM です。　
（実はBGIは、PGMのほかにIlluminaのマシンも使って読んでいるのですがそれはさておき）
Ion Torrent社のホームページへ行くと大きく宣伝しているのがわかるでしょう。

今はまだ数少ない、PGMのデータ（この大腸菌のリード）はここから入手できます。
Sffフォーマット：　http://lifetech-it.hosted.jivesoftware.com/docs/DOC-1621
Fastqフォーマット：　http://lifetech-it.hosted.jivesoftware.com/docs/DOC-1516

もしかするとユーザー登録が必要かもしれません。
また、アプリケーションノートやビデオなどは、Ion Torrentを良く知るために大変役に立ちますのでお勧めです。
念のため申しますと、私は Ion Torrent 社及びライフテクノロジーズ社と、特別な関係、があるわけではありませんのでご安心下さい。　

ひとつ断わりをいれておきますが、Ion Torrent社のサイトから落とせるデータは、BGIが読んだデータではありません。
ストレインが違います。
現在5つのストレインでゲノムが読まれていて

• TY2482 (BGI in collaboration with University Medical Centre Hamburg-Eppendorf)
• LB226692 (Life Tech in-house in collaboration with University of Muenster)
• H112180280 (Health Protection Agency, Colindale, UK)
• 2 isolates, unnamed (Gottingen Genomics Lab, Germany)

BGIはTY2482を読んでいます。

Ion Torrentのサイトから落とせるデータは、LB226692のデータです。

BGIのデータを使いたい方は、

NCBI SRA からは、SRX067313　で検索すると出てきます。

BGIのサイトからは、ftp://ftp.genomics.org.cn/pub/Ecoli_TY-2482
フォーマットはfastqです。

解析、de novo assembly については、彼ら（BGI）のワークフローが参考になるかと思います。
バクテリアゲノムアセンブリを初めてされる方は見ておいて損は無いでしょう。
http://climb.genomics.cn/Ecoli_TY-2482

見てわかる通り、BGIのチームは、IlluminaとIon Torrentの両方で読んで、アセンブルしました。
Ion Torrent PGM データのアセンブリに、Newbler を使ったのですね。
Roche 454 と同じ sff フォーマットなのでこれが使えるのでしょう。
PGMは、sff を出した後、fastqでも出力してくれます。　これならvelvetなどフリーツールでも使えますね。
現在、PGM に付属するソフトには、残念ながら独自のアセンブリツールが無いので、今のところPGMを使って de novo assembly をしようとすると、アセンブリソフトを別に求める必要があります。


さて、Ion Torrent のサイトから入手できる、LB226692のデータに話を戻しますね。
先のサイトから、sff のファイルをダウンロードしてきました。
8ファイル（8ラン分）あります。
私はRocheの機械を持っていませんし、アカデミアではなく企業の人間なので、タダでNewblerは使えません。
fastqを落としてアセンブってもいいがせっかく sff のフォーマットがあるのでできればsffのままアセンブりたい。　
さて、、、どのアセンブリツールを使おうか。

手元にちょうどCLC Genomics Workbench があるのでそれでやってみることにします。
最初に言い訳がましくなりますが、このPGMのデータだけを使ってアセンブリしても、1本にはなりません。
BGIは、7ランのPGMのデータ、200x以上のIlluminaのシングルリード、さらにIlluminaペアエンド、を使って読んでやっとドラフトゲノムを完成させています。
私がこれからお見せしようとしているのは、PGMのsffデータだけでアセンブルしたら、どれくらいのContigができるのか、というものです。
CLC のGenomics Workbenchには、sff フォーマットがインポートでき、アセンブルも簡単にできます。

取り込んだ後の様子

アセンブルした後のマッピングした結果のContigリスト

Contigをクリックしたときの様子

データ量：　8ラン分のデータ、平均98bpの117万リード、合計114Mbp

アセンブルに要した時間：　1分10秒
マシンの性能：　64Bit Linux, メモリ24Mb

Contigについて（200bp以上の長さのみ）

2,700本のContigができました。　最長16,000 bpでした。
ふーん、こんなもんかな。
こんどはちゃんと配列の前処理をやってからランしようかと、思った次第です。

と、ここで終わってしまっては結構さみしいので、次に、BGIのTY2482のリードデータを使って、
Ion Torrent 7ラン分　＋　Illumina 200x分　のハイブリッドアセンブリをやってみました。
フォーマットはfastq です。　結果をまとめると、

データ量：　合計　20.5 M reads, 1.8 G bp
ランタイム：　15分
Contigについて
512本のContig、最長224,799 bp、N50 = 60,618 bp！

実はBGIはペアエンドでも読んでいるのですが、データは公開されていません。　
ゆえにこれは、Ion Torrent PGM ＋Illumina GAII シングリリード　のハイブリッドアセンブリの結果、です。
ショートリード（GAII）とロングリード（PGM）がうまく互いを補いながら、長くつないでいるのでしょうね。
ちなみにCLC Genomics Workbench のパラメータはデフォルトです。
このソフトについてはこちらを参照
http://www.w-fusion.com/J/CLC_wb.html　

・・・・・・今回は少し長文になりました。



大腸菌についての参考資料
広島県保健環境センター研究報告，No.12，p1－12，2004
http://www.pref.hiroshima.lg.jp/hec/press/pdf/kenkyuhoukoku12/01.pdf