OrthoFinder の使い方

### OrthoFinder の使い方
(2016.02.20, 2016.03.09−10, El Capitan)

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

また簡単なだけではなく、[OrthoFinder の論文](http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4531804/) によると、
[OrthoBench](http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3193375/) を使用した評価で、最も良く使用される OrthoMCL をはじめ、
様々な予測手法よりも高性能となったそうです。

---
#### <i class="fa fa-cube" style="font-size:1em;"></i> 準備

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
```

---
#### <i class="fa fa-cube" style="font-size:1em;"></i> OrthoFinder のダウンロード

[OrthoFinder](https://github.com/davidemms/OrthoFinder/releases/download/v0.4/OrthoFinder-0.4.tar.gz) の最新版を ~/tools にダウンロードします。

tar コマンドで、解凍します。

```
tar xvzf OrthoFinder-0.4.tar.gz 
```

---
#### <i class="fa fa-cube" style="font-size:1em;"></i> テストデータセット

解凍をすると、すぐ直下のフォルダにテストデータセットが用意されています。
`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
```

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

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

---
#### <i class="fa fa-cube" style="font-size:1em;"></i> テスト解析

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

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

```
python2 orthofinder.py -f ExampleDataset -t 4
```

`ExampleDataset` フォルダの直下に `Results_Date` という結果が作成されるはずです。
.csv ファイルなどにオルソグループの情報が記述されていました。

---
#### <i class="fa fa-cube" style="font-size:1em;"></i> Phytozome からダウンロードしたプロテオームデータファイルに拡張子を追加（2016.03.09）

[Phytozome v11.0](https://phytozome.jgi.doe.gov/pz/portal.html) からダウンロードしたマルチファスタファイルは、
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
```

---
#### <i class="fa fa-cube" style="font-size:1em;"></i> 実際に使用してみる (2016.03.09)

[Phytozome v11.0](https://phytozome.jgi.doe.gov/pz/portal.html) から、例えば、適当に６つの植物種のプロテオームデータをダウンロードして、
`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` のようになっているのが結果が格納されているフォルダになります。

---
#### <i class="fa fa-cube" style="font-size:1em;"></i> 解析結果のファイルフォーマットを調べる（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 が記述されている
- ２列目からは、そのオルソグループに所属する各生物種の遺伝子名がタブ区切りで記述されている__
となっていました。

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

---
#### <i class="fa fa-cube" style="font-size:1em;"></i> シロイヌナズナのオルソログからアノテーションを付ける（2016.03.10）

最も機能解明が進んでいるシロイヌナズナ遺伝子から、機能未知が多いトマト遺伝子にアノテーションします。
シロイヌナズナ遺伝子のアノテーションを得るため、[EnsemblPlants](http://plants.ensembl.org/index.html) の 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` という名前で保存しました。
以下のプログラムでトマト遺伝子にシロイヌナズナのアノテーションを付与しました。

```python
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` に結果が出力されるはずです。