OrthoFinder の使い方
(2016.02.20, 2016.03.09−10, El Capitan)
OrthoFinder というオルソグループ（オルソロググループ）予測手法の使い方を調べました。
以前使用したことがある OrthoMCL に比べて非常に簡単に使用することができました。
また簡単なだけではなく、OrthoFinder の論文 によると、
OrthoBench を使用した評価で、最も良く使用される OrthoMCL をはじめ、
様々な予測手法よりも高性能となったそうです。
準備
python2.7, scipy, BLAST+, MCL が必要です。
どのような方法で用意しても良いと思いますが、１例を紹介します。
python and scipy
まず、home brew を使用できるようにします。
/usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
python2 をインストールし、pip で scipy をインストールします。
brew install python
pip intall scipy
BLAST+
BLAST+ をダウンロードします。
Mac用の最新版
(ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.3.0+-universal-macosx.tar.gz)
を ~/tools にダウンロードしました。
tar コマンドで、解凍します。
cd ~/tools
tar xvzf ncbi-blast-2.3.0+-universal-macosx.tar.gz
.zshenv に以下を追記し、source .zshenv とします。
# for BLAST+
export PATH="$HOME/tools/ncbi-blast-2.3.0+/bin:$PATH"
MCL
MCL クラスタリングを使用できるようにします。
brew tap homebrew/science
brew install mcl
OrthoFinder のダウンロード
OrthoFinder の最新版を ~/tools にダウンロードします。
tar コマンドで、解凍します。
tar xvzf OrthoFinder-0.4.tar.gz
テストデータセット
解凍をすると、すぐ直下のフォルダにテストデータセットが用意されています。
cd ~/tools/OrthoFinder-0.4 として、
ls ExampleDataset とすると、どのようなデータが入っているか調べることができます。
Mycoplasma_agalactiae_5632_FP671138.faa Mycoplasma_genitalium_uid97_L43967.faa
Mycoplasma_gallisepticum_uid409_AE015450.faa Mycoplasma_hyopneumoniae_AE017243.faa
１生物につき、１つのプロテオームデータ（すべてのタンパク質配列）が、
マルチファスタファイル形式で入ったディレクトリが用意されています。
４生物のプロテオームデータがあることがわかります。
自分にあった解析を行うためには、適当なディレクトリを一つ作り、同様にして、
そこに自分が解析したい生物種のプロテオームデータを配置するだけで良いそうです。
テスト解析
python orthofinder.py -f の後に、ファスタファイルの入ったディレクトリを指定するだけで、
後は、OrthoFinder が自動でオルソグループ作成を行ってくれました。
以下は、既に用意された ExampleDataset を使用して、４並列で解析するコマンドです。
python2 orthofinder.py -f ExampleDataset -t 4
ExampleDataset フォルダの直下に Results_Date という結果が作成されるはずです。
.csv ファイルなどにオルソグループの情報が記述されていました。
Phytozome からダウンロードしたプロテオームデータファイルに拡張子を追加（2016.03.09）
Phytozome v11.0 からダウンロードしたマルチファスタファイルは、
OrthoFinder で認識される以下の拡張子を持たないため、自分で拡張子を付ける必要がありました。
fastaExtensions = {"fa", "faa", "fasta", "fas"}
ダウンロードリンクのファイル名は、.fa.gz となっているのですが、
ダウンロードすると何故かファイル名が、 .gz だけになっていました。。
例えば、Athaliana_167_TAIR10.protein_primaryTranscriptOnly.fa.gz をダウンロードして、
解凍すると Athaliana_167_TAIR10.protein_primaryTranscriptOnly という名前のファイルができました。
そこで、Athaliana_167_TAIR10.protein_primaryTranscriptOnly のあるフォルダに移動して、
以下のコマンドで .fa という拡張子を追加しました。
brew install rename
rename 's/Only/Only.fa/' *Only
実際に使用してみる (2016.03.09)
Phytozome v11.0 から、例えば、適当に６つの植物種のプロテオームデータをダウンロードして、
Dataset という名前のフォルダに入れ、ls Dataset と入力すると以下のような感じに出力されるようにします。
Athaliana_167_TAIR10.protein_primaryTranscriptOnly.fa Osativa_323_v7.0.protein_primaryTranscriptOnly.fa
Zmays_284_5b+.protein_primaryTranscriptOnly.fa BrapaFPsc_277_v1.3.protein_primaryTranscriptOnly.fa
Gmax_275_Wm82.a2.v1.protein_primaryTranscriptOnly.fa Slycopersicum_225_iTAGv2.3.protein_primaryTranscriptOnly.fa
ここでは、シロイヌナズナ、ライス、トウモロコシ、フィールドマスタード、大豆、トマト のプロテオームデータ
を使用しています。
テスト解析の時と同様に、python2 orthofinder.py -f Dataset -t 4 と入力すれば解析がスタートします。
自分のパソコンでは、だいたい 24h くらいで解析が終了しました。
ls Dataset とすると、結果のフォルダが作成されていることを確認することができます。
Athaliana_167_TAIR10.protein_primaryTranscriptOnly.fa Osativa_323_v7.0.protein_primaryTranscriptOnly.fa
Results_Mar09 Zmays_284_5b+.protein_primaryTranscriptOnly.fa
BrapaFPsc_277_v1.3.protein_primaryTranscriptOnly.fa Gmax_275_Wm82.a2.v1.protein_primaryTranscriptOnly.fa
Slycopersicum_225_iTAGv2.3.protein_primaryTranscriptOnly.fa
Results_Mar09 のようになっているのが結果が格納されているフォルダになります。
解析結果のファイルフォーマットを調べる（2016.03.10）
ここでは、解析をスタートしたのが３月９日なので、Results_Mar09 フォルダに結果が格納されています。
ls Results_Mar09 とすると以下のように表示されるもののうち、
OrthologousGroups.csv にオルソグループの情報が出力されていました。
OrthologousGroups.csv OrthologousGroups.txt OrthologousGroups_UnassignedGenes.csv WorkingDirectory
OrthologousGroups.csv の１行目を head -1 Results_Mar09/OrthologousGroups.csv として調べます。
Athaliana_167_TAIR10.protein_primaryTranscriptOnly.fa BrapaFPsc_277_v1.3.protein_primaryTranscriptOnly.fa Gmax_275_Wm82.a2.v1.protein_primaryTranscriptOnly.fa Osativa_323_v7.0.protein_primaryTranscriptOnly.fa Slycopersicum_225_iTAGv2.3.protein_primaryTranscriptOnly.fa Zmays_284_5b+.protein_primaryTranscriptOnly.fa
上記から、以下のような内容であることがわかります。
- 最初の１行目は、各々の生物種のプロテオームデータファイル名がタブ区切りで記述されている
２行目からは、１行につき１つのオルソグループが出力されていました。
例えば、OG0007911 という ID のオルソグループが記述された行は、以下のような内容になっていました。
OG0007911 AT4G31780.2 Brara.K01636.1.p, Brara.A00599.1.p Glyma.14G104100.1.p, Glyma.17G221900.1.p LOC_Os09g25580.1 Solyc08g006640.2.1 GRMZM2G142873_P01
オルソグループが記述された２行目以降のファイルフォーマットは、
- １列目には、オルソグループの ID が記述されている
- ２列目からは、そのオルソグループに所属する各生物種の遺伝子名がタブ区切りで記述されている
となっていました。
オルソグループには、１つの生物種由来の複数の遺伝子が含まれる場合があります。このような場合、
- １生物種から複数遺伝子が同一オルソグループに所属する場合、それらの遺伝子はカンマ区切りで記述
となっていました。
シロイヌナズナのオルソログからアノテーションを付ける（2016.03.10）
最も機能解明が進んでいるシロイヌナズナ遺伝子から、機能未知が多いトマト遺伝子にアノテーションします。
シロイヌナズナ遺伝子のアノテーションを得るため、EnsemblPlants の Biomart を利用して、
以下のような情報をダウンロードします。
Gene stable ID Gene name Source of gene name Gene description
AT3G18710 PUB29 UniProtKB Gene Name U-box domain-containing protein 29 [Source:UniProtKB/Swiss-Prot;Acc:Q9LSA6]
AT4G25880 APUM6 UniProtKB Gene Name Pumilio homolog 6, chloroplastic [Source:UniProtKB/Swiss-Prot;Acc:Q9C5E7]
AT1G71695 PER12 UniProtKB Gene Name Peroxidase 12 [Source:UniProtKB/Swiss-Prot;Acc:Q96520]
AT1G32045 transposable element gene [Source:TAIR;Acc:AT1G32045]
AT5G41480 EMB9 TAIR Gene Name Folylpolyglutamate synthetase family protein [Source:TAIR;Acc:AT5G41480]
...
...
...
ここでは、ダウンロードしたファイルを ath_agi_name_source_desc.tsv という名前で保存しました。
以下のプログラムでトマト遺伝子にシロイヌナズナのアノテーションを付与しました。
def get_ath_annot(file_in):
annot = {}
with open(file_in, 'rb') as inf:
for line in inf:
agi, name, source, desc = line.rstrip('\n').split('\t')
if name:
annot[agi] = (name + ' (' + agi + ') ' + '[Source:' + source + ']', desc)
else:
annot[agi] = (agi, desc)
return annot
def annot_with_ath_annot(annot, ath_col_idx, col_idx, file_in, file_out):
with open(file_in, 'rb') as inf, open(file_out, "wb") as outf:
inf.readline()
for line in inf:
agis = line.split('\t')[ath_col_idx].split(', ')
solyc_ids = line.split('\t')[col_idx].split(', ')
if solyc_ids[0] and agis[0]:
for solyc_id in solyc_ids:
for agi in agis:
outf.write(solyc_id[0:16] + '\t' + annot[agi[0:9]][0] + '\t' + annot[agi[0:9]][1] + '\n')
if __name__ == "__main__":
annot = get_ath_annot('ath_agi_name_source_desc.tsv')
annot_with_ath_annot(annot, 1, 5, 'OrthologousGroups.csv', 'tomato_gene_with_ath_annot.tsv')
annot_with_ath_annot.py という名前で保存して、python annot_with_ath_annot.py とすると、
tomato_gene_with_ath_annot.tsv に結果が出力されるはずです。