最近Orthofinder2開始陸續(xù)更新，文章已經(jīng)投發(fā)到biorkix上，在網(wǎng)上搜了一圈，關(guān)于Orthofinder的使用的中文教程幾乎是空白的，這周就借此機(jī)會(huì)和大家一起來(lái)學(xué)習(xí)一下這款軟件。

Orthofinder介紹

OrthoFinder是比較基因組學(xué)中的實(shí)用的，運(yùn)行快速，準(zhǔn)確的全面的工具。它的主要功能是，找到了正交群和直系同源物，推斷出所有正交群的根基因樹，并識(shí)別那些基因樹中的所有基因重復(fù)事件。它還為所分析的物種推斷出有根的物種樹，并將基因重復(fù)事件從基因樹比對(duì)到物種樹的分支中。

另外，OrthoFinder還為比較基因組分析提供全面的統(tǒng)計(jì)數(shù)據(jù)。相比OrthoMCL， OrthoFinder使用簡(jiǎn)單安裝都很簡(jiǎn)單，安裝只需要用過(guò)conda就行，運(yùn)行只需一組FASTA格式的蛋白質(zhì)序列文件（每個(gè)物種一個(gè)）。

基礎(chǔ)知識(shí)介紹

為了讓大家更好的理解使用Orthofinder，在正式介紹其使用的流程前，先來(lái)回顧一下與之相關(guān)的生物學(xué)知識(shí)。

直系同源物和旁系同源物

Orthologs(直系同源物)是在兩個(gè)物種的最后共同祖先（LCA）中來(lái)自單個(gè)基因的一對(duì)基因。直系同源物是同源性基因，是物種形成事件的結(jié)果。Paralogs（旁系同源物)是同源基因，是重復(fù)事件的結(jié)果。下圖就可以看到，不同物種間的alpha-chain gene互為Orthologs(直系同源物)。這時(shí)候可以引用一個(gè)新名詞orthogroup （正交群）就用來(lái)形容自一組物種的LCA中的單個(gè)基因的基因組（在圖中就是alpha chain gene）。然后同一物種間alpha 和beta chain gene互為Paralogs(旁系同源物)。最后所有這些關(guān)系都可以由OrthoFinder來(lái)識(shí)別。

為什么要研究正交群？

1. 正交群中的所有基因都來(lái)自單個(gè)祖先基因。因此，正交群中的所有基因都有類似的序列和功能。由于基因重復(fù)和丟失在進(jìn)化中經(jīng)常發(fā)生，一對(duì)一的直系同源物很少見(jiàn)，通過(guò)分析orhtogroup所有直系同源的情況（一對(duì)一，多對(duì)一，多對(duì)多），我們可以分析數(shù)據(jù)的所有情況。
2. 直系同源物是由系統(tǒng)發(fā)育定義的。它不能通過(guò)氨基酸含量，密碼子偏好，GC含量或其他序列相似性測(cè)量定義。在沒(méi)有系統(tǒng)發(fā)育的情況下使用這些分?jǐn)?shù)來(lái)定義，直向同源物的方法只能提供猜測(cè)。從序列相似性提供粗略的類比猜測(cè)直系同源，類似于從氣味中猜測(cè)顏色。因此，真正識(shí)別直向同源物的唯一方法是通過(guò)分析系統(tǒng)發(fā)育樹。確保直系分配正確的唯一方法是對(duì)所考慮的物種的最后共同祖先的單個(gè)基因下降的所有基因進(jìn)行系統(tǒng)發(fā)育重建。這組基因是正交群。因此，定義直系同源的唯一方法是分析正交群。

安裝

廢話說(shuō)了一大篇，是時(shí)候開始安裝軟件了。在這里推薦大家使用conda來(lái)安裝，簡(jiǎn)單明了，不用擔(dān)心其他Dependencies的安裝

conda install orthofinder

接著準(zhǔn)備好測(cè)試數(shù)據(jù)，為了展示的方便，我將原始測(cè)試數(shù)據(jù)重新命名：

##重命名
cd /work1/ExampleData
mv Mycoplasma_agalactiae.faa Maga.faa
mv Mycoplasma_gallisepticum.faa Mgal.faa
mv Mycoplasma_genitalium.faa Mgen.faa
mv Mycoplasma_hyopneumoniae.faa Mhyo.faa

##創(chuàng)建新的數(shù)據(jù)目錄
mkdir Mycoplasma/
cp *faa Mycoplasma/

正式運(yùn)行Orthofinder，相當(dāng)簡(jiǎn)單的操作,-f輸入目錄，里面包含你需要運(yùn)行的蛋白質(zhì)fasta文件， -t 所用到的CPU數(shù)目?；镜挠梅ň腿缦?，更多的可以去manual中查看。

orthofinder -f Mycoplasma/ -t 2

生成的結(jié)果會(huì)存儲(chǔ)于Results_XXX文件中，現(xiàn)在簡(jiǎn)單看看里面有啥。
正交群以兩種不同的格式返回，一種是制表符分隔的表格格式（* .csv），另一種格式與orthoMCL（* .txt）的輸出格式相同。還應(yīng)該有一個(gè)名為WorkingDirectory的目錄，其中包含運(yùn)算過(guò)程的中間文件，例如blast結(jié)果。

