top of page
このホームページは
.com
を使って作成されました。あなたも無料で作ってみませんか？今すぐはじめる

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%程度あれば良いと見なしていいのかな。

​

​

bottom of page