Trinity

2018 年 11 月 2日 改訂

このページはまだ未完成です.

Trinity はトランスクリプトーム解析専用の k-mer アセンブラです.ここでは, Trinity とその他のプログラムを用いて NGS で得られた トランスクリプトームデータを de novo assembly (参照配列を用いずに配列張り合わせ) する手順を紹介します.得られたリードが 150bp より長い場合は,Newbler の方が良いみたいです.

SGE (Sun Grid Engine) スケジューラーで解析を行う場合は,こちら のジョブスクリプトを参照してください.


1. Trinity

ここでは Trinity を用いて,参照配列を用いない配列張り合わせ (de novo assembly) を行います.Trinity はfastq ファイルを読み込み,一つの fasta file (すべてのトランスクリプトーム配列を含む) を算出します.

解析を行う前に,Solexa QA 2 や condetri などによって (上記の手順),quality が低いリードを除外してからアセンブルを実行した方が良いらしいです.

こちらは英語の資料ですが,Trinity の全体像,つまり,Inchworm[シャクトリ虫], Chrysalis[サナギ], Butterfly[チョウ] の役割,をつかむには良いです (とくに HOW IT WORKS というビデオ).日本語の資料としてはこちら (1, 2) を参考にしました.


インストール
こちらをご覧ください.


テストデータの解析

ダウンロード & コンパイルして得られたファイル,

trinityrnaseq_r2013-02-25/sample_data/test_Trinity_Assembly

に入ってください.そこで,

sh runMe.sh

と入力すれば,解析が始まるはずです.trinity_out_dir にアウトファイルが保存されます.Trinity.fasta が最終的に必要なファスタファイルです.
実際にこのスクリプト内部で解析を行っている行は,

../../Trinity.pl --seqType fq --JM 2G --left reads.left.fq --right reads.right.fq --SS_lib_type RF --CPU 4 --no_cleanup --monitoring

です.もう一行,

../../util/RSEM_util/run_RSEM_align_n_estimate.pl --transcripts trinity_out_dir/Trinity.fasta --seqType fq --left reads.left.fq --right reads.right.fq --SS_lib_type RF

がありますが,これは発現量解析などで使うコマンドのようです.こちらをご覧下さい.私は使わないです.

Slurm の job file は以下です.

#!/bin/bash
#SBATCH --job-name=TriTrimed
#SBATCH --mail-user="jun.inoue@oist.jp"
#SBATCH --partition=compute
#SBATCH --mem=100G
#SBATCH --cpus-per-task=20
#SBATCH --ntasks=1 # 1 task

software_dir='/work/SatohU/0inoue/software'
module load cmake/3.11.4
export PATH=$software_dir/salmon-0.11.3-linux_x86_64/bin:$PATH
export PATH=$software_dir/bowtie2-2.3.4.3-linux-x86_64:$PATH
export PATH=$software_dir/jellyfish-2.2.10/bin:$PATH
export PATH=$software_dir/trinityrnaseq-Trinity-v2.8.4:$PATH
export PATH=$software_dir/samtools-1.9/bin:$PATH

TRINITY_HOME='/work/SatohU/0inoue/software/trinityrnaseq-Trinity-v2.8.4'

${TRINITY_HOME}/Trinity --seqType fq --max_memory 100G \
--left paired_output_sub1.fq \
--right paired_output_sub2.fq\
--CPU 20 \
--output trinity_out1


--max_memory
スケジューラーで割り当てたメモリー数を書きます.

Trinity.pl がうまく走ると,以下のようなスクリーンアウトが得られるはずです.

[jun-inoue:test_Trinity_Assembly]$ perl /home/j/jun-inoue/bin/trinityrnaseq_r2013-02-25/Trinity.pl --seqType fq --JM 2G --left reads.left.fq --right reads.right.fq --SS_lib_type RF --CPU 4 --no_cleanup --monitoring
Use of uninitialized value in pattern match (m//) at /home/j/jun-inoue/bin/trinityrnaseq_r2013-02-25/Trinity.pl line 474.
WARNING, --monitoring can only be used on linux. Turning it off.

Current settings:
core file size (blocks, -c) 0
data seg size (kbytes, -d) unlimited
scheduling priority (-e) 0
file size (blocks, -f) unlimited
pending signals (-i) 399360
max locked memory (kbytes, -l) unlimited
max memory size (kbytes, -m) unlimited
open files (-n) 1024
pipe size (512 bytes, -p) 8
POSIX message queues (bytes, -q) 819200
real-time priority (-r) 0
stack size (kbytes, -s) unlimited
cpu time (seconds, -t) unlimited
max user processes (-u) 399360
virtual memory (kbytes, -v) unlimited
file locks (-x) unlimited

Paired mode requires bowtie. Found bowtie at: /usr/bin/bowtie

-since butterfly will eventually be run, lets test for proper execution of java
#######################################
Running Java Tests
CMD: java -Xmx64m -jar /imports/home/j/jun-inoue/bin/trinityrnaseq_r2013-02-25/util/ExitTester.jar 0
CMD finished (0 seconds)
CMD: java -Xmx64m -jar /imports/home/j/jun-inoue/bin/trinityrnaseq_r2013-02-25/util/ExitTester.jar 1

-we properly captured the java failure status, as needed. Looking good.
Java tests succeeded.
###################################

..........

---------------------------------------------------------------
-------------------- Butterfly --------------------------------
-- (Reconstruct transcripts from reads and de Bruijn graphs) --
---------------------------------------------------------------

CMD: /imports/home/j/jun-inoue/bin/trinity........
Number of Commands: 52
succeeded(52) 100% completed.

All commands completed successfully. :-)

