CN108763872A - 一种分析预测癌症突变影响lir模体功能的方法 - Google Patents

Abstract

本发明公开了一种分析预测癌症突变影响LIR模体功能的方法，包括以下步骤：(1)利用LIR模体数据训练预测模型；(2)利用模型对癌症突变前后的肽段数据分别打分；(3)利用模型对候选肽段的肽段数据进行打分，得出打分结果与该候选肽段对应的肽段数据其属于真实LIR模体或非真实LIR模体的贝叶斯概率之间的关系，接着将癌症突变前后对应的肽段分别作为候选肽段，并代入步骤(2)中的打分结果，得到癌症突变前后对应的肽段它们的贝叶斯概率；(4)计算该癌症突变对于LIR功能的可选影响方式的概率。本发明通过利用LIR模体数据、以及癌症突变前后对应的肽段数据，系统性的研究癌症突变对LIR模体功能的影响，能够预测出癌症突变对于LIR功能的可选影响方式的概率。

Description

一种分析预测癌症突变影响LIR模体功能的方法

技术领域

本发明属于自噬领域，更具体地，涉及一种分析预测癌症突变影响LIR模体功能的方法，利用该方法可以分析预测已有的癌症突变对LIR模体进而对自噬通路的影响，当然，随着未来新癌症突变的发现，该方法也可以用于判断这些新癌症突变对LIR模体进而对自噬通路的影响。

背景技术

细胞自噬是细胞内清除损坏组件和蛋白质或在饥饿条件下通过降解自身一部分蛋白质以渡过难关的一种机制。在正常生理条件下，细胞通过自噬清除受损的细胞器或蛋白质聚集物，自噬的受阻会引起多种疾病。在自噬的核心机制中，自噬体的形成一个复杂而保守的过程，其中，LC3(酵母中ATG8的人类同源)对于自噬体的形成起着至关重要的作用，其关系到自噬体膜的弯曲、自噬底物的募集的过程，是自噬过程中唯一一直定位于自噬体膜的蛋白质。除了影响到自噬体膜的弯曲外，LC3决定了需要募集到自噬体内并降解的细胞内容物，其通过一个保守的模体，即LIR(LC3 interaction region，亦称AIM，ATG8interaction motif)模体来识别和特异性的选择底物。研究表明，几乎所有LIR模体都符合[W/Y/F]XX[L/I/V]的序列模式。这一序列模式的改变或是周围氨基酸微环境的改变，都将有可能影响到LC3与底物的结合进而影响到细胞自噬的过程。影响LIR功能的一大因素是癌症突变。目前，主流提供癌症突变信息的数据库有TCGA，ICGC，ClinVar，COSMIC等等。对于LIR，目前尚缺乏系统性的研究。2014年，Ioanna Kalvari等人开发了在线服务平台iLIR并于2016年更新(Ioanna Kalvari et.al.,Autopahgy,Volume 10,2014-Issue 5)。iLIR使用收集的LIR数据，结合PSSM算法，在八个模式生物中鉴定了潜在的LIR模体。然而，其并没有对鉴定的数据进行进一步分析与验证，并且，PSSM算法的预测准确性有待进一步提高。而系统性的研究癌症突变对LIR模体功能影响的工作更为缺乏。

发明内容

针对现有技术的以上缺陷或改进需求，本发明的目的在于提供一种分析预测癌症突变影响LIR模体功能的方法，通过利用LIR模体数据、以及癌症突变前后对应的肽段数据，并通过对该分析预测方法的整体流程步骤设置、关键步骤(如预测模型的训练步骤、贝叶斯概率的计算步骤、可选影响方式的概率计算步骤等)的具体数据处理及计算过程等进行改进和进一步优选，系统性的研究癌症突变对LIR模体功能的影响，进而可以预测出癌症突变对于LIR功能的若干种可选影响方式的概率，填补了现有技术的空白；并且由于本发明是基于逻辑回归的训练方式对预测模型进行训练，具有训练速度快，分析预测效率高等优点。

为实现上述目的，按照本发明，提供了一种分析预测癌症突变影响LIR模体功能的方法，其特征在于，包括以下步骤：

