实验目的:比较基因组学,从基因组获取单拷贝基因,并构建有根的进化树,并估计分化时间,最终画出有分化时间的进化树
实验流程:
orthofinder2获取单拷贝基因——muscle多序列比对——Gblocks寻找保守位点——iqtree2构建进化树——MCMCTREE计算分歧时间——
(1)单拷贝基因的查找 Orthofinder;
(2)多序列比对 Muscle / MAFFT / ClustalW / T-coffee, 本实验使用Muscle,准确度最高。 mafft的准确性真的低
(3)调取保守区域,并收尾连接,形成supergene ,使用软件Gblocks
(4)进化树构建 RaxML,PhyML或Mrbayes,本实验使用iqtree2
构建树的教程:https://www.yuque.com/wusheng/gw7a9p/mcc73y
(5)分化时间分析 divergence time,本实验使用mcmctree 是PAML 软件包的一个功能。分化时间进化树可视化使用MCMCTreeR
(6)基因扩张收缩分析 CAFE
(7)基因正选择codeML PAML中一个程序
- 提取所有基因组每个基因的最长的转录本
使用getLongProteins.py脚本获取,最后一个参数是指定相同的基因后的分隔符
getLongProteins.py input.pep.fa out.pep.fa .
checkFa.py *.pep.fa
检测fa文件序列中是否有其他字符,如果有的话会输出其他的字符,需要提前删除这些其他字符。
2.使用Orthofinder进行分析
把上面处理后的蛋白文件放入OrthoFinder/primary
文件夹中,或者是其他名称也行,注意修改下面的-f后文件夹名称。一般尽量把蛋白序列命名为abbr.pep.fa这种格式。
第一次运行使用
orthofinder -t 16 -a 16 -f OrthoFinder/primary
-t 16 -a 16是指定使用的cpu进程数量为16
-f是指定输入的蛋白文件的位置
如果你分析完成了,想再添加新的物种,只需要新建一个文件夹run3,把新加的物种蛋白序列放入run3。再运行下面的命令即可
orthofinder -t 16 -a 16 -b OrthoFinder/Results_Jul08 -f run3
-b参数是指定之前分析的输出的文件夹
-f是指定新增的蛋白序列的文件夹
如果是想删除一个物种,同时新增一些物种,命令还是下面这样
orthofinder -t 16 -a 16 -b OrthoFinder/Results_Jul08 -f run3
但是注意,需要进入 OrthoFinder/Results_Jul08/WorkingDirectory这个文件夹中,找到SpeciesIDs.txt这个文件,在需要删除的物种前加上#
符号即可。新增的物种蛋白序列还是放入run3
orthofinder输出的结果目录文件结构说明
Species_Tree
目录 里面是物种的进化树
SpeciesTree_rooted.txt 基于所有蛋白序列构建的进化树(可能不是真实的进化关系,因为涉及到染色体倒位等结构变异)
SpeciesTree_rooted_node_labels.txt 和上面一样的树,有label
Single_Copy_Orthologue_Sequences
里面是单拷贝直系同源基因的fa文件,每个fa文件里是每一个物种的一条序列。
Orthogroups
目录
Orthogroups_SingleCopyOrthologues.txt 单拷贝基因的OG
开头的编号的列表文件。总行数就是每个物种里单拷贝基因的数量。
Comparative_Genomics_Statistics
目录
OrthologuesStats_one-to-one.tsv 这个文件里可以看到每个基因组和其他基因组比较,1对1的直系同源基因的数量。
如果singlecopy的基因数量比较少,这时候需要使用参数过滤,例如:在70%的物种里是单拷贝的基因的列表,这样能够多一些。使用getSingleCopy.py Results_XXXXX/Orthogroups/Orthogroups.GeneCount.tsv all 0.7
即可获取对应的基因家族的名称。
从Results_XXXXX/Orthogroup_Sequence里找到对应家族的序列。
3.使用singlecopy2tree.bash
脚本完成比对和保守位点提取
bash singlecopy2tree.bash OrthoFinder/Results_Jul08/Single_Copy_Orthologue_Sequences orthofinder_inputpath
- 使用iqtree2构建单拷贝基因的进化树
##使用iqtree2构建进化树(此处需要16个cpu,其实少一点也可以,速度很快)
iqtree2 -s all.fa -T 16 -B 1000 -m MFP --bnni -o ${outgroup}
5.进化树进化时间估计
从http://www.timetree.org/查询两个物种之间的分化时间。输入拉丁学名,即可查询到。
一般情况下,使用的是CI后面的两个值,对应图中的3.9和53.1。这个例子不太好,差的有点太大了。最好选择两个物种之间估计的分化时间差距不大的。注意不是两个物种分化时间不大,是CI后的括号内的值的差距小。
PAML使用的是CI后面的值,R8S使用的是Median Time值。
PAML 的mcmctree用于计算分化时间(速度比较慢,但是比较准)
R8S也可以用于计算分化时间(速度快,准确率没有mcmctree高)
#R8S安装方法
wget -q https://sourceforge.net/projects/r8s/files/r8s1.81.tar.gz \
&& tar -zxvf r8s1.81.tar.gz \
&& cd r8s1.81/src \
&& sed -i 's|/usr/include/sys/errno.h||' Makefile.linux \
&& sed -i 's/continuousML.o //' Makefile.linux \
&& sed -i 's/continuousML.o:/#continuousML.o:/' Makefile.linux \
&& make -f Makefile.linux
OrthoFinder/Results_Jul08/Orthogroups/Orthogroups.GeneCount.tsv:每一行对应一个基因家族在该基因组对应的基因数目;
获取所有的直系同源基因的ID列表
cat name.list|tr ">" "\t" >SingleGeneID.tsv
awk '
{
for (i=1; i<=NF; i++) {
if (NR == 1) {
transposed[i] = $i
} else {
transposed[i] = transposed[i] "\t" $i
}
}
}
END {
for (i=1; i<=NF; i++) {
print transposed[i]
}
}' SingleGeneID.tsv > Specs.singlecopy.id
tail -1 Specs.singlecopy.id >header
sed -i '$ d' Specs.singlecopy.id
cat header Specs.singlecopy.id > Species.singlecopy.tsv
rm -rf SingleGeneID.tsv Specs.singlecopy.id header
最终Species.singlecopy.tsv是所有的单拷贝基因的列表,第一行是每个物种的名称,每一列是对应的基因的ID。每一行之间是单拷贝直系同源基因。