ex801 Functional annotation using similarity-based methods

演習801: ホモロジー検索を用いた機能アノテーション

ここでは、近縁種のモデル生物で、良質のアノテーションがついたゲノムデータが利用可能であるケースを想定して、ホモロジー検索等を用いてアノテーションを行ってみる。使用するデータは、酵母 Saccharomyces eubayanus のトランスクリプトームデータ( GEO accesion: GSE133146 )である。この酵母は、ラガービールの生産に使われる酵母 S. pastorianus の祖先種の一つで、S. pastorianus はこの酵母とS. cerevisiaeのハイブリッドによって成立したとされており、その起源に迫ることが元の研究の目的となっている。この研究ではハイブリッドした祖先種に近いと考えられるS. eubayanusのヒマラヤ株を対象として、そのゲノムシーケンスも行っているが、ここでは公開されている別株のゲノムを用いてマッピングし、主に S. cerevisiae との比較に基づいてアノテーションを行うという想定で解析を行ってみよう。なお、ここではアノテーションまでを行い、それに基づくエンリッチメント解析は演習803で行う。

データ

bias5で作業する。 ディレクトリ~/gitc/data/IU/yeastに、以下のファイルがある。

file	contents	remarks
seub_genome.fa	S. eubayanus ゲノム配列	入力データ
stringtie_merged.gtf	stringtieによるアセンブル結果	入力データ
topTags.gene.txt	EdgeRの結果（トランスクリプトレベル）	入力データ／演習803で用いる
topTags.gene.txt	EdgeRの結果（遺伝子レベル）	入力データ／演習803で用いる
seub_genes.pep	S. eubayanus 遺伝子翻訳配列	中間結果
seub_genes6.pep	S. eubayanus 遺伝子翻訳配列縮小版	中間結果（一部）
scer_prot.fa	S. cerevisiae 遺伝子翻訳配列	公的データから取得
blastout.tab	S. eub x S. cer BLAST結果	中間結果
seub_genes6.iprscan.tsv	seub_genes6に対する InterProScan結果	出力結果（一部）
seub_genes.emapper.annotations	seub_genesに対するEggNOG-mapper結果	出力結果／演習803で用いる

状況としては、ゲノム配列に対してRNA-seqリードをマッピングしてゲノムベースでアセンブルし、トランスクリプトごとに頻度をカウントして、EdgeRで有意な発現変動を示すトランスクリプト／遺伝子を抽出するところまでが終わっていることを想定している。すなわち、初期状態として最初の４つのファイルが存在しており、これからアノテーションをおこなっていく。途中の過程をスキップできるように、いくつか中間結果のファイルも置いてある。これらを上書きしないように、別のディレクトリを作成して、そこにファイルをコピーして作業することを勧める。

% cd ~/gitc/data/IU/
% mkdir ex1
% cp yeast/* ex1
% cd ex1

Step 1: トランスクリプト配列の作成

stringtieでは、トランスクリプトの座標を記録したGTFファイルを作成するが、配列ファイルは作成しないので、まずこのGTFファイルとゲノム配列からトランスクリプト配列を作成する必要がある。これはgffreadコマンドで行う。

% gffread stringtie_merged.gtf -g seub_genome.fa -w seub_transcripts.fa

Step 2: コード領域(CDS)の推定と翻訳配列の作成　(講義ではスキップする予定)

作成したトランスクリプト配列からCDSを推定する。これは、TransDecoderを用いて行うが、長いORFを抽出する段階と、そこからCDSを推定する段階の２段階で行う。

% TransDecoder.LongOrfs -t seub_transcripts.fa
% TransDecoder.Predict -t seub_transcripts.fa --single_best_only --cpu 8

デフォルトでは、一つのトランスクリプト中に重複しない複数のCDSが同定された際は、それらを複数のCDSとして出力するが、後処理がやや面倒になるので、ここではそういった場合には一つのCDSのみを出力するオプションを指定している。

出力結果として、seub_transcripts.fa.trandecoder.pep （アミノ酸配列）のほか、ファイル名のサフィックスが .cds （塩基配列）および　.gff3（トランスクリプト中のCDSの座標）の３つのファイルが作成される。アミノ酸配列ファイルをこの後の解析に用いるが、遺伝子名が、元のGTFファイル中のtranscript_idから少し変わっている。

>DI49_1142.p1 GENE.DI49_1142~~DI49_1142.p1  ORF type:complete len:236 (-),score=50.10 DI49_1142:234-941(-)
MLPLIASRNRRPISLTIRKLFRTMSIVKGKPEEAKIVEARHVKDTSDCKWIGLQKIIYKD
PNGNEREWDSAVRTTRNSGGVDGIGILTILKYKDGKPDEILLQKQFRPPVEGVCIEMPAG
LIDAGEDVDTAALRELKEETGYKGKIISKSPTVFNDPGFTNTNLCLVTVEVDMSLPENQK
PVTQLEDNEFIECFSVELHKFPDEMVKLDQQGYKLDARVQNVAQGILMAKQYNIQ*

これは後々問題になるので、元のtranscript_idに名前を戻してしまおう（先に--single_best_onlyを指定して一対一対応がつくようにしているので、これが可能になっている）。これは、seqkit の文字列置換コマンド(replace)を用いて以下のようにして行える。

seqkit replace -p '\.p[0-9] ' -r ' ' seub_genes.fa.transdecoder.pep | sed 's/\*//' > seub_genes.pep

