ex701 - nibb-gitc/gitc2025mar-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 アセンブルする。

Set up

Software

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

  • 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

準備

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

$mkdir ex701
$cd ex701

Trinity が正常に動くか確認。

$ Trinity --version
Trinity version: Trinity-v2.15.1

もしTrinityが起動しない場合は、RCCS環境では、modulefilesというシステムでソフトウェアを管理しているので、 module load Trinityrnaseq/2.15.1 と入力して、Trinityのver. 2.15.1 をロードする。

Input readsの準備

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

$ ln -s ~/gitc/data/EX/case1/IlluminaReads/L1_R1.fastq
$ ln -s ~/gitc/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
   7269   22895  359356

$ 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(%)  GC(%)
trinity_out_dir.Trinity.fasta  FASTA   DNA      7,269  3,010,808      201    414.2    6,372  238  301  469        0  453       0       0  44.84

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

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

結果例

################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':	7092
Total trinity transcripts:	7269
Percent GC: 44.84

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

	Contig N10: 1255
	Contig N20: 914
	Contig N30: 714
	Contig N40: 564
	Contig N50: 453

	Median contig length: 301
	Average contig: 414.20
	Total assembled bases: 3010808


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

	Contig N10: 1236
	Contig N20: 891
	Contig N30: 689
	Contig N40: 549
	Contig N50: 438

	Median contig length: 299
	Average contig: 407.41
	Total assembled bases: 2889367

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

Revision History

2023-2-25

  • minor update (Trinity v2.15.1 on bias5)

2022-8-28

  • minor update (Trinity v2.14.0)

2022-2-26

  • minor update (Trinity v2.13.2)

2021-9-16

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

2021-9-11

  • v2.13.1 に対応。

2021-3-X