利用 orthofinder 寻找单拷贝基因构建系统发育树

# 1.orthofinder 介绍

OrthoFinder 是一种快速、准确和全面的比较基因组学分析工具。它可以找到直系和正群，为所有的正群推断基因树，并为所分析的物种推断一个有根的物种树。OrthoFinder 还为比较基因组分析提供全面的统计数据。OrthoFinder 使用简单，只需运行一组 FASTA 格式的蛋白质序列文件（每个物种一个）。

# 2. 基础知识介绍

Orthologue（直系同源基因）指的是来自两个物种的基因。Orthologue 是由两个物种的最后共同祖先 (LCA) 上的单个基因进化而来的成对基因 (图 1A 和 B)。正群是同源概念在物种群中的自然延伸。一个 Orthogroup（正交群）是由一个物种的 LCA 中的单个基因进化而来的一组基因 (图 1A)。当观察基因树时，一个邻位群体中基因的第一次分化是一个物种形成事件，对同源基因来说也是如此。

作为基因复制事件的结果，当观察直系同源基因和正交群时，可能会有来自同一物种的多个基因。在这个例子中 (图 1A 和 B)，人类和老鼠的 HuA 基因是鸡中 ChA1 和 ChA2 的同源基因。再看一下正交群，我们发现有两个鸡的基因 (图 1A)，但是只有一个来自老鼠和人类的基因。一些作者将 ChA1 和 ChA2 基因作为 HuA 的共同源基因，以强调存在多个同源基因的事实。由于基因重复和丢失在进化中经常发生，一对一的直系同源物很少见，通过分析正交群所有直系同源的情况（一对一，多对一，多对多），我们可以分析数据的所有情况。

paralogues (旁系同源基因) 是指在基因复制事件中从单个基因中分离出来的成对基因，鸡的两个基因 ChA1 和 ChA2 是旁系同源基因 (图 1C)。来自不同物种的两个基因如果在基因重复事件中彼此分离，也可能是同源的。由于基因树中所有的分支事件要么是物种形成事件 (产生直系同源基因)，要么是重复事件 (产生旁系同源基因)，因此同一正交群中任何不是直系同源基因的基因必然是旁系同源基因。

直系同源物是同源性基因，是物种形成事件的结果。Paralogs（旁系同源物) 是同源基因，是重复事件的结果。可以看到 (图 2)，不同物种间的 α-chain gene 互为 Orthologs (直系同源物)。正交群用来形容自一组物种的 LCA 中的单个基因的基因组（α-chain gene）。然后同一物种间 α 和 β chain gene 互为 Paralogs (旁系同源物)。最后所有这些关系都可以由 OrthoFinder 来识别。

# 3. 安装

在这里推荐大家使用 conda 来安装，简单明了，不用担心其他 Dependencies 的安装

$conda install orthofinder

# 4. 运行 orthofinder

运行 Orthofinder，相当简单的操作，"-f" 输入目录，里面包含你需要运行的蛋白质 fasta 文件，“-t” 所用到的 CPU 数目，“-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 &

生成的结果会存储于 Orthofinder/Results_XXX 文件中，现在简单看看里面有啥。

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 里的单拷贝基因构建系统发育树

WorkingDirectory 其中包含运算过程的中间文件，例如 blast 结果，如果我们想去掉某一物种，在 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. 利用单拷贝基因构建系统发育树

Orthofinder 的 Single_Copy_Orthologue_Sequences 下存放着单拷贝同源基因的序列，我们可以利用这些序列构建系统发育树

# 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 添加至环境变量

单拷贝的序列较多，所以使用 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

添加至环境变量

$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 设置序列的类型，可选的值是 p,d,c 分别代表 protein， DNA， Codons 。
-b1= default:（ 序列条数的 50% + 1 ）, 设定保守性位点必须有 >= 该值的序列数。该参数后接一个 integer 数，默认下比序列条数的 50% 大 1.
-b2= default: 序列条数的 85%, 确定保守位点的侧翼位点时，其位点必须有 >= 该值的序列数。该值必须要比 -b1 的值要大。
-b3= default: 8, 最大连续非保守位点的长度。
-b4= default: 10, 保守位点区块的最小长度。该值必须 >=2 。
-b5= default: n, 设置允许含有 Gap 位点。可选的值有 n,h,a 分别代表 None, With Half, All 。 当为 h 时，表示
- e= default: .2, 设置输出结果的后缀。

# 5.3. 序列合并

比对好之后我们需要将所有文件合并，在合并之前，需要对每个序列文件进行排序，并将多行序列转换为一行序列，此时用到的工具是 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

此时我们已经把所有的单拷贝序列合并，接下来使用 notepad++ 把所有的 ID 更改

# 5.4. 选择合适的替代模型

修改结束后选择合适的氨基酸替代模型准备构建系统发育树，在这里我使用 prottest 预测合适的模型，prottest 使用 phy 格式文件，所以用 python 脚本先将 fa 文件转换为 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)

从 [] 官网上](https://www.softpedia.com/dyn-postdownload.php/cc45406e35260b47bfa4132e67f8c446/5f13ff63/286c7/4/1) 下载 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 中可以看到最佳模型

# 5.5. 构建系统发育树

使用 raxml 构建系统发育树

$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

添加至环境变量

选择 JTT+I+G+F 模型构建发育树，外群选择江南卷柏和小立碗藓

$nohup raxmlHPC-PTHREADS-SSE3 -T 16 -f a -x 123 -p 123 -N 1000 -m PROTGAMMAJTT -k -O -o Smoellendorffii,Ppatens /

       -n all.tre -s all.fa &

运行结束后使用 figtree 打开 RAxML_bipartitions.all.tre，对进化树进行修改。

至此，利用单拷贝基因构建系统进化已经完成

# 参考链接

https://www.jianshu.com/p/16e0bbb2ba19

http://www.chenlianfu.com/?p=2217

https://github.com/stamatak/standard-RAxML

https://github.com/davidemms/OrthoFinder

转载请注明周小钊的博客 >> 利用 orthofinder 寻找单拷贝基因构建系统发育树