cd Mycoplasma/Results_Nov10/
ll -thrl
-rw-r----- 1 hhu hhu  12K Nov 10 16:50 Orthogroups.GeneCount.csv
-rw-r----- 1 hhu hhu  59K Nov 10 16:50 Orthogroups.csv
-rw-r----- 1 hhu hhu  89K Nov 10 16:50 Orthogroups.txt
-rw-r----- 1 hhu hhu  110 Nov 10 16:50 Orthogroups_SpeciesOverlaps.csv
-rw-r----- 1 hhu hhu  33K Nov 10 16:50 Orthogroups_UnassignedGenes.csv
drwxr-x--- 6 hhu hhu 4.0K Nov 10 16:50 Orthologues_Nov10
-rw-r----- 1 hhu hhu 2.5K Nov 10 16:50 SingleCopyOrthogroups.txt
-rw-r----- 1 hhu hhu 1.5K Nov 10 16:50 Statistics_Overall.csv
-rw-r----- 1 hhu hhu 2.5K Nov 10 16:50 Statistics_PerSpecies.csv
drwxr-x--- 2 hhu hhu 4.0K Nov 10 16:50 WorkingDirectory

使用R簡(jiǎn)單分析Orthofinder的結(jié)果

安裝加載以下R包

install.packages("micropan")
install.packages("phangorn")

library(micropan) 
library(phangorn)

讀取Orthofinder中正交群的文件。結(jié)果文件將正交群存儲(chǔ)為制表符分隔表，其中每一行是正交群，其名為OG，后跟7位數(shù)字（例如OG0000123）。表中的每列對(duì)應(yīng)一個(gè)物種，并且與fasta文件具有相同的名稱。表中的單元格包含序列ID。由于正交組可以包含來(lái)自單個(gè)物種的幾個(gè)基因（即，旁系同源物），因此表中的單個(gè)cell可以包含由逗號(hào)分隔的幾個(gè)序列ID。下面是讀取結(jié)果文件的例子：

# Define orthofinder results directory (may vary)
resultsDir <- "orthofinder/Mycoplasma/Results_Oct22"

grpsTbl <- read.table( file.path(resultsDir,"OrthologousGroups.csv"),
                       header=T, sep = "\t", row.names = 1, stringsAsFactors=F)

# Split the comma-separated paralogs
grps <- strsplit( as.matrix(grpsTbl),", ")
# grps is now a list of character vectors..
# Add dimensions to grps:
dim(grps) <- dim(grpsTbl)
# Add row and column names:
dimnames(grps) <- dimnames(grpsTbl)

獲取特定正交群的信息：

grps["OG0000018", ]```

可以清晰的看到，在對(duì)應(yīng)的fa文件中，在該特定的正交群OG000018下所包含的fasta序列，其結(jié)果如下：

## $Maga.faa
## character(0)
## 
## $Mgal.faa
## [1] "Mgal_AAP56913.2" "Mgal_AAP56907.2" "Mgal_AAP56908.2"
## 
## $Mgen.faa
## [1] "Mgen_AAC71531.2" "Mgen_AAC71529.1" "Mgen_AAC71563.2"
## 
## $Mhyo.faa
## character(0)

接著讓我們來(lái)看開正交群的數(shù)目，為了獲得特定正交群的序列，需要在每個(gè)正交群中構(gòu)建每個(gè)物種的基因數(shù)量的矩陣：

grpSize <- apply(grps,2,sapply,length)

我們可以在正交群OG0000018中檢查每個(gè)物種中的旁系同源物的數(shù)量：

grpSize["OG0000018", ]

## Maga.faa Mgal.faa Mgen.faa Mhyo.faa 
##        0        3        3        0

我們還可以繪制每個(gè)正交群中存在的基因數(shù)量：

plot(table(apply(grpSize>0,1,sum)),ylab="no. orthogroups",xlab="no. genomes represented")

另外一寫有趣的事我們可以查看，每個(gè)物種中只有一個(gè)直系同源的正交群（1：1：1：1正交群）：

sum(apply(grpSize==1,1,all)) # number of 1:1:1:1 groups
## [1] 254

正好只存在于兩個(gè)物種（1：1：0：0正交群）的正交群：

sum(apply(grpSize==1,1,sum)==2 & apply(grpSize > 1,1,sum)==0)
## [1] 155

讓我們找出這個(gè)數(shù)字背后所對(duì)應(yīng)的物種：

# get the 1:1:0:0 groups
is1100 <- apply(grpSize==1,1,sum)==2 & apply(grpSize > 1,1,sum)==0

spcs <- sub("\\.faa","",colnames(grps)) # get species names

# Count the combination of species-pairs in the 1:1:0:0 groups
table(apply(grpSize[is1100, ], 1, function(grpSizeRow){
  paste(spcs[ grpSizeRow==1 ], collapse="&")
}))
## 
## Maga&Mgal Maga&Mgen Maga&Mhyo Mgal&Mgen Mgal&Mhyo Mgen&Mhyo 
##        15         2        57        70         8         3

基本的用法和分析就講解到這里，當(dāng)然這只是最基本的東西。其結(jié)果涉及到的下游分析還有很多很多，例如構(gòu)建直系同源基因樹等等的知識(shí)，大家可以自行進(jìn)行探討。