ex701 - nibb-gitc/gitc2021sep-rnaseq GitHub Wiki

ex701: de novo RNA-seq assembly using Trinity

Trinity でRNA-seq readsをde novo assembling. 実戦演習1 case1 でも使っているシーケンスデータを使って演習する。L1_R1.fastq　と　D1_R1.fastq　のRNA-seq reads をde novo アセンブルする。データの中身の説明については実戦演習2 case2 のwiki pageを参照の事。

Set up

本練習問題は、bias5のLinux環境で解析する。

Software

本解析に必要なソフトウェアは以下の通り。bias5にインストール済みである。

  • Trinity

Data

実戦演習1のイルミナシーケンスデータを使う。MiSeqによる 76-base long の　Single-end データ。

  • ~/gitc/data/EX/case1/IlluminaReads/L1_R1.fastq
  • ~/gitc/data/EX/case1/IlluminaReads/D1_R1.fastq

準備

bias5サーバへ ssh ログインする。

(example)

$ ssh [email protected]

(course01の部分は各自配布されたアカウントに置き換える）。

ex701ディレクトリを新しく作成しその下で解析を行おう。

$mkdir ex701
$cd ex701

Trinity が正常に動くか確認。

$ Trinity --version
Trinity version: Trinity-v2.13.1

Input readsの準備

扱いやすいよう、L1_R1.fastq　と　D1_R1.fastq　をワーキングディレクトリにコピーもしくはシンボリックリンクしておく。

$ ln -s ../data/EX/case1/IlluminaReads/L1_R1.fastq
$ ln -s ../data/EX/case1/IlluminaReads/D1_R1.fastq

Run Trinity

Trinity --seqType fq --single L1_R1.fastq,D1_R1.fastq --CPU 4 --max_memory 10G > run_Trinity.log 2>&1 &
  • --CPU, --max_memory は実行するコンピュータの環境によって変更する。--CPUの引数の値を増やせば、計算に使うCPUが増え計算時間は短縮できる。ただし、今回の演習では、多数の受講生が同時に解析するので、決して4以上には増やさない事。ちなみにこの条件で、計算に約1時間かかる。

（共有の大型計算機を使う際は、本来はジョブキューイングシステムを使うべきである。今回はトレーニングコース用に特別にログインノードで直接実行することを許可する。）

上のコマンドでは、標準出力とエラー出力をまとめて　run_Trinity.log　に保存するように指定している。(> run_Trinity.log 2>&1 の部分)。より詳しく知りたい人は、標準出力、標準エラー、リダイレクトのキーワードでUNIX/Linuxのコマンドラインの基礎の書籍やオンラインリソースを調べて欲しい。一番最後の"&"はバックグラウンドで実行するためのおまじない。

実行中にログファイルの記録を追いかけながら見る為には tail -f を使う。

$ tail -f run_Trinity.log

問題なく完了すると、ログファイルの末尾に、

#############################################################################
Finished.  Final Trinity assemblies are written to /mnt/gpfsA/home/courset1/prep_shige/ex701/trinity_out_dir.Trinity.fasta
#############################################################################

と書かれているはず。つまり、最終結果は、「trinity_out_dir.Trinity.fasta」

trinity_out_dir　ディレクトリ以下には様々な中間ファイルが記録されている。

その他の重要なファイル

  1. trinity_out_dir.Trinity.fasta.gene_trans_map

This is the gene-to-transcript mapping file.

(example)

TRINITY_DN84_c0_g1      TRINITY_DN84_c0_g1_i1
TRINITY_DN84_c0_g1      TRINITY_DN84_c0_g1_i2
TRINITY_DN70_c0_g1      TRINITY_DN70_c0_g1_i1
TRINITY_DN70_c0_g1      TRINITY_DN70_c0_g1_i2
TRINITY_DN70_c1_g1      TRINITY_DN70_c1_g1_i1
TRINITY_DN86_c0_g1      TRINITY_DN86_c0_g1_i1
TRINITY_DN97_c0_g1      TRINITY_DN97_c0_g1_i1
TRINITY_DN97_c0_g2      TRINITY_DN97_c0_g2_i1
TRINITY_DN46_c0_g1      TRINITY_DN46_c0_g1_i1
TRINITY_DN28_c0_g1      TRINITY_DN28_c0_g1_i1
TRINITY_DN28_c0_g1      TRINITY_DN28_c0_g1_i3
TRINITY_DN28_c0_g1      TRINITY_DN28_c0_g1_i4
TRINITY_DN28_c0_g1      TRINITY_DN28_c0_g1_i5
...

遺伝子単位でカウントをまとめる際に必要となる。

Trinity のnaming ruleについて：Trinity_DN[num]_c[num]_g[num]_i[num]. Trinity_DN[num]_c[num]_g[num]までが共通な名前は同じ遺伝子の異なるトランスクリプト（splicing variants）と解釈。

Inspect results

アセンブル結果、"trinity_out_dir.Trinity.fasta"　を検証しよう。

アセンブリの善し悪しを評価する為にはいくつかの指標がある。基本的なものでは、コンティグの数、総塩基数、平均長、N50など。論文では、遺伝子カタログの完全性の検証としてBUSCOなどが使われることが多い。

  1. コンティグの数がいくつあるか、UNIXのコマンドラインで調べよう。(hint: grep, wc を使う)。seqkit statsも使える。
$ grep "^>" trinity_out_dir.Trinity.fasta |wc
   7277   22963  360102

$ seqkit stats -a trinity_out_dir.Trinity.fasta
file                           format  type  num_seqs    sum_len  min_len  avg_len  max_len   Q1   Q2   Q3  sum_gap  N50  Q20(%)  Q30(%)
trinity_out_dir.Trinity.fasta  FASTA   DNA      7,277  3,013,917      201    414.2    6,372  238  302  470        0  453       0       0

contig数7413, N50=453bp であることがわかった。

  1. Trinityソフトウェアに含まれる　TrinityStats.pl でassembly 諸statisticsを計算。
$ /bio/package/Trinityrnaseq/2.13.1/util/TrinityStats.pl trinity_out_dir.Trinity.fasta

結果例


################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':	7095
Total trinity transcripts:	7277
Percent GC: 44.85

########################################
Stats based on ALL transcript contigs:
########################################

	Contig N10: 1246
	Contig N20: 911
	Contig N30: 715
	Contig N40: 565
	Contig N50: 453

	Median contig length: 302
	Average contig: 414.17
	Total assembled bases: 3013917


#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

	Contig N10: 1234
	Contig N20: 891
	Contig N30: 690
	Contig N40: 550
	Contig N50: 438

	Median contig length: 299
	Average contig: 407.58
	Total assembled bases: 2891788

Trinityは複数のisoformが存在する場合別々のエントリーとして出力する。TrinityStats.pl　はgene levelと、(isoformw を区別した) transcript levelの集計結果を解析して表示してくれる。なお、gene <=> isoform の関係は、trinity_out_dir.Trinity.fasta.gene_trans_map　に保存されている。

Revision History

2021-9-16

  • .gene_trans_map　についての情報を追加。

2021-9-11

  • v2.13.1 に対応。

2021-3-X