我觉得处理任何事情，不仅要做到“知其然”，还要做到“知其所以然”。在做事情时，如果仅知道低头蛮干，不知道抬头看路；只知道跟着“牵牛绳”的方向，在鞭子的抽打下前进或停止，那么最终等待我们的就是被淘汰。
01. 定义
祖先重建（或特征映射）是指从个体或种群的测量特征推断其共同祖先的特征，它是系统发育、个体、种群或物种与祖先进化关系的重建和研究的重要应用。在演化生物学的研究中，祖先重建可以用来重建数百万年前生活的生物的不同类型的祖先特征状态。这些状态包括遗传序列（祖先序列重建）、蛋白质的氨基酸序列、基因组组成、生物表型以及祖先种群或物种的地理范围（祖先范围重建）。由于现代遗传序列本质上是古代序列的变异，因此进入古代序列可能会发现其他变异和可能由这些序列产生的生物体。除了遗传序列外，生物学家可能会试图追踪一个性状到另一个性状的演化历史，例如，鳍向四肢的转变。
02. 方法原理
祖先重建常常依赖于一个足够真实的演化统计模型来准确地推断祖先状态，这些模型利用已经通过系统发育等方法获得的遗传信息来确定演化历史以及发生演化事件的时间。然而，无论模型多么接近实际的演化历史，准确重建祖先状态的能力都会随着该祖先与其后代之间的演化时间的增加而减弱。祖先重建的方法主要有三种：最大简约法（maximum parsimony）、最大似然法（maximum likelihood）以及贝叶斯法（bayesian）。其中，最大简约法认为所有演化事件的可能性相同；而最大似然法认为某些类别演化事件具有不同的可能性；贝叶斯推断则将事件的条件概率与系统发育树的可能性以及与该发育树相关联的不确定性的数量联系起来。
03. 方法差异
最大简约法主要提供快速和简单的方法来推断祖先的状态，然而，由于演化过程过于简单化（例如，没有考虑分支长度和演化时间），简约法的准确性有限。最大似然法是基于性状演化的概率模型。在各种条件下，利用理论论证和模拟研究，最大似然方法的性能优于简约法。仿真结果表明，由于考虑了进化时间和分支长度，即使是最简单的模型也能产生比简约更精确的结果，并且对系统发育的不确定性具有稳健性。贝叶斯方法是利用马尔科夫链蒙特卡洛方法（MCMC）来推断祖先特征状态、系统发育树和模型参数的联合后验分布，这涉及到描述序列演化的复杂概率模型、分子钟、人口统计学以及研究特征的演化。由于贝叶斯方法具有丰富的选择性以及灵活性而备受欢迎，例如著名的系统发育分析软件BEAST。然而，MCMC具有较高的计算成本，对于大数据集，无法实现所有树、参数和性状分布的联合推理。
04. 性状演化模型
许多模型已经发展起来用以估计从现生后代到祖先的离散和连续型性状的演变。这些模型假设一个性状通过时间的演化可以被建模为一个随机过程。
a. 连续型性状（如，体重、脑重等）
连续型性状演化模型被用来理解许多有趣的演化现象，包括：性状演化的速率，性状演化的模式，这些性状如何随着时间或谱系之间的变化以及理解导致这种变化的生物和非生物因素。
布朗运动模型（Brownian Motion models），最简单的连续性状演化模型假定性状在布朗运动模型下演化。在此模型下，沿分支预期的性状变化量为零，性状变化的方差与时间和速率参数成正比。它可以估算某性状整体的演化速率、检测演化速度随时间的变化、检测不同谱系进化速率的变化、检测与状态相关的演化速率以及检测性状相关性演化（如，性状X与性状Y在演化中是否具有相关性）。
Ornstein-Uhlenbeck 模型主要估计演化最优值以及演化最优的变化检测。
b. 离散型性状（如，花的颜色、形状）
对于这种性状，通常是马尔科夫链模型（Markov chain）-马尔可夫链是描述一系列可能事件的随机模型，其中每个事件的概率仅取决于在前一事件中达到的状态，在连续时间内，它被称为马尔可夫过程。
05. 软件工具
有许多软件包可以执行祖先状态重建，一般来说，这些软件包是通过相关领域的科学家的努力开发和维护的，并在免费软件许可下发布。下表展示了一些比较经典的软件工具：
| Name | Methods | Platform | ! Character Types | Continuous (C) or Discrete Characters (D) | Software License | |
|---|---|---|---|---|---|---|
| PAML | Maximum Likelihood | Unix, Mac, Win | PHYLIP, NEXUS, FASTA | Nucleotide, Protein | D | GNU General Public License, version 3 |
| BEAST2 | Bayesian | Unix, Mac, Win | NEXUS, BEAST XML | Nucleotide, Protein, Geographic | C, D | GNU Lesser General Public License |
| phytools | Maximum Likelihood | Unix, Mac, Win | newick, nexus | Qualitative and quantitative traits | C, D | GNU General Public License |
| APE | Maximum Likelihood | Unix, Mac, Win | NEXUS, FASTA, CLUSTAL | Nucleotide, Protein | C, D | GNU General Public License |
| Diversitree | Maximum Likelihood | Unix, Mac, Win | NEXUS | Qualitative and quantitative traits, Geographic | C, D | GNU General Public License, version 2 |
| HyPhy | Maximum Likelihood | Unix, Mac, Win | MEGA, NEXUS, FASTA, PHYLIP | Nucleotide, Protein (customizable) | D | GNU Free Documentation License1.3 |
| BayesTraits | Bayesian | Unix, Mac, Win | TSV or space delimited table. Rows are species, columns are traits. | Qualitative and quantitative traits | C, D | Creative Commons Attribution License |
| Lagrange | Maximum Likelihood | Linux, Mac, Win | TSV/CSV of species regions. Rows are species and columns are geographic regions | Geographic | - | GNU General Public License, version 2 |
| Mesquite | Parsimony, Maximum Likelihood | Unix, Mac, Win | Fasta, NBRF, Genbank, PHYLIP, CLUSTAL, TSV | Nucleotide, Protein, Geographic | C, D | Creative Commons Attribution 3.0 License |
| Phylomapper | Maximum Likelihood, Bayesian (as of version 2) | Unix, Mac, Win | NEXUS | Geographic, Ecological niche | C, D | - |
| Ancestors | Maximum Likelihood | Web | Fasta | Nucleotide (indels) | D | - |
| Phyrex | Maximum Parsimony | Linux | Fasta | Gene expression | C, D | Proprietary |
| SIMMAP | Stochastic Mapping | Mac | XML-like format | Nucleotide, qualitative traits | D | Proprietary |
| MrBayes | Bayesian | Unix, Mac, Win | NEXUS | Nucleotide, Protein | D | GNU General Public License |
| PARANA | Maximum Parsimony | Unix, Mac, Win | Newick | Biological networks | D | Apache License |
| PHAST(PREQUEL) | Maximum Likelihood | Unix, Mac, Win | Multiple Alignment | Nucleotide | D | BSD License |
| RASP | Maximum Likelihood, Bayesian | Unix, Mac, Win | Newick | Geographic | D | - |
| VIP | Maximum Parsimony | Linux, Win | Newick | Geographic | D (grid) | GPL Creative Commons |
| FastML | Maximum Likelihood | Web, Unix | Fasta | Nucleotide, Protein | D | Copyright |
| MLGO | Maximum likelihood | Web | Custom | Gene order permutation | D | GNU |
| BADGER | Bayesian | Unix, Mac, Win | Custom | Gene order permutation | D | GNU GPL version 2 |
| COUNT | Maximum Parsimony, maximum likelihood | Unix, Mac, Win | Tab-delimited text file of rows for taxa and count data in columns. | Count (numerical) data (e.g., homolog family size) | D | BSD |
| MEGA | Maximum parsimony, maximum likelihood. | Mac, Win | MEGA | Nucleotide, Protein | D | Proprietary |
| ANGES | Local Parsimony | Unix | Custom | Genome maps | D | GNU General Public License, version 3 |
| EREM | Maximum likelihood. | Win, Unix, Matlab module | Custom text format for model parameters, tree, observed character values. | Binary | D | None specified, although site indicates software is freely available. |
06. 最新网页版软件
最近有两款软件分别发布于MBE期刊（PastML）与BMC evolutionary biology期刊（pastview），网页版的优点就是使用方便、容易出图。一般只需要准备一颗系统发育树以及对应的性状矩阵文件，即可快速完成祖先状态的重建。
07. 软件实操
在上述中，我已经列举了许多软件，接下来主要是选择R包phytools来进行演示如何重构祖先状态。
a. 连续型数据
install.packages("phytools")
library(phytools)
anole.tree <- read.tree("http://www.phytools.org/eqg2015/data/anole.tre") # 树文件
plotTree(anole.tree,type="fan",ftype="i")
svl <- read.csv("http://www.phytools.org/eqg2015/data/svl.csv",row.names=1) #数据矩阵
svl<-as.matrix(svl)[,1] #把这个变成一个向量
svl
# 估算祖先状态，并计算方差与95%置信度水平区间
fit <- fastAnc(anole.tree,svl,vars=TRUE,CI=TRUE)
fit
fit$CI[1,]
range(svl)
obj<-contMap(anole.tree,svl,plot=FALSE)
plot(obj,type="fan",legend=0.7*max(nodeHeights(anole.tree)),fsize=c(0.7,0.9))
phenogram(anole.tree,svl,fsize=0.6,spread.costs=c(1,0))
## 祖先状态的属性模拟
## simulate a tree & some data
tree <- pbtree(n=26,scale=1,tip.label=LETTERS)
## simulate with ancestral states
x <- fastBM(tree,internal=TRUE)
## ancestral states
a <- x[as.character(1:tree$Nnode+Ntip(tree))]
## tip data
x <- x[tree$tip.label]
# 使用Fastanc估计祖先状态
fit <- fastAnc(tree,x,CI=TRUE)
fit
# 比较真实状态与估计状态
plot(a,fit$ace,xlab="true states",ylab="estimated states")
lines(range(c(x,a)),range(c(x,a)),lty="dashed",col="red") ## 1:1 line
b. 离散型数据
data(anoletree)
x <- getStates(anoletree,"tips")
tree <- anoletree
rm(anoletree)
tree
x
plotTree(tree,type="fan",fsize=0.8,ftype="i")
cols <- setNames(palette()[1:length(unique(x))],sort(unique(x)))
tiplabels(pie=to.matrix(x,sort(unique(x))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],y=-max(nodeHeights(tree)),fsize=0.8)
# 重构祖先状态
fitER <- ace(x,tree,model="ER",type="discrete")
fitER
round(fitER$lik.anc,3)
plotTree(tree,type="fan",fsize=0.8,ftype="i")
nodelabels(node=1:tree$Nnode+Ntip(tree),
pie=fitER$lik.anc,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(x,sort(unique(x))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=-max(nodeHeights(tree)),fsize=0.8)
参考文献：
- RevBayes: Bayesian Phylogenetic Inference Using Graphical Models and an Interactive Model-Specification Language.
- A Fast Likelihood Method to Reconstruct and Visualize Ancestral Scenarios.
- PastView: a user-friendly interface to explore evolutionary scenarios.
- Wikipedia：Ancestral reconstruction
网友评论