(1)利用LIR模体数据训练预测模型；该模型用于以LIR模体数据作为参照物评价某一肽段数据与这些参照物之间的相似度分数；

(2)利用所述步骤(1)得到的预测模型，基于生物的正常蛋白质序列和该生物的癌症突变的肽段数据，对突变前后的肽段数据与所述LIR模体数据之间的相似度分别打分；

(3)利用所述步骤(1)得到的预测模型对候选肽段的肽段数据进行打分，并计算得出打分结果与该候选肽段对应的肽段数据其属于真实LIR模体或非真实LIR模体的贝叶斯概率之间的关系；接着，将癌症突变前后对应的肽段分别作为候选肽段，并将所述步骤(2)得到的对应打分结果代入，从而得到所述癌症突变前后对应的肽段它们的贝叶斯概率；

(4)接着，利用所述步骤(3)得出的所述癌症突变前后对应的肽段它们的贝叶斯概率进一步计算该癌症突变对于LIR功能的若干种可选影响方式的概率，进而判断得出该癌症突变对于LIR功能的最终影响方式。

作为本发明的进一步优选，所述步骤(1)中，所述LIR模体数据为多个LIR模体数据，任意一个LIR模体数据包括14个氨基酸，这14个氨基酸由核心LIR模体的4个氨基酸、以及该核心LIR模体左侧的6个氨基酸和右侧的4个氨基酸构成。

作为本发明的进一步优选，所述步骤(2)中，所述生物的正常蛋白质序列为人类全蛋白质序列，所述生物的癌症突变的肽段数据为人类癌症突变的肽段数据。

作为本发明的进一步优选，所述步骤(3)是先以LIR模体数据作为阳性数据，以与所述阳性数据中任意一个LIR模体数据出现在同一条蛋白上且符合[W/Y/F]XX[L/I/V]序列模式但未被所述阳性数据涵盖的肽段数据作为阴性数据，通过对阳性数据打分所得分值、以及阴性数据打分所得分值分别进行拟合正态分布，从而分别得到阳性数据和阴性数据分值分布的概率密度函数pdf_p和pdf_n，从而进一步在所述步骤(3)中基于阳性数据分值分布的概率密度函数pdf_p和阴性数据分值分布的概率密度函数pdf_n得出候选肽段的打分结果与该候选肽段对应的肽段数据其属于真实LIR模体或非真实LIR模体的贝叶斯概率之间的关系；

所述步骤(3)中，记所述候选肽段为pep，所述预测模型对该候选肽段的肽段数据进行打分给出的分值为S_pep，则该候选肽段其预测结果为真实LIR模体的概率p(True|S_pep)满足：

该候选肽段其预测结果为非真实LIR模体的概率p(False|S_pep)满足：

其中，pdf_p(S_pep)、pdf_n(S_pep)为将S_pep分别代入所述概率密度函数pdf_p和所述概率密度函数pdf_n对应得到的值，p(True)和p(False)均为预先设定的值；优选的，所述p(True)为所述步骤(1)得到的预测模型利用留一法检测得到的ROC曲线的曲线下面积，所述p(False)为1。

作为本发明的进一步优选，所述步骤(4)中，所述LIR功能的若干种可选影响方式包括突变前属于真实LIR模体而突变后属于非真实LIR模体、以及突变前属于非真实LIR模体而突变后属于真实LIR模体两种情况；所述最终影响方式对应为所述若干种可选影响方式中的概率最大值者。

作为本发明的进一步优选，所述步骤(4)中，记癌症突变前的正常蛋白质序列肽段为W，癌症突变后的癌症突变肽段为M，将W和M分别作为候选肽段为pep得到p(True|S_W)、p(False|S_W)、p(True|S_M)、p(False|S_M)，则突变前属于真实LIR模体而突变后属于非真实LIR模体的概率p(loss)满足：

p(loss)＝p(True|S_W)*p(False|S_M)；

突变前属于非真实LIR模体而突变后属于真实LIR模体的概率p(gain)满足：

p(gain)＝p(False|S_W)*p(True|S_M)。

作为本发明的进一步优选，所述步骤(1)得到的预测模型利用留一法检测得到的ROC曲线的曲线下面积为0.96。

