系统发育比较分析—R

Python020

系统发育比较分析—R,第1张

系统发育树 是研究物种进化历史必不可少的信息,我们可以利用它得到一些重要历史线索,如:

首先,安装系统发育分析所需的软件包

其实,此处的树文件就是一个字符串列表(列表还可以是数字型)。

接下来,主要是看一下这些对象是如何存储在变量中的:

我们可以看到,树中所有分类单元之间的关系信息都包含在每条边的起始节点和结束节点中,共享共同起始节点编号的边是直接共同祖先的后代。

我们还可以测试树文件,编辑以及去除一些末端枝:

布朗模型:

(1). 连续型性状进化的一个模型;

(2). 性状随着时间不断变化;

(3). 经过一段时间后,预期的状态将服从正态分布;

(1). 有时称为维纳过程;

(2). 连续时间随机过程;

(3). 描述连续型性状的“随机演化”;

(1). BM可以用来描述由大量独立的弱力组合而产生的运动;

(2). 加入许多小的自变量后,无论原始分布如何(中心极限定理),都会产生正态分布;

演化过程近似布朗运动

Ⅰ. 遗传漂变

Ⅱ. 随机改变

Ⅲ. 相对于考虑的时间间隔弱的自然选择

Ⅳ. 随时间随机变化的自然选择

系统发育独立差(Phylogenetic Independent Contrast, PIC)是去除性状分析中物种系统发育关系的一种方法,是美国进化生物学家 J. Felsenstein于1985年提出的。 参考链接

当时,有动物学家检验了若干哺乳动物脑容量和体重的相关性。然而,所用统计方法都假设各个种的性状是相互独立的,采集自同一个正态分布的总体。由于不同种之间存在系统发育上的联系,因此不能直接套用常规统计模型,所以所得结果也就不可信了。在Felsenstein之前,一些学者已经意识到这个问题,但是都没能给出理想的解决方案。

Felsenstain认真思考了这个问题,并且认为:研究多个不同性状的相关性时,必须考虑系统发育关系。他提出:假设物种性状的进化服从布朗运动,那么用系统发育关系上最近的分类单元性状的差值,再经过枝长加权,就可以得到一组新的数据。这组新的数据去掉了进化信息的干扰,因此可以用常规统计方法检验。上述方法新获得的数据就称为系统发育独立差 (The Phylogenetic Independent Contrast, PIC)。

举例来说,如果要研究哺乳动物脑容量和体重的关系,就要先获得这些种之间的进化关系,也就是经过分子钟校订的进化树。然后再计算脑容量的PIC和体重的PIC,再用脑容量的PIC和体重的PIC建立线性模型,进行相关性检验。

现在让我们来模拟一下布朗运动在树枝上的移动。坚持离散时间,我们首先需要一个离散时间的系统演化,我们可以用pbtree在phytools中进行模拟。

现在,为了在树的所有分支上进行模拟,我们只需要分别在每个分支上模拟,然后根据每个分支的父分支的结束状态“移动”每个子分支。我们可以这样做,记住,因为布朗演化的结果在每个时间步骤是独立于所有其他的时间步骤。

在本教程中,我们将学习应用Felsenstein的独立对比方法来估计性状之间的进化相关性。Felsenstein发现一个以前已经认识到的问题,但这一问题没有得到充分的认识:即从统计分析的角度来看,物种数据不能被视为独立的数据点。

为了学习对比方法,我们首先需要学习一些用r语言拟合线性模型的基础知识。

这个模拟非常简单地表明,当满足标准统计假设时,我们的参数估计是无偏的,并且第一类误差在标称水平或附近。

系统发育数据的困难在于,对于y的观测不再是独立的和相同的分布:

从这个例子中我们可以看出,对于系统发育来说,诱发I类型的错误并不难。这是因为紧密相关的分类群具有高度相似的表型。换句话说,它们并不是关于树上x和y的演化过程的独立数据点.

PGLS还可以使用多个预测变量:

另一种检测系统发育信号的方法是Pagel的lambda。Lambda是一种相对于内部分支延伸尖端分支的树变换,使树越来越像一个完整的多段切除法。如果我们估计的λ=0,则推断这些性状没有系统发育信号。λ=1对应于布朗运动模型,0<lambda<1介于两者之间。

