ex701 - nibb-gitc/gitc2021mar-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.11.0

（注）最新版はv2.12.0ですが、今回はv2.11.0 を使います。

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

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

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

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

$ tail -f run_trinitity.log

Result => trinity_out_dir ディレクトリ以下。"Trinity.fasta" が最終的なアセンブル結果。

Inspect results

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

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

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

$ 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,413  3,088,520      201    416.6    6,372  238  302  473        0  460       0       0

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

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

結果例

################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':	7180
Total trinity transcripts:	7413
Percent GC: 44.84

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

	Contig N10: 1271
	Contig N20: 923
	Contig N30: 720
	Contig N40: 573
	Contig N50: 460

	Median contig length: 302
	Average contig: 416.64
	Total assembled bases: 3088520


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

	Contig N10: 1238
	Contig N20: 896
	Contig N30: 692
	Contig N40: 553
	Contig N50: 439

	Median contig length: 298
	Average contig: 408.05
	Total assembled bases: 2929774

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