ex801 - nibb-gitc/gitc2025mar-rnaseq GitHub Wiki

ex801 Functional annotation using similarity-based methods

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

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

データ

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

file	contents	remarks
seub_genome.fa	S. eubayanus ゲノム配列	入力データ
stringtie_merged.gtf	stringtieによるアセンブル結果	入力データ
topTags.transcript.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 ex801
% cp yeast/* ex801
% cd ex801

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から少し変わっている（配列名の後に.p1がついている。1はこのトランスクリプト配列からとれたCDSの番号で、2以上になることもある）。

>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_transcripts.fa.transdecoder.pep | sed 's/\*//' > seub_genes.pep

-p で取り除く部分の正規表現を指定し、それを -r で指定した空白に置き換えている。なお、元の配列ファイルはアミノ酸配列の末尾にストップコドンの存在を示す*が入っているが、ソフトウェアによってはこれが問題になることがあるので、併せて sed コマンドを用いてこれを除去している。加工後は、以下のようになる。

>DI49_1142 GENE.DI49_1142~~DI49_1142  ORF type:complete len:236 (-),score=50.10 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: DIAMONDを用いたホモロジー検索（講義ではこちらを実施する）

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

なお、DIAMONDの最新版では、BLAST用のインデックスをそのまま使った検索も可能になり、BLAST用インデックスが存在する場合はそちらを優先するようである。そこで、以下ではStep 3で作成したBLAST用データベースと名前がかち合わないように、 scer2という名前でDIAMOND用データベースを作成している。

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

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

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

% sort -s -k 1,1 -u diamondout.tab > diamond_top.tab

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

% sort -k 2,2 -k 4,4g -k 5,5nr diamondout.tab | sort -s -k 2,2 -u > diamond_top_rev.tab

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

% sort diamond_top.tab diamond_top_rev.tab | uniq -d > diamond_bbh.tab

（追加課題）Step 3で作成したBLASTによる検索結果 blastout.tab を用いて、同様にベストヒット、および双方向ベストヒットを抽出せよ。DIAMOND検索の時と比べて出力カラムが異なっているので、カラム位置がずれている点に注意せよ。

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

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

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

結果はseub_genes6で始まるいくつかのファイルとして出力される。このうち、タブ区切りでアノテーションを記載したseub_genes6.tsvは、上記ディレクトリ上に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

(註）出力のカラム構成はバージョンによって異なるので、何番目のカラムを抜き出すかはヘッダを見て確認すること。最新のver 2.1.12では、GOのカラムは10番目になっている。

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

（以下は、分子科学研究所計算科学研究センターのシステム(RCCS)を使って行う実習手順である。参考までに掲載するが、今回の実習では行わない）

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

% perl ~/gitc/data/IU/bin/split_seq.pl -BLOCK_SIZE=100000 seub_genes.pep

分割された配列ファイルは、query_seub_genesというディレクトリに格納される(この場合、30個程度のファイルに分割される)。併せて、これらの配列をクエリとして、PBSというジョブ管理システムを用いてクラスター上で並列に実行するための、qsub_blast.sh というファイルが作成される。このファイルの末尾にあるコマンド群のうち、blastp をコメントアウトし（先頭に#をつける）、emapperのコメントを外す（先頭の#を削除する）。（同様に、他のコマンドもコメントを外すことにより実行できるが、blast, diamondについては、検索対象データベース名として、この上で設定されるDB変数に適切な名前（ここではscer）を入れる必要がある）

また、計算科学センターでは、ジョブを走らせる際に、使用する計算リソースを決められたオプションに従って指定することが規則として決められている。これは、先頭の #PBS で始まる行に -l オプションとして指定する。ここでは、ジョブ当たり10コアを使い（ncpus, ompthreads）、最大1時間の実行時間（walltime）を指定している。また、eggNOG-mapper を使うために、モジュールのロードを行う必要がある。

これらを反映させるために、qsub_blast.sh を以下の様に修正する（#<<<< がついた行を修正する。コピペする場合は、#PBS行末尾のコメントは削除すること）。

#!/bin/sh
#PBS -J 1-30
#PBS -N blastjob
#PBS -S /bin/sh
#PBS -V
#PBS -l select=1:ncpus=10:mpiprocs=1:ompthreads=10   #<<<<　書き換え(リソース指定、CPUコア数とメモリ）
#PBS -l walltime=1:00:00                           #<<<<　書き換え(リソース指定、最大実行時間）

source /apl/bio/etc/bio.sh          #<<<<　追加(生物学関連アプリを使うための環境設定）
module load eggNOG-mapper           #<<<<　追加(モジュールのロード）

cd $PBS_O_WORKDIR
DB=scer                             #<<<<　データベースを検索する場合はここを変更
SEQDIR=query_seub_genes
RESULT_OUT_DIR=output_seub_genes
QUERY_OUT_DIR=query_seub_genes

INFILE=$QUERY_OUT_DIR/seub_genes.$PBS_ARRAY_INDEX
OUTFILE=seub_genes.$PBS_ARRAY_INDEX

if [ ! -d $RESULT_OUT_DIR ]; then
    mkdir $RESULT_OUT_DIR
fi
#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

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

% jsub qsub_blast.sh

実行状況は、jobinfoで確認できる。ジョブの情報が表示されなくなると、実行は終わっている。あまりに早く終了する場合は、失敗している可能性が高いので、カレントディレクトリに作成される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 -m no_search --override

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