公開日: 最終更新日:

Trinity で複数のRNA-seqのデータをまとめてアセンブル + 出力データの見方

最近更新を続けている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での出力結果のヘッダ行の見方である。