通过本发明所构思的以上技术方案，与现有技术相比，由于同时利用LIR模体数据、以及癌症突变前后对应的肽段数据，能够系统性的研究癌症突变对LIR模体功能的影响，进而预测出癌症突变对于LIR功能的若干种可选影响方式的概率。

本发明通过对阳性数据、阴性数据基于以下原则：同一条蛋白上最多只出现一个LIR模体数据，当一条蛋白上已经存在一个实验验证的LIR模体数据时，其他位于该蛋白上且满足[W/Y/F]XX[L/I/V]序列模式的肽段序列但是未被实验证实为类LIR模体的数据均属于非LIR模体数据(实验验证的LIR模体数据拥有最高优先级，已经被实验验证的LIR模体数据必然属于阳性数据)。

本发明可通过人工收集实验验证的数据，并进一步计算分析预测癌症突变对LIR模体功能的影响。本发明中的预测模型可以基于GPS(Group-based prediction system)系列算法进行训练，通过控制训练流程，尤其是基于逻辑回归的训练方式进行训练，能够大幅提高训练速度，为优化调参提供更多的操作空间。

附图说明

图1是本发明分析预测方法的概括流程示意图。

图2是本发明实施例1中分析预测方法的具体流程示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明中的分析预测方法概述起来，包括5个部分，分别为：模型训练、肽段打分、概率估算、结果筛选、数据整合。如图1所示，各个部分之间的关联如下：

①通过⑤提供的已知LIR模体数据训练预测模型；

②通过⑤提供的人类全蛋白质序列和人类癌症突变数据的肽段为突变前后的肽段数据打分；

③基于②的打分结果估算候选肽段突变前后的贝叶斯概率，并基于此概率估算突变对于LIR功能影响方式和相应的概率；

④基于③提供的概率筛选结果。

实施例1

如图2所示，本发明中的分析预测方法具体包括以下分析流程：

S1为训练模块提供数据；

S2训练预测模型；

S3基于上一步训练的模型构建待选肽段打分器；

打分器可用于打分评价候选肽段与已知阳性肽段之间的整体相似性(已知阳性肽段数据即LIR模体数据，需要利用LIR模体数据来训练该预测模型)，该打分器所基于的打分规则可参考相关现有技术，例如可以采用GPS算法等。

以GPS算法为例，GPS算法基于相似的序列具有相似的功能的假设，通过计算候选序列与阳性数据集之间的平均相似度来评价候选序列是阳性数据的可能性。以两条长度为m的肽段Pep1和Pep2为例，假设第i(0<i<m+1)个位置上的氨基酸分别为A1和A2，查询氨基酸相似性矩阵BLOSUM62得到A1和A2的相似度为S_i,则Pep1和Pep2的相似度为：

上述公式默认假设了每个位置的氨基酸重要性相同，但实际情况下并非如此，故假设每个位置的重要性由权重数组W给出，W_i为第i个位置的重要性，则上述公式演变为：

因此，对于一条候选肽段PepX与N条(记为Pep₁-Pep_N)阳性肽段之间的平均相似度计算公式可以如下：

得到的Avg_Sim即可作为打分结果。

S4整合S9和S10提供的数据；

S5通过S4得到突变前肽段；

S6通过S4得到突变后的肽段；

S7由S3对S5打分得到；

S8由S3对S6打分得到；

S9从公共数据库整合得到；

S10从公共数据库整合得到；

S11设计从预测分值估算贝叶斯概率的算法；

S12通过S11对S7的分值估算得到；

S13通过S11对S8的分值估算得到；

S14通过S12和S13的概率估算突变对LIR功能不同类型的改变的概率；

S15筛选结果，对于通过分析高可信度的概率结果了解突变影响LIR模体功能的方式。

进一步的，LIR模体的预测可以基于本发明发明人此前开发的GPS(Group-basedprediction system)系列算法。

