case4 - nibb-gitc/gitc2025mar-rnaseq GitHub Wiki

Quick Annotation of Trinity Assembly Using BLAST

ダイコン(Raphanus sativus var. hortensis (long and thick root radish)) のRNA-seq データをTrinityでアセンブルした。Trinityで得られたダイコンのコンティグはそれぞれがどのような遺伝子をコードしているだろうか？ここでは、ダイコンのコンティグを対象に、BLASTを使った簡易アノテーションを学ぶ。

BLASTによる相同性検索はおおまかなアノテーションを行なうのに便利な手法である。ここでは、シロイヌナズナのタンパク質データベースを検索することにより、各コンティグがシロイヌナズナのどのタンパク質に対応するかを調べる。ダイコンとシロイヌナズナはどちらもアブラナ科に属するので、進化距離が近く配列も比較似ていると考えられ、配列によるアノテーションは効果的であると期待できる。

Data

必要なデータは、"~/gitc/data/EX/case4/" ディレクトリにある。

  • Trinity_Daikon.fasta -- ダイコンのIllumina RNAseq からTrinityによってde novo assembly したコンティグ配列

（参考）このTrinity_Daikon.fasta は、以下の論文で用いられているデータの一部（accession#: DRR010353, DRR014773）を使って、今回のコースのために新たに生成したものである。

Mitsui, Y. et al. (2015). The radish genome and comprehensive gene expression profile of tuberous root formation and development. Scientific reports 5(1), 10835. https://dx.doi.org/10.1038/srep1083

Software

  • NCBI BLAST+
  • PythonとRubyのスクリプトを実行する

Setup

case4 ディレクトリを作成し、その下で作業しよう。

$ mkdir case4
$ cd case4

Trinity アセンブルファイルをカレントディレクトリにコピーしておこう。

$ cp ~/gitc/data/EX/case4/Trinity_Daikon.fasta ./

less などで、Trinity_Daikon.fasta の中身を確認しよう。エントリー数を確認しておくとさらに良いだろう。

Build BLAST DB

国際コンソーシアムの運営するシロイヌナズナデータベース TAIRから、シロイヌナズナのタンパク質アミノ酸配列セットをダウンロードする。

ダウンロード

（このファイルは、~/gitc/data/EX/case4/ ディレクトリにもコピーしてある）。

ダウンロードしたファイルを''TAIR10.pep''の名前に変更。

$mv TAIR10_pep_20110103_representative_gene_model_updated TAIR10.pep

BLAST DBをビルド。

$ makeblastdb -in TAIR10.pep -dbtype prot -parse_seqids

Similarity search with BLAST

$ blastx -query Trinity_Daikon.fasta -db TAIR10.pep -num_threads 4 \
       -evalue 1.0e-8 -outfmt 6 > blastx_results.txt

計算に時間がかかるため、参考のため、あらかじめ計算済みの結果をdataディレクトリに保存しておいた。~/gitc/data/EX/case4/blastx_results.txt

less などで blastx_results.txt の中身を確認しておこう。