CMD finished (16 seconds)
CMD: /imports/home/j/jun-inoue/bin/trinityrnaseq_r2013-02-25/util/print_butterfly_assemblies.pl /work/........
CMD finished (0 seconds)

###################################################################
Butterfly assemblies are written to /work/SinclairU/inoue/tr
###################################################################

[jun-inoue:test_Trinity_Assembly]$



リードの方向

こちらのページに,まとめられています.
(以下,作成中です)
--SS_lib_type オプションで設定します.わからない時はこのオプションは使わない方が良いです.ライブラリの作り方がわかっていない場合におかしな設定を行うと,結果が無茶苦茶になるそうです.以下は mouse SRA のページを例としています.

青のラインがある部分で,NEBnext Ultra RNA Library Prep Kit を使った,と書いてあります.この製品試薬のページから判定します.


アウトファイル
trinity_out_dir というディレクトリが自動的に作成され,ここに結果が保存されます.このなかにある,

Trinity.fasta

というファイルに,アッセンブルの結果得られた contig が fasta 形式で保存されています.Trinity.fasta が出力されない場合は,何か解析に問題があったことになります.Butterfly はメモリをたくさん使うようで,私の場合は Butterfly が終了せず Trinity.fasta が得られないことがよくあります.
 他にたくさんのアウトファイルが得られますが,普通は使いません.

#####

TrinityStats.pl を使うと,得られた結果の内容を知ることができます.

[jun-inoue:TEST]$ perl $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta
Total trinity transcripts: 90
Total trinity components: 47
Contig N50: 4954



データ量

リード数が 30M〜60M が常識的なデータ量だそうです.これぐらいなら,Trinity でそのまま解析して良いです.これより多い場合は,リード数を間引くなどの処理をしないと,結果の質が悪くなることがあるそうです.良いアッセンブルを行うのに 100M 必要,という人もいますが,30M で十分だそうです.


エラー
ヘッダー:ヘッダーを変更するように,エラーメッセージが出る.

[cluster:Symsagittifera-roscoffensis]$ cat slurm-4622056.out
Trinity-v2.8.4
....
CMD: seqtk-trinity seq -A /work/SatohU/0inoue/Symsagittifera-roscoffensis/paired_output_sub2.fq >> right.fa
Error, not recognizing read name formatting: [SRR5760179.40838091]

If your data come from SRA, be sure to dump the fastq file like so:

SRA_TOOLKIT/fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files file.sra

Thread 1 terminated abnormally: Error, cmd: seqtk-trinity seq -A /work/SatohU/0inoue/Symsagittifera-roscoffensis/paired_output_sub1.fq >> left.fa died with ret 512 at /work/SatohU/0inoue/software/trinityrnaseq-Trinity-v2.8.4/util/insilico_read_normalization.pl line 762.

この例題は SRA データを解析中に出たものです.bold の行は,ヘッダーに問題があるから,変更する必要がある,というメッセージが出ています. NCBI が配っている SRA tool kit に入っている fastq-dump を用いて SRA ファイルを二つの fastq ファイルに分割する際に,上記 bold のオプションをつけるように支持されています.fastq-dump の使い方はこちらをご覧ください.

bowtie:bowtie が見つからないというメッセージが出ました.

[i2:test_Trinity_Assembly]$ sh runMe.sh
....
Error, cannot find path to bowtie, which is now needed as part of Chrysalis' read scaffolding step at ../../Trinity.pl line 622.

こちらのサイトにあるダウンロードサイトから bowtie-1.0.0-macos-i386.zip をダウンロードしました (uname -p で CPU の種類を確認).解凍して得られた bowtie を ~/bin にコピーしたら問題なく動きました.bowtie-build もない場合は,~/bin にコピーしましょう.



データの前処理は必要?
Solexa QA 2 と filterPCRdupl.pl によって前処理を行った場合,Trinity の結果がどの程度少なくなるか Trinity のテストデータ (reads.left.fqreads.right.fq) を用いて比較してみました.解析には TrinityStats.pl を使いました.

