|
写这篇文章的原因： 前段时间注意到一篇给CorMut提意见的文章，一直没有时间回应，终于在假期用几天时间把CorMut的代码逐行研读，现在做一个简单的说明。所提问题主要为三个：第一，序列处理时删除gap造成密码子错位。CorMut依照序列文件第一条的参考序列进行去除gap，不是去除所有的gap，去除的也是跟参考序列对应位置的gap，不可能造成密码子错位的发生。 第二，考虑到序列中可能存在混合碱基，所以有一个去除混合碱基的步骤，这是按照核苷酸为单位进行，遍历到哪个核苷酸返回那个位置涉及的密码子然后与相应参考序列位置的密码子比较，看是否引起氨基酸突变，最终选择引起突变的核苷酸替换混合碱基。数数也是没有问题的。 第三，计算LOD值时会用到二项分布，原版本是直接用分布的公式进行计算，这样在数值较大时的确会造成误差（大样本时choose函数会产生Inf等异常值），所以使用二项分布的dbinom会避免大样本时异常值的产生，但是1000以下的小样本对结果是没有影响的，不可能产生NaN。这的确是当时没有考虑到的问题。在新版本中已经进行了更新。 经过对代码逐行检查，并没有发现其他原则性的问题。同时在这个版本对部分内容进行了优化，对一些有分母的公式进行了一个平滑处理，避免了分母是0时造成异常值的情况；对画图功能增加了选择网络呈现形式的参数，用户可以自由设置网络以何种方式呈现。
就像其他软件一样，CorMut肯定时有bug的，过去我们也使用了大量数据进行测试，让CorMut变的更加完美。CorMut始终是开放的，我们会积极的接受任务批评和建议，遇到问题的话，欢迎与我进行交流。
CorMut开发始末： 最早CorMut使用的主流模型是用在HIV耐药位点识别上。当时并没有现成的工具来解决这件事情，当时研究组正在开展HIV耐药位点鉴定的工作，因此我们便首先根据这个算法编写了初步的脚本，使用这个算法对自己的数据进行分析，并对耐药位点实验验证，发现了几个全新的耐药位点。并对多套数据进行分析发现了与生物学相一致的结果。后来萌生了编写一个软件包的想法，然后就有了CorMut。原有的算法是ka/ks和条件ka/ks，单单这个肯定是不够的，后来增加了互信息，杰卡德系数等算法，对CorMut进行了较大的扩充，并增加了绘制网络图的功能。2014年博士毕业之后，参加工作，研究方向发生了较大转变，因为种种原因，CorMut也基本停止了更新。最近意识到CorMut不能丢掉，对CorMut的维护也是自我学习的一个过程。首先，在这里公布一个比较全面的使用指南，供参考。
CorMut使用指南： 相关突变的背景介绍： 在遗传学中，Ka/Ks的比值是指每个非同义位点非同义取代的数目与每个同义位点同义取代数目的比值。Ka/Ks比值可以被用作度量作用于蛋白质编码基因选择压力的指标。Ka/Ks等于1表示中性选择，也就是说观察到的非同义突变与同义突变的比值与期望的随机突变模型的比值相匹配。因此，氨基酸改变既不被选择也不被淘汰。Ka/Ks大于1表明氨基酸的改变是进化偏好的，也就是说这些突变能够更强组织的适应性。这种不寻常的状态可能反映了基因功能的改变或强迫机体去适应的环境条件改变。例如，高度变异的病毒，如HIV和HCV，对于新的抗病毒药物的突变可能在接受这些药物治疗的感染者体内是受到阳性选择的。相关突变（CorrelatedMutation）是进化生物学中的一个基本概念，在一个可编码蛋白质的基因上，氨基酸突变会受到蛋白质功能的限制。基于蛋白质功能的限制，一个突变的发生可能引发另外一个补偿突变的发生，因此，在特定条件下特定的突变模式就会形成，理解这些突变模式对于了解特定条件下基因的功能的改变至关重要。 CorMut组合了Ka/Ks比值和相关突变分析。软件包中所使用的计算单个位点Ka/Ks的方法是由Chen等提出的。CorMut提供了计算单个位点或特异氨基酸Ka/Ks并发现他们之间相关突变的函数。CorMut提供了三种方法来发现相关突变，包括条件选择压力、互信息和Jaccard系数。计算主要有两个步骤组成：第一，发现阳性选择位点；第二，计算阳性选择位点之间的突变关联。 值得注意的是第一个步骤是可选择的。同时，CorMut可以通过构建相关突变网络方便的比较两种条件下的相关突变关系。 CorMut原理介绍： 方法： 1. 单个密码子和单个氨基酸取代Ka/Ks的计算 单个密码子或单个特异的氨基酸取代的Ka/Ks是基于Chen的算法确定的。单个密码子或单个特异的氨基酸取代(X2Y)的Ka/Ks由下列公式计算： 分子中，NY和NS分别表示某个氨基酸发生有义突变和无义突变的数目；然后使用分母中的随机突变模型（也就是不存在选择压力情况下的模型）对该值进行校正，也就是式子的分母部分。在随机突变模型中，ft和fv分别表示转化和颠换的频率。它们的计算公式为：ft=Nt/ntS和fv=Nv/nvS（其中S表示样本总数；Nt和Nv表示所有样本中RT区实际发生转换和颠换的数量；nt和nv表示在所检测的序列中可能发生转换和颠换的数量（分别等于样本的数量和两倍的样本数量）。分母中，nY,t和nSt分别表示序列中既发生颠换又发生无义突变的数目和既发生颠换又发生有义突变的数目，nY,v和nSv分别表示序列中既发生转换又发生无义突变的数目和既发生转换又发生有义突变的数目。同时，对突变使用对数比值比（logodds ratio，LOD）进行评估。LOD的计算公式如下：
N表示密码子的突变总数；q计算公式如下： 如果LOD大于2，阳性选择被认为是显著的。 2. 条件Ka/Ks的计算 条件选择比值可以定义为氨基酸Y突变时氨基酸Z突变的Ka/Ks与氨基酸Y没有突变时氨基酸Z突变的Ka/Ks的比值，计算公式为： 公式中NZaYa表示同时存在氨基酸突变Y和Z的样本总数； NZsYa表示Z同义突变， Y氨基酸突变的样本总数；NZaYo和NZsYo分别表示在不发生氨基酸Y突变的情况下，氨基酸Z非同义突变和同义突变的样本总数。 LOD计算公式如下： 3. 互信息 互信息(MutualInformation)是信息论里提到的一种的信息度量的方法，被用来度量两个事件集合之间的相关性[37]。已有文献报道互信息被用作度量残基取代之间的相关性。也就是说，在一个具有N个密码子的蛋白质的多序列多序列比对中N列中的每一列被认为是一个离散变量Xi(1<><>。当计算密码子之间的互信息时，Xi会计算20种氨基酸的每一种的概率，而计算两个密码子特定氨基酸的互信息时，Xi只计算两个部分的概率，也就特异的氨基酸取代或者其他取代。互信息描述两个变量Xi和Xj之间的依赖关系，通过关联Xi和Xj两列的相关指标来度量两者依赖的强弱。
4. Jaccard 系数 Jaccard系数可以度量两个变量之间的相似性，并且已经广泛应用于度量相关突变[38-40]。对于一个突变对X和Y，Jaccard系数计算的公式为Nxy/(Nxy+Nx0+Ny0)，在式中，Nxy代表X和Y都发生突变的序列的条数，Nx0代表只有X突变而Y没有突变的序列的数目，Ny0只有Y突变而X没有突变的序列的数目。Jaccard系数被认为可以有效的避免X和Y都没有发生突变序列条数较多引起的的稀有突变结果的夸大计算。
CorMut包含的主要函数及使用说明： ckaksCodon, ckaksAA, miCodon, miAA, jiAA设置了参数 'setPosition'，可以指定感兴趣的位置进行计算，避免其他位点影响结果的展示，降低运算的复杂度。 CorMut实战： 1. 多序列比对文件的处理 seqFormat函数能够用普通碱基取代混合碱基，基于参考序列删除序列比对中的空缺。由于一个混合碱基可能与多个普通碱基相对应，seqFormat会随机选择能够引发氨基酸突变的普通碱基。这里，未治疗和治疗的HIVPR序列被用作为例子来演示函数的使用： library(CorMut) examplefile=system.file('extdata','PI_treatment.aln',package='CorMut') examplefile02=system.file('extdata','PI_treatment_naive.aln',package='CorMut') example=seqFormat(examplefile) example02=seqFormat(examplefile02) 2. 计算单个密码子位置的kaks result=kaksCodon(example) fresult=filterSites(result) head(fresult)
3. 计算单个氨基酸突变的Ka/Ks result=kaksAA(example) fresult=filterSites(result) head(fresult) 4. 计算密码子之间的条件Ka/Ks（条件选择压力） result=ckaksCodon(example) plot(result) fresult=filterSites(result) head(fresult) 5. 计算氨基酸突变之间的条件Ka/Ks（条件选择压力） result=ckaksAA(example) plot(result) fresult=filterSites(result) head(fresult)
6. 计算密码子之间的信息熵 包含两个矩阵的列表将会被返回，他们分别是互信息矩阵和p值矩阵。矩阵的第i行第j列元素代表第i个密码子和第j密码子之间互信息的值或p值。 result=miCodon(example) plot(result) fresult=filterSites(result) head(fresult) 7. 计算氨基酸突变之间的信息熵 包含两个矩阵的列表将会被返回，他们分别是互信息矩阵和p值矩阵。矩阵的第i行第j列元素代表第i个氨基酸突变和第j个氨基酸突变之间互信息的值或p值。 result=miAA(example) plot(result) fresult=filterSites(result) head(fresult)
8. 计算氨基酸突变之间的Jaccard系数 包含两个矩阵的列表将会被返回，他们分别是Jaccard系数矩阵和p值矩阵。矩阵的第i行第j列元素代表第i个氨基酸突变和第j氨基酸突变之间Jaccard系数的值或p值。 result=jiAA(example) plot(result) fresult=filterSites(result) head(fresult) 9. 通过构建相关突变网络比较两种条件下的相关突变 biCompare函数能够比较两种条件下的相关突变关系，比如HIV未治疗和治疗情况下相关突变关系的改变. biCompare 的结果能够通过plot函数实现可视化。只有阳性选择密码子或氨基酸突变被展现在图上，蓝色圆圈代表第一种条件下独特的阳性选择密码子或氨基酸，也就是说这些位点在第二种条件下是受到中性选择或阴性选择.而红色圆圈代表第二种条件下出现的独特的阳性选择的密码子或氨基酸。可以通过控制Plot的参数plotUnchanged展现那些两种条件下均为阳性选择的位点，plotUnchanged为FALSE时可以实现这一点。 举例: biexample=biCkaksAA(example02,example) #biexample=biCkaksCodon (example02,example) #biexample=biMICodon (example02,example) #biexample=biMIAA (example02,example) #biexample=biJIAA (example02,example) result=biCompare(biexample) plot(result)
|
|
|
