2020/05/04 公開

Transdecoder による遺伝子予測とアノテーション


Transdecoder は、Trinity でアセンブルされたcontigから、アミノ酸配列をコードしている領域を推定できるソフトウエアである。hmmerと組み合わせて、どの遺伝子ファミリーに属する配列であるか推定することができるので、遺伝子のアノテーションにも役立つ。

このマニュアルでは、コマンド入力をするところという意味で「:>」という記号を行頭につける。この記号は除いて、後のコマンドを改行せずに一行で入力する。
 

用語解説


contig (コンティグ): アセンブルした結果出来上がった、ひと繋がりの塩基配列。

hmmer (ハンマー or ハマー): 隠れマルコフモデルを用いて類似する配列を探索するソフト。アミノ酸ファミリーのデータベースであるPfamから、関連するドメイン構造を検索する際に有効なソフト。本マニュアルの最初にある環境構築で、Anacondaからインストールしてあるはずである。
:> conda list | grep hmmer
として確認する。インストールされていない場合は、
:> conda install hmmer

seqkit(セックキット): シークエンスデータを取り扱うのに便利なソフト。BiorubyやBiopythonなどを用いても同様のことが可能ではあるが、個人的にはseqkitが最も作業手順が少なくて便利だと思う。
 

準備


uniprot のデータをダウンロード
https://www.uniprot.org

にアクセスして、

Download latest release
から
Reviewed (Swiss-Prot) の fasta をダウンロードする
 

データベースフォルダーを作る

:> cd ~/Analysis
:> mkdir -p /DB/BLASTDB


データベースの場所を設定する

:> cd ~/
:> echo 'export BLASTDB=/Users/アカウント名/Analysis/DB/BLASTDB' >> .zshrc



uniprot_sprot.fasta をBLASTDBに移動して、BLAST検索用のデータベースを構築

:> cd ~/Analysis/DB/BLASTDB
:> makeblastdb -in ./uniprot_sprot.fasta -dbtype prot -out uniprot_sprot



Pfam のデータをダウンロード
https://pfam.xfam.org

Pfam-A.hmm.gz

をダウンロード。この記事を制作している時点では、32が最新版。

データベースフォルダーを作る

:> cd ~/Analysis/DB
:> mkdir Pfam32


ダウンロードした Pfam-A.hmm.gz を Pfam32 に移動して展開

HMMER で検索できるようにデータベースを構築

:> cd ~Analysis//DB/Pfam32
:> hmmpress ./Pfam-A.hmm

 
ここまでの作業で、BLASTによるuniprot検索、hmmerによるPfam検索の準備が整った。
BLASTDB内にnr等のデータを用意すれば、下記と同じ方法で検索が可能になる。
 

解析0


Trinityの結果内に、解析用のフォルダーを作る

例: :> mkdir Transdecoder


作成したフォルダーに移動する

:> cd Transdecoder

 

解析1


最長のORFを抜き出す

:> TransDecoder.LongOrfs -t ../Trinity.fasta

 

解析2


Pfamで遺伝子ドメインを推測

:> hmmscan --cpu 8 -E 1e-10 --domE 1e-10 -o pfam.out --domtblout pfam.out_domtblout ~Analysis/DB/Pfam32/Pfam-A.hmm ./Trinity.fasta.transdecoder_dir/longest_orfs.pep

 

--cpuの数は、自分が使っているMacのCPU数に合わせて変更する

 

解析3


Blastで遺伝子名を推測

:> blastp -query ./Trinity.fasta.transdecoder_dir/longest_orfs.pep -db uniprot_sprot -max_target_seqs 1 -outfmt 6 -evalue 1e-10 -num_threads 8 -out blastp.outfmt6

 

--num_threadsの数は、自分が使っているMacのCPU数に合わせて変更する

 

解析4


PfamとBlastの結果を参照して遺伝子を予測

:> TransDecoder.Predict -t ../Trinity.fasta --retain_pfam_hits ./pfam.out_domtblout --retain_blastp_hits ./blastp.outfmt6 --cpu 0

 

--cpuの数を0にすることで、使える全てのCPUを使ってくれる



これで、タンパク質をコードしている範囲(.cds)とアミノ酸配列に変換したFastA(.pep)ができる。
 

解析5


遺伝子ファミリーの推定

:> hmmscan --cpu 8 -E 1e-10 --domE 1e-10 -o pfam_Final.out --domtblout pfam_Final.out_domtblout ~Analysis/DB/Pfam32/Pfam-A.hmm ./Trinity.fasta.transdecoder.pep

 

--cpuの数は、自分が使っているMacのCPU数に合わせて変更する

 

解析6


ターゲット遺伝子のアミノ酸配列を抜き出す
(昆虫味覚受容体ファミリー; 7tm_7 の遺伝子を抜き出したい場合)
(下記の7tm_7のところを変更すれば、自身が興味のある遺伝子ファミリーの配列が取り出せる)

:> mkdir Genes
:> cd Genes
:> grep 7tm_7 ../pfam_Final.out_domtblout | awk '{print $4}' > 7tm_7_list.txt
:> for file in `cat ./7tm_7_list.txt`;do seqkit grep -p $file ../Trinity.fasta.transdecoder.pep;done > 7tm_7_pep.fasta

 
組織間の発現量の比較を行わず、新規遺伝子を探索するのであればここまでで完了。
 

コマンドの解説


grep: 指定したファイルからキーワード検索を行う

grep キーワード ファイル


awk スペース等で区切られた文字列から、任意の部分を取り出す

awk '{print $4}'

スペース区切りの文字列から、4つ目のワードを表示
$の数字を変えて、必要な文字列を取り出せる


seqkit grep -p: パターン一致で配列データを検索する

使用例: seqkit grep -p 検索する文字列 検索対象ファイル

目次

  1. 環境構築
  2. リードのクオリティチェックとアセンブル
  3. Transdecoder による遺伝子予測とアノテーション
  4. リード数のカウント
  5. カウントマトリクスの作成
  6. リードカウントのQC解析
  7. DE解析
  8. DEG取り出し
  9. おまけ1: Transdecoderの自動化
  10. おまけ2: 発現量解析の自動化

Q&A

皆様からいただいた質問と回答など、こちらへ掲載します。

質問やご意見はこちらへ

※ボタンのリンク先は、皆で語り合う生命誌研究館の掲示板「みんなの広場」です。ご投稿いただいたご質問やご意見、館員からのお返事は公開されますのでご了承ください。