macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

Trinotateを使ってde novo transcriptome のアセンブリ配列にアノテーションをつける

2020 9/16 pythonのバージョンを指定して導入, configファイルについて追記

2020 9/16 コメント追加, signalPとTMHMM、rnammerの初期設定追加

2020 9/27 わかりにくい表現を修正

 

Trinotarteは、Trinityののde novoトランスクリプトームアセンブリを機能的にアノテーションするためのプロトコルおよび支援ソフトウェア。

 

以前紹介したTrinotarteの説明が分かりづらかったので、簡潔にまとめ直します。


TrinotateはデフォルトではTrinityのアセンブリを使ってアノテーション付けを行います。そのため、別のde novoアセンブラを使った場合は少しだけ工夫が必要になります。下に書いた手順を確認して下さい、

 

wiki

TrinotateWeb: Graphical Interface for Navigating Trinotate Annotations and Expression Analyses

https://github.com/Trinotate/Trinotate.github.io/wiki/TrinotateWeb

 

 

インストール

macだと一部のバイナリ配布プログラムが動作しないのでlinuxを使う。ubuntu18.04LTSにて、condaでpython3.7の仮想環境を作ってテストした。

本体 Github 

#bioconda(link) 依存が多いツールは仮想環境を作って導入するのが無難
conda create -n trinotate python=3.7 -y
conda activate trinotate
conda install -c bioconda -y trinotate

#transdecoderは別途インストールが必要
conda install -c bioconda -y transdecoder

#最近のバージョンではhmmerとblastは入っているはずだが、無いなら導入
conda install -c bioconda -y hmmer
conda install -c bioconda -y blast

Trinotateのヘルプ

autoTrinotate.pl

$ autoTrinotate.pl

 

##############################################################################

#

# Required:

#

#  --Trinotate_sqlite <string>                Trinotate.sqlite boilerplate database

#

#  --transcripts <string>                     transcripts.fasta

#

#  --gene_to_trans_map <string>               gene-to-transcript mapping file

#

#  --conf <string>                            config file

#

#  --CPU <int>                                number of threads to use.

#

#

##############################################################################

 

 

 RNAMMER、signalP、TMHMMはcondaでは入らない。それぞれダウンロードしてパスを通す。

まずHPからRNAMMER、signalP、TMHMMのプログラムをダウンロードする(ダウンロードを選ぶと記載したメールアドレスにダウンロードリンクが届く)。signalPはv5が出ているが、v4.1 の方をダウンロードする。

RNAMMER、signalP、TMHMM

https://services.healthtech.dtu.dk/software.php

 

signalPはダウンロード後初期設定が必要。signalp-4.1/signalPを開いてGENERAL SETTINGSのパスを修正する。ここでは/data/signalp-4.1/に実行スクリプトsignalpがあるので、以下のようにした。

f:id:kazumaxneo:20200916202557p:plain

FASTA.pmが見えないと言われたので
export PERL5LIB=/data/signalp-4.1/lib/で対処した。

 

TMHMMも初期設定が必要。中に入っているREADMEと TMHMM2.0.htmlに書いてあるが、tmhmmformat.plとtmhmmの1行目のperlパスを環境に合わせて修正する。また、decodeanhmm.Linux〜という実行形式ファイルが複数入っているので、動作するバージョンをdecodeanhmmにリネームする。ここではdecodeanhmm.Linux_x86_64をコピーしてリネーム。

> cp decodeanhmm.Linux_x86_64 decodeanhmm

 

rnammerも初期設定が必要。中に入っているrnammer-1.2.readmeに書いてあるが、$INSTALL_PATHのパスを環境に合わせて修正する。

f:id:kazumaxneo:20200917122740p:plain

他にも色々入れる必要がある。マニュアルの4ページ目によると、

http://arcsb2017.weebly.com/uploads/4/6/1/9/46199061/trinotate_install.pdf

hmmer2 を入れる(link)。

cd hmmer-2.3 
./configure
make
make check
make install

rnammerのhmmsearchパスを修正する。導入したhmmsearch2のパスを記載する。

f:id:kazumaxneo:20200917151700p:plain

環境変数$PATHのhmmsearchもhmmsearch2シンボリックリンクに書き換えた。

mv /root/miniconda3/envs/trinotate/bin/hmmsearch /root/miniconda3/envs/trinotate/bin/hmmsearch3

