![[Share]Twitter](/common/img/lab/sns_tw.png){kind=small}
![[Share]Facebook](/common/img/lab/sns_fb.png){kind=small}
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 検索する文字列 検索対象ファイル
目次
- 環境構築
- リードのクオリティチェックとアセンブル
- Transdecoder による遺伝子予測とアノテーション
- リード数のカウント
- カウントマトリクスの作成
- リードカウントのQC解析
- DE解析
- DEG取り出し
- おまけ1: Transdecoderの自動化
- おまけ2: 発現量解析の自動化
- おまけ3: Anacondaコマンドの使い方一覧