1. 选择分析
1.1. omega(ω)检测选择
dNdS(KAKS),以及 omega(ω) 的定义,以及 ω 用于检测选择压可以参考博客一文说清KAKS、dNdS、DnDs,以及ω:https://yanzhongsino.github.io/2023/10/09/evolution_selection_dNdS_intro
- omega(ω) 的含义
- omega(ω) = dN/dS(KA/KS)
- omega(ω) 可用于推断作用于蛋白质编码基因的自然选择的方向和程度
- 当 ω > 1 时,意味着正选择 (positive selection);
- 当 0 < ω < 1 意味着负选择 (negative selection,也叫净化选择或纯化选择 purifying selection);
- 当 ω = 1 则表示中性进化(neutral selection),即不受选择。
- 背景知识
- 通常,由于未改变编码氨基酸,大多数同义突变都是中性的,即没有受到正选择或者负选择;由于改变了编码氨基酸,大多数非同义突变受到负选择,没有保留下来。
- 所以,在现存类群中,观察到的同义突变数量比非同义突变数量会高得多,通常dN < dS。计算现存类群间 ω 时,结果通常在0.05-0.40之间,只有很少数的基因(e.g. MHC)发现 ω > 1。
- 可能的误差
- 在基因的不同位点或进化的不同时期,正选择和负选择的组合可能会相互抵消。
- 由此产生的平均值可能会掩盖其中一种选择的存在,并降低另一种选择的似然程度。
- 统计分析
- 有必要进行统计分析,以确定 ω 的结果是否与 1 有显著差异,或是否由于数据集有限而出现任何明显差异。
- 近似法的适当统计检验包括用正态近似值对 dN - dS 进行近似,并确定 0 是否位于近似值的中心区域。
- 可以使用更复杂的似然法来分析最大似然分析的结果,通过进行卡方检验来区分空模型和观测结果。
1.2. 选择分析的目的
- 绝大多数情况下 ω < 1,因此选择分析通常是为了检测到较少见的正选择(ω >1)
- 在一个基因的所有位点或系统发育上的所有谱系(分支)上的平均 ω 值很少大于1。因此,用平均 ω 值检测正选择是非常保守的(很可能检测不出来)。
- 检测系统发育树上的特定分支或基因的单个位点是否经历正选择更有价值。所以选择分析的重点通常是找到系统发育树上受到正选择的分支,以及基因上受到正选择的位点。
- 另外,结合系统发育树,还可以分析 ω 在进化树上的变化,虽然进化树上的 ω 都小于1,但某一支比其他所有支的 ω 都要小,这也说明这一支受到负选择的放松。
1.3. 选择分析的情况
- 基因在两个物种的选择分析,较简单。
- 直接比较这两个物种的密码子序列,计算dN/dS,即 ω 值,通过 ω 值判断选择的方向和程度。若 ω > 1,即表明该基因在物种进化过程中,即由其祖先物种分化成这两个物种时,基因受到了正选择。
- 基因在多物种的选择分析,如果仍然按照两个物种的方法,结果可能不好解释。
- 因为该基因可能在某一类群中序列很相似,其两两比较时,ω <= 1;而在另外一类群中两两比较时,很多时候 ω > 1。最后软件可以从总体上给一个ω值,但该值不可以拿来简单地评价该基因是否受到了正选择。
- 所以,对多个物种进行正选择分析时,没法直接评价该基因是否受到了正选择。正选择只有在进行两两序列比较的时候,才能计算 ω 值,从而得到结果。
- 基因在多物种的选择分析,目的则是:比较某个分枝上祖先节点和后裔节点(可以理解成,对无根树上某分枝两侧的两组物种进行比较,依然属于两两比较),从而计算该分枝的 ω 值。
- 在实际数据中,基因在不同的进化分枝上具有不同的 ω 值,在序列不同的位点也具有不同的 ω 值
- 可以同时分析目标分枝的 ω 值和序列位点的 ω 值,从而判断哪一支的哪个基因受到正选择。
- 目标分枝两侧的物种数量较多时,可以对序列上的每个位点进行 ω 值分析,从而鉴定出正选择位点。
1.4. 选择分析的软件
基于计算 ω 值 (dN/dS) 来进行选择分析的软件有很多,最经典的是PAML的codeML,本文接下来便介绍这个程序。
PAML(Phylogenetic Analysis by Maximum Likelihood)是Yang Zihen lab开发的使用最大似然法进行DNA和蛋白序列的系统发育分析的程序包,包括许多程序,其中codeML程序整合了许多选择分析的模型,用于系统发育树上的选择分析。
2. 使用codeML进行选择分析的流程
- 配置零假设模型(null model)和备选模型(正选择模型)的参数(配置文件codeml.ctl);
- 运行codeML程序分别对两个模型进行分析,获得各自的似然值 (lnL);
- 通过统计检验(如卡方检验)比较两个似然值(lnL),自由度为两个模型之间自由参数数量的差异,并计算p值,用来判断两个模型间是否存在显著差异;
- 根据统计检验结果解读,是同意备选模型 (存在正选择) 还是拒绝备选模型。
- 如果是位点模型或分支位点模型,还可以进一步根据BEB方法判断受到正选择的位点的位置。
3. codeML进行选择分析的四种模型
四种模型的主要区别在于假设的不同,是否允许位点间的 ω 值和分支间的 ω 值不同,是基于先验知识选择的模型。
3.1. 位点模型 (site model)
3.1.1. 位点模型介绍
- 位点模型的作用
- 主要用于检测基因中的正选择位点。
- 位点模型的主要假设
- 主要假设数据集中不同密码子位点受的选择压是来自统计分布的随机变量,因此允许密码子位点之间的 ω 不同。
- 不考虑不同支系间受的选择压力差异,假设进化树中各分枝的 ω 值是一致的。
- 所以,使用位点模型能在整体水平上检测基因的正选择位点,而不能表明基因在某个进化分枝上是否受到正选择压。
- 正选择定义为存在某些密码子的 ω >1. 执行似然比检验(LRT)与null模型的一般密码子比较,不允许任何密码子ω>1。
- 位点模型使用的模型
- M0(单一比率),即:One-ratio model,假设所有位点具有相同的 ω 值;
model = 0, NSites = 0
- M1a(近中性),假设仅有保守位点(0<ω <1)和中性位点( ω =1)而没有正选择位点 (ω >1)存在,这两类位点的比率分别为p0和p1,其对应的ω 值分别为ω0、ω1;
model = 0, NSites = 1
- M2a(正选择),该模型在M1基础上增加了第三类ω值,即假设除了保守位点和中性位点外,还存在处于正选择压力下的位点 (ω >1),这三类位点的比率分别为p0、p1和p2,其对应的ω 值分别为ω0、ω1和ω2;
model = 0, NSites = 2
- M3(离散),假设所有的位点ω 值呈简单的离散分布趋势;
model = 0, NSites = 3
- M7(beta),假设所有位点的 ω 属于矩阵(0, 1)并呈beta分布;
model = 0, NSites = 7
- M8(beta & ω ) ,该模型在M7基础上增加另一类ω 值(ω >1);
model = 0, NSites = 8
- M8a(beta & ω =1),与M8模型类似,但将ω 值固定为1(ω =0);
model = 0, NSites = 8
- 常选的成对模型(前两对被广泛使用)
- M1a(选择约束放松)vs. M2a(正选择)
- M7 (beta) (选择约束放松)vs. M8 (beta&ω) (正选择)
- M3(分散比率) vs. M0(负选择)
- M8* (正选择)vs. M8a(选择约束放松)
3.1.2. 位点模型的分析过程
- 选择零假设模型(null model)和备选模型(正选择模型)两种模型
- 有两对位点模型特别有效,M1a(选择约束放松)vs. M2a(正选择),以及 M7 (beta) (选择约束放松)vs. M8 (beta&ω) (正选择):
- (1) 模型M1a是null model,认为所有位点的 ω 值 < 1 或 = 1 两类;
- (2) 模型M2a是正选择模型,存在 ω <1、=1或> 1的位点。
- 使用codeML分别分析两个模型,获得各自的似然值 (lnL)
- 通过统计检验(如卡方检验)比较两个似然值(lnL),并计算p值,用来判断两个模型间是否存在显著差异
- 似然率检验 (LRT)统计量、或两倍于两个比较模型之间的对数似然差(2Δ),可用于卡方检验。
- 自由度为两个模型之间自由参数数量的差异。例如,M1a有2个自由参数,M2a有4个,因此自由度为2,M7-M8的比较也有2个自由度。
- 若p值 < 0.05,则否定null model,认为存在正选择位点。
- 同时,推荐采用比较模型M7和M8,其p值结果比上一种比较方法更宽松,能检测到更多的正选择基因。
- M7 模型是一个零假设模型,所有位点在 0 ≤ ω ≤ 1 的区间内,与M1a不同的是把 ω 分为10类。不允许 ω > 1 的站点。
- 备选模型 M8 模型增加了 ω > 1 的第 11 类位点。对每个位点进行测试,以确定其所属类别。
- 如果LRT倾向于M2a或M8,我们可能会问:基因中的哪些位点处于正选择之下?
- 利用贝叶斯定理计算每个 ω > 1的候选正选择位点的后验概率来回答这个问题;
- 在codeML中实现了两种方法:
- (1) Naïve Empirical Bayes (NEB)方法(Nielsen and Yang 1998)通过使用模型中参数的最大似然估计(MLEs)计算每个位点来自不同位点类别的后验概率,而不考虑其抽样误差或不确定性。
- (2) 贝叶斯经验贝叶斯(Bayes Empirical Bayes, BEB)方法是对NEB方法的改进,并适应了MLEs中的不确定性(Yang et al. 2005)。虽然codeML报告两种方法的结果,但应该使用BEB而不是NEB。
3.2. 分支位点模型 (branch-site model)
3.2.1. 分支位点模型介绍
- 分支位点模型的作用
- 主要用于检测基因在某个进化枝上是否存在的正选择位点。
- 分支位点模型的主要假设
- 主要假设不同氨基酸位点的和不同支系间受的选择压力均存在差异(既考虑位点间也考虑支系间的 ω 值存在差异)。
- 接受正选择测试的分支称为前景分支,而树上的所有其他分支都是背景分支。
- 背景分支有两类位点,保守位点0 < ω0 < 1 ,中性位点 ω1 = 1;前景分支,正选择存在的位点 ω2 > 1 。通过比较零假设模型(ω2 = 1)和备择假设模型(ω2 > 1)的差异显著性来判断是否存在正选择。
- 分支位点模型使用的模型和主要参数
- Model A null:
model=2, NSites=2, ncatG=ignored, fix_omega = 1, omega = 1
,设定 ω 为固定值1。 - Model A :
model=2, NSites=2, ncatG=ignored, fix_omega = 0, omega = 1.5
, 估算 ω 是否大于1。 - Model B :
model=2, NSites=3, ncatG=ignored
- Model C :
model=3, NSites=2, ncatG=ignored
- Model D :
model=3, NSites=3, ncatG=2 or 3
- 常选的成对模型
- Model A (正选择) vs. Model A null (中性进化)
3.2.2. 分支位点模型的分析过程
- 选择零假设模型(null model)和备选模型(正选择模型)两种模型
- 正选择模型为Model A,将 ω 值分成<1、=1、>1的三类,这和site model中的一样;
- 零假设模型和正选择类似,只是将 ω 固定成1,作为null model(Model A null)。
- 使用codeML分别分析两个模型,获得各自的似然值 (lnL)
- 比较两种模型的似然差异,利用卡方检验(自由度为2)算p值(chi2命令算出的值除以2)。
- 若p值< 0.05,则否定null model,认为存在正选择位点。
- 与位点模型一样,通过BEB方法识别前景支上可能处于正选择的密码子位点。
- 通过BEB方法计算正选择位点的后验概率,若存在概率值 > 0.95正选择位点,则表示基因在目标分枝上受到正选择。
- notes
- 分支位点模型并不给出分枝上的 ω 值。这表明虽然考虑了目标分枝上具有不同的 ω 值,但仍然以分析位点上的 ω 为主。
- 值得注意的是,在分支位点模型下可能检测到正选择位点,但在目标分枝上的 ω 值仍然可能低于1。可能软件作者基于这点考虑,就没有给出目标分枝上的 ω 值,以免影响一些人对正选择结果的判断。
- 要注意应该预先指定前景支。如果在没有先验生物学假设的情况下使用相同的数据集对树的多个分支进行正选择测试,则可能需要对多个测试进行校正。
- Bonferroni修正可能过于保守,而Rom程序(Rom 1990)的功效略高,是首选(Anisimova and Yang 2007)。人们还可以使用控制假阳性(FDR)的程序,它是所有被拒绝的零假设中真零的预期比例,或者是所有阳性测试结果中假阳性结果的比例(Benjamini和Hochberg 1995,2000)。需要注意的是,如果这些序列非常分散或存在严重的模型违规,则多重测试校正可能是不可靠的。
3.3. 分支模型 (branch model)
- 分支模型的作用
- 主要用于检测在某个分枝上,其 ω 值是否显著高于背景分枝,即基因在目标分枝上进化速度是否加快。
- 分支模型的主要假设
- 分支模型认为基因序列上所有位点的 ω 值是一致的;但系统发育树上不同分支的 ω 值是不同的。
- 分支模型的分析过程
- 对两种模型进行比较:
- (1) 第一种模型为null model,所有分枝具有相同的 ω 值;
- (2) 第二种备择模型认为目标分枝具有一个 ω 值,其它所有分枝具有一个相同的 ω 值。
- 比较两种模型的似然差异,利用卡方检验(自由度为1)算p值。若p值 <= 0.05,且目标分枝上的 ω 值高于背景值,则认为该基因为快速进化基因。
- 一般情况下,该方法计算得到的p值会低于分支位点模型的结果。
- 分支模型使用的模型
- 分支模型主要用于对系统发育树中不同支系 ω 值差异性进行界定,主要有三个模型:
- One-ratio model :
model = 0, NSites = 0
假设系统发育树中所有支系的 ω 值相等; - Free-ratio model :
model = 1, NSites = 0
假设系统发育树中所有支系的 ω 值不相等; - Two-ratio model :
model = 2, NSites = 0
假设前景枝和背景枝的ω 值不同;
- 常选的成对模型
- one ratio(负选择) vs. free ratio(自由比率)
- one ratio(负选择) vs. two ratio (正选择)
- two ratio natural(中性进化) vs. two ratio (正选择)
3.4. 进化支模型 (clade model)
与分支位点模型类型 (branch-site model),能同时检测多个进化支 (clade),该模型并没有将背景支的dN/dS值约束在(0,1)。有CmC和 CmD 两种模型,主要参数如下:
- 进化支模型使用的模型
- M2a_model:
model = 0, NSsites = 22
- Clade Model C (CmC):
model = 3, NSsites = 2, ncatG=2 or 3
- Clade Model D (CmD):
model = 3, NSsites = 3, ncatG=ignored
- 常选的成对模型
- Clade Model C vs. M2a_model
4. 用codeML的模型执行选择分析的常见场景
根据不同的分析目的,可以选择不同的模型进行选择分析,常见场景包括以下几种。以后还可以每种应用场景写一篇博客。
- 计算平均选择压 (平均 ω 值 )
- 在所有位点和分支均为 1 个相同的 ω 的模型(m0)下,计算 ω 作为对基因的平均选择压力的度量。
- 用位点模型 (site model) 检测编码序列中的哪些位点处于正选择
- 用分支模型(branch model) 检测系统发育树上一个或多个分支是否处于正选择或负选择是否放松
- 用分支位点模型 (branch-site model) 检测特定分支发生正选择的位点集。
- codeML的最新实现,可以使用由数千个基因比对组成的基因组数据集进行正选择检测。
5. 用codeML进行选择分析的实操
linux系统下使用PAML的codeML,windows系统下可以使用PAML的图形界面版本PAMLX,或者基于PAML的可视化软件EasyCodeML。
5.1. 安装PAML
5.1.1. 方法一:conda安装【方便快捷】
conda install -c bioconda paml
5.1.2. 方法二:下载github库并编译【开发者手册推荐】
- 下载最新版本的PAML:
git clone git@github.com:abacus-gene/paml.git
- 用git clone是复制了PAML文件夹,进入paml/src目录,编译
1 | Move to the cloned repository. Make sure you have cloned it in the correct location of your local disk. |
- 可以把包含所有执行文件的paml/bin/目录放到环境变量中
5.2. 运行codeML
5.2.1. 必需的输入文件
输入的序列比对文件和树文件的格式可以参考https://github.com/abacus-gene/paml-tutorial/tree/main/positive-selection/00_data 中的 Mx_aln.phy Mx_root.tree,Mx_unroot.tre。
- 序列比对文件 (cds_aln.phy)
- phylip格式,加上样品数量和序列长度两个数字组成的首行。
- 比对好的编码序列 (cds) 文件,碱基数量是3的倍数(codon模式比对)
- 建议删除gaps和难以align的区域,但要按密码子删除(删除和保留都是3的倍数个碱基)
- 树文件 (tree.newick)
- newick格式,并在首行前加上一行,包括两个数字,一个是物种/样本数量,一个是树的数量(一般是1),空格隔开。
- 只需要拓扑结构,删除枝长、支持率等信息。
- 如果包含多棵树,需要一个头文件。
- 配置文件 (codeml.ctl)
- codeml.ctl 的例子在安装的程序包paml/examples/下有,可以复制到分析目录下,修改后使用
- codeml.ctl 的参数需要根据使用的模式,运行的目的等具体情况进行修改
5.2.2. 配置文件的参数解释
1 | * 输入输出参数: |
5.2.3. 运行命令
codeml codeml.ctl
5.2.4. 结果文件
- 主要结果文件在codeml.ctl设置的outfile中,例子中默认是mlc。
- 主要关注mlc文件中的两个信息
- lnL是似然值(likeilhood value)的自然对数,之所以是负数,是因为计算出似然值是一个非常小的小数,如果不取对数,结果显示就是0,难以使用。
- omega (dN/dS)
1 | lnL(ntime: 27 np: 29): -29984.121043 +0.000000 |
- 似然比检验(likelihood-ratio test, LRT)
- 通常在零假设模型和备择假设模型下各自得到lnL值和omega值,还需要对两个似然值进行似然比检验(LRT)。
- LRT是指根据两个竞争的统计模型的似然值的比值,评估两者的拟合度,其中一个是最大化整个参数空间,另一个则是做一些限制。如果限制条件(零假设)被观测数据所支持,那么两者的似然值的差异不会超过抽样误差。
- 因此,LRT检验的是,比值是不是和1有显著区别,或者说比值的自然对数和0有显著区别。
- 使用卡方检验(chi2 test)来分析LRT的统计量是否有显著性,即零假设模型和备择假设模型下的似然值的差异是否显著。
6. 用codeML的模型执行选择分析的常见问题
6.1. 有根树和无根树的选择
- 有根树的祖先节点是二叉树,而无根树为三叉树。
- 零模型在所有测试中都是M0,它假设所有分支都是相同的ω。除检验4中的备择假设外,所有模型都使用无根树,在备择假设中,树的根周围的两个分支被赋予不同的ω比,因此需要有根树。
- 在使用codeml时,如果没有指定有根树参数却使用了有根树作为输入,那么在输出结果中将会得到这样的报错信息:”This is a rooted tree. Please check!”。
- 对于大多数模型,即使使用有根树,其该模型似然值仍然是正确的,但是root周围两个分支的长度不稳定,因为它们的和是估计值。对于其他模型,似然估计和参数估计都是不正确的。因此,分析时确实应该注意到这一信息,并尽可能使用一棵无根树。
6.2. 是否删除多序列比对中的gap与模糊字符
- 在多序列比对过程中,对齐gap是极其困难的,paml软件包无法处理gap。因此,我们可以通过设置cleandata = 1来去除gap;此外,还可以将gap当作为模糊字符进行处理。但是,这都不是最好的解决办法,这两种策略都低估了序列差异。
- 除了一个或两个序列之外,大多数序列都有序列信息的位点也许应该保留,而除了一个或两个序列之外,所有序列都有对齐间隙的位点最好被移除。
- 因此,选择合适的多序列比对软件以及过滤软件,在codeml分析前就去除gap和歧义序列,然后参数cleandata = 0是更好的选择。
6.3. 当前景支的dN/dS值显著大于背景支时,该如何解释
- 如果检测前景支时,其dN/dS > 1时,我们可以认为它受到正选择作用。
- 但是,如果其dN/dS <1但大于背景支时,就不能认为它是受正选择作用的,选择压力约束放松可能是较为合理的解释。
- 此外,在Ohta’s的微有害突变假说下,净化选择在大种群中比在小种群中更有效,因此不同谱系的种群规模的差异提供了另一个相容的假设。如果氨基酸的变化是稍微有害的,我们预计在大群体中它们从群体中移除的速度会比在小群体中更高。
- 因此,即使两系在选择压力或基因功能上没有差异,我们也期望在一个大群体中看到一个较小的dN/dS比值,例如,许多核基因的dN/dS比值在啮齿类动物中低于灵长类或偶蹄类。
6.4. 不同的模型鉴定出不同的位点,该选择相信哪一种模型?
- 通常,如果一个位点在一种模型下出现在选择列表中,那么在另一种模型下也会有相当大的概率。
- 确定位点的问题很困难,而且容易出错。因此,我们通常会认为后验概率大于95%或者99%的位点,是比较可信的。
6.5. 如何标记前景支
- 如果标记某一个节点,则可以使用”#”;如果是标记某一类群(节点和节点后的所有类群),则可以使用”$”。
- 符号#的优先级高于$,树顶端的进化枝标签优先于接近根处的祖先节点的进化枝标签。
- 位点模型(Site models)和free ratio模型分析不需要标记前景支。但其他三个与枝相关的模型均需要提前进行前景枝的标记。枝模型(Branch model)和进化枝模型(Clade model)可以在单次分析中标记多个前景枝(或进化枝),但枝位点模型(Branch-site model)只能一次标记一个前景枝。
6.6. dN/dS代表进化速率而非突变率
- 较高的dN/dS值可以解释为正选择或者快速进化。尽管突变本身也处于选择压力之中(大多数为净化选择),但不可以解释为“基因A增加了突变的选择压力”,因为突变是随机发生的。原则上讲,突变率会同时影响dN、dS,但通常dN/dS不受突变率影响。
- dN/dS是一种进化速率,但它不是突变速率,因为同义和非同义替代率具有不同的选择约束水平。 选择压力测试的基本原理是:它假设同义替换是中性进行,也就是说,它们大多是在遗传漂移下进化的。如果这是真实情况,那么dS可以用作(中性)突变率的替代。然而,非同义替代率总是处于净化选择压力,在正选择下程度较小。
- 因此,dN/dS是中性偏离的度量。所以,dN > dS,既dN/dS > 1,则受到正选择;如果dN小于dS,则dN/dS < 1,则是净化选择。选择压力测试的关键就是它通过特定基因的同义替换的“中性”进化速率归一化非同义替换的速率。
6.7. 基因树or物种树
- 在任何情况下,使用代表真实演化历史的基因树都是最好的。
- 但是,有时可能无法轻易判断是否符合真实演化历史,那么你可以选择物种树作为代替。
- 基因组水平分析,那么推荐使用物种树。你可以可以使用基因树和物种树进行数据稳健性测试。
6.8. 祖先节点有正选择信号而整个枝系却没有检测到如何解释
- 如果以哺乳动物祖先为前景,假设该基因在共同祖先中存在适应性,可能是由于获得了新的功能,但随后该基因在净化选择下得到了保守进化。如果你把整个clade作为前景,那么假设该基因在整个哺乳动物中的所有分支都承受着持续的压力来改变或多样化,如果该基因涉及到防御或免疫,则可能是这种情况。
- 对于到底是检测祖先,还是整个枝系,取决于生物学问题。例如,溶菌酶在所有colobine猴子中都应该具有相同的功能,因此蛋白质在clade内被期望受到选择性约束,但在colobine clade的分支上,该酶显然获得了一个新的功能,其正选择驱动了氨基酸的变化。有了这个假设,你就应该把分支的祖先贴上clade的标签,而不是那些在clade中的分支。
6.9. Clade模型判别背景支是否同样受正选择
- 经常遇到审稿人的审稿意见是这样的:前景分支上正选择的基因的重要支持并不意味着在背景分支上没有正选择,这些基因可能仍在许多(即使不是全部)背景分支上处于正选择状态。
- 为了检验原假设(基因仅在前景支受到正选择)是否正确,可以进一步通过Clade模型检验,因为进化枝模型允许在不限制背景dN/dS小于1的约束下估计前景支与背景支dS/dN的比率。
6.10. 基因受到正选择但无正选择位点
- 一种可能性是,对整个基因有积极选择的证据,但每个单独位点的信息或证据太弱。 可以查看rst文件,该文件具有所有站点的后验概率,以查看是否是这种情况,mlc文件只列出后验概率高于0.5的文件。
6.11. 有正选择位点但是通过序列比对却找不到
- 可能是因为codeml在删除gap或模糊字符的列,之后重新编号位点了(cleandata =1)。
6.12. dN/dS(ω)值过大,结果是否可信?
- 遇到这种极大值的dN/dS时,比如ω = 999,首先要确保你的序列是否正确;其次,该位置的dn,ds是否远小于0.0001,枝长是否过小。
- 显然,高度相似的序列和非常发散的序列都不是信息丰富的,很难指定确切的值。
- 为了避免此类问题的出现,可以先通过M0模型获得枝长,再将具有分支长度的进化树应用到codeml,并在ctl中设置FIX_blength=2。
6.13. Branch-site检测适应性趋同进化
- 如图所示,红色分支代表表型趋同的进化分支,如果你想通过分支位点模型来检测适应性趋同进化,应该将所有红色分支统一设置为前景分支。当然,前提是假设所有前景支均具有相同的位点受到正选择。 对于背景支是否同样存在类似的适应性趋同,可以通过clade进化支模型进行检测。
6.14. 若p2(正选择位点概率)接近于0
- p0/ω0,p1/ω1,p2=(1-p0-p1)/ω2:在正选择分析的备择假设结果文件中,通常会获得三个p值,其中,p0代表处于净化选择下的位点概率;p1代表处于中性进化下的位点概率;p2则代表处于正选择下的位点概率。
6.15. 分支模型检测正选择
- 可以使用two ratio备择假设(fix_omega = 0 omega = 1)与two ratio零假设(fix_omega = 1 omega = 1)进行统计假设检验。
6.16. 选择约束放松
- codeml检测选择约束放松可以通过2步完成:首先,确定dN/dS显著增加的情况(由于正选择或者选择约束放松);然后,过滤掉显著正选择的情况。
6.17. 不同模型下的起始ω
- 针对于CladeC和CladeD模型,一般要设置几个不同的起始ω,测试其lnL值是否稳定(ω=0.001,ω=0.01,ω=0.1,ω=0.5,ω=1,ω=1.5),最终一般会选取lnL较大的值作为最终值。
6.18. 枝长会影响PAML结果
- 正常分析,应该先用M0估计树的枝长,Kappa值,然后用跑出来的树作为初始树,并设置fix_blength = 2。
6.19. 基因家族选择压力分析
- https://groups.google.com/g/pamlsoftware/c/bsfGWYHPCgA
- https://groups.google.com/g/pamlsoftware/c/KqhrPFYnWMI
6.20. CladeC检测正选择
- CladeC经常用于检测,不同分支的分化选择压力,但有时候检测到前景枝dN/dS>1,此时我们需要进一步利用CladeC的零假设进一步测试(fix_omega=1,omega=1),或者利用branch-site进一步验证一致性。
6.21. Free ratio检测结果是否可以作为正选择的证据
- free ratio模型估计通常会造成较大的抽样误差,例如,较短的分支通常会具有较大的dS/dN。所以,一般对于free ratio出现dN/dS > 999或者dN、dS < 0.0005结果不建议采用。
6.22. CladeC模型检测多组支系选择压力差异
- 在dataset中,有三个clade:A, B, C,那么在clade 1和clade 2之间是否存在显著差异?
- 假设,前景支的标注为
((A1,A2)$1 , (B1,B2)$2 , (C1,C2)$0);
- 首先,为了测试clades之间的显着差异,可以比较CladeC与M2a_rel。M2a_rel假设$1、$2、$0等都是在相同的选择压力下进化的,这个测试应该有两个自由度。
- 其次,为了测试clade A和B之间的显着差异,同时允许clade C是不同的,您可以比较使用上面提供的树运行CMC与使用更简单的树运行CMC的契合程度。在这种情况下,更简单的树会将clade A和B分配给同一个组,这个测试应该有一个自由度。
6.23. 多个前景支的支位点模型检测
当数据集中有多个前景支时:1)进行多个测试,然后在每个测试中设置一个感兴趣的分支作为前景分支;2)只进行一次测试,在此测试中中将所有感兴趣的分支设置为前景分支。那么,此时又会出现另一个问题就是进行多个测试时其他感兴趣分支是否应该被去除?这也许应该取决于具体的生物学问题。
6.24. p值校正
如果校正之后没有显著结果,我们可以选择adjP排序。
参考:https://groups.google.com/g/pamlsoftware/c/vty5QrRCUCk
6.25. 分支位点模型 (branch-site model) 与位点模型 (site model) 的比较
- 位点模型 (site model) 的一个相当大的缺点是,ω 是以位点所有密码子的平均值计算的。因此,位点模型不适用于 ω 在不同品系之间也有差异的数据。
- 分支位点模型 (branch-site model) 也是寻找正选择的位点,但会预先指定可能出现不同 ω 率的支系。
- 相比于位点模型的优点是考虑了不同的分枝具有不同的选择压,即具有不同的 ω 值。
- 分支位点模型将不同支系序列被先验地划分为可能发生正选择的前景支系群和只发生净化选择或中性进化的背景支系群。
- 分支位点模型让目标分枝具有一个不同的 ω 值,并没有让所有分枝的 ω 值独立进行计算(理论上这样是最好的)。这样算法很复杂,程序运行非常非常消耗时间。但其实也没必要这样做,因为正选择分析其实是两条序列比较后,分析dN/dS,再找正选择位点,其分析结果就应该是某个分枝上基因是否受到正选择,在序列某个位点上受到正选择。
6.26. 正选择基因的判断
- Q:若在目标分枝上,其 ω 值小于1,但是却能找到正选择位点。即该基因在该分枝上的dN/dS < 1,但是在某些位点上,dN/dS > 1。那么该基因是否属于正选择基因?
- A:属于。之所以为正选择基因,主要是因为基因的个别位点或多个位点存在正选择。当只有个别位点受到正选择压时,而其它多个位点存在纯化选择时,可能导致整体上的 ω 值小于1。此时,该基因也应该是属于正选择基因。
7. reference
- PAML User Guide:https://github.com/abacus-gene/paml/blob/master/doc/pamlDOC.pdf
- PAML 中文用户手册:https://blog.sciencenet.cn/blog-3433349-1241310.html
- PAML开发者网站:http://abacus.gene.ucl.ac.uk/software/#phylogenetic-analysis-by-maximum-likelihood-paml
- 陈连福博客-基因正选择分析原理:http://www.chenlianfu.com/?p=3084%E5%92%8Chttp://blog.sciencenet.cn/blog-460481-1163040.html
- 视频-使用codeml的分支模型(branch test)检测进化树上某一支的选择压力:https://www.bilibili.com/video/av10469605/?from=search&seid=8859495264060497812&vd_source=86d4997d664e6bb659b1d0f0dcd15267
- wiki-Ka/Ks ratio:https://en.wikipedia.org/wiki/Ka/Ks_ratio
- paml检测正选择的初学者指南:https://wap.sciencenet.cn/blog-3434047-1389140.html?mobile=1
- PAML选择压力分析:https://yuzhenpeng.github.io/2019/12/24/paml/
- 欢迎关注微信公众号:生信技工
- 公众号主要分享生信分析、生信软件、基因组学、转录组学、植物进化、生物学概念等相关内容,包括生物信息学工具的基本原理、操作步骤和学习心得。