case4 - nibb-gitc/gitc2021mar-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から、シロイヌナズナのタンパク質アミノ酸配列セットをダウンロードする。

ダウンロード

ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_blastsets/TAIR10_pep_20110103_representative_gene_model_updated

（このファイルは、~/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コマンドを駆使する（内山先生の講義資料のp11を参照）。

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%

Populate annotations

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

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

ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_functional_descriptions

（このファイルは、~/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

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