教程 | 柚子的祖先啥时候出来的？物种分歧时间推断
介绍
还是接上之前Orthofinder分析后，很多老师和同学都会选择做时间分歧分析。
"时间分歧" 简单来说，就是确定物种分化的大概时间区域。
分歧时间区域计算原理：基于已知的时间节点（已被报道的分支时间节点/具有化石时间）和分子钟计算，推算系统发育树的其他节点发生的时间。
简单补充分子钟原理：任意两个物种，自从分化成两个物种后，他们之间的遗传差异的进化速率应该与分化的时间保持相对稳定。举个例子：人与青蛙的分化点肯定比人与猴子的早，那么人与猴子的序列差异理论上会比人与青蛙的序列差异要小。
目前，采用估算物种分歧时间的程序常见有：mcmctree、r8s..
这次推送主要分享mcmctree的使用方法~
运行
PAML安装
mcmctree主要是在PAML软件包。
wget http://abacus.gene.ucl.ac.uk/software/paml4.9j.tgztar zxf paml4.9j.tgzcd paml4.9jrm bin/* && cd srcmake -f Makefilecp baseml basemlg chi2 codeml evolver infinitesites mcmctree pamp yn00 ../bincd ..export PATH=$PATH:`pwd`/binconda安装~
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 .#直系同源数据集对应的物种基因IDcp $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.phy2. 准备树文件
在进行分歧时间预测，我们至少知道一个节点的大概分歧时间
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=2mcmctree查看FigTree.tre
#NEXUSBEGIN 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