なお、元の配列ファイルはアミノ酸配列の末尾にストップコドンの存在を示す*が入っているが、ソフトウェアによってはこれが問題になることがあるので、併せて sed コマンドを用いてこれを除去している。加工後は、以下のようになる。

>DI49_1142 Gene.6155::DI49_1142::g.6155  ORF type:complete len:236 (-) DI49_1142:234-941(-)
MLPLIASRNRRPISLTIRKLFRTMSIVKGKPEEAKIVEARHVKDTSDCKWIGLQKIIYKD
PNGNEREWDSAVRTTRNSGGVDGIGILTILKYKDGKPDEILLQKQFRPPVEGVCIEMPAG
LIDAGEDVDTAALRELKEETGYKGKIISKSPTVFNDPGFTNTNLCLVTVEVDMSLPENQK
PVTQLEDNEFIECFSVELHKFPDEMVKLDQQGYKLDARVQNVAQGILMAKQYNIQ

Step 3: BLASTを用いたホモロジー検索　(講義ではスキップする予定)

前ステップで作成したアミノ酸配列 seub_genes.pep を用いて、S. cerevisiaeゲノムのアミノ酸配列 scer_prot.fa と、BLASTによる総当たりのホモロジー検索を行う(これらの配列ファイルは、上記データディレクトリ上に置いてある)。

まず、scer_prot.fa をもとに、BLAST検索用データベースを作成する。

% makeblastdb -in scer_prot.fa -dbtype prot -parse_seqids -out scer

結果として、scerで始まる複数のファイルが作成される。これを用いてBLAST検索を実行する。

% blastp -query seub_genes.pep -db scer -evalue 0.001 -outfmt "6 std stitle" -max_target_seqs 10 -num_threads 8 > blastout.tab

Step 4: BLAST結果からベストヒット関係の抽出

前ステップで作成したBLAST検索結果から、各クエリ配列につきスコアが最高のヒット（ベストヒット）ひとつだけを抽出する。これは、元の検索結果が、予めクエリ配列ごとにスコア順に並んでいることを前提として、sortコマンドのstable option (-s)とunique option (-u) を用いることで、元の順序を維持しつつクエリ配列ごとに最初の一つのヒットのみを出力することで実現している。

% sort -s -k 1,1 -u blastout.tab > blast_top.tab

オーソログ同定においては、ペアワイズのゲノム比較において、一方向のベストヒットだけでなく、双方向のベストヒットを確認することで、その精度を高めることができる。これを行うため、まず逆方向のベストヒットを抽出する。データベース配列(S. cerevisiae)ごとのスコア順に並べ替え、前コマンドと同様にしてユニークな配列を抽出する。

% sort -s -k 2,2 -k 11,11g -k 12,12nr blastout.tab | sort -s -k 2,2 -u > blast_top_rev.tab

その後、両方向のベストヒットを合わせてソートし、重複した行を出力する(uniq -d)。これにより、双方向のベストヒットが抽出される。

% sort blast_top.tab blast_top_rev.tab | uniq -d > blast_bbh.tab

Step 5: DIAMONDを用いたホモロジー検索（オプショナルだが、講義ではこちらを実施する）

DIAMONDは、BLASTと比べて精度はやや落ちるが、圧倒的に高速であることから、大規模な検索において広く用いられているツールである。使い方はBLASTとよく似ているが、各コマンドはdiamondのサブコマンドとして呼び出すこと、およびオプションの指定がハイフン２つになるところが異なるので注意する。２番目のコマンドは長いので横スクロールして全体を確認すること。なお、diamondではデフォルトでevalueの閾値が0.001なので、そのままでよければ--evalue オプションは省略可能である。

% diamond makedb --in scer_prot.fa --db scer
% diamond blastp --query seub_genes.pep --db scer --max-target-seqs 10 --outfmt 6 qseqid sseqid pident evalue bitscore stitle --threads 4 --out diamondout.tab