前処理をやらない場合

Total trinity transcripts: 90
Total trinity components: 47
Contig N50: 4954

前処理をやった場合

Total trinity transcripts: 76
Total trinity components: 41
Contig N50: 3739

クオリティについては,そのうち検証してみます.


2. CD-HIT

CD-HIT は,Trinity などで行ったアッセンブル後に得られた配列群から,似たような配列を取り除くソフトウェアです.こちらの左にあるコラムの下の方からダウンロードができます.マニュアルはこちらです.

類似の配列をたくさん集めたいなど,場合によっては CD-HIT による解析は行わない方が良いです.よく検討してください.

CD-HIT の解析は割と時間がかかるので,クラスタで行った方が良いでしょう.

Trinity の例題 (reads.left.fqreads.right.fq) で得られたファイル (Trinity.fasta) を解析したのですが,ファイルが小さすぎたようで (Total trinity transcripts: 90) うまく走りませんでした.


コンパイル
ダウンロード & 解凍して得られたファイルに入り,

make

と入力してください.
multi-threading に対応させるには,

make openmp=yes

と入力してください.


トランスクリプトームデータの解析
cd-hit-est を使います.

cd-hit-est -r 1 -T 1 -M 1000 -n 5 -c 0.8 -i INFILE.fa -o OUTFILE.fs >& log.txt &

-r 1
両方のストランドを比較します.

-T 1

0 だとすべてのスレッドを使います.TOMBO で解析を行うときは 1 や 4 にしましょう.

-M 1000

メモリ 100Mb を意味します.1000Mb で十分だそうです.

-n 5
ワードサイズ.なくても (デフォルトでも) 大丈夫だそうです.

-c 0.8

アライメントの類似性

>&

エラーと出力ファイルを同じファイルに書き込む.

&

qsub する場合はいりません.



ログファイル

解析がうまく走り出したら,以下のようなログファイルが作成されます.

[jun-inoue:clearnDir]$ cat log.txt
================================================================
Program: CD-HIT, V4.6, Jul 16 2013, 13:59:58
Command:
/apps/SinclairU/cd-hit-v4.6.1-2012-08-27_withoutMULTI/cd-hit-est
-T 4 -M 1000 -n 5 -c 0.8 -i
trinity_out_dir/Trinity.fasta -o eyeClearn.fas

Started: Tue Jul 16 14:59:53 2013
================================================================
Output
----------------------------------------------------------------
Option -T is ignored: multi-threading with OpenMP is NOT enabled!
total seq: 136837
longest and shortest : 12149 and 201
Total letters: 75184757
Sequences have been sorted

Approximated minimal memory consumption:
Sequence : 92M
Buffer : 1 X 16M = 16M
Table : 1 X 2M = 2M
Miscellaneous : 1M
Total : 112M

Table limit with the given memory limit:
Max number of representatives: 1080875
Max number of word counting entries: 110945009

comparing sequences from 0 to 136837


アウトファイル
以下 2 種類のアウトファイルが作成されます.

*.fas
得られた結果です.

*.fas.clstr
おそらく得られた結果の長さなどが書いてあります.詳しくはわかりません.


コマンド集

grep ^@ <fastqファイル名> | wc -l

fastq ファイル内部のリード数をカウントします.@ で始まる行数を調べています.

grep ^\> <fastaファイル名> | wc -l

fasta ファイル内部のレコード数をカウントします.> で始まる行数を調べています.
注意: > の前にバックスラッシュを入れることでエスケープします.これをしないと,> がファイルを作成するなどリダイレクションと判断され,インファイル内部の情報が消えてしまいます.

リンク

FASTAX-Toolkit

PCR duplicate を取り除いたり,read の qulity trim, fastaq から fasta の変換などを行うソフトとして有名です.

Newbler

Roche 社が出しているアセンブラーです.

velvet-oases

トランスクリプトーム解析用のソフト.有名なアセンブラー.解析を行う前に,quality が低い塩基は除外してからアセンブルした方が良いらしいです.

非モデル生物の RNA-seq 解析

次世代シーケンサーを用いた de novo トランスクリプトーム解析

 

トランスクリプトームデータ解析シリーズ

次回は「4. TransDecoder による転写配列の推定」 のページです.アセンブルして得られた配列の翻訳される方向や 1st codon poistion を推定します.
1. SRA データのダウンロードと変換
2. FastQC による fastq データの検証
3. fastq データの精製
4. Trinity によるアッセンブル.
5. TransDecoder による転写配列の推定
6. ORTHOSCOPE によるオーソログ推定.

このページは主に OIST に所属する研究者の方から教えていただいた情報をもとに作成しています.皆さんのご協力に感謝します.