Comparative Genomics：物种进化树

比较基因组分析-构建物种进化树

比较基因组分析：构建带有分化时间的物种进化树

1 准备蛋白质序列输入文件

1.1 激活环境


xxxxxxxxxx
1
1
conda activate comp_genome

1.2 安装软件

Gblocks、muscle、modeltest、raxml-ng、PAML、mcmctree


1
1
conda install <software>

1.3 多序列比对和保守区域提取：muscle和Gblocks

输入文件：Orthofinder的运行结果文件夹中Single_Copy_Orthologue_Sequences的112条单拷贝基因

输入文件目录

新建muscle_And_gblocks.sh运行文件，先进行多序列比对获得保守结构域，再将保守结构域进行提取。


10
1
cd /home/hztext/dai/Comp_Genome/Nine_ortho_result/Single_Copy_Orthologue_Sequences/
2
for i in *.fa
3
do
4
 file=${i%.fa*}
5
 muscle -in "${file}.fa" -out "${file}.aln" ## 多序列比对
6
 Gblocks "${file}.aln" -b4=5 -b5=h -t=p -e=.gb ## 提取保守序列
7
 seqkit seq "${file}.aln.gb" -w 0 > "${file}.new.aln" ## 多行序列并为一行
8
 awk '{if (NR%2==1) print substr ($1, 1, 4); else print $0}' "${file}.new.aln" > "${file}.final.aln" ## ID 修饰
9
done
10
bash /home/hztext/dai/Comp_Genome/Demo_deal/reminder.sh "您的多序列比对等结果已输出"

1.4 串联蛋白质序列


1
cd /home/hztext/dai/Comp_Genome/Nine_ortho_result/Single_Copy_Orthologue_Sequences/
2
#串联文件
3
seqkit concat *.final.aln > all.fa
4
#并去除文件中的空格（因为不同aa之间存在空格，像ncbi的格式）
5
sed -i "s/\ //g" all.fa
6
bash /home/hztext/dai/Comp_Genome/Demo_deal/reminder.sh "您的单拷贝序列已串联完毕"

final.aln文件：

final.aln文件

all.fa文件：将串联文件中的序列ID改为物种名简写，检查每个物种的氨基酸数量是否一致（必须一致）。

串联文件

1.5 fa格式转换为phy格式

新建trans_fa_To_phy.py命令文件。


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

新建do_trans_fa_To_phy.sh命令文件执行trans_fa_To_phy.py程序。


4
1
cd /home/hztext/dai/Comp_Genome/Nine_ortho_result/Single_Copy_Orthologue_Sequences/
2
python /home/hztext/dai/Comp_Genome/Demo_deal/trans_fa_To_phy.py
3
bash /home/hztext/dai/Comp_Genome/Demo_deal/reminder.sh "您转phy格式已完成"
4

格式转化结束后，在每个物种的名后添加8个空格

phy格式

2 物种进化树构建

2.1 利用modeltest-ng预测最佳的模型

新建modeltest.sh文件并执行


3
1
cd /home/hztext/dai/Comp_Genome/Nine_ortho_result/Single_Copy_Orthologue_Sequences/
2
modeltest-ng -i all.fa -d aa -o modeltest -p 45
3
bash /home/hztext/dai/Comp_Genome/Demo_deal/reminder.sh "您的建树模型预测已完成"

-i：蛋白质序列文件→包含所有物种比对后的保守结构域串联文件

-d：aa代表输入的文件是氨基酸序列（amino acid）

-o：输出文件前缀

-p：线程数

建树模型预测的结果文件

modeltest.out文件中的BIC（贝叶斯信息准则）最佳模型

最佳建树模型

2.2 利用raxml-ng来构建物种进化树


3
1
cd /home/hztext/dai/Comp_Genome/Nine_ortho_result/Single_Copy_Orthologue_Sequences/
2
raxml-ng --all --msa all.fa --model JTT+I+G4+F --bs-trees 100
3
bash /home/hztext/dai/Comp_Genome/Demo_deal/reminder.sh "您的raxml-ng建树已完成"

–msa：输入蛋白质序列文件；

–model：最佳建树模型。

raxml-ng结果文件

all.fa.raxml.support是构建的物种进化树文件，如下


xxxxxxxxxx
1
1
(((Cma:0.008612,Ccl:0.007167)100:0.099212,((TLG:0.012426,Tci:0.010765)100:0.005563,Tsi:0.020005)100:0.063715)100:0.038323,((Atr:0.575439,Vvi:0.159512)100:0.033988,Egr:0.287300)100:0.074467,Aya:0.132865);

3 构建带有进化时间的物种进化树

3.1 利用PAML估计物种的进化时间

基于已知物种间的分化时间，估计其余物种的进化时间，可以从已发表的基因组文章、TimeTree网站获取。

TimeTree网站

3.2 物种进化树输入文件预处理

下载all.fa.raxml.support文件（无根树）到windows，使用iTOL网站查看进化树，修改为有根树并指定根分支（外类群中进化最早的）。

iTOL网站

根据进化树中分支的长度，可以发现Atr的分支长度最长，即相比于其他物种的进化时间更早。

物种进化树

为了后续进化树的展示效果，我们将Atr分支设置为根分支（root），并忽略分支长度。

iTOL设置根分支

导出修改后的进化树文件，上传至服务器DiffTimetree文件夹。

导出修改后的进化树

新建tree_deal.r命令文件。

树文件：