ln -s /data/hmmer-2.3/src/hmmsearch /root/miniconda3/envs/trinotate/bin/hmmsearch

core-rnammerの”--cpu 1”を消す。2箇所ある。

f:id:kazumaxneo:20200917125340p:plain

f:id:kazumaxneo:20200917125343p:plain

rnammerのテストランはしたほうがいい。

cd rnammer-1.2.src/example/
rnammer -S bac -m lsu,ssu,tsu -gff - < ecoli.fsa

出力

##gff-version2

##source-version RNAmmer-1.2

##date 2020-09-17

##Type DNA

# seqname           source                      feature     start      end   score   +/-  frame  attribute

# ---------------------------------------------------------------------------------------------------------

ecoli_section RNAmmer-1.2 rRNA 21067 21181 86.3 + . 5s_rRNA

ecoli_section RNAmmer-1.2 rRNA 16177 17706 1950.6 + . 16s_rRNA

ecoli_section RNAmmer-1.2 rRNA 18068 20969 3754.0 + . 23s_rRNA

# ---------------------------------------------------------------------------------------------------------

O.K! 

もしXML::Simpleが無いなら、perl lib/にXML::Simpleをダウンロードして以下コマンドでビルドする。cpanmで入れると、gccのバージョン指定ミスなのかコンパイル時のチェックでこける。

https://metacpan.org/pod/XML::Simple

cd XML-Simple-2.25/
perl Makefile.PL
make
make install

 エラーがないか改めてテストする。

rnammer -S bac -m lsu,ssu,tsu -gff - < rnammer-1.2.src/example/

#自分のデータでもやったほうがいい。
rnammer -S bac -m lsu,ssu,tsu -gff - < tour_RNA_aasembly.fasta

 出力される.gffでrRNAが検出されているか確認する。ないとは思うが、TrinotateはrRNAがないデータだと、とrnammerのステップでエラー扱いになって止まる。

 

 

 データベースの準備

1、sqliteデータベースその他をダウンロードする。

Build_Trinotate_Boilerplate_SQLite_db.pl Trinotate

ダウンロードされた。

f:id:kazumaxneo:20200131193613p:plain

2、ダウンロードされたuniprot_sprot.pepを元に blastデータベースを作成する。

makeblastdb -in uniprot_sprot.pep -dbtype prot

 3、続いてpfamデータベースを作成する。

#解凍
gunzip Pfam-A.hmm.gz

#index作成
hmmpress Pfam-A.hmm

f:id:kazumaxneo:20200131193817p:plain

 

レポジトリのスクリプトが必要。また、configファイル のテンプレートが置いてあるので、これをベースに修正していく。

git clone https://github.com/Trinotate/Trinotate.git
cp Trinotate/auto/conf.txt .

conf.txt

f:id:kazumaxneo:20200916121356p:plain

 

/dataにダウンロードしたプログラムsignalp-4.1、tmhmm-2.0c、rnammer-1.2.src、trinityrnaseqのレポジトリ、上のデータベースを配置して、以下のように修正した。condaで入れたツールはcondaのパスを記載。

f:id:kazumaxneo:20200917120133p:plain



 

 

実行方法

1、Trinityアセンブリ配列以外にgene-to-transcript mapping fileが必要。Trinityrnaseqのユーティリティで作成する(下記参照)。

git clone https://github.com/trinityrnaseq/trinityrnaseq.git

#trinityのアセンブリ結果のtrinity.fastaを指定する
perl trinityrnaseq/util/support_scripts/get_Trinity_gene_to_trans_map.pl trinity.fasta > gene_trans_map

> head gene_trans_map

$ head gene_trans_map

TRINITY_DN144_c0_g1 TRINITY_DN144_c0_g1_i1

TRINITY_DN179_c0_g1 TRINITY_DN179_c0_g1_i1

TRINITY_DN159_c0_g1 TRINITY_DN159_c0_g1_i1

TRINITY_DN159_c0_g2 TRINITY_DN159_c0_g2_i1

TRINITY_DN153_c0_g1 TRINITY_DN153_c0_g1_i1

TRINITY_DN153_c0_g2 TRINITY_DN153_c0_g2_i1

TRINITY_DN130_c0_g1 TRINITY_DN130_c0_g1_i1

TRINITY_DN106_c0_g1 TRINITY_DN106_c0_g1_i1

TRINITY_DN106_c0_g2 TRINITY_DN106_c0_g2_i1

