[发明专利]一种基于基因突变频率的癌症驱动基因的筛选方法

摘要

本发明公开了一种基于基因突变频率的癌症驱动基因的筛选方法，属于癌症医学领域。该方法包含如下步骤：(1)肿瘤基因突变数据获取；(2)突变数据预处理；(3)筛选每个基因的邻近基因；(4)计算每种突变的背景突变数据；(5)根据突变分值使用假设检验筛选突变基因。本发明方法不仅利用现代高通量测序技术及DNA数据处理软件，使用经典的聚类算法和统计方法，而且与影响基因突变的生物因素相结合使癌症驱动基因的筛选更加精确，对新型抗癌药物的研发和癌症临床诊疗都具有重要意义。

主权项

一种基于基因突变频率的癌症驱动基因的筛选方法，其特征在于，步骤如下：(1)肿瘤基因突变数据获取：对多名患同种癌症的患者的肿瘤细胞和正常细胞的DNA进行高通量测序，对测序得到的DNA序列与标准基因hg19进行比对，得到肿瘤细胞DNA和正常细胞DNA的基因突变位点，取肿瘤细胞DNA的特有突变位点，对突变位点进行注释，得到突变的基因名，突变类型；最后将这些数据整理成数据集：突变数据，覆盖区域，协变量，表格如下：表1 突变数据表1中包含的信息有每个突变位点所在对应的基因、病人编号、突变影响及突变类别；突变影响为突变对蛋白质合成影响，包含silent、nonsilent、noncoding三种，silent：即同义突变，nonsilent：可以导致蛋白质发生改变的突变，noncoding：发生在非编码区的突变；突变类别包含7种突变类别：1.CpG发生转换，2.CpG发生颠换，3.CpG外的C:G发生转换，4.CpG外的C:G发生颠换，5.A:T发生转换，6.A:T发生颠换，7.null+indel突变，包含无义突变、插入删除突变和剪接位点发生突变；表2 覆盖区域表2中包含的信息有每个突变位点的基因、突变影响、突变类别及病人编号，其中基因、突变影响、突变类别具体信息同表1；L1、L2、……是病人编号，病人编号对应的信息是病人基因里可能发生同一影响相应类别突变的所有的碱基数；表3 协变量基因的表达水平、复制时间及染色体状态的数据从NCBI数据库中获得；(2)数据预处理：对初始肿瘤基因突变数据进行整理；将突变数据整理成三个3维矩阵，3个维度g为基因，c为突变类型，p为病人编号；矩阵统计的是每个病人、每个基因内发生每种突变影响、每种突变类型的基因突变总个数；将覆盖区域表整理成三个3维矩阵，3个维度g为基因，c为突变类型，p为病人编号；矩阵统计的是每个病人，每个基因内能发生每种突变影响、每种突变类型的碱基总个数；其中c维度在所有突变类别的基础上增加一列nc+1，统计的为所有突变类型的突变个数总和；将协变量表整理成矩阵Vv,g，其中v为协变量，g为基因；将Vv,g进行标准化得到Zv,g，即用公式(1)将各个协变量数据转化成均值为0、方差为1的数据Zv,g=Vv,g-1ngΣi=1ngVv,i1ng-1Σj=1ng(Vv,i-1ngΣi=1ngVv,i)2---(1)]]>其中，ng为基因总数；i,j为选中的一个基因；Vv,i为基因i的协变量值；(3)筛选每个基因的邻近基因：在协变量差距不大情况下，基因内碱基的突变概率大致相同的基因视为该基因的邻近基因；筛选步骤如下：1)首先使用K‑means算法对Vv,g进行聚类，同时使用轮廓系数法确定聚类的类别数，得到每个基因的类别，轮廓系数计算方法如下：si＝(bi‑ai)/max(bi,ai)   (2)ai用于量化簇内凝聚度：对第i个基因gi，计算gi与其同一个簇内的所有其他元素距离的平均值；bi用于量化簇之间分离度：选取gi外的一个簇b，计算gi与b中所有点的平均距离，遍历所有其他簇，找到gi与其他所有簇的最近的平均距离，记作bi；计算所有基因g的轮廓系数，求出平均值即为当前聚类的整体轮廓系数，挑选最大的轮廓系数对应的聚类类别数K；2)然后在基因所属的类里使用假设检验算法选出每个基因的邻近基因，其中零假设为基因i为基因g的邻近基因，数据采用背景突变数据；统计背景突变数据即发生在非编码区和非同义突变区域内的突变，计算方法如下：ngbkgd=Σp=1np(ng,nc+1,psilent+ng,nc+1,pnoncoding)---(3)]]>Ngbkgd=Σp=1np(Ng,nc+1,psilent+Ng,nc+1,pnoncoding)---(4)]]>同一基因内，由于碱基处在同一环境下，每个碱基发生突变的概率相同，所以基因内N个碱基发生n个突变的概率分布属于二项分布；若基因i是基因g的邻近基因，即它们的突变属于同一个二项分布，则基因i和基因g的突变数据服从beta‑二项分布；以此做假设检验，零假设为基因i是基因g的邻近基因，p值为公式(5)中的Qi,g，在同类协变量的基因中筛选出p值大于0.05的基因作为基因g的邻近基因Zg；Qi,g=2min(Qi,gleft,1-Qi,gleft)---(5)]]>Qi,gleft=HC(nibkgd,Nibkgd,ngbkgd,Ngbkgd)---(6)]]>HC为beta‑二项分布H的连加，具体的计算如下：HC(n1,N1,n2,N2)=Σn=0n1H(n,N1,n2,N2)---(7)]]>H(n1,N1,n2,N2)=N1n1B(n1+α,N1-n1+β)B(α,β)=Γ(N1+1)Γ(N2+1)Γ(n1+n2+1)Γ(N1+N2-n1-n2+1)Γ(n1+1)Γ(n2+1)Γ(N1-n1+1)Γ(N2-n2+1)Γ(N1+N2+2)---(8)]]>其中，α＝n2+1，β＝N2‑n2+1，Γ为gamma函数；筛选出每个基因的邻近基因，统计基因和所有邻近基因的背景突变碱基数xg和每种突变所在区域碱基总数Xg；xg=ngbkgd+Σi∈Zgnibkgd---(9)]]>Xg=Ngbkgd+Σi∈ZgNibkgd---(10)]]>公式(9)和(10)中，i∈Zg基因i是基因g的邻近基因；(4)平均突变数据，计算每个突变位点的背景突变数据，统计所有样本中每个病人，每种突变类别的总突变数，根据突变频率计算每个基因、每个病人、每种突变类型对应的背景突变数据；忽略突变影响，统计每个基因、每个病人、每种突变类型的突变数据及区域碱基总数ng,c,ptotal=ng,c,psilent+ng,c,pnonsilent+ng,c,pnoncoding---(11)]]>Ng,c,ptotal=Ng,c,psilent+Ng,c,pnonsilent+Ng,c,pnoncoding---(12)]]>统计每种突变类别的突变数据及区域碱基总数nctotal=Σg=1ngΣp=1npnc,ptotal---(13)]]>Nctotal=Σg=1ngΣp=1npNc,ptotal---(14)]]>统计所有突变位点的突变数据及区域碱基总数noveralltotal=nnc+1total---(15)]]>Noveralltotal=Nnc+1total---(16)]]>统计每个病人的突变数据及区域碱基总数nptotal=nnc+1,ptotal---(17)]]>Nptotal=Nnc+1,ptotal---(18)]]>计算每个基因、每个病人、每种突变类型对应的背景突变数据xg,c,p及背景区域碱基总数Xg,c,p：xg,c,p=xgnctotalnoveralltotalnptotal1npnoveralltotal---(19)]]>Xg,c,p=XgNctotalNoveralltotalNptotal1npNoveralltotal---(20)]]>(5)筛选驱动基因：基于突变概率及突变类型设计一种突变分值，计算样本突变数据每个基因的总分值，使用假设检验算出每个基因是驱动基因的p值，算出对应的错误发现率，根据错误发现率筛选出驱动基因；在基因内部，nonsilent区域和所在的背景突变区域的突变数据服从beta‑二项分布，所以非同义突变区域里有0个、1个碱基发生突变的概率为大于等于2个发生突变的概率为计算公式如下：Pg,c,p(0)=H(0,Ng,c,pnonsilent,xg,c,p,Xg,c,p)---(21)]]>Pg,c,p(1)=H(1,Ng,c,pnonsilent,xg,c,p,Xg,c,p)---(22)]]>Pg,c,p(2+)=1-Pg,c,p(0)-Pg,c,p(1)---(23)]]>对每个病人的同一基因内1个碱基发生突变的概率按突变类型进行降序排序，选取前两类d1、d2；公式(24)计算每个病人每种基因选取类别组合的概率值；Pg,p(d1,d2)=Πd=1ncPg,d,p(0)ifd1=0,d2=0Pg,d1,p(1)Πd=1d1-1Pg,d,p(0)Πd=d1+1ncPg,d,p(0)ifd1>0,d2=0Pg,d1,p(1)(Pg,d2,p(1)+Pg,d2,p(2+))Πd=d2+1d1-1Pg,d,p(0)Πd=d1+1ncPg,d,p(0)ifd1>0,0<d2<d1Pg,d1,p(2+)Πd=d1+1ncPg,d,p(0)ifd1>0,d2=d10ifd2>d1---(24)]]>在非同义突变区域内的每种突变都会影响到基因表达，即影响到蛋白质的合成；但null和indel突变对蛋白质合成的影响最大，本筛选方法对每种突变组合方式赋予一定的分值；数值的具体计算方式如公式(25)，为凸显null和indel突变对基因表达的影响Snull值的设定如公式(26)：Sg,p(d1,d2)=0ifd1=0,d2=0Snull-log10Pg,d1,p(1)ifd1>0,d2=0Snull-log10Pg,d1,p(1)-log10Pg,d2,p(1)ifd1>0,0<d2<d1Snull-log10Pg,d1,p(2+)ifd1>0,d2=d10ifd2>d1---(25)]]>Snull=0ifd1<nc3ifd1=nc---(26)]]>利用卷积计算每个基因的每个分值对应的概率Pg,p(S=x)=⊗d1=0nc⊗d2=0ncPg,p(d1,d2)δ(x-Sg,p(d1,d2))---(27)]]>Pg(S=x)=⊗p=1npPg,p(S=x)---(28)]]>使用公式(29)计算样本突变数据的每个基因的分值，其中Emin为最小效应值取1.25，目的为降低使用背景突变率的不确定性，得到每个基因的Sgobs=Σp=1npSg,pDg,pEmin---(29)]]>公式(29)中，Dg,p为按样本突变数据每个病人每个基因选取的两种突变类型矩阵；使用假设检验计算每个基因是驱动基因的p值，其中零假设是：基因g是癌症的驱动基因，则基因g的p值计算如公式(30)、(31)：Pg(S≥Sgobs)=∫Sgobs∞Pg(S=x)dx---(30)]]>Pg(S≥Sgobs)=1-∫0SgobsPg(S=x)dx---(31)]]>在假设检验中，由于p值只能控制发生第一类错误的概率，因此最终使用错误发现率(FDR)筛选基因；错误发现率计算方法：对所有p值进行升序排序：P(1)≤P(2)≤…≤P(m)，使用公式(32)计算每个p值对应的FDR值，当基因的fdr≤0.1时，认为该基因为驱动基因；fdr=mip---(32).]]>