一种突变候选基因的筛选方法
阅读说明：本技术 一种突变候选基因的筛选方法 (Screening method of mutant candidate gene ) 是由 刘珍 刘志岩 王海宁 于 2021-08-11 设计创作，主要内容包括：本发明涉及基因组测序分析技术领域,具体涉及一种突变候选基因的筛选方法,通过使用已有公共正常样本的数据构建正常突变分类模型,进行模拟突变发生,而不是仅仅使用已知突变位点以及总体的突变频率,避免直接使用数据结论造成的结果偏差。同时当研究样本数量很少或不存在对照样本时,无法构建突变分类模型,也可直接使用正常突变分类模型,进行突变模拟和计算。该方法可快速并准确地从大量的非同义突变数据中,寻找并筛选出具有重要意义的主效基因,以便进行下一步研究验证。(The invention relates to the technical field of genome sequencing analysis, in particular to a screening method of mutant candidate genes. Meanwhile, when the number of research samples is small or no control sample exists, a mutation classification model cannot be constructed, and a normal mutation classification model can be directly used for mutation simulation and calculation. The method can quickly and accurately find and screen out major genes with important significance from a large amount of non-synonymous mutation data so as to carry out the next research and verification.)
技术领域
本发明涉及基因组测序分析技术领域，尤其涉及一种突变候选基因的筛选方法。
背景技术
随着二代测序技术的发展与成熟，测序价格持续降低，二代测序已经越来越多的应用于人类疾病的研究，其中以DNA测序最为广泛。但随着测序深度的增加，越来越多的非同义突变在研究中被发现，在大量的非同义突变中如何有效的筛选出有意义的突变和基因成为研究的难点。
目前解决的办法主要有：1)通过现有的已知的公共数据库(dbSNP、gnomAD等)筛选过滤，但由于公共数据库人群数量限制，只能去除少部分非同义突变，虽然突变数量能够减少，但无法准确的从剩余突变中筛选出更有意义的基因；2)采用对照研究去除非同义突变，这种方法仅适用于具有对照组织的疾病(比如肿瘤)，无法适用于自身免疫性疾病以及更多未知的复杂性疾病，此外即使能够采用对照去除，但对于非特异性肿瘤以及当研究样本数量较大时，仍有大量突变数据，难以准确筛选出主效基因；3)根据目前已有的研究计算的背景非同义突变率，然后对基因进行统计分析，找出主效基因，但这种方法虽然能够评估背景突变频率，但由于人群以及不同疾病存在的巨大差异性，导致直接使用背景突变频率进行计算差异性较大，从而使结果偏离。
发明内容
为解决上述技术问题，本发明的目的在于提供一种突变候选基因的筛选方法，所述突变候选基因的筛选方法可快速并准确地从大量的非同义突变数据中，寻找并筛选出具有重要意义的主效基因，以便进行下一步研究验证。
为达到上述技术效果，本发明采用了以下技术方案：
一种突变候选基因的筛选方法，包括步骤S1-S4：
S1：以公共数据库中的正常健康人群样本作为对照样本，以基因为单位，根据已知的非同义突变位点，结合两侧碱基类型进行分类，构建突变分类模型，作为正常突变分类模型数据。
S2:根据所需研究疾病的对照样本数据作为研究样本，以基因为单位，将发现的突变位点，结合两侧碱基类型进行分类，与正常突变分类模型数据结合，作为最终正常突变分类模型数据。
S3：对每个基因，遍历每个编码区碱基位点，根据最终正常突变分类模型数据随机预测其是否发生突变以及突变类型；统计基因上非同义突变数量，计算模拟环境下真实疾病样本群体特异性基因的背景非同义突变频率。进行多次模拟后，获取多次模拟下基因的背景非同义突变频率。
S4：将多次模拟下背景非同义突变频率以及实际基因发生的非同义突变数量，运用二项分布计算是否为主效基因。
进一步地，所述步骤S1中的数据库为1000g数据库、gnomAD数据库、EXAC数据库、esp6500数据库中的至少一种。优选为1000g数据库。与现有技术相比，本发明的有益效果为：
第一方面，本发明提供的一种突变候选基因的筛选方法使用已有公共正常样本的数据构建正常突变分类模型，进行模拟突变发生，而不是仅仅使用已知突变位点以及总体的突变频率，避免直接使用数据结论造成的结果偏差。同时当研究样本数量很少或不存在对照样本时，无法构建突变分类模型，也可直接使用正常突变分类模型，进行突变模拟和计算。
第二方面，本发明提供的一种突变候选基因的筛选方法由于算法利用并结合研究样本的真实突变(采用对照样本的非同义突变)，与正常突变分类模型合并，构建最终正常突变分类模型，模拟出来的突变数据能够准备反映特定疾病环境下患病人群真实突变情况。
第三方面，本发明提供的一种突变候选基因的筛选方法通过构建突变分类模型，将突变碱基位点，并结合两侧碱基位点，作为一个单位进行汇总统计，能够结合实际序列环境反映突变的实际发生环境。
第四方面，本发明提供的一种突变候选基因的筛选方法可以进行多次模拟，根据多次模拟计算背景非同义突变频率，避免样本数据量较小或偶然模拟造成的实验偏差，使结果更稳定可信。
附图说明
图1为本发明提供的最终正常突变分类模型数据的构建流程图；
图2为本发明提供的突变分类模型模拟流程图；
具体实施方式
下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例仅用于更加清楚地说明本发明的技术方案，因此只作为示例，而不能以此来限制本发明的保护范围。
如图1-2所示，本发明提供的一种突变候选基因的筛选方法包括步骤S1-S4，其中：
步骤S1:使用1000G正常健康人群样本数据，利用已知的非同义突变，结合两侧碱基，将每个突变点根据不同的改变进行分类统计，构建正常突变分类模型。
具体为：假设某基因编码区某一位点C(左侧碱基)-G(位点)-G(右侧碱基)在1000G正常人类样本中此位点发生非同义突变，其中发现带有G>A的改变样本人数为13人，发现带有G>T的改变样本数为8人，发现带有G>C的改变样本数为4人，未发生改变的样本数为5人，根据碱基可能的四种改变方式统计此位点分类信息，即：CGG>CAG13(人)，CGG>CTG8(人)，CGG>CCG4(人)，CGG>CGG5(人)。
将所有正常样本发生的非同义突变，根据上述规则，进行分类统计汇总。如果一个基因上有两个位于不同位置的碱基位点，但碱基相同，其两侧碱基也相同，比如基因TP53上两位点：
位点1：chr17-7579707G(左侧碱基)-T(位点)-T(右侧碱基)，
位点1突变分类信息列表[GTT>GAT3，GTT>GTT17,GTT>GCT10，GTT>GGT0]，
位点2：chr17-7579710G(左侧碱基)-T(位点)-T(右侧碱基)，
位点2突变分类信息列表[GTT>GAT1,GTT>GTT24,GTT>GCT0，GTT>GGT5]。
将其相加，最终突变分类信息列表为[GTT>GAT3+1,GTT>GTT17+24,GTT>GCT10+0，GTT>GGT0+5]。
步骤S2：对于研究样本，利用对照样本中的非同义突变，构建突变分类模型，并与正常突变分类模型相加，增加模型的稳健性。
具体为：假设某基因编码区某一位点C(左侧碱基)-G(位点)-G(右侧碱基)在研究对照样本中发现非同义突变，其中发现带有G>A的改变样本人数为8人，发现带有G>T的改变样本数为2人，发现带有G>C的改变样本数为0人，未发生改变的样本数为20人。则此位点根据碱基可能有四种改变方式进行统计，即：CGG>CAG8(人),CGG>CTG2(人),CGG>CCG0(人),CGG>CCG20(人)。
将所有对照样本发生的非同义突变，根据上述规则，进行分类统计汇总。如果一个基因上有两个位于不同位置的碱基位点，碱基相同，其两侧碱基也相同，比如基因TP53上两位点：
位点1：chr17-7579707G(左侧碱基)-T(位点)-T(右侧碱基)，
位点1突变分类信息列表[GTT>GAT3，GTT>GTT17,GTT>GCT10，GTT>GGT0]，
位点2：chr17-7579710G(左侧碱基)-T(位点)-T(右侧碱基)，
位点2突变分类信息列表[GTT>GAT1,GTT>GTT24,GTT>GCT0，GTT>GGT5]。
将其相加，最终突变分类信息列表为[GTT>GAT3+1,GTT>GTT17+24,GTT>GCT10+0，GTT>GGT0+5]。
将研究样本的突变分类模型数据与1000G正常人类样本的正常突变分类数据模型进行合并，即：
CGG>CAG13+8(人),CGG>CTG8+2(人),CGG>CCG4+0(人),CGG>CGG5+20(人)。
步骤S3：利用每个基因上的突变分类模型，将基因编码区每个碱基进行遍历，随机获取突变改变(突变分类模型中次数越大说明取到的可能性越大，但并不是百分之百取到)，预测是否发生突变。
具体为：设某基因编码区长度为5bp，即含有5个碱基-c(内含子)AGTCAg(内含子)，遍历每个碱基，即：将一个碱基以及其两侧碱基作为一个单位，单位列表为[cAG,AGT,GTC,TCA,CAg]。
根据构建的正常突变分类模型，包含一个碱基以及两侧碱基组成的一个单位，其发生四种碱基改变的可能性，比如：
[cAG>cAG10,cAG>cTG3,cAG>cCG1,cAG>cGG1,
AGT>AAT0,AGT>ATT3,AGT>ACT2,AGT>AGT10,
GTC>GAC1,GTC>GAC9,GTC>GCC2,GTC>GGC3,
TCA>TAA3,TCA>TTA3,TCA>TCA8,TCA>TGA1,
CAg>CAg12,CAg>CTg1,CAg>CCg0,CAg>CGg2,]
将遍历的单位在正常突变分类模型中进行搜索，寻找其碱基完全一致的单位，比如第一个cAG找到[cAG>cAG10,cAG>cTG3,cAG>cCG1,cAG>cGG1]。
将[cAG>cAG10,cAG>cTG3,cAG>cCG1,cAG>cGG1]作为一个列表，其中10、3、1、1代表此种改变在列表中出现的次数；随机从列表中进行取值,比如取到cGG。
2)如果在基因的编码区中cAG>cGG为非同义突变类型，记为1,其他类型记为0；统计预测模拟下的基因背景非同义突变频率；进行多次模拟，计算每次模拟下的背景非同义突变频率。
根据上述参考基因编码区长度为5bp，即基因编码碱基总长度＝5；根据上述遍历模拟突变发生，5个碱基共有2个非同义突变发生，即：
基因背景非同义突变频率＝2/5；
根据上述流程，可以进行多次突变模拟，每次计算背景非同义突变频率。
步骤S4：根据多次模拟计算的背景非同义突变频率，同时计算实际基因的非同义突变频率，即基因实际发生的非同义突变数除以基因编码区长度，运用二项分布计算P值，从而判定基因是否为主效基因。
根据上述模拟，5次模拟计算的背景非同义突变频率列表结果为[0.3,0.1,0.5,0.2,0.2]，基因实际发生的背景突变频率(假设为0.4)。
通过实际频率与模拟计算频率比较，发现有1次模拟背景突变频率大于0.4，进行二项分布计算:
其中，k为实际背景突变频率大于模拟计算的频率的次数，示例中为1次；
n为模拟次数，示例中为5次；
p为模拟情况下实际背景突变频率大于模拟的背景突变频率发生的概率，固定概率为0.5。
如果计算的P值小于0.05，则提示此基因为重要或主效基因，即实际发生的非同义突变数是一个选择性(比如疾病的发生)造成的而非随机产生。
将所有计算的基因P值排序，P值越小则提示基因在疾病的发生中作用越大。
以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。本发明未详细描述的技术、形状、构造部分均为公知技术。