TRINITY_DN154_c0_g1 TRINITY_DN154_c0_g1_i1

見ればわかるが、ヘッダーを少し変えただけ。このスクリプトに頼らずとも、配列のheaderを少し改良するだけで作成できる(*1)。

 

2、Trinotateを実行する。trinity.fastaと1で作成したgene_trans_map、configファイルを最低限指定する必要がある。

#trinityのアセンブリ結果のtrinity.fastagene_trans_mapを指定する
autoTrinotate.pl --Trinotate_sqlite Trinotate.sqlite --transcripts Trinity.fasta --gene_to_trans_map gene_trans_map --conf config.txt --CPU 12

 

 

ラン中にどこかのステップでエラーが起きたら、こちら

http://arcsb2017.weebly.com/uploads/4/6/1/9/46199061/trinotate_install.pdf

を確認してエラーが起きたステップのコマンドを実行し、エラーが起きるか調べてください。自分はsignalPのステップでエラーが起きた時、上記リンク先のコマンドを打ってオプションが認識されないことでsignalPのバージョンが違うことに気づきました。

TrinotateレポジトリのTrinotate/sample_data/data/Trinity.fastaはexampleデータですが、これを使ってうまくランするかテストすると、私の環境ではrnammerのランでエラーになるため使えませんでした(rnammerは正常に動いており、実際、rnammerだけ動かすと他のデータではrRNAを予測できても、このファイルではrRNAを予測できない)。注意して下さい。

追記

ラン最後のsqliteデータベースを作るところでTrinotateコマンドとextract_GO_assignments_from_Trinotate_xls.pl, annotation_importerがないとエラーが出てしまった。Trinotateのパスを調べてエラーの通りのパスにシンボリックリンクを置くことで対処した。

ln -s /root/miniconda3/envs/trinotate/bin/Trinotate /root/miniconda3/envs/trinotate/bin/..//

mkdir /root/miniconda3/envs/trinotate/bin/..//util
ln -s /root/miniconda3/envs/trinotate/bin/extract_GO_assignments_from_Trinotate_xls.pl /root/miniconda3/envs/trinotate/bin/..//util/

mkdir /root/miniconda3/envs/trinotate/bin/..//util/annotation_importer
ln -s /root/miniconda3/envs/trinotate/bin/import_transcript_names.pl /root/miniconda3/envs/trinotate/bin/..//util/annotation_importer/

Trinotateはチェックポイントを作っています。どこかのステップでランに失敗した場合、そのステップのコマンドが動作するようにして、失敗した同じコマンドを打てば、そのステップから再開でき、無駄な時間を取られずに済みます。

追記

Error, cannot find obo/go-basic.obo.gz

が出たら、ここを参照して修復する。

cd /root/miniconda3/pkgs/perl-extutils-makemaker-7.36-pl526_1/lib/5.26.2/x86_64-linux-thread-multi/
mkdir obo
cd obo/
cp <your_path>/Trinotate/PerlLib/obo/go-basic.obo.gz .
gzip -dv go-basic.obo.gz

追記

入力配列の数が多すぎると最後のデータベース作成でエラーになることがあるようです。50000配列(61MB)ではエラーを起こしましたが、半分ずつに分けてランすると最後までエラーを吐かずに動かすことができるようになりました。 

引用

The De Novo Transcriptome and Its Functional Annotation in the Seed Beetle Callosobruchus maculatus

Ahmed Sayadi, Elina Immonen, Helen Bayram, and Göran Arnqvist*

PLoS One. 2016; 11(7): e0158565. Published online 2016 Jul 21.

 

関連 


 

コメント

dockerイメージもdockerhubで探せば見つかります。導入に苦戦している方は探してみてください。

 

*1

https://github.com/Trinotate/Trinotate.github.io/issues/7

例えば

>TRINITY_DN144_c0_g1_i1 len=231 path=[227:0-162 228:163-230] [-1, 227, 228, -2]

というヘッダー名なら、

TRINITY_DN144_c0_g1 TRINITY_DN144_c0_g1_i1

にする。もう少しわかりやすく言うと、ヘッダーの赤字の部分

TRINITY_DN144_c0_g1_i1

をtab区切りで二回書く。2回目は"_i1"をつける。すなわち

TRINITY_DN144_c0_g1 TRINITY_DN144_c0_g1_i1

にする。これを全配列に対して実行するスクリプトを書けばよい。