教程 | 柚子的祖先啥时候出来的？物种分歧推断

介绍

还是接上之前Orthofinder分析后，很多老师和同学都会选择做时间分歧分析。

Orthofinder下篇

Orthofinder上篇

“时间分歧” 简单来说，就是确定物种分化的大概时间区域。

分歧时间区域计算原理：基于已知的时间节点（已被报道的分支时间节点/具有化石时间）和分子钟计算，推算系统发育树的其他节点发生的时间。

简单补充分子钟原理：任意两个物种，自从分化成两个物种后，他们之间的遗传差异的进化速率应该与分化的时间保持相对稳定。举个例子：人与青蛙的分化点肯定比人与猴子的早，那么人与猴子的序列差异理论上会比人与青蛙的序列差异要小。

目前，采用估算物种分歧时间的程序常见有：mcmctree、r8s…

这次推送主要分享mcmctree的使用方法~

运行

PAML安装

mcmctree主要是在PAML软件包。

wget http://abacus.gene.ucl.ac.uk/software/paml4.9j.tgz
tar zxf paml4.9j.tgz
cd paml4.9j
rm bin/*  && cd src
make -f Makefile
cp baseml basemlg chi2 codeml evolver infinitesites mcmctree pamp yn00 ../bin
cd ..
export PATH=$PATH:`pwd`/bin

conda安装~

conda install -c bioconda paml

两个方法下载的paml的软件包都是4.9版本的~

主要用到paml软件包的mcmc，其中学习手册的介绍：Bayesian estimation of species divergence times incorporating uncertainties in fossil calibrations (mcmctree).

paml的学习手册：http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf

准备文件

我们在上次推文已经完成了orthofinder的分析，所以从其分析结果获mcmctree运行所需的文件。

主要需要两个文件：

1.单拷贝直系同源基因数据集合的基因ID

2. 物种的系统发育树文件

mkdir mcmctreeAnalysis && cd mcmctreeAnalysis

#定义我Orthofinder结果路径
OrthofinderPATH=$mypath/OrthoFinder/Results_Jun12

#复制单拷贝直系同源的数据集合列表到工作目录
cp $OrthofinderPATH/Orthogroups/Orthogroups_SingleCopyOrthologues.txt .
#直系同源数据集对应的物种基因ID
cp $OrthofinderPATH/Orthogroups/Orthogroups.txt .
#复制物种系统发育树文件作惨啊空
cp $OrthofinderPATH/Species_Tree/SpeciesTree_rooted.txt .

1. 准备序列比对文件

#提取各物种单拷贝基因的CDS序列
#使用的六个物种，所以我也准备好四个物种的最长CDS序列；
#Cpapaya.fasta
#Csinensis.fasta
#Vvinifera.fasta
#Fvesca.fasta
#Dzibethinus.fasta
#Tcacao.fasta

cat *.fasta > all.cds.fasta #因为我每个物种的CDS都以fasta结尾，所以直接合并
#可选,为了让基因ID可以简洁一点
cut -f1  -d " " all.cds.fasta >  tmp.fasta && mv tmp.fasta all.cds.fasta 

mkdir CDSfile #存放单个文件cds序列
for i in `cat Orthogroups_SingleCopyOrthologues.txt`; do grep $i Orthogroups.txt|sed 's/ /\n/g'  |sed '1d' > CDSfile/${i}.id ;done

for i in `cat Orthogroups_SingleCopyOrthologues.txt`;do seqkit grep -f CDSfile/${i}.id all.cds.fasta > CDSfile/${i}.cds ;done 

#为了让所有id都全部转换成以物种作id, 最好保证每个物种名字不超过7个字符（phy格式）
for i in `cat Orthogroups_SingleCopyOrthologues.txt`;do sed 's/Cs[0-9g.]\+/Csi/' CDSfile/${i}.cds|sed 's/orange[0-9t.]\+/Csi/' | sed 's/evm[a-z_0-9.]\+/Cpa/' |sed 's/GSVIVT[0-9]\+/Vvi/' |sed 's/LOC[0-9]\+/Dzi/' |sed 's/FvH4_[0-9g.t]\+/Fve/'|sed 's/Thecc.[0-9G.v]\+/Tca/' > CDSfile/${i}.rename.cds;done

#查看了一下一共多少个文件，事实上我物种数量不多，所以相对的共同单拷贝基因数量会多；如果增加了物种数量，单拷贝基因就会少。
ls CDSfile/*.id|wc
# 1207

#比对,形成phy格式（如果单拷贝基因选得比较多，比对时间会稍微有点长）
for i in `cat Orthogroups_SingleCopyOrthologues.txt`;do muscle -in  CDSfile/${i}.rename.cds -physout   CDSfile/${i}.phy;done 
cat CDSfile/*.phy > all.singlecopy.phy

2. 准备树文件

在进行分歧时间预测，我们至少知道一个节点的大概分歧时间

TIMETREE数据库：http://timetree.org/ 是一个记录多个植物分歧进化时间的公共数据库。而且提供相关文献，我们可以进一步进行确定

这边我们暂定葡萄与拟南芥/柑橘的分歧点约为大于 100

cat  SpeciesTree_rooted.txt
#(Vvinifera:0.0882857,(Fvesca:0.224125,(Csinensis:0.181856,((Tcacao:0.0563349,Dzibethinus:0.0762371)1:0.105748,Cpapaya:0.200521)1:0.0199089)1:0.0346594)1:0.0882857);

参考orthofinder的系统进化关系

vi treefile
#treefile
#记得修改和序列比对文件对应的物种名称
6 1
(Vvi,(Fve,(Csi((Dzi,Tca),Cpa)))'L(1.0,.1,.2)');

稍微解释一下如何编写这个物种树，借鉴Orthofinder的构建好的树结构，补充了葡萄分歧时间大于100MYB。

也可以直接使用 ‘L(1.0)’、’>1.0’

如果是在某个时间段可以：‘B(1.0,1.17)’、’>1.0<1.17’

如果是小于某个时间节点可以：‘U(1.17)’、’<1.17’

而括号后的两个参数分为 p=0.1, c=0.2, 由于物种数量比较小而且距离还相对较远，增加两个内置函数更严格限制已知的分歧点。

更详细请参考PAML User Guide的 第9节mcmctree的Fossil calibration 部分。

mcmctree运行

#在安装路径寻找mcmctree配置文件
cp $pamlPath/mcmctree.ctl .
##如果conda安装的没有配置文件，vim mcmctree.ctl

简单说一下，最主要填写的项：

seqfile = all.singlecopy.phy

treefile = treefile

outfile = out *输出前缀

ndata = 30

建议将随机采取的比对序列情况设置大一点；因为我这里六个物种单拷贝基因数量太多了，运行起来耗时多，所以我仅采用30个情况进行时间分歧计算。

usedata

设置是否利用多序列比对的数据：0，表示不使用多序列比对数据，则不会进行likelihood计算；1，表示使用多序列比对数据使用Felsenstein（1981）的似然计算，也是MCMC常使用的参数，不过计算速度较慢; 2，进行approximation likelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的in.BV文件。该文件是使用usedata = 3参数生成的out.BV文件重命名而来的。

clock = 3

model = 4

0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 说明书列举了四种情况，其中JC69是适合比较简单的关系，HKY85是相对比较远源的关系（更复杂的模型，配合alpha=0.5使用）。

alpha = 0.5

首先设置usedate=3

#第一次运行
mcmctree

获得 out.BV文件

usedate=2

mv out.BV in.BV
#修改mcmctree.ctl 将usedate=2
mcmctree

查看FigTree.tre

#NEXUS
BEGIN TREES;

        UTREE 1 = (Vvi: 1.017597, (Fve: 1.012650, (Csi: 0.888532, ((Dzi: 0.267560, Tca: 0.267560) [&95%={0.241975, 0.297312}]: 0.517603, Cpa: 0.785163) [&95%={0.736974, 0.849596}]: 0.103369) [&95%={0.839088, 0.957735}]: 0.124118) [&95%={0.963919, 1.08662}]: 0.004947) [&95%={0.968735, 1.09186}];

在Figtree可以可视化，后期导出SVP。

至于如何更美化这个分歧树… 就各自修行。

mcmctree.ctl配置

seed = -1
       seqfile = all.singlecopy.phy
      treefile = treefile
       outfile = out

         ndata = 20
       seqtype = 0  * 0: nucleotides; 1:codons; 2:AAs
       usedata = 2    * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
         clock = 3    * 1: global clock; 2: independent rates; 3: correlated rates
       RootAge = '<1.0'  * safe constraint on root age, used if no fossil for root.

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

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

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

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

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

         print = 1
        burnin = 2000
      sampfreq = 10
       nsample = 20000