最近更新を続けているEnTAPでde novo assemble するRNA-seqデータにアノテーションをつけるプロジェクトの下準備編。
今回はTrinityでde novo アセンブルを行う。
2020/07/09 追記: 出力結果の見方を追記
目的と実行環境
複数のRNA-seqデータをマッピングするためのContigを作成したい。さらに今回はEnTAPを用いてそのContigsに対してGOなどのアノテーションを行いたい。
今回は非モデル生物のため、Trinity を用いて新たにアセンブルを行うこととした。
実行環境
- AMD Ryzen 9 3900X (12コア24スレッド)
- 64GB RAM
- 1TB SSD
- Ubuntu 20.04
手法
Trinity のインストール
https://bioconda.github.io/recipes/trinity/README.html
$ conda install trinity
これでいけるらしい。らしい、というのも随分前にインストールしたのでどうやって入れたか覚えていない……
とりあえず現段階ではv2.6.6が入っていた。最新版はv2.11.0なのでそちらをインストールしておくとよいだろう。
実行
実行記録を残すためにも簡単なシェルスクリプトを残した。
Trinity --seqType fq --left ./c04/paired_output_1.fq.gz,.c05/paired_output_1.fq.gz,./c06/paired_output_1.fq.gz,\ --right ./c04/paired_output_2.fq.gz,./c05/paired_output_2.fq.gz,./c06/paired_output_2.fq.gz,\ --CPU 24 \ --max_memory 60G \ --output trinity_out
これをTrinity.shとして保存。だいぶ省いているが、全部で24個のペアエンドのFASTQがあった。
あとは
$ sh Trinity.sh
で実行される。実行時間はおよそ10時間程度。
つまづいた点
躓いた点としては、カンマの後に半角スペースを入れる癖があるのだが、それを入れてしまうとプログラムが上手く動作しない。別のものとして認識してしまうようだ。
あとは今回インプットにペアエンドの合計24ファイルのFASTQを入れた。gzで圧縮してあるが、1つあたり1.2GB程度ある。そのせいでTrinityの一時ファイルのファイルサイズが300GBを超え、ディスク容量が足りなくなり結果途中で止まってしまった。
そのためある程度ディスク容量には余裕をもたせておきたい。
ちなみに一時ファイルのため、一定以上段階が進むと削除されるのでご安心を。
最終的には1GBも行かないFASTAファイルが出力される。
次はこれを用いてアセンブルされた配列にアノテーションを行う。
追記: Trinity 出力データの解釈
新しい記事にするのも微妙な長さなので、ここにTrinityの出力データ、Trinity.fastaについての解釈を置いておく。
参考: 公式ドキュメント
以下のようなマルチFASTAフォーマットで出力される。
>TRINITY_DN100000_c0_g1_i1 len=344 path=[1:0-123 2:124-279 4:280-335] ATATATATAGAGAGAGATGTCCCCCCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA…… >……
TRINITY_DN100000_c0_g1_i1 len=344 path=[1:0-123 2:124-279 4:280-335]
のうち、DN100000_c0はTrinityリードクラスタと呼び、g1は遺伝子、i1はアイソフォームの1ということを意味する。
遺伝子g1はそのリードクラスタ内で一意に割り当てられる番号である。つまり、ある「遺伝子」に着目したとき、取得すべき配列はTRINITY_DN100000_c0_g1であるということだ。
TRINITY_DN100000_c0_g1_i1はアイソフォームi1をコードする遺伝子TRINITY_DN100000_c0_g1という意味になるようだ。
path=[1:0-123 2:124-279 4:280-335]についてだが、1，2，4などの:の前についている数字はノードの番号を指す。これは異なるアイソフォーム間でも共有されている番号であり、異なる配列間での共通配列を見つけ出すのに有効である。:のあとの数字、0-123のようなものは配列のどこからどこまでを指すかの範囲である。ここでは0番目の塩基(人間からすると1番目)から123番目の塩基までを指す。
例えば以下のような遺伝子があったとしよう。
>TRINITY_DN1000_c7_g3_i3 len=336 path=[1:0-123 2:124-279 4:280-335] ATTGATATATATATTTTTTTTTTTTT…… >TRINITY_DN1000_c7_g3_i2 len=302 path=[0:0-57 2:58-213 3:214-301] CACACACACAACACACA…… >TRINITY_DN1000_c7_g3_i1 len=270 path=[0:0-57 2:58-213 4:214-269] CACACACACAACACACA……
こういう配列があったとき、この3つの配列はTRINITY_DN1000_c7_g3遺伝子のアイソフォーム1から3が存在するということになる。
これらのアイソフォームの配列間を比較するときに、共通している配列があり、それがpath=[]の中に書かれているノード番号である。i1, i2ではノード0，2を共有しており、実際その配列を見ると同じである。また、1～3においてもノード2を共有している、といったようなことがわかる。
以上がTrinityでの出力結果のヘッダ行の見方である。
