新的系列系列教程已上傳，大家可以去看新的教程

1.orthofinder介紹

OrthoFinder是一種快速铲汪、準確和全面的比較基因組學分析工具勤婚。它可以找到直系和正群，為所有的正群推斷基因樹口锭，并為所分析的物種推斷一個有根的物種樹类嗤。OrthoFinder還為比較基因組分析提供全面的統(tǒng)計數(shù)據(jù)祈餐。OrthoFinder使用簡單瀑踢，只需運行一組FASTA格式的蛋白質(zhì)序列文件（每個物種一個）扳还。

2.基礎(chǔ)知識介紹

Orthologue（直系同源基因）指的是來自兩個物種的基因。Orthologue是由兩個物種的最后共同祖先(LCA)上的單個基因進化而來的成對基因(圖1A和B)橱夭。正群是同源概念在物種群中的自然延伸氨距。一個Orthogroup（正交群）是由一個物種的LCA中的單個基因進化而來的一組基因(圖1A)。當觀察基因樹時棘劣，一個鄰位群體中基因的第一次分化是一個物種形成事件俏让，對同源基因來說也是如此。

作為基因復(fù)制事件的結(jié)果，當觀察直系同源基因和正交群時舆驶，可能會有來自同一物種的多個基因橱健。在這個例子中(圖1A和B)而钞，人類和老鼠的HuA基因是雞中ChA1和ChA2的同源基因沙廉。再看一下正交群，我們發(fā)現(xiàn)有兩個雞的基因(圖1A)臼节，但是只有一個來自老鼠和人類的基因撬陵。一些作者將ChA1和ChA2基因作為HuA的共同源基因，以強調(diào)存在多個同源基因的事實网缝。由于基因重復(fù)和丟失在進化中經(jīng)常發(fā)生巨税，一對一的直系同源物很少見，通過分析正交群所有直系同源的情況（一對一粉臊，多對一草添，多對多），我們可以分析數(shù)據(jù)的所有情況扼仲。

paralogues (旁系同源基因)是指在基因復(fù)制事件中從單個基因中分離出來的成對基因远寸，雞的兩個基因ChA1和ChA2是旁系同源基因(圖1C)。來自不同物種的兩個基因如果在基因重復(fù)事件中彼此分離屠凶，也可能是同源的驰后。由于基因樹中所有的分支事件要么是物種形成事件(產(chǎn)生直系同源基因)，要么是重復(fù)事件(產(chǎn)生旁系同源基因)矗愧，因此同一正交群中任何不是直系同源基因的基因必然是旁系同源基因灶芝。

圖1

直系同源物是同源性基因，是物種形成事件的結(jié)果唉韭。Paralogs（旁系同源物)是同源基因夜涕，是重復(fù)事件的結(jié)果∈舴撸可以看到(圖2)钠乏，不同物種間的α-chain gene互為Orthologs(直系同源物)。正交群用來形容自一組物種的LCA中的單個基因的基因組（α-chain gene）春塌。然后同一物種間α 和β chain gene互為Paralogs(旁系同源物)晓避。最后所有這些關(guān)系都可以由OrthoFinder來識別。

圖2

3.安裝

在這里推薦大家使用conda來安裝只壳，簡單明了俏拱，不用擔心其他Dependencies的安裝

$conda install orthofinder

4.運行orthofinder

運行Orthofinder，相當簡單的操作, "-f"輸入目錄吼句，里面包含你需要運行的蛋白質(zhì)fasta文件蜡豹，“-t”所用到的CPU數(shù)目，“-S”選擇的比對模式（默認diamod旬痹，可選blast）。

我使用的物種有九個：大豆（G.max）驹愚、菠蘿（A.comosus）、柑橘（C.sinensis）劣纲、無油樟（A.trichopoda）逢捺、擬南芥（A.thaliana）、黃瓜（C.sativus）癞季、水稻（O.sativa）劫瞳、小立碗蘚（P.patens）、江南卷柏（S.moellendorffii）

將所有物種的蛋白文件放到Angiospermae文件夾下

$nohup orthofinder -t 16 -f Angiospermae/ -S blast &

生成的結(jié)果會存儲于Orthofinder/Results_XXX文件中绷柒，現(xiàn)在簡單看看里面有啥志于。

Citation.txt                     Orthologues
Comparative_Genomics_Statistics  Phylogenetically_Misplaced_Genes
Gene_Duplication_Events          Putative_Xenologs
Gene_Trees                       Resolved_Gene_Trees
Log.txt                          Single_Copy_Orthologue_Sequences
Orthogroups                      Species_Tree
Orthogroup_Sequences             WorkingDirectory

我們主要使用Orthoogroups查看正交群的基因和使用 Single_Copy_Orthologue_Sequences里的單拷貝基因構(gòu)建系統(tǒng)發(fā)育樹

WorkingDirectory其中包含運算過程的中間文件，例如blast結(jié)果废睦，如果我們想去掉某一物種，在SpeciesIDs.txt中將該物種注釋掉

$vi SpeciesIDs.txt
#0: Acomosus.pro.fa
1: Athaliana.pro.fa
2: Atrichopoda.pro.fa
3: Csativus.pro.fa
4: Csinensis.pro.fa
5: Gmax.pro.fa
6: Osativa.pro.fa
7: Ppatens.pro.fa
8: Smoellendorffii.pro.fa
~                          

$nohup orthofinder -b WorkingDirectory

如果我們想增加額外的物種進行分析，可以按照如下方式運行

$nohup orthofinder -b WorkingDirectory -f new_fasta_directory