推断祖先状态一直是系统发育比较生物学中一个重要的目标,本教程主要介绍如何在BM以及OU模型下重构祖先性状。第一个例子就是蜥蜴祖先体型大小的状态重建,这是一个连续型变量。

现在我们可以估计祖先的状态,还将计算方差&每个节点的95%置信区间:

写在最前面:

首先需要说一下,本文的bootstrap和jackknife都算是蒙特卡罗方法(Monte Carlo method)的一种。应用广泛的的MCMC链(马尔可夫链蒙特卡洛方法Markov chain Monte Carlo)也是蒙特卡罗与马尔可夫链的结合。简单来说,蒙特卡罗方法就是从已知样本的分布中随机抽取新的样本集进行评估,然后放回,再次抽取的方法。根据具体方法的不同,抽取样本集的手段也不同。

bootstrap抽样方法将观测到的样本视为一个有限的总体,是唯一的信息来源,从中有放回的随机抽样来评估总体特征,以及对抽样总体进行推断统计。bootstrap 也分参数bootstrap和非参数bootstrap,前者的分布已完全知道。但在生信领域一般没有这种情况。所以下面讨论的是非参数bootstrap。

直接上例子:

假设现在有bootstrap包中的law数据集如下,

现在我们要计算LSAT成绩(美国法学入学考试)和GPA之间的相关系数。但因为样本量太少了,所以我们使用bootstrap重复抽样评估其标准误。

200次循环抽样后,计算得se.R标准误为0.1474629

得到如下的图:

1e6次循环抽样后,计算得se.R标准误为0.1333802

得到如下的图:

如果用bootstrap包的bootstrap函数会快一些:

bootstrap函数的用法: bootstrap(抽取样本范围,重复次数,进行bootstrap的函数,bootstrap的数据集)

偏差定义为bootstrap结果(多个数值)与原数据统计结果(单个数值)的均值:

得到bias大约为0.001817608,比较小

换一个包,boot包

这里用了三种方法计算置信区间:basic、正态和百分数。样本相关系数分布接近正态,则正态置信区间接近百分数区间。此外还有“Better Bootstrap Confivendence Interval” 更好的bootstrap置信区间,称为BCa区间,使用偏差和偏度对百分数置信区间进行矫正。设置type="bca"即可。

简单的说,bootstrap是从原有真实样本中有放回地抽取n个。jacknife就是每次都抽取n-1个样本,也就是每次只剔除一个原样本。

同样地,如果以bootstrap包中的law数据进行演示:

Jackknife计算的bias为-0.006473623。 这里jackknife的偏差公式相比于bootstrap有一个(n-1)系数,推导就不写了。

标准误se为0.1425186,与bootstrap得出的比较接近。

当统计量不太平滑的时候,Jacknife有很大误差。比如说对中位数进行统计,其变化很大。在进行Jacknife之后最好再跑一次bootstrap,看看是否相差很大。

居然还能这么嵌套着玩,针对每次bootstrap形成的数列向量计算jackknife的标准差,这样可以看出bootstrap若干次取样之间的差异。

算出来分别为0.1344824和0.08545141。后者较小,表面bootstrap取样之间的variance较小。

简单来说就是一种数据分割检验的方法,将数据分割为K份,称为"K-fold"交叉检验,每次第i个子集作为测试集来评估模型,其余的用来构建模型。Admixture使用的就是这个原理。Jackknife也属于Cross Validation的应用之一。

现在我创建一个这样的alignment:

这棵树长这样,符合遗传距离:

进行bootstrap:

phylogeny的bootstrap是对每一个节点都进行bootstrap取样并建树,比如说在9号节点,查看其bootstrap子集建的树符合系统发育关系((human2,human4,human3)(human8,human1,human6,human7,human5))的百分比(不管内部怎么样,先看这个节点)。发现Node1支持率是100(1000次都符合)。而后移到下一个节点,并且只看节点内部的分支支持率是多少。

其实原理都比较简单,计算bootstrap也会有专门的软件。

参考资料:

1)中科大张伟平教授课件

2) https://ecomorph.wordpress.com/2014/10/09/phylogenetic-trees-in-r-4/