本发明可以采用三种GPS模型的训练方法，即，初始的GPS训练算法(例如，可参考Liu Zexian et.al.,PLoS One,doi:10.1371/journal.pone.0034370)、基于模拟退火算法(例如，可参考Deng Wankun et.al.,Brief Bioinform,doi:10.1093/bib/bbw041)和基于逻辑回归的算法(例如，可参考Boyd,C.R et.al.,The Journal of trauma.doi:10.1097/00005373-198704000-00005)。

1.GPS训练算法

GPS训练算法，如前所述，基于相似的序列具有相似的功能的假设，通过计算候选序列与阳性数据集之间的平均相似度来评价候选序列是阳性数据的可能性。对于两条长度为m的肽段Pep1和Pep2，假设第i(0<i<m+1)个位置上的氨基酸分别为A1和A2，查询氨基酸相似性矩阵BLOSUM62得到A1和A2的相似度为S_i,则Pep1和Pep2的相似度为：

上述公式默认假设了每个位置的氨基酸重要性相同，但实际情况下并非如此，故假设每个位置的重要性由权重数组W给出，W_i为第i个位置的重要性，则上述公式演变为：

因此，对于一条候选肽段PepX与N条(记为Pep₁-Pep_N)阳性肽段之间的平均相似度计算公式可以如下：

因此，想得到一个好的预测结果，GPS模型需要合理设置两个参数，即位置权重和氨基酸相似性矩阵。默认情况下，位置权重为各个位置相同均为1，而氨基酸相似性矩阵则设置为BLOSUM62。GPS算法最初的训练策略为随机突变：

对于位置权重而言，随机抽取一个位置，对其随机加或减一个固定数值X，通过留一法检测模型性能，若预测性能提升则保留这个更改，若不提升则回到突变前的模型。对于氨基酸相似性矩阵而言，策略也是类似的。

2.基于模拟退火的训练方式

基于模拟退火算法的GPS模型训练过程与上述训练过程最大区别在于：训练过程中若出现性能下降(即当随机突变一个权重或氨基酸相似性矩阵中的一个值之后，留一法检验模型预测性能下降)，有一定概率会保留此结果并在此结果基础上继续训练。模拟退火训练过程参数设置主要有三个：突变总次数上限，退火速度和初始温度.接受性能下降的结果的概率由以下公式给出：

p＝e^delta/T0

其中，delta为性能下降的幅度，一般用最优模型的Sn减去当前模型的Sn，T0为当前温度，受初始温度和退火速度设置的影响。模拟退火算法能一定程度上解决局部最优的问题，但是对于训练速度反而会降低，在GPS模型的训练中，一般设置突变总次数为20,000次，初始温度设置为100,000，退火速度设置为T0＝T0*Math.exp(-0.001*DenStep)，其中，DenStep为已进行的退火次数；当连续退火超过50次并且无性能提升时，温度置为0.25485142899，退火次数置为160。

3.基于逻辑回归的训练方式

通常的训练算法是：

(1)根据位置权重和打分矩阵，计算每个peptide和positive数据的分值，取平均

(2)随机突变矩阵中分值(+1或-1)，保留性能提高的突变；以肽段长度为15为例：

N₊代表阳性数据个数，

P_j代表给定一条肽段上第j个位置氨基酸，

T_ij代表阳性训练集中第i条肽段第j个位置上氨基酸，

W_j代表第j个位置上权重，

M代表打分矩阵

权重W未知时，对以上公式改写下可得下式：

我们令：

则

进一步写成：

因此，可以利用逻辑回归计算回归系数W：

即，将相似性分值通过logistic函数转化为预测结果为真的拟概率表示(P>0.5则为真)，再通过逻辑回归拟合回归系数W。类似的，也可以将氨基酸相似性矩阵转化为回归系数求解。

本发明通过改进训练流程，尤其可以将基于模拟退火算法的训练方式改进为基于逻辑回归的训练方式，大幅提高了训练速度，为优化调参提供了更多的操作空间，在此，逻辑回归可以使用Python sklearn模块进行训练，具体参数设置为：Cs＝[0.1,0.01,0.001,0.000001,1,10,1000],penalty＝'l1',cv＝10,solver＝'liblinear',refit＝True,n_jobs＝1。可以使用从核心LIR模体(4个氨基酸)左侧6个氨基酸到其右侧4个氨基酸，总长14个氨基酸的肽段于LIR模体的预测分析(即，核心肽段上游6个氨基酸，下游4个氨基酸，假设某核心肽段序列为WABL，那么对应的用于预测的氨基酸序列应该是XXXXXXWABLXXXX)。