（追加課題）上記を実行すると、diamondout.tabという結果ファイルができる。この結果を用いて、前項と同様にベストヒット、および双方向ベストヒットを抽出せよ。BLAST検索の時と比べて出力カラムを減らしているので、カラム位置がずれている点に注意せよ。

Step 6: InterProScan を用いたモチーフ／ドメイン検索（オプショナル。講義ではスキップする予定）

モチーフ／ドメイン検索は、機能推定のためのより詳細な手がかりをうる手段として、ホモロジー検索と併せて実施されることが多い。InterProScanは、一つのコマンドで多岐にわたるデータベースを一度に検索できるツールとして広く用いられている。検索には時間がかかるため、試してみる場合は、クエリ配列のごく一部を抽出した縮小版seub_genes6.pepを用いること。なお、InterProScan はbias5上ではモジュールで管理されているので、パスが通っていない場合はmodule addコマンドでロードする。

% module add interproscan
% interproscan.sh -i seub_genes6.pep -b seub_genes6 -goterms -pa --cpu 4

結果はseub.genes6.iprscanで始まるいくつかのファイルとして出力される。このうち、seub.genes6.iprscan.tsvは上記ディレクトリ上に置いてある。

Step 7: EggNOG mapperを用いたオーソログ検索

EggNOG mapperは、あらかじめ作成したオーソロググループと系統樹の情報を用いてクエリ配列のオーソログを同定し、それに基づいてアノテーションづけするツールである。検索エンジンにはDIAMONDを用いているが、データベースが大きいために検索には時間がかかる。試してみる場合は、縮小版のseub_genes6.pepを用いること。

% emapper.py -i seub_genes6.pep -o seub_genes6 -m diamond --cpu 6

配列全体に対して実行した結果は seub_genes.emapper.annotaionsとして置いてあるので、これを用いてテキスト記述およびGOリストを抽出する。

% cut -f1,8 seub_genes.emapper.annotations | grep -v ^# > seub_genes.emapper.tit
% cut -f1,13 seub_genes.emapper.annotations | grep -v ^# > seub_genes.emapper.go

（発展課題）クラスター計算機を用いた大量検索の高速実行

クエリを分割して並列に実行することで検索速度を上げることができる。実際に配列全体に対して実行してみたい人は、配列を分割してBIASのクラスターシステムを用いて実行してみよう。配列を分割するためのコマンドsplit_seq.plが用意してあるので、これを実行する。-BLOCK_SIZEオプションで、分割の単位となる配列長の総和を指定する。

% split_seq.pl -BLOCK_SIZE=10000 seub_genes.pep

分割された配列ファイルは、query_seub_genesというディレクトリに格納される(この場合、300個程度のファイルに分割される)。併せて、qsub_blast.sh というファイルが作成されるので、末尾にあるコマンド群のうち、blastp をコメントアウトし、emapperのコメントを外す。（同様に、他のコマンドもコメントを外すことにより実行できるが、blast, diamondについては、この上で設定されるDB変数に適切な名前（ここではscer）を入れる必要がある）

#blastp -db $DB -query $INFILE -out $RESULT_OUT_DIR/blastp_$OUTFILE -num_threads $NCPUS -outfmt 6 -evalue 0.001
#diamond blastp --db $DB --query $INFILE --out $RESULT_OUT_DIR/diamond_$OUTFILE --threads $NCPUS 
emapper.py -m diamond --no_annot --no_file_comments -i $INFILE -o $RESULT_OUT_DIR/emapper_$OUTFILE --cpu $NCPUS
#interproscan.sh -goterms -pa -i $INFILE -b $RESULT_OUT_DIR/iprscan_$OUTFILE -cpu $NCPUS

以下のコマンドで実行する。

% qsub qsub_blast.sh

実行状況は、qstat -u ユーザ名で確認できる（-uをつけないと全員分が表示される）。何も表示されなくなると、実行は終わっている。あまりに早く終了する場合は、失敗している可能性が高いので、カレントディレクトリに作成されるblastjob.e#### というファイルに実行したコマンドの標準エラー出力が記録されているので、これを参照して対処すること。

実行が無事に終了すると、分割したクエリ配列ごとの検索結果がoutput_seub_genesディレクトリの下に格納される。この下のファイルをcat コマンドでまとめてカレントディレクトリにコピーする。EggNOG mapperの場合、この段階ではseed ortholog の検索のみが終わっている(オプションに--no_annotを指定したため)ので、この結果をもとに再度emapperを実行して、アノテーションづけを行う。

% cat output_seub_genes/*.seed_orthologs > input.emapper.seed_orthologs
% emapper.py --annotate_hits_table input.emapper.seed_orthologs --no_file_comments -o seub_genes --cpu 10

結果はseub_genes.emapper.annotationsというファイルとして作成される。