3-1. TRINITYを用いたde novoアセンブル
2021.06.10 初稿
Trinityを用いたde novo assemble
以下はスパコンの中のTrinityでアセンブルする際のシェルスクリプト。
-------------------
#!/bin/bash
#$ -S /bin/sh
#$ -l s_vmem=164G,mem_req=164G
#$ -cwd
#$ -o ~/results_sh_eando
#$ -e ~/results_sh_eando
ulimit -s unlimited
echo “pwd: $(pwd)”
echo HOME: $HOME
echo USER: $USER
echo HOSTNAME: $HOSTNAME
echo starting at
date
singularity exec -e /usr/local/biotools/t/trinity\:2.8.5--hdbcaa40_0 Trinity --full_cleanup --seqType fq --max_memory 200G --CPU 12 --left ~/RNAseq/sample_1.fq.gz --right ~/RNAseq/sample_2.fq.gz --output trinity_out --no_version_check
echo ending at
date
-------------------
このTrinityのバージョンでもアセンブルは可能。
しかし、のちのステップでリードカウントをしたりする際にRSEMとの連携がなく、最新版では使えるTrinity内のスクリプトが使えなくなっている。
そこで最新版のTrinityをダウンロードする。
ーーー（homeにダウンロード後）
$ wget https://data.broadinstitute.org/Trinity/TRINITY_SINGULARITY/trinityrnaseq.v2.11.0.simg
$ export LANG=C
$ singularity exec trinityrnaseq.v2.11.0.simg Trinity
ーーー
上のTrinityスクリプトのうち、実行コマンドを以下のように変える。
$ TRINITY_HOME=/usr/local/bin/trinityrnaseq
$ singularity exec -e ${HOME}/trinityrnaseq.v2.11.0.simg Trinity \
--seqType fq \
--CPU 12 \
--max_memory 200G \
--left ${HOME}/QC_out/qc_left.fastq.gz \
--right ${HOME}/QC_out/qc_right.fastq.gz \
--output trinity_out \
--full_cleanup \
--no_version_check
2つのリードファイルだけをアセンブルするときには２ファイルだけを指定したらよいが、例えばサンプル1~5まで全部使ってアセンブルする場合にはもう一捻り必要。
ペアエンドのleftとrightをそれぞれ一つのファイルにまとめる。
$ ls ./QC_out/qc_*_1* | sed -e “s/\n/,/g”> left.txt
$ cat left.txt > sample_all_left.fq.gz #これだとテキストファイルが映されるだけで、fq.gzを一つにまとめてくれない。
$ cat qc_sampleA*_1* qc_sampleB*_1* qc_sampleC*_1* qc_sampleD*_1* qc_sampleE*_1* > sample_all_right.fq.gz
上記のようにconcatenateしなくても、―left, ―rightの後ろにコンマで繋いでもOK。
Trinityを走らせる時はなんどもメモリが足りないと怒られて途中で止まっていることがしばしばありました。
Memoryを最大の200Gまで指定してあげてCPUを12で指定するとうまく走りました。
スパコン上でのCPU指定がどれだけ意味があるのか私はよくわかっていないのですが、メモリ不足のためにTrinity解析が停止することをどうにか回避するためのおまじないとして大事だと思ってます。
TrinityStatsを用いてde novo アセンブルの中身をチェック
まずはTrinityStatsがスパコン 上のどこにあるのかSingularity shellを用いてチェックする。
[~]$ singularity shell ./trinityrnaseq.v2.11.0.simg
Singularity> which Trinity
/usr/local/bin/trinityrnaseq/Trinity
Singularity> which TrinityStats.pl
Singularity> which align_and_estimate_abundance.pl
whichコマンドでは.plファイルを探すことができなかった。lsコマンドでutilディレクトリを除いてみる。
Singularity> ls -l /usr/local/bin/trinityrnaseq/util/
total 96
drwxr-xr-x 2 root root 352 Jul 1 2020 PBS
drwxr-xr-x 2 root root 84 Jul 1 2020 R
-rwxr-xr-x 1 root root 3866 Jul 1 2020 TrinityStats.pl
-rwxr-xr-x 1 root root 11279 Jul 1 2020 abundance_estimates_to_matrix.pl
-rwxr-xr-x 1 root root 29559 Jul 1 2020 align_and_estimate_abundance.pl
-rwxr-xr-x 1 root root 8024 Jul 1 2020 analyze_blastPlus_topHit_coverage.pl
-rwxr-xr-x 1 root root 9540 Jul 1 2020 filter_low_expr_transcripts.pl
-rwxr-xr-x 1 root root 28966 Jul 1 2020 insilico_read_normalization.pl
drwxr-xr-x 14 root root 6338 Jul 1 2020 misc
-rwxr-xr-x 1 root root 1058 Jul 1 2020 retrieve_sequences_from_fasta.pl
-rwxr-xr-x 1 root root 3823 Jul 1 2020 sift_bam_max_cov.pl
drwxr-xr-x 3 root root 2095 Jul 1 2020 support_scripts
TrinityStatsは/usr/local/bin/trinityrnaseq/util/にあることがわかった。
-------------------
TRINITY_HOME=/usr/local/bin/trinityrnaseq
singularity exec -e ${HOME}/trinityrnaseq.v2.11.0.simg \
${TRINITY_HOME}/util/TrinityStats.pl \
./trinity_out/trinity_out.Trinity.fasta
-------------------
結果は以下のように出てくる。
################################
## Counts of transcripts, etc.
################################
Total trinity 'genes': 122650
Total trinity transcripts: 239267
Percent GC: 40.68
########################################
Stats based on ALL transcript contigs:
########################################
Contig N10: 4271
Contig N20: 3328
Contig N30: 2764
Contig N40: 2326
Contig N50: 1949
Median contig length: 628
Average contig: 1100.31
Total assembled bases: 263268226
#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################
Contig N10: 3885
Contig N20: 2902
Contig N30: 2282
Contig N40: 1758
Contig N50: 1250
Median contig length: 407
Average contig: 749.74
Total assembled bases: 91956102
参考にしたウェブページ
・実際のTrinityの動かし方：https://github.com/trinityrnaseq/trinityrnaseq/blob/master/Trinity
de novo assembleした結果の評価
いくつか評価の仕方がある。
●BUSCOによる評価●
一つ目はBUSCO (Benchmarking Universal Single-Copy Orthologs) というどの生物にも存在している遺伝子セットが、de novoでアセンブルしたtranscriptsの中にどれだけ含まれているかを評価する方法。ゲノムであれば95%以上、トランスクリプトームであれば80%以上あればうまくつながっていると見なせる。
遺伝研スパコンでBUSCOを使うには>32GBメモリが必要なので、いつも通りqloginするとメモリ不足で動かない。以下のようにメモリ指定してqloginする(いつも通りqloginし、メモリ指定したjobをqsubで投げるのもあり)。
```
$ qlogin -l s_vmem=64G -l mem_req=64G
```
以下のコマンドでどのようなデータベースが使えるか調べられる。
$ singularity exec --no-mount tmp /usr/local/biotools/b/busco:5.2.2--pyhdfd78af_0 busco --list-datasets
QueryにはTrinity.fastaもしくはそれをTransDecoderで翻訳したfastaファイルを用いる。（queryにはゲノム、トランスクリプトーム、タンパク質が選べる）
$ singularity exec --no-mount tmp /usr/local/biotools/b/busco:5.2.2--pyhdfd78af_0 busco -m protein -i input.pep -o busco_pep -l arthropoda_odb10 -c 20
すると以下のような結果が得られる。
2021-06-10 14:06:12 INFO:busco.BuscoRunner
-----------------------------------------------
|Results from dataset arthropoda_odb10 |
--------------------------------------------------
|C:70.9%[S:69.6%,D:1.3%],F:13.8%,M:15.3%,n:1013 |
|718 Complete BUSCOs (C) |
|705 Complete and single-copy BUSCOs (S) |
|13 Complete and duplicated BUSCOs (D) |
|140 Fragmented BUSCOs (F) |
|155 Missing BUSCOs (M) |
|1013 Total BUSCO groups searched |
--------------------------------------------------
2021-06-10 14:06:12 INFO:busco.BuscoRunner BUSCO analysis done. Total running time: 78 seconds
ちなみに私が遺伝研スパコンでBUSCOをかけようとした日はbusco v5.2.2が使えず、v4.1.4など前のバージョンを使っても動かず、次の日になったらなぜかv5.2.2が使えるようになっていました。不思議…
また、Trinity.fastaをまるまるqueryにするとcoredumpして止まってしまったので、seqkitでqueryを10こに分割してBUSCOをかけました。
---
singularity exec -e /usr/local/biotools/s/seqkit:2.0.0--h9ee0642_0 seqkit split -p 10 longest_orfs.pep
---
以下が、10このqueryを一気に走らせるスクリプト。
---
"busco.sh"
#!/bin/bash
#$ -S /bin/sh
#$ -t 01-10
#$ -l medium
#$ -l s_vmem=64G,mem_req=64G
#$ -cwd
#$ -o ~/results_sh_eando
#$ -e ~/results_sh_eando
echo "pwd: $(pwd)"
echo HOME: $HOME
echo USER: $USER
echo JOB_ID: $JOB_ID
echo JOB_NAME: $JOB_NAME
echo HOSTNAME: $HOSTNAME
echo SGE_TASK_ID: $SGE_TASK_ID
echo SGE_TASK_FIRST: $SGE_TASK_FIRST
echo SGE_TASK_LAST: $SGE_TASK_LAST
echo SGE_TASK_STEPSIZE:
echo starting at
date
singularity exec --no-mount tmp /usr/local/biotools/b/busco:5.2.2--pyhdfd78af_0 \
busco -m prot \
-i longest_orfs.pep.part_00${SGE_TASK_ID}.fa \
-o busco_pep_0${SGE_TASK_ID} \
-l insecta_odb10
echo ending at
date
---
Busco_pep_0*のディレクトリ以下にrun_insecta_odb10/full_table.tsvというテーブルができている。
このテーブル10こをcatを用いて1つに結合して、ダウンロード。
2列目のstatusがDuplicatedかCompleteだったら、その遺伝子が存在していたということになるので、結合したfull_table.tsvを用いてDuplicatedかCompleteのstatusを持つ遺伝子数（=Busco id）を数える。冗長な遺伝子（=Busco id）もあるし、Trinity_IDも重複を除く。insecta_odb10は全体で1367個の遺伝子が準備されているので、得られたcompleteとduplicatedの遺伝子数を1367で割れば、BUSCO scoreを算出できる。
●リードをアセンブル結果にマッピングする●
アセンブルした配列に対してリードをマッピングし、どの程度のマップ率になるかで、アセンブル結果を評価することができる。
80%のマッピング率が最低ライン、モデル生物のillumina ハイクオリティリードの場合98%くらいになることもあるらしい。90%程度あれば良いと見なしていいのかな。