计算分子进化-搞懂PAML的正选择分析

xuzhougeng 2023-03-28 原文

许多基因组的文章里都会提到使用PAML进行基因的正选择分析(positive selection), 网上也有一些教程介绍如何用PAML进行分析。无论是文章，还是教程，大多只介绍了过程，读完之后，是能够做相应的分析了，但却不知道为什么要这样子做，这篇教程就做这一方面的补充。

我们应该知道组成蛋白序列的氨基酸对应的核苷酸突变可以分为两类，同义置换(synonymous substitution)和非同义置换(nonsynonymous substitution),通过计算非同一置换速率和同义置换速率的比值, omega=dN/dS, 我们可以衡量蛋白的选择压力。如果选择不影响物种适应环境，那么比值两者的速率应该相等，因此比值为1，如果非同义突变会降低物种的适应性，那么dN < dS, 因此比值小于1，如果非同义突变让提高物种的适应性，那么dN > dS, 则比值大于1. 于是乎，非同义突变率显著高于同义突变即为蛋白质适应性进化的证据。

接下来，我们以书中被子植物光敏色素(phy)的适应性进化作为例子结合理论来讲解，数据方式如下

# PAML的安装和配置不再赘述
# codeml的配置文件
wget http://abacus.gene.ucl.ac.uk/ziheng/data/phyACF.codeml.ctl 
# 15个物种的codon联配结果, 包含gap
wget http://abacus.gene.ucl.ac.uk/ziheng/data/phyACF.txt
# 15个物种的系统发育树
wget http://abacus.gene.ucl.ac.uk/ziheng/data/phyACF.trees

下载的配置文件phyACF.codeml.ctl信息显示如下

      seqfile = phyACF.txt
     treefile = phyACF.trees

      outfile = mlc   * main result file name
        noisy = 3 * 0,1,2,3,9: how much rubbish on the screen
      verbose = 0  * 1: detailed output, 0: concise output
      runmode = 0

      seqtype = 1  * 1:codons; 2:AAs; 3:codons-->AAs
    CodonFreq = 2  * 0:1/61 each, 1:F1X4, 2:F3X4, 3:Fcodon
       aaDist = 0  * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a

        model = 2
      NSsites = 2

        icode = 0  * 0:universal code; 1:mammalian mt; 2-10:see below

    fix_kappa = 0  * 1: kappa fixed, 0: kappa to be estimated
        kappa = 5  * initial or fixed kappa
    fix_omega = 0  * 1: omega or omega_1 fixed, 0: estimate
        omega = 0.1

    fix_alpha = 1  * 0: estimate gamma shape parameter; 1: fix it at alpha
        alpha = .0  * initial or fixed alpha, 0:infinity (constant rate)
        ncatG = 3  * # of categories in dG of NSsites models

        clock = 0   * 0:no clock, 1:global clock; 2:local clock; 3:TipDate
        getSE = 0  * 0: don't want them, 1: want S.E.s of estimates
 RateAncestor = 1  * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)

   Small_Diff = .5e-6
*    cleandata = 1  * remove sites with ambiguity data (1:yes, 0:no)?
*        ndata = 3
       method = 0  * 0: simultaneous; 1: one branch at a time

虽然参数很多，但正选择分析上, 我们所需要修改的参数是 model, NSsites, fix_omega, omega 这四项，来决定codeml的分析模式。

首先，修改配置文件中的 model = 0 和 NSsites = 0. 此时, codeml会计算全局omega.

使用codeml phyACF.codeml.ctl 运行，输出结果在当前目录的mlc文件中。里面信息很多，我们重点关注如下几项

...
lnL(ntime: 27  np: 29): -29984.121043      +0.000000
...
Detailed output identifying parameters

kappa (ts/tv) =  1.98351

omega (dN/dS) =  0.08975

lnL是似然值(likeilhood value)的自然对数，之所以是负数，是因为计算出似然值是一个非常小的小数，如果不取对数，结果显示就是0，难以使用。

从结果来看，我们算出的全局omega非常小，约为0.09。这很容易理解，因为我们都知道一个蛋白序列，保守的位点肯定远远多于不保守的位点，那么平均下来，整体的值就会很小。因此正选择通常分析的是系统发育关系的特定谱系或者是蛋白质的某几个位点。