x
1
(Atr:0.2877195,(Vvi:0.159512,(Egr:0.287300,(Aya:0.132865,((Ccl:0.007167,Cma:0.008612)100:0.099212,(Tsi:0.020005,(Tci:0.010765,TLG:0.012426)100:0.005563)100:0.063715)100:0.038323)100:0.074467)100:0.033988):0.2877195);

tree_deal.r：


6
1
if(!requireNamespace("ape", quietly = TRUE)) install.packages("ape")
2
library(ape)
3
MLtree <- read.tree("Tree_reroot.txt") #读取单拷贝进化树
4
MLtree$edge.length <- NULL #删除分支长度信息
5
MLtree$node.label <- NULL
6
write.tree(MLtree,"Tree_reroot.tree") #输出只有物种名的进化树

新建do_tree_deal.sh命令文件执行tree_deal.r，处理进化树文件。


2
1
cd /home/hztext/dai/Comp_Genome/Nine_ortho_result/DiffTimetree/
2
Rscript /home/hztext/dai/Comp_Genome/Demo_deal/tree_deal.r

修改上述进化树，增加物种分化时间（只添加确定的2-3个），转为以下格式的进化树输入文件：


x
1
9 1
2
(Atr,(Vvi,(Egr,(Aya,((Ccl,Cma)'>1.5<5.7',(Tsi,(Tci,TLG))))))'>179.9<205.0');

9代表树种数量，1固定不变。

3.3 利用PAML的mcmctree构建带有分化时间的进化树

准备输入文件all.phy、Tree_reroot.tree、mcmctree.ctl

mcmctree.ctl：（修改seqfile、 treefile的文件位置，usedata设为3）


x
1
       seed = -1
2
       seqfile = ./all.phy
3
      treefile = ./Tree_reroot.tree
4
       outfile = out
5

6
         ndata = 1 ## 必须注意为 1 
7
       seqtype = 2  * 0: nucleotides; 1:codons; 2:AAs ## 序列格式，蛋白就选择 AA
8
       usedata = 3    * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV ## 此步骤先设为 3
9
         clock = 3    * 1: global clock; 2: independent rates; 3: correlated rates
10
       RootAge = <450.0  * safe constraint on root age, used if no fossil for root.
11

12
         model = 0    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
13
         alpha = 0    * alpha for gamma rates at sites
14
         ncatG = 5    * No. categories in discrete gamma
15

16
     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?
17

18
       BDparas = 1 1 0    * birth, death, sampling
19
   kappa_gamma = 6 2      * gamma prior for kappa
20
   alpha_gamma = 1 1      * gamma prior for alpha
21

22
   rgene_gamma = 2 2   * gamma prior for overall rates for genes
23
  sigma2_gamma = 1 10    * gamma prior for sigma^2     (for clock=2 or 3)
24

25
      finetune = 1: 0.1  0.1  0.1  0.01 .5  * auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr
26

27
         print = 1
28
        burnin = 20000
29
      sampfreq = 2
30
       nsample = 100000

命令行输入以下命令执行。


1
1
mcmctree mcmctree.ctl

获得的结果文件如下：

mcmctree第一次运行结果

利用第一次运行结果中的out.BV作为输入文件，进行mcmctree的第二次建树。

新建Second_run_mcmctree.sh运行文件。


xxxxxxxxxx
1
10
9
1
cd /home/hztext/dai/Comp_Genome/Nine_ortho_result/DiffTimetree
2
#新建第二次运行的目录
3
mkdir -p ../SecondRun
4
#拷贝第一次运行结果out.BV，重命名为in.BV
5
cp ./out.BV ../SecondRun/in.BV
6
#拷贝比对文件、树文件和mcmctree的参数文件
7
cp ./all.phy ../SecondRun/all.phy
8
cp ./Tree_reroot.tree ../SecondRun/Tree_reroot.tree
9
cp ./mcmctree.ctl ../SecondRun/mcmctree.ctl

准备mcmctree.ctl：（usedata设为2）


xxxxxxxxxx
1
31
1
       seed = -1
2
       seqfile = ./all.phy
3
      treefile = ./Tree_reroot.tree
4
       outfile = out
5

6
         ndata = 1 ## 必须注意为 1 
7
       seqtype = 2  * 0: nucleotides; 1:codons; 2:AAs ## 序列格式，蛋白就选择 AA
8
       usedata = 2    * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV ## 此步骤先设为 3
9
         clock = 3    * 1: global clock; 2: independent rates; 3: correlated rates
10
       RootAge = <450.0  * safe constraint on root age, used if no fossil for root.
11

12
         model = 0    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
13
         alpha = 0    * alpha for gamma rates at sites
14
         ncatG = 5    * No. categories in discrete gamma
15

16
     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?
17

18
       BDparas = 1 1 0    * birth, death, sampling
19
   kappa_gamma = 6 2      * gamma prior for kappa
20
   alpha_gamma = 1 1      * gamma prior for alpha
21

22
   rgene_gamma = 2 2   * gamma prior for overall rates for genes
23
  sigma2_gamma = 1 10    * gamma prior for sigma^2     (for clock=2 or 3)
24

25
      finetune = 1: 0.1  0.1  0.1  0.01 .5  * auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr
26

27
         print = 1
28
        burnin = 20000
29
      sampfreq = 2
30
       nsample = 100000

命令行输入以下命令执行。


1
1
mcmctree mcmctree.ctl

获得如下结果文件，FigTree.tre是我们需要的最终进化树文件。

mcmctree获得的物种进化树

4 FigTree软件美化物种进化树

FigTree下载链接：https://github.com/rambaut/figtree/releases

大家先自行摸索，后续再出美化教程哈哈哈哈哈哈哈哈！！！