bioinformatics

OrthoFinder の使い方

More than 1 year has passed since last update.

OrthoFinder の使い方

(2016.02.20, 2016.03.09−10, El Capitan)

OrthoFinder というオルソグループ(オルソロググループ)予測手法の使い方を調べました。
以前使用したことがある OrthoMCL に比べて非常に簡単に使用することができました。

また簡単なだけではなく、OrthoFinder の論文 によると、
OrthoBench を使用した評価で、最も良く使用される OrthoMCL をはじめ、
様々な予測手法よりも高性能となったそうです。


準備

python2.7, scipy, BLAST+, MCL が必要です。
どのような方法で用意しても良いと思いますが、1例を紹介します。

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

1生物につき、1つのプロテオームデータ(すべてのタンパク質配列)が、
マルチファスタファイル形式で入ったディレクトリが用意されています。
4生物のプロテオームデータがあることがわかります。

自分にあった解析を行うためには、適当なディレクトリを一つ作り、同様にして、
そこに自分が解析したい生物種のプロテオームデータを配置するだけで良いそうです。


テスト解析

python orthofinder.py -f の後に、ファスタファイルの入ったディレクトリを指定するだけで、
後は、OrthoFinder が自動でオルソグループ作成を行ってくれました。

以下は、既に用意された ExampleDataset を使用して、4並列で解析するコマンドです。

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 から、例えば、適当に6つの植物種のプロテオームデータをダウンロードして、
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)

ここでは、解析をスタートしたのが3月9日なので、Results_Mar09 フォルダに結果が格納されています。
ls Results_Mar09 とすると以下のように表示されるもののうち、
OrthologousGroups.csv にオルソグループの情報が出力されていました。

OrthologousGroups.csv     OrthologousGroups.txt     OrthologousGroups_UnassignedGenes.csv     WorkingDirectory

OrthologousGroups.csv の1行目を 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

上記から、以下のような内容であることがわかります。
- 最初の1行目は、各々の生物種のプロテオームデータファイル名がタブ区切りで記述されている

2行目からは、1行につき1つのオルソグループが出力されていました。
例えば、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

オルソグループが記述された2行目以降のファイルフォーマットは、
- 1列目には、オルソグループの ID が記述されている
- 2列目からは、そのオルソグループに所属する各生物種の遺伝子名がタブ区切りで記述されている

となっていました。

オルソグループには、1つの生物種由来の複数の遺伝子が含まれる場合があります。このような場合、
- 1生物種から複数遺伝子が同一オルソグループに所属する場合、それらの遺伝子はカンマ区切りで記述
となっていました。


シロイヌナズナのオルソログからアノテーションを付ける(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 に結果が出力されるはずです。