接下来，我们修改 model = 2 和 NSsites = 0, 此时codeml会分析我们提供系统发育树中某个分支(foreground)相对于其他分支(background)是否处于正选择。问题来了，codeml如何判断哪个是foreground，哪些是background呢？此时需要看下 phyACF.trees.

1


((C.Sorg:1.414715, (F.Tom:1.174355, C.Arab:1.734907):0.401510):7.949045 #1, (Oat3:0.217161, (A.Rice:0.255094, (A.Zea:0.084488, A.Sorg:0.041302):0.239315):0.038530):1.329584, ((A1.Pea1:0.304507, A.Soy:0.370169):0.329161, ((A.Pars:0.912283, (A.Tob:0.154040, (A.Tom:0.051204, A.Pot:0.054514):0.100673):0.456387):0.182221, (A.Zuc:0.729785, A.Arab:0.792669):0.136483):0.130833):0.547791);

我们可以使用iTOL这个网页工具对这个树进行展示，树形如下。

iTOL tree

不难从图中发现，有一个明显和其他格格不入的分支，其中包含C. Sorg(高粱,PhyC)和 C.Arab(拟南芥,PhyC), F.Tom(番茄,PhyF)。我们想要检验这一支，是不是受到了正选择。为了强调这一支，我们需要在树结构对应的位置上加上#1，表示foreground, 而这在我们下载的.tree文件中已经有了，你可以检查下。

我们运行codeml phyACF.codeml.ctl，关注输出文件mlc里的如下内容

...
lnL(ntime: 27  np: 30): -29983.513876      +0.000000
...
Detailed output identifying parameters

kappa (ts/tv) =  1.99551

w (dN/dS) for branches:  0.08998 0.03881

在dN/dS这一行有两个omega值，第一个是background, 第二个是foreground。虽然后者比前者大，但是也没有超过1。另外，似然比检验(LRT)也显示，这两个模式没有明显差异。

这里，我们提到了似然比检验(likelihood-ratio test, LRT), 这个概念很重要，我们需要稍稍展开说明下。它指的是，根据两个竞争的统计模型的似然值的比值，评估两者的拟合度，其中一个是最大化整个参数空间，另一个则是做一些限制。如果限制条件（零假设）被观测数据所支持，那么两者的似然值的差异不会超过抽样误差。因此，LRT检验的是，比值是不是和1有显著区别，或者说比值的自然对数和0有显著区别。

通常似然比检验的统计量表现为两个对数似然值的差值

$\lambda_{\mathrm{LR}}=-2\left[\ell\left(\theta_{0}\right)-\ell(\hat{\theta})\right]$

这个统计量，我们可以使用卡方检验(chi2)来分析它的显著性。

# 用R算出
abs(-2*(-29983.513876 - -29984.121043))
# PAML计算chi2, 自由度设置见最后的补充
chi2 1 1.27
df =  1  prob = 0.270475480 = 2.705e-01

对比我们model=2和model=0, 显著性0.27 > 0.05, 不足以拒绝零假设，即我们检测分支的omega相对于全局的omege没有明显差异。

接下来，我们来介绍在所有教程都会提到的branch-site model. 也就是设置model=2 NSsites=2. 它会分析目标分支里的位点是不是受到了正选择。

我们需要建立两个假设，分别是零假设和备择假设

零假设: 检验的分支里的位点不受选择，我们设置参数fix_omega=1, omega = 1
备择假设: 检验的分支里的位点受到正选择，我们设置fix_omege=0, omege = 1.1

在零假设时，输出结果的内容如下，记录lnL= -29704.738847

lnL(ntime: 27  np: 31): -29704.738847      +0.000000
...
Detailed output identifying parameters

kappa (ts/tv) =  2.16177

MLEs of dN/dS (w) for site classes (K=4)

site class             0        1       2a       2b
proportion       0.77433  0.07298  0.13953  0.01315
background w     0.07767  1.00000  0.07767  1.00000
foreground w     0.07767  1.00000  1.00000  1.00000
...

在备择假设时，输出结果的内容如下，记录lnL= -29694.784206

lnL(ntime: 27  np: 31): -29694.784206      +0.000000

...
Detailed output identifying parameters

kappa (ts/tv) =  2.18201

MLEs of dN/dS (w) for site classes (K=4)

site class             0        1       2a       2b
proportion       0.81323  0.07539  0.10194  0.00945
background w     0.07958  1.00000  0.07958  1.00000
foreground w     0.07958  1.00000 17.59530 17.59530

计算统计量

abs(-2 *(-29694.784206 - -29704.738847 ) )
# 19.90928

卡方检验(关于自由度是1, 见最后的补充)

chi2 1 19.9
df =  1  prob = 0.000008160 = 8.160e-06

p远远小于0.01, 为正选择提供了强有力的证据。此时检查mlc输出文件的如下内容，我们还能够确定被选择的位点有哪些。 *表示显著性。

Bayes Empirical Bayes (BEB) analysis (Yang, Wong & Nielsen 2005. Mol. Biol. Evol. 22:1107-1118)
Positive sites for foreground lineages Prob(w>1):
    17 T 0.539
    33 G 0.634
    34 D 0.672
    35 S 0.750
    43 E 0.621
    55 R 0.982*
    61 I 0.643
    66 H 0.924
    71 K 0.701
   102 T 0.974*
   104 V 0.932
   105 S 0.994**
   115 D 0.875
   117 P 0.990*
   120 G 0.819
   130 T 0.976*
....

最后补充下为什么自由度是1, 一方面这是PAML的文档中提到的。。

image.png

另一方面，根据定义，自由度指的是统计量计算时不受限制的变量个数。在我们检验分支的时候，只有一个foreground w是需要自由，所以df=1; 在我们对比branch-sites model时，备择假设相对零假设，只有一个w2是自由(fix_omega = 0), 所以df=1

以上是我在学习正选择分析时的整理结果，由于数学功底太弱，有些疑问我还没有结果，比如说似然值为什么那么小？ML模型参数是如何确定的？输出结果中Bayes Empirical Bayes是怎么运算的？这些还需要不断的学习。

正选进化 code section

有关计算分子进化-搞懂PAML的正选择分析的更多相关文章

ruby-on-rails - 使用一系列等级计算字母等级 - 2
这里是Ruby新手。完成一些练习后碰壁了。练习:计算一系列成绩的字母等级创建一个方法get_grade来接受测试分数数组。数组中的每个分数应介于0和100之间，其中100是最大分数。计算平均分并将字母等级作为字符串返回，即“A”、“B”、“C”、“D”、“E”或“F”。我一直返回错误:avg.rb:1:syntaxerror,unexpectedtLBRACK,expecting')'defget_grade([100,90,80])^avg.rb:1:syntaxerror,unexpected')',expecting$end这是我目前所拥有的。我想坚持使用下面的方法或.join，
ruby - Rails 3 的 RGB 颜色选择器 - 2
状态:我正在构建一个应用程序，其中需要一个可供用户选择颜色的字段，该字段将包含RGB颜色代码字符串。我已经测试了一个看起来很漂亮但效果不佳的。它是“挑剔的颜色”，并托管在此存储库中:https://github.com/Astorsoft/picky-color.在这里我打开一个关于它的一些问题的问题。问题:请建议我在Rails3应用程序中使用一些颜色选择器。最佳答案也许页面上的列表jQueryUIDevelopment:ColorPicker为您提供开箱即用的产品。原因是jQuery现在包含在Rails3应用程序中，因此使用基
报告回顾丨模型进化狂飙，DetectGPT能否识别最新模型生成结果？ - 2
导读语言模型给我们的生产生活带来了极大便利，但同时不少人也利用他们从事作弊工作。如何规避这些难辨真伪的文字所产生的负面影响也成为一大难题。在3月9日智源Live第33期活动「DetectGPT：判断文本是否为机器生成的工具」中，主讲人Eric为我们讲解了DetectGPT工作背后的思路——一种基于概率曲率检测的用于检测模型生成文本的工具，它可以帮助我们更好地分辨文章的来源和可信度，对保护信息真实、防止欺诈等方面具有重要意义。本次报告主要围绕其功能，实现和效果等展开。（文末点击“阅读原文”，查看活动回放。）Ericmitchell斯坦福大学计算机系四年级博士生，由ChelseaFinn和Chri
计算机毕业设计ssm+vue基本微信小程序的小学生兴趣延时班预约小程序 - 2
项目介绍随着我国经济迅速发展,人们对手机的需求越来越大,各种手机软件也都在被广泛应用,但是对于手机进行数据信息管理,对于手机的各种软件也是备受用户的喜爱小学生兴趣延时班预约小程序的设计与开发被用户普遍使用,为方便用户能够可以随时进行小学生兴趣延时班预约小程序的设计与开发的数据信息管理,特开发了小程序的设计与开发的管理系统。小学生兴趣延时班预约小程序的设计与开发的开发利用现有的成熟技术参考,以源代码为模板,分析功能调整与小学生兴趣延时班预约小程序的设计与开发的实际需求相结合,讨论了小学生兴趣延时班预约小程序的设计与开发的使用。开发环境开发说明：前端使用微信微信小程序开发工具：后端使用ssm：VU
ruby - 如何计算 Liquid 中的变量 +1 - 2
我对如何计算通过{%assignvar=0%}赋值的变量加一完全感到困惑。这应该是最简单的任务。到目前为止，这是我尝试过的:{%assignamount=0%}{%forvariantinproduct.variants%}{%assignamount=amount+1%}{%endfor%}Amount:{{amount}}结果总是0。也许我忽略了一些明显的东西。也许有更好的方法。我想要存档的只是获取运行的迭代次数。最佳答案因为{{incrementamount}}将输出您的变量值并且不会影响{%assign%}定义的变量，我
ruby - 使用 Ruby，计算 n x m 数组的每一列中有多少个 true 的简单方法是什么？ - 2
给定一个nxmbool数组:[[true,true,false],[false,true,true],[false,true,true]]有什么简单的方法可以返回“该列中有多少个true？”结果应该是[1,3,2] 最佳答案使用转置得到一个数组，其中每个子数组代表一列，然后将每一列映射到其中的true数:arr.transpose.map{|subarr|subarr.count(true)}这是一个带有inject的版本，应该在1.8.6上运行，没有任何依赖:arr.transpose.map{|subarr|subarr.in
ruby - 我正在学习编程并选择了 Ruby。我应该升级到 Ruby 1.9 吗？ - 2
我完全不是程序员，正在学习使用Ruby和Rails框架进行编程。我目前正在使用Ruby1.8.7和Rails3.0.3，但我想知道我是否应该升级到Ruby1.9，因为我真的没有任何升级的“遗留”成本。缺点是什么？我是否会遇到与普通gem的兼容性问题，或者甚至其他我不太了解甚至无法预料的问题？最佳答案你应该升级。不要坚持从1.8.7开始。如果您发现不支持1.9.2的gem，请避免使用它们(因为它们很可能不被维护)。如果您对gem是否兼容1.9.2有任何疑问，您可以在以下位置查看:http://www.railsplugins.or
ruby-on-rails - Rails 单选按钮 - 模型中多列的一种选择 - 2
我希望用户从一个模型的三个选项中选择一个。即我有一个模型视频，可以被评为正面/负面/未知目前我有三列bool值(pos/neg/unknown)。这是处理这种情况的最佳方式吗？为此，表单应该是什么样的？目前我有类似的东西但显然它允许多项选择，而我试图将它限制为只有一个..怎么办？最佳答案如果要使用字符串列，让我们说rating。然后在你的表单中:#...#...它只允许一个选择编辑完全相同但使用radio_button_tag: 关于ruby-on-rails-Rails单选按钮-模
ruby-on-rails - CarrierWave - PDF - 只选择第一页 - 2
我的Rails应用程序中安装了carrierwave。但是，当用户上传多页pdf时，我只希望应用程序获取文档中的第一页并将其转换为jpeg。这可能吗？用什么命令？这是我的uploader。#encoding:utf-8classImageUploader[200,300]##defscale(width,height)##dosomething#end#Createdifferentversionsofyouruploadedfiles:version:thumbdoprocess:resize_to_fill=>[150,210]process:convert=>:jpgdefful
arrays - 计算数组中的匹配元素 - 2
给定两个大小相等的数组，如何找到不考虑位置的匹配元素的数量？例如:[0,0,5]和[0,5,5]将返回2的匹配项，因为有一个0和一个5共同；[1,0,0,3]和[0,0,1,4]将返回3的匹配项，因为0有两场，1有一场；[1,2,2,3]和[1,2,3,4]将返回3的匹配项。我尝试了很多想法，但它们都变得相当粗糙和令人费解。我猜想有一些不错的Ruby习惯用法，或者可能是一个正则表达式，可以很好地回答这个解决方案。最佳答案您可以使用count完成它:a.count{|e|index=b.index(e)andb.delete_at

计算分子进化-搞懂PAML的正选择分析

有关计算分子进化-搞懂PAML的正选择分析的更多相关文章

随机推荐