一つのqueryに対し複数のヒットが記録されているのがわかる。下流の解析を簡単にするために、トップヒットのみを抽出する。 そのためにsortコマンドを駆使する（内山先生の講義資料のp14[BLASTタブ区切り出力からのベストヒット抽出]を参照）[https://github.com/nibb-gitc/gitc2024jul-rnaseq/blob/main/textbook/8_GO%E8%A7%A3%E6%9E%90.pdf]。

sort -k1,1 -s -u blastx_results.txt > blastx_results_tophit.txt

〜〜〜参考〜〜〜

同じことを実現する（トップヒットのみを抽出する）、python scriptも示しておく。

imput_file = "blastx_results.txt"

prev_id = ""
f = open(imput_file)

for line in f:
    ary = line.split("\t")
    curr_id = ary[0]
    if curr_id == prev_id:
        pass
    else:
        print(line.rstrip("\n") )
    prev_id = curr_id

f.close()

これを、"select_first_hit.py" というファイル名で保存しよう。（dataフォルダにも置いてある）。

$ python select_first_hit.py > blastx_results_tophit.txt

=> blastx_results_tophit.txt

〜〜〜参考終わり〜〜〜

hit率を計算してみよう。inputのTrinity_Daikon.fastaにはいくつのエントリーがあって、そのうちいくつがシロイヌナズナのタンパク質にヒットしただろうか。

$ seqkit stat -a Trinity_Daikon.fasta
file                  format  type  num_seqs     sum_len  min_len  avg_len  max_len   Q1     Q2     Q3  sum_gap    N50  Q20(%)  Q30(%)
Trinity_Daikon.fasta  FASTA   DNA     23,167  35,138,830      500  1,516.8   16,191  818  1,269  1,888        0  1,801       0       0

$ wc blastx_results_tophit.txt
  21173  254076 1664743 blastx_results_tophit.txt

21173 / 23167 = 0.9139

Ans. 91.4%

Similarity search with diamond

上記のBLAST検索はかなり時間を要する。diamondを使えば、検索精度をほとんど犠牲にすることなく超高速に検索することが可能である。ここではdiamondを使ってみよう。

diamond makedb --in TAIR10.pep --db TAIR10.dmnd --threads 8

diamond blastx --query Trinity_Daikon.fasta \
  --db TAIR10.dmnd --evalue 1.0e-8 --max-target-seqs 10 \
  --outfmt 6 --threads 8 --out diamond_results.txt

驚くほど高速に結果が得られる。

出力フォーマットは--outfmtオプションでBLASTのformat 6 と同じ形式を指定している。従って得られたテーブルは上記で説明したのと全く同じ方法で処理できる。

sort -k1,1 -s -u diamond_results.txt > diamond_results_tophit.txt
wc diamond_results_tophit.txt
  21011  252132 1639624 diamond_results_tophit.txt

map率はblastxと比べてどうだろうか？

21011 / 23167 = 0.9069

diamondは90.7%。blastxの91.4% に比べ僅かに下回るだけで、高速性によって検索の感度はほとんど犠牲になっていないといえる。

Populate annotations

blastx_results_tophit.txt （およびdiamond_results_tophit.txt）は、ヒットしたAdabidopsisのIDが記録されているのみなので、具体的になんという名前のどのような機能の遺伝子なのかはいちいちデータベースを検索しなければわからない。これらの情報を追記したいというニーズは大きいだろう。

モデル生物の場合、大半の遺伝子に詳細なアノテーションがついているので、それらの情報を取り入れることを考える。シロイヌナズナの場合、各遺伝子のfunctional annotationは以下のファイルにまとめられており、TAIRのウェブサイトからダウンロードすることができる。

（このファイルは、~/gitc/data/EX/case4/ ディレクトリにもコピーしてある）。

blastx_results_tophit.txt に TAIR10_functional_descriptions のアノテーション情報を付加する短い Ruby scriptを用意した。

#=== conf ===
result_blast = "blastx_results_tophit.txt"
gene_annotation_file = "TAIR10_functional_descriptions"
#===

## load gene annotation
data_annot = {}
File.open(gene_annotation_file).each_with_index do |l, i|
  next if i == 0
  a = l.chomp.split(/\t/)
  gene = a[0]
  data_annot[gene] = a
end

## load blast tophit
File.open(result_blast).each do |l|
  a = l.chomp.split(/\t/)
  id = a[0]
  hitid = a[1]
  out = [a, data_annot[hitid]].flatten.join("\t")
  puts out
end

これを　add_TAIRdesc_to_blast6.rb　という名前で保存し、以下のように実行する。

$ ruby add_TAIRdesc_to_blast6.rb > blastx_results_tophit_annot.txt

スクリプトの中の、

result_blast = "blastx_results_tophit.txt"

result_blast = "diamond_results_tophit.txt"

に変更すれば、diamondの結果も処理できる。

タブ区切りテキストになっているので、MS Excelなどのスプレッドシートで開くと見易いだろう。

Links

BLASTは汎用性が高くそして奥の深い配列解析ツールである。私たちは基生研ゲノムインフォマティクストレーニングコースの一つとして、BLAST自由自在を不定期に開催している。以下に最近のコースの資料を公開しているので是非参考にしていただきたい。

Revision History

2024-7-31

  • 動作確認 on bias5 (NCBI BLAST: 2.13.0; diamond 2.0.15)

2023-2-25

  • 動作確認 on bias5
  • NCBI BLAST: 2.13.0; diamond 2.0.15

2022-8-28

  • 動作確認。軽微な修正。(NCBI BLAST: 2.13.0; diamond: 2.0.12)

2021-9-12

  • TAIR10のファイルのリンクをアップデート
  • diamond での解析を追加。

2021-3-XX