使用留一法(LOO，leave one out)对模型性能检测，其ROC曲线下面积为0.96，模型稳定且预测性能令人满意。作为比较，我们也实现了了iLIR所采用的PSSM算法并通过留一法对其性能进行检验，结果显示，PSSM算法ROC曲线下面积为0.87。

对于训练好的GPS预测模型，分别对阳性数据和阴性数据打分，得到阳性数据与阴性数据的预测分值，经过拟合正态分布，得到阳性和阴性数据分值分布的概率密度函数(pdf_p和pdf_n)。阳性数据指的是训练模型时所使用的收集的实验验证的LIR模体数据，阴性数据是指与收集的阳性数据位于同一条蛋白上且符合[W/Y/F]XX[L/I/V]序列模式(W/Y/F/L/I/V均代表指定氨基酸，X则代表任一氨基酸)但是未被实验证实为LIR模体的数据。本发明中阳性数据、阴性数据的来源基于以下原则：同一条蛋白上最多只出现一个LIR模体数据，当一条蛋白上已经存在一个实验验证的LIR模体数据时，其他位于该蛋白上且满足[W/Y/F]XX[L/I/V]序列模式的肽段序列但是未被实验证实为类LIR模体的数据均属于非LIR模体数据(实验验证的LIR模体数据拥有最高优先级，已经被实验验证的LIR模体数据必然属于阳性数据)。

对于候选肽段Pep，假设其由GPS算法给出的相似性分值为Spep其预测结果为真实LIR模体的概率(即，预测结果为真的概率，也就是候选肽段是真实LIR模体的概率)由以下公式给出：

其中p(True)由ROC曲线的曲线下面积给出，而p(False)则设置为1。

同样的，Pep预测结果为非真实LIR模体的概率(即，预测结果为假的概率，也就是候选肽段为不真实LIR模体的概率)由以下公式给出：

对于发生在LIR模体序列区间的癌症突变，考虑两种情况：突变前与LC3/ATG8结合而突变后无法与LC3/ATG8结合(loss)；突变前不与LC3/ATG8结合而突变后与LC3/ATG8结合(gain)；真实LIR模体会与LC3/ATG8结合，非真实LIR模体则无法与LC3/ATG8结合。对于某条肽段，假设其突变前的相似性分值为SW而突变后相似性分值为SM，则该突变对于影响属于loss和gain的概率分别由以下公式给出：

p(loss)＝p(True|S_W)*p(False|S_M)

p(gain)＝p(False|S_W)*p(True|S_M)

通过上述公式，可以有效的估算癌症突变对于LIR模体功能的影响的概率和方式。

本发明中的所采用的癌症突变信息例如可取自TCGA，ICGC，ClinVar，COSMIC数据库，例如可以预测某一种或多种甚至是全部已知癌症突变对LIR模体功能的影响；当然，也不排除对其他可能新发现的癌症突变其对LIR模体功能影响的预测。LIR模体数据例如可以取自iLIR平台收集的LIR数据。在预测模型训练过程中，本发明可以利用由多个LIR模体数据(即，多个实验验证的LIR模体数据，当然也可以收集所有实验验证的LIR模体数据)构成的LIR数据集训练预测模型，本发明训练预测模型所采用的LIR数据集可以是已知LIR数据集，当然，随着未来新LIR的发现，也可以将新发现的LIR更新至该LIR数据集中。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种分析预测癌症突变影响LIR模体功能的方法，其特征在于，包括以下步骤：

(1)利用LIR模体数据训练预测模型；该模型用于以LIR模体数据作为参照物评价某一肽段数据与这些参照物之间的相似度分数；

(2)利用所述步骤(1)得到的预测模型，基于生物的正常蛋白质序列和该生物的癌症突变的肽段数据，对突变前后的肽段数据与所述LIR模体数据之间的相似度分别打分；

