本发明属于癌症医学领域，涉及一种基于基因突变频率的癌症驱动基因的筛选方法。
背景技术：
癌症是由体细胞突变引起。致癌因子诱发正常细胞进行变异，细胞中某些特定的基因位点发生突变时，会极大的促进癌症细胞的产生和发展，这些特殊的基因称为驱动基因。驱动基因突变在癌症发展过程中起着关键性的作用，是决定癌症的最主要原因。
在肿瘤细胞中并不是检测到的所有突变基因都是驱动基因。与驱动基因对应的是乘客基因，乘客基因突变对癌症的产生与发展的影响非常小。因此如何有效的识别癌症驱动基因仍是当前癌症研究中亟待解决的难题。癌症驱动基因的筛选对癌症靶向治疗药物的研发，癌症的预防、早期检测和诊断、分期分型以及康复治疗都有重要作用。
借助高通量基因组测序技术，可以快速准确的找到癌症基因序列中碱基组成与结构的异常变异，将肿瘤中的基因突变转化成数据，使用数据处理技术筛选癌症驱动基因。已有的使用癌症突变数据筛选癌症驱动基因相关研究中，多是基于突变频率或突变功能的方法。现有的癌症驱动基因相关专利，申请号：cn201510111810.9公开了“一种基于生物网络的癌症驱动基因的筛选方法”其中使用的是构建癌症生物分子网络筛选驱动基因。此方法是通过比较正常细胞与肿瘤细胞之间蛋白质的差异筛选突变基因。但蛋白质的合成不仅与基因相关，还与离子环境，ph值等细胞环境相关，因此仅通过比较蛋白质的差异无法精确的反应基因的具体突变情况和基因表达情况。基于突变频率的方法最先用于驱动基因研究。但基因突变频率不仅与致癌因子有关，基因的表达水平、复制时间、染色体状态等基因的固有特征也能导致各基因突变频率不同。因此在筛选癌症驱动基因时需考虑这些因素。申请号：cn201310284338.x公开了“一种检测非小细胞肺癌驱动基因突变谱的方法及试剂盒与应用”其中使用的是检测突变基因频率的方法判断驱动基因。此方法仅仅依赖基因的突变率，并未考虑基因的一些固有属性同时可以影响到基因突变频率。
技术实现要素：
本发明要解决的技术问题是提供一种基于基因突变频率的癌症驱动基因的筛选方法，使用聚类算法对影响基因突变协变量进行聚类分析，然后使用统计算法对聚类后的基因突变数据进行处理，筛选出癌症驱动基因。
本发明的技术方案：
基于基因突变频率的癌症驱动基因的筛选方法，步骤如下：
(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)数据预处理：对初始肿瘤基因突变数据进行整理。将突变数据整理成
将覆盖区域表整理成
将协变量表整理成矩阵vv,g，其中v为协变量，g为基因。将vv,g进行标准化得到zv,g，即用公式(1)将各个协变量数据转化成均值为0、方差为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中所有点的平均距离，遍历所有其他簇，找到最近的这个平均距离，记作bi。
计算所有基因g的轮廓系数，求出平均值即为当前聚类的整体轮廓系数，挑选最大的轮廓系数对应的聚类类别数k。
2)然后在基因所属的类里使用假设检验算法选出每个基因的邻近基因。其中零假设为基因i为基因g的邻近基因，数据采用背景突变数据。
统计背景突变的数据
同一基因内，由于碱基处在同一环境下，每个碱基发生突变的概率相同，所以基因内n个碱基发生n个突变的概率分布属于二项分布。若基因i是基因g的邻近基因，即它们的突变属于同一个二项分布，则基因i和基因g的突变数据服从beta-二项分布。以此做假设检验，零假设为基因i是基因g的邻近基因，p值的大小为公式(5)中的qi,g，在同类协变量的基因中筛选出p值大于0.05的基因作为基因g的邻近基因zg。
hc为beta-二项分布h的连加，具体的计算如下：
其中α＝n2+1，β＝n2-n2+1。γ为gamma函数。
筛选出每个基因的邻近基因，统计基因和所有邻近基因的背景突变碱基数xg和每种突变所在区域碱基总数xg。
公式(9)(10),i∈zg基因i是基因g的邻近基因
(4)平均突变数据，计算每个突变位点的背景突变数据。统计所有样本中每个病人，每种突变类别的总突变数，根据突变频率计算每个基因、每个病人、每种突变类型对应的背景突变数据。
忽略突变影响，统计每个基因、每个病人、每种突变类型的突变数据
统计每种突变类别的突变数据
统计所有突变位点的突变数据
统计每个病人的突变数据
计算每个基因、每个病人、每种突变类型对应的背景突变数据xg,c,p及背景区域碱基总数xg,c,p：
(5)筛选驱动基因：基于突变概率及突变类型设计一种突变分值，计算样本突变数据每个基因的总分值，使用假设检验算出每个基因是驱动基因的p值，算出对应的错误发现率，根据错误发现率筛选出驱动基因。
在基因内部，nonsilent区域和所在的背景突变区域的突变数据服从beta-二项分布，所以非同义突变区域里有0个，1个碱基发生突变的概率为
对每个病人的同一基因内1个碱基发生突变的概率
在非同义突变区域内的每种突变都会影响到基因表达，即影响到蛋白质的合成。但null和indel突变对蛋白质合成的影响最大，本文对每种突变组合方式赋予一定的分值。数值的具体计算方式如公式(25)，为凸显null和indel突变对基因表达的影响snull值的设定如公式(26)：
利用卷积计算每个基因的每个分值对应的概率
使用公式(29)计算样本突变数据的每个基因的分值，其中emin为最小效应值取1.25，目的为降低使用背景突变率的不确定性，得到每个基因的
公式(29)中，dg,p为按样本突变数据每个病人每个基因选取的两种突变类型矩阵。
使用假设检验计算每个基因是驱动基因的p值，其中零假设是：基因g是癌症的驱动基因，则基因g的p值计算如公式(30)、(31)：
在假设检验中，由于p值只能控制发生第一类错误的概率，因此最终使用错误发现率(fdr)筛选基因。
错误发现率计算方法：对所有p值进行升序排序：p(1)≤p(2)≤…≤p(m)，使用公式(32)计算每个p值对应的fdr值，当基因的fdr≤0.1时，认为该基因为驱动基因。
所属步骤(1)中的对细胞dna进行高通量测序，具体包括但不限于illumina的测序仪。
所属步骤(1)中的对dna序列处理的工作环境为linux系统具体包括但不限于ubuntu系统。
所属步骤(1)中的突变检测工具具体包括但不限于mutect2或varcan2。
所属步骤(1)中的位点注释工具具体包括但不限于annovar、oncotator。
所属步骤(2)中的整理数据工具具体包括但不限于r。
所属步骤(3)中的背景突变计算中，若测序时只对外显子组进行测序导致非编码区的突变数据过少失去统计意义，可以忽略非编码区突变数据，仅考虑同义突变数据。
本发明所述方法中，从得到病人的肿瘤和正常细胞dna，进行一系列数据处理，最终根据得到的每个基因是驱动基因的错误发现率筛选驱动基因。在此步骤中考虑了影响基因突变的自然因素，借助聚类算法对基因进行分类，考虑了相同环境下每种碱基的突变概率不同对基因表达的影响不同，因此对突变分类时分成7类，最后使用统计算法筛选驱动基因。此方法考虑了多种基因突变及基因表达的影响因素，最终得到的错误发现率也能精确地表达基因是否是癌症的驱动基因。
本发明的有益效果在于本发明对测序平台没有限制，对数据处理软件没有限制，有从初始细胞到最终得到驱动基因的详细过程，使用灵活方便，可移植性强。
本发明方法不仅和现代高通量测序技术及dna数据处理软件紧密结合，而且使用经典的聚类算法和统计方法，与影响基因突变的生物因素相结合使癌症驱动基因的筛选更加精确，对新型抗癌药物研发和癌症临床诊疗中都具有重要意义。
具体实施方式
以下结合发明内容详细说明本发明的具体实施方式:
下面详细介绍了肺癌的驱动基因筛选过程，其中包含数据提取过程和数据处理过程，最终筛选出驱动基因。本发明通过以下案例进一步说明本发明的用途和使用方法，但本发明并不受限于此。
1.实验材料与环境配置
177位肺癌患者的肿瘤细胞和正常细胞dna。
软件系统：linux系统
数据处理需求软件：bwa、samtools、picard、gatk、annovar、oncotator下载于工具官方网站。
r：开源计算平台，下载于r官方网站
2.实验方法
(1)肿瘤基因突变数据获取：对177位肺癌患者的肿瘤细胞和正常细胞的dna进行高通量测序，其中每种样本数据包含2个dna序列，包含多个短读，数据格式为fastq格式。初始数据包含3个信息：短读的id及测序仪信息、短读具体序列、碱基的测序质量。在linux系统下对初始序列进行处理。详细数据处理步骤如下：
(1.1)dna序列比对
应用bwa将测序短读与人类基因组hgl9进行比对，得到每段短读所在染色体上的具体位置。生成文件为sam文件。
(1.2)比对结果处理
使用picard软件对sam文件中短读按染色体号、同一染色体按位点顺序从小到大进行排序，然后标记由pcr过量扩增出来的完全相同的序列。因为过量扩增的短读并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由pcr扩增所形成的重复序列。生成文件为bam文件。
(1.3)突变检测
使用gatk对bam文件与hg19和dbsnp_138.hg19进行比对，得到插入删除位点和单核苷酸突变位点。生成具有突变结果的vcf文件。突变结果包含突变所在的染色体位点和突变前后的碱基类别。
(1.4)突变注释
使用oncotator和annovar对突变结果进行注释，得到突变的类型和影响和突变位点所在的基因名。
最后将这些数据整理成数据集突变数据、覆盖区域、协变量。表格如下：
表1突变数据
表2覆盖区域
表3协变量
(2)数据预处理：对初始突变数据进行整理。将突变数据整理成
(3)使用协变量和背景突变数据筛选每个基因的邻近基因：首先使用k-means算法对vv,g进行聚类，同时使用公式(2)计算所有基因g的轮廓系数，求出平均值即为当前聚类的整体轮廓系数。最终聚类结果的类别数为100，得到每个基因的类别。轮廓系数计算方法如下：
si＝(bi-ai)/max(bi,ai)(2)
然后在基因所属的类里筛选出碱基突变概率相同使用假设检验算法选出每个基因的邻近基因。
统计背景突变的数据：
由于基因i和基因g的突变数据服从beta-二项分布。以此做假设检验，零假设为基因i是基因g的邻近基因，p值的大小为公式(5)中的qi,g，具体计算方式如下列公式。最终在同类协变量的基因中筛选出p值大于0.05的基因作为基因g的邻近基因zg。
hc计算如下：
其中α＝n2+1，β＝n2-n2+1。γ为gamma函数。
筛选出每个基因的邻近基因，统计此基因和所有邻近基因的突变碱基数xg和每种突变所在区域碱基总数xg。
(4)使用统计方法，将每人每个基因每种突变类型的邻近基因的n值和n值进行平均，消除个体间的差异。
统计每种突变类别的突变数据
统计所有突变位点
统计每个病人的突变数据
求平均值，忽略个体间的差异
(5)筛选驱动基因：计算突变分值，使用假设检验算出每个基因是驱动基因的p值、错误发现率，根据错误发现率筛选出驱动基因。
首先计算
对不同突变类型相同病人和基因的
计算每种突变组合方式的突变分值。数值的具体计算方式如公式(25)，snull值的设定如公式(26)：
计算每个基因的每个分值对应的概率
使用公式(29)计算样本突变数据的每个基因的分值，其中emin取1.25，得到每个基因的
公式(29)中，dg,p为按照样本突变数据每个基因每个病人选取的两种突变类型矩阵。
使用假设检验计算每个基因是驱动基因的p值，其中零假设是：基因g是癌症的驱动基因，则基因g的p值计算如公式(30)、(31)：
计算错误发现率：对所有p值进行升序排序：p(1)≤p(2)≤…≤p(m)，使用公式(32)计算每个p值对应的错误发现率值，当基因的fdr≤0.1时，认为该基因为驱动基因。
3实验结果
最终得到错误发现率值低于0.1的基因为tp53，keap1，pten,cdkn2a,mll2，nfe2l2，rb1，fbxw7，hla-a,notch1，pik3ca,or2g6。通过查阅相关文献，这些基因变异会对极大提高肺癌细胞的生存能力，与细胞的生殖规律和科学主流认知相符合。