8. 翻訳とblastとアノテーション
2021.06.09 初稿
TransDecoderによる翻訳
フィルタリング後のfasta配列を用いて、TransDecoderによりORFを予測し、アミノ酸に翻訳する。
---------------------------------------------
#!/bin/bash
#$ -S /bin/sh
#$ -l s_vmem=64G,mem_req=64G
#$ -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/transdecoder:5.5.0--pl526_2 \
TransDecoder.LongOrfs -m 100 -t filtered_Trinity.fasta
echo ending at
date
---------------------------------------------
longest_orfs.pepというファイルがORFとして予測された配列を翻訳したもの。.p1, .p2...というように一つのアイソフォームに対してたくさんの翻訳されたバージョンができるので最長配列だけを取り出す。
blastをかける
TransDecoderで翻訳した遺伝子（ORF）がどのような遺伝子をコードしているのかを知る必要がある。そのために行うのがblast。
自分の材料の近縁種でゲノムが読まれているもの、もしくは材料の分類群でモデル生物と言われているもの、これまでに読まれた全生物、などどの生物のゲノムをreferenceとするかは目的によって異なる。
また、blastp、blastx、tblastxのどれを行うかも目的によって異なる。
ここではよく行われているblastxのやり方を提示する。
まずは遺伝研スパコン内にblastがあるか確認。
$ ls /usr/local/seq /usr/local/biotools/b/blast*
多数のバージョンのblastがあることがわかる。
Trinityの時もそうだったが、最新のものだとうまくいかないことがあるので２番目に新しいものを使うことにする。
その後、「遺伝研スパコン blast」とGoogleで検索。
以下のページ（https://sc.ddbj.nig.ac.jp/ja/guide/database#info1）によるとNCBI-nr、UniProt、Pfam、Genbank、EMBLデータベースのデータを全てスパコンから利用できる。
blastデータベースのパスは、/usr/local/biotools/b/blast/ncbi/nr.000(数字).*(アルファベット)
以下で中身を確認。
$ ls /usr/local/biotools/b/blast/ncbi/nr*
ここでたくさん出てくる「*.phd」「*.phi」「*.phr」「*.pin」というのはタンパク質配列のfastaファイルにインデックスをつけたもので、blastxやblastpを行う場合に使用できる。塩基の場合は「*.nhd」「*.nhi」「*.nhr」「*.nin」となる。
スパコンにおいて、blastのオプション確認をするときは以下。
$ singularity exec --bind /usr/local/seq /usr/local/biotools/b/blast:2.9.0--pl526h979a64d_3 blastx --help
最低限必要なコマンドは以下。-outfmtはデフォルトだと7になっている。
blastx \
-max_target_seqs 1 \
-query ${HOME}/RNA-seq/Trinity_filt.fa \
-db reference name \
-outfmt 6 \
-evalue 1e-06 \
-out out_blast.txt
このコマンドを遺伝研スパコンで使える仕様にする。
--bind /usr/local/seqをつけることでスパコン内のデータベースにアクセスできるようになる
（参照：https://sc.ddbj.nig.ac.jp/ja/guide/database#info1）。
データベース名は/usr/local/seq/の後、「利用可能DB一覧」に示されたパスを参照する。
例えばnr databaseをreferenceに使う場合は「/usr/local/seq/blast/ncbi/nr」と指定すれば良い。Uniprotの場合は/usr/local/seq/flat/uniprot。
※ちなみにこのスクリプトをランするとNCBIの遺伝子IDは得られるが遺伝子名や遺伝子機能（gene description）は得られず、後から自分で付け足す必要がある。NCBI Batch Entrez（https://www.ncbi.nlm.nih.gov/sites/batchentrez）というページで、得られたprotein IDを入れることでgene descriptionを得ることができる。
singularity exec --bind /usr/local/seq /usr/local/biotools/b/blast:2.9.0--pl526h979a64d_3 blastx \
-query ${HOME}/RNA-seq/Trinity_filt.fa \
-db /usr/local/seq/blast/ncbi/nr \
-evalue 1e-06 \
-outfmt 6 \
-out out_nr.txt
NCBIのデータベースではなく、遺伝研が独自に作成したデータベースを使用することもできる。
DDBJのBLASTについての説明ページ（https://www.ddbj.nig.ac.jp/services/blast.html）によると、検索ページ（http://ddbj.nig.ac.jp/blast/）においてはblastするDIVISIONを選ぶことが可能となっている。スパコンからはこのDIVISIONで使用されているデータベースにアクセスすることができ、それは /usr/local/resources/trad/ddbjnew/unified-all/blastdb/に存在する。例えば私が、DIVISIONの中のPlant (PLN)をreferenceとして使用したい場合には以下のコマンドを使う。
singularity exec --bind /usr/local/seq /usr/local/biotools/b/blast:2.9.0--pl526h979a64d_3 \
tblastx -query ${HOME}/blast/my_file.fa \
-db /usr/local/seq/blast/ddbj-unified-all/pln \
-evalue 1e-04 \
-outfmt 6 \
-out out_blast.txt
自分でreferenceとなるデータベースを作成することももちろん可能。まず、対象となる生物のデータを入れるディレクトリを作成。
$ mkdir arabi_ref
NCBI等のサイトで欲しい生物種の塩基配列もしくはタンパク質配列のfastaファイルをダウンロードする。例えばシロイヌナズナのゲノム配列をダウンロードするとする（もちろんNCBIではなくTAIRを使ってもいい）。
NCBIでgenomeのカテゴリを指定し、「arabidopsis thaliana」を検索。（添付画像）
RefSeqをクリックし、GCF_000001735.4_TAIR10.1/ をクリック。
たくさんの形式が異なるファイル名が現れるが、使用するのは「GCF_000001735.4_TAIR10.1_genomic.fna.gz」もしくは「GCF_000001735.4_TAIR10.1_protein.faa.gz」。
欲しいファイル名にカーソルを合わせ、右クリックして「リンクをコピー」。そうするとhttpsアドレスが取得できる。
それをwgetで取得。
ファイルが圧縮されているので解凍
$ gunzip GCF_000001735.4_TAIR10.1_genomic.fna.gz
$ ls -l
でファイルサイズを確認。
ファイル名が長いのでmvで変更。
$ mv GCF_000001735.4_TAIR10.1_genomic.fna tair_genomic.fna
以下のコマンドでblastのためのreference databaseを作成
$ makeblastdb -in tair_genomic.fna -out tair_genome -dbtype nucl -hash_index
Referenceが塩基の場合はtblastxを行うので以下のようなスクリプトとなる。
-------------
singularity exec --bind /usr/local/seq /usr/local/biotools/b/blast:2.9.0--pl526h979a64d_3 \
tblastx -query ${HOME}/RNA-seq/Trinity_filt.fa \
-max_target_seqs 1 \
-db ${HOME}/arabi_ref/tair_genome \
-evalue 1e-05 \
-outfmt 6 \
-out out_blast.txt
-------------
同様にNCBIもしくはTAIRのサイトからタンパク質配列のfastaファイルをダウンロードし、blastxもしくはblastpをかけることも可能。
Trinity.fastaファイルは非常に重たいため、いくつかに分割し、さらにアレイジョブとして走らせると早く終わる。
例えば50,000 transcriptがあるTrinity.fastaを10,000ずつ5つのファイルに分ける場合、以下により最初の10,000行をTrinity_temp1.fastaに出力。
head -n 10000 Trinity.fasta > Trinity_temp1.fasta
続いて、次の10,000行をTrinity_temp2.fastaに出力。
head -n 20000 Trinity.fasta | tail -n 10000 > Trinity_temp2.fasta
というパターンを繰り返して5ファイルに分割できる。手順をまとめると以下。
head -n 10000 Trinity.fasta > Trinity_temp1.fasta
head -n 20000 Trinity.fasta | tail -n 10000 > Trinity_temp2.fasta
head -n 30000 Trinity.fasta | tail -n 10000 > Trinity_temp3.fasta
head -n 40000 Trinity.fasta | tail -n 10000 > Trinity_temp4.fasta
head -n 50000 Trinity.fasta | tail -n 10000 > Trinity_temp5.fasta
※seqkitを使うともっと早くできます。コマンドは以下。
---
singularity exec -e /usr/local/biotools/s/seqkit:2.0.0--h9ee0642_0 \ seqkit split -p 5 Trinity.fasta
---
スパコンでアレイジョブを投げるスクリプトは以下。
-------------------
#!/bin/bash
#$ -S /bin/sh
#$ -t 1-8
#$ -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 starting at
date
singularity exec --bind /usr/local/seq /usr/local/biotools/b/blast:2.9.0--pl526h979a64d_3 \
blastx -query ${HOME}/RNA-seq/split_${SGE_TASK_ID}.fa \
-max_target_seqs 1 \
-db /usr/local/seq/blast/ncbi/nr \
-evalue 1e-05 \
-outfmt 6 \
-out out_${SGE_TASK_ID}.txt
echo ending at
date
-------------------
上記ではnrデータベースをreferenceにした場合のものだが、nrでかけたときは3回やって3回とも途中でcoredumpエラーが出て止まっていたので途中で断念...。メモリ量が不足していたかも。UniProtでのblastに切り替え。
アノテーション付け
Trinityでde novo assembleをした後だけでなく、Hisat2でマッピングが終わった後も遺伝子アノテーションをつける必要がある。
ゲノムが読まれている場合、gtfファイルもしくはアミノ酸配列のfastaファイルに遺伝子情報も載っている。
このアノテーションを自分のファイルに結合する。
また、ゲノムが読まれているとはいえGOがついていない種もある（例えば、私がこれまで使用してきたO. longistaminataやC. campestrisなど）。この後の「10. GO解析」を行うにはGOがついているモデル植物の遺伝子IDが必要になる（私の例で言えばO. sativaやA. thaliana）。そこでモデル植物の遺伝子IDも紐付けることが必要となる。
自分の材料に対してモデル植物の遺伝子IDを紐付けるにはいくつか手法はあると思うが、一番手っ取り早いのはblastをおこなった後にトップヒットだけ取ってくることだろう。
Blastを自前でやって後から遺伝子名や遺伝子機能をくっつける作業をするなんてめんどくさい。全部一気に行ってくれるソフトはないのか？
と探していたところ、見つけたのがこちら→ EggNOG-mapper (http://eggnog-mapper.embl.de)
eggNOG-mapperにTransDecoderで翻訳したlong_orfs.faをfa.gzの形に圧縮して投げるだけでOK。数十分で遺伝子IDおよび遺伝子機能、GOのついた結果が返ってくる。