(3)利用所述步骤(1)得到的预测模型对候选肽段的肽段数据进行打分，并计算得出打分结果与该候选肽段对应的肽段数据其属于真实LIR模体或非真实LIR模体的贝叶斯概率之间的关系；接着，将癌症突变前后对应的肽段分别作为候选肽段，并将所述步骤(2)得到的对应打分结果代入，从而得到所述癌症突变前后对应的肽段它们的贝叶斯概率；

(4)接着，利用所述步骤(3)得出的所述癌症突变前后对应的肽段它们的贝叶斯概率进一步计算该癌症突变对于LIR功能的若干种可选影响方式的概率，进而判断得出该癌症突变对于LIR功能的最终影响方式。

2.如权利要求1所述分析预测癌症突变影响LIR模体功能的方法，其特征在于，所述步骤(1)中，所述LIR模体数据为多个LIR模体数据，任意一个LIR模体数据包括14个氨基酸，这14个氨基酸由核心LIR模体的4个氨基酸、以及该核心LIR模体左侧的6个氨基酸和右侧的4个氨基酸构成。

3.如权利要求1所述分析预测癌症突变影响LIR模体功能的方法，其特征在于，所述步骤(2)中，所述生物的正常蛋白质序列为人类全蛋白质序列，所述生物的癌症突变的肽段数据为人类癌症突变的肽段数据。

4.如权利要求1所述分析预测癌症突变影响LIR模体功能的方法，其特征在于，所述步骤(3)是先以LIR模体数据作为阳性数据，以与所述阳性数据中任意一个LIR模体数据出现在同一条蛋白上且符合[W/Y/F]XX[L/I/V]序列模式但未被所述阳性数据涵盖的肽段数据作为阴性数据，通过对阳性数据打分所得分值、以及阴性数据打分所得分值分别进行拟合正态分布，从而分别得到阳性数据和阴性数据分值分布的概率密度函数pdf_p和pdf_n，从而进一步在所述步骤(3)中基于阳性数据分值分布的概率密度函数pdf_p和阴性数据分值分布的概率密度函数pdf_n得出候选肽段的打分结果与该候选肽段对应的肽段数据其属于真实LIR模体或非真实LIR模体的贝叶斯概率之间的关系；

所述步骤(3)中，记所述候选肽段为pep，所述预测模型对该候选肽段的肽段数据进行打分给出的分值为S_pep，则该候选肽段其预测结果为真实LIR模体的概率p(True|S_pep)满足：

该候选肽段其预测结果为非真实LIR模体的概率p(False|S_pep)满足：

其中，pdf_p(S_pep)、pdf_n(S_pep)为将S_pep分别代入所述概率密度函数pdf_p和所述概率密度函数pdf_n对应得到的值，p(True)和p(False)均为预先设定的值；优选的，所述p(True)为所述步骤(1)得到的预测模型利用留一法检测得到的ROC曲线的曲线下面积，所述p(False)为1。

5.如权利要求1所述分析预测癌症突变影响LIR模体功能的方法，其特征在于，所述步骤(4)中，所述LIR功能的若干种可选影响方式包括突变前属于真实LIR模体而突变后属于非真实LIR模体、以及突变前属于非真实LIR模体而突变后属于真实LIR模体两种情况；所述最终影响方式对应为所述若干种可选影响方式中的概率最大值者。

6.如权利要求5所述分析预测癌症突变影响LIR模体功能的方法，其特征在于，所述步骤(4)中，记癌症突变前的正常蛋白质序列肽段为W，癌症突变后的癌症突变肽段为M，将W和M分别作为候选肽段为pep得到p(True|S_W)、p(False|S_W)、p(True|S_M)、p(False|S_M)，则突变前属于真实LIR模体而突变后属于非真实LIR模体的概率p(loss)满足：

p(loss)＝p(True|S_W)*p(False|S_M)；

突变前属于非真实LIR模体而突变后属于真实LIR模体的概率p(gain)满足：

p(gain)＝p(False|S_W)*p(True|S_M)。

7.如权利要求1所述分析预测癌症突变影响LIR模体功能的方法，其特征在于，所述步骤(1)得到的预测模型利用留一法检测得到的ROC曲线的曲线下面积为0.96。