5.利用單拷貝基因構(gòu)建系統(tǒng)發(fā)育樹

Orthofinder的 Single_Copy_Orthologue_Sequences下存放著單拷貝同源基因的序列钥组，我們可以利用這些序列構(gòu)建系統(tǒng)發(fā)育樹

5.1.多序列比對

多序列比對推薦使用muscle

$wget http://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz
$tar xzvf muscle3.8.31_i86linux64.tar.gz
$chmod +x muscle
將muscle添加至環(huán)境變量

單拷貝的序列較多，所以使用shell批量運行

$vi bash.sh

for i in *.fa
do muscle -in $i -out $i.1
done

$sh bash.sh

5.2.提取保守序列

*.1文件是比對好的序列文件今瀑，接下來使用Gblocks提取保守序列

$wget http://molevol.cmima.csic.es/castresana/Gblocks/Gblocks_Linux64_0.91b.tar.Z
$tar Zxf Gblocks_Linux64_0.91b.tar.Z
添加至環(huán)境變量
$vi bash2.sh
for i in *.1
do Gblocks $i -b4=5 -b5=h -t=p -e=.2
done
$sh bash2.sh

-t= default:p設(shè)置序列的類型程梦，可選的值是 p,d,c 分別代表 protein， DNA橘荠， Codons 屿附。
-b1= default:（ 序列條數(shù)的 50% + 1 ）,設(shè)定保守性位點必須有 >= 該值的序列數(shù)。該參數(shù)后接一個 integer 數(shù)哥童，默認下比序列條數(shù)的 50% 大 1.
-b2= default: 序列條數(shù)的 85%,確定保守位點的側(cè)翼位點時贮懈，其位點必須有 >= 該值的序列數(shù)匀泊。該值必須要比 -b1 的值要大。
-b3= default: 8,最大連續(xù)非保守位點的長度朵你。
-b4= default: 10,保守位點區(qū)塊的最小長度各聘。該值必須 >=2 。
-b5= default: n, 設(shè)置允許含有 Gap 位點〈舐觯可選的值有 n,h,a 分別代表 None, With Half, All 搞监。 當為 h 時，表示
-e= default: .2, 設(shè)置輸出結(jié)果的后綴镰矿。

5.3.序列合并

比對好之后我們需要將所有文件合并琐驴，在合并之前，需要對每個序列文件進行排序衡怀，并將多行序列轉(zhuǎn)換為一行序列棍矛，此時用到的工具是seqkit

$conda install seqkit
$vi bash3.sh
for i in *.2
do seqkit sort $i >$i.3
seqkit seq $i.3 -w 0 > $i.3.4
done
$sh bash3.sh
$mkdir new
$mv *.4 new/
$cd new
$paste -d " " *.4 > all.fa
$sed -i "s\ \\g" all.fa

此時我們已經(jīng)把所有的單拷貝序列合并安疗，接下來使用notepad++把所有的ID更改

修改前

修改后

5.4.選擇合適的替代模型

修改結(jié)束后選擇合適的氨基酸替代模型準備構(gòu)建系統(tǒng)發(fā)育樹抛杨，在這里我使用prottest預(yù)測合適的模型，prottest使用phy格式文件荐类，所以用python腳本先將fa文件轉(zhuǎn)換為phy文件怖现，

import re
with open('all.fa', 'r') as fin:
    sequences = [(m.group(1), ''.join(m.group(2).split()))
    for m in re.finditer(r'(?m)^>([^ \n]+)[^\n]*([^>]*)', fin.read())]
with open('all.phy', 'w') as fout:
    fout.write('%d %d\n' % (len(sequences), len(sequences[0][1])))
    for item in sequences:
        fout.write('%-20s %s\n' % item)

從官網(wǎng)下載prottest

$ tar zxf prottest-3.4-20140123.tar.gz
$java -jar /opt/biosoft/prottest-3.4-20140123/prottest-3.4.jar -i all.phy -all-distributions -F -AIC -BIC -tc 0.5 -threads 24 -o prottest.out

在prottest.out中可以看到最佳模型

prottest.out

5.5.構(gòu)建系統(tǒng)發(fā)育樹

使用raxml構(gòu)建系統(tǒng)發(fā)育樹

$wget https://github.com/stamatak/standard-RAxML/archive/master.zip
$unzip master.zip
$cd standard-RAxML 
$make -f Makefile.SSE3.gcc
$make -f Makefile.SSE3.PTHREADS.gcc
$make -f Makefile.SSE3.MPI.gcc
$make -f Makefile.SSE3.HYBRID.gcc
添加至環(huán)境變量

選擇JTT+I+G+F模型構(gòu)建發(fā)育樹，外群選擇江南卷柏和小立碗蘚

$nohup raxmlHPC-PTHREADS-SSE3 -T 16 -f a -x 123 -p 123 -N 1000 -m PROTGAMMAIJTTF -k -O -o Smoellendorffii,Ppatens /
       -n all.tre -s all.fa &

運行結(jié)束后使用figtree打開RAxML_bipartitions.all.tre玉罐，對進化樹進行修改屈嗤。

tree

至此，利用單拷貝基因構(gòu)建系統(tǒng)進化已經(jīng)完成

參考鏈接

http://www.reibang.com/p/16e0bbb2ba19
http://www.chenlianfu.com/?p=2217
https://github.com/stamatak/standard-RAxML
https://github.com/davidemms/OrthoFinder

轉(zhuǎn)載請注明周小釗的博客>>利用orthofinder尋找單拷貝基因構(gòu)建系統(tǒng)發(fā)育樹