本发明涉及生物学技术领域，具体涉及一种基于外显子测序技术的人单基因遗传疾病检测生物信息分析方法。
背景技术：
按照国内医学遗传学界的传统，遗传性疾病一般可分为染色体病、单基因病、线粒体遗传病、多基因病和体细胞遗传病五大类。其中，单基因遗传病是由单一基因突变导致的，是造成出生缺陷的主要原因之一，其遗传规律符合孟德尔遗传方式，所以也称孟德尔式遗传病。目前，已经报道的单基因病(或性状)超过12000种(www.omin.org)，而ensembl网站上所统计的人类蛋白质编码基因为20469个(www.ensembl.org)。在单基因遗传病的致病基因定位和克隆研究中，近二十多年所采用的连锁分析和定位克隆技术已被公认为是最有效、最准确的方法之一。然而，自2009年9月，一种全新的寻找疾病致病基因或易感基因的方法-全外显子组测序(wholeexomesequencing，wes)开始崭露头角，发挥了越来越多的重要作用，并且已经被用于临床基因诊断。
外显子(exon)是真核生物基因组dna序列的一部分，可被转录成为成熟的mrna，并在翻译过程中被表达为蛋白质。外显子组(exome)是一个物种基因组中的全部外显子的总和，人类的外显子组序列只占全基因组序列的1％。外显子组测序(exomesequencing)技术是一种利用目标序列捕获技术将全基因组中的所有外显子区域dna序列捕获并富集后再进行高通量测序获得蛋白质编码区域dna的基因组分析方法。
利用外显子测序技术对人单基因遗传病进行检测，并对测序结果进行生物信息分析，可以检测出单基因遗传病的突变位点，有助于对受检者的患病风险进行评估。单基因遗传病的发病率高，而且对人类健康造成了较大威胁，不仅对患者的健康造成了严重危害，也给家庭和社会带来了沉重的精神和社会负担。因此，利用外显子测序技术对人单基因遗传病的检测数据进行准确全面的分析具有重要的社会意义。目前，还没有形成系统的利用外显子测序技术对人单基因遗传病的检测数据进行分析的方法。
技术实现要素：
本发明的目的在于为克服上述现有技术的不足之处而提供一种基于外显子测序技术的人单基因遗传疾病检测生物信息分析方法。
为实现上述目的，本发明采取的技术方案为：
一种基于外显子测序技术的人单基因遗传疾病检测生物信息分析方法，包括以下步骤：
s1.突变有害性位点筛选；
s2.对人单基因遗传疾病检测的样本关系进行筛选；
s3.基因功能筛选。
进一步地，所述步骤s1中突变有害性位点筛选，是基于突变位点的筛选方法筛选出有害性的突变位点，并对筛选得到的突变位点进行注释分析。
进一步地，所述步骤s2中对人单基因遗传疾病检测样本关系进行筛选主要包含：常染色体隐形遗传模型筛选、常染色体显性遗传模型筛选、新生突变基因筛选及共有突变基因筛选，其中，常染色体隐形遗传和常染色体显性遗传筛选需要提供家系模型，新生突变筛选需要提供双亲的信息。
进一步地，所述步骤s3中基因功能筛选的具体方法为：在对样本关系筛选后，将筛选到的突变基因并集进行go和kegg富集分析，以发现疾病与基因功能之间的关联。
本发明的有益效果：(1)本发明能够快速清晰地对人外显子测序数据进行有害突变位点的筛选；(2)本发明为了全面进行基因筛选，将所检测到的snp以及indel等基因组变异和外部数据库进行注释分析，以确定变异的基因组位置、变异频率、蛋白有害性、基因型杂合性以及所在的功能通路等信息；(3)本发明利用常染色体隐形遗传模型筛选、常染色体显性遗传模型筛选、新生突变筛选及共有突变基因筛选，从而确定候选基因，并对候选基因进行go和kegg富集分析，为确定疾病相关的候选基因提供了强有力的证据。本发明检测数据生物信息分析步骤简洁清晰、分析内容全面，对遗传疾病的筛选和预防具有重要的指导意义。
附图说明
图1为本发明一种基于外显子测序技术的人单基因遗传疾病检测生物信息分析方法框架图。
图2为本发明实施例中筛选出的候选基因分子功能富集结果图，其中description：goterm的功能描述；generatio：注释到某个goterm的差异基因数以及占总差异基因数目的百分比；bgratio：注释到某个goterm的基因数以及占总基因数目的百分比；pvalue：p值；p-adjust：p值多重检验校正后的值。
图3为本发明实施例中筛选出的候选基因的go分类柱状图。
图4为本发明实施例中筛选出的候选基因的kegg分析热图。
具体实施方式
以下通过特定的具体实例并结合附图说明本发明的实施方式，本领域技术人员可由本说明书所揭示的内容轻易地了解本发明的其它优点与功效。本发明亦可通过其它不同的具体实例加以施行或应用，本说明书中的各项细节亦可基于不同观点与应用，在不背离本发明的精神下进行各种修饰与变更。
实施例1
对肾发育不全的人样本(其中包含3个正常样本和1个患病样本)外显子测序数据进行分析，参考基因组为hg19，其中正常样本分别为：s1881-16、s1882-16及s0405-16，患病样本为：s0401-16。家系信息为：s1881-16、s1882-16、s0405-16及s0401-16。如图1所示，一种基于外显子测序技术的人单基因遗传病检测生物信息分析方法，具体步骤如下：
步骤s1，突变有害性位点筛选：具体的筛选步骤如下：
1.过滤千人基因组数据库(人群中频率大于0.01)的变异位点；
2.过滤数据库exac、esp6500中频率大于0.01的变异位点；
3.保留突变在外显子区(exonic)或剪接位点区(splicing)的变异；
4.过滤同义突变，保留非同义突变；
5.突变保守性筛选。依据sift、polyphen、mutationtaster、cadd这4个软件，要求这4个软件中，有分值的软件中至少有一半支持该位点可能有害，该位点才能被保留。siftscore＜0.05，polyphen和mutationtaster大于0.85被认为是保守的(有害的)。
得到的样本突变位点筛选的统计结果如表1所示：
表1.突变筛选统计表
其中：total：所有样本的变异总数；g1000：过滤千人基因数据库中频率大于0.001所剩的变异数目；exac：使用上步结果，过滤exac频率大于0.01所剩的变异数目；esp6500：使用上步结果，过滤esp6500频率大于0.01所剩的变异数目；exonic/splicing：使用上步结果，保留外显子或剪切位点区域的变异数目；synonymous：使用上步结果，过滤同义突变所剩下变异数目；deleterious：使用上步结果，保留有害(保守)变异(此步骤对indel无效)。
接着，对上述筛选得到的变异位点进行注释分析：利用annovar软件将所检测到的单核苷酸多态性(snp)以及插入缺失标记(indel)等基因组变异和外部数据库进行注释分析，以确定变异的基因组位置、变异频率、蛋白有害性、基因型杂合性以及所在的功能通路等信息。
步骤s2，对人单基因遗传疾病检测的样本关系筛选，具体包括：
s2.1常染色体隐形遗传模型(需提供家系模型)
家系：是指某一家族的成员数目、亲属关系以及有关遗传形状或遗传病在该家系中分布情况。在此，成员主要包括父亲、母亲以及二者所生的子代，子代的数量为一个或两个以上。
常染色体隐形遗传致病主要包含两种情况：基因纯合变异和复杂杂合变异。
1.纯合模式候选位点筛选的具体标准如下：1)保留家系所有患者共有的纯合突变；2)此突变在家系正常人为杂合突变或没有突变位点。如：双亲均为aa或者aa+aa型，则筛选后代患者中aa基因型。
2.复杂杂合模式筛选是指如果单基因病在家系中为隐形模式遗传，则保留患者和正常人都不为纯合突变的位点，且要求一个基因在患者中至少有两个杂合突变位点，且患者此基因上的突变位点分布不能与正常人(此基因)的突变位点分布一样，也不能是任何一个正常人(此基因)突变位点的子集。例如：正常的亲本基因型分别为aa1和aa2，则筛选患者中的a1a2基因型为复合杂合基因型。常染色体隐形遗传模型筛选得到的结果如表2所示：
表2常染色体隐形遗传模型筛选结果
表2、3、6中，chr：染色体；start:snp/indel起始位点；end:snp/indel终止位点；ref:参考碱基；mut:突变碱基；variation_type:突变类型；structure_type：所在基因位置类型；structure_gene：所在基因位置名称；function_type：突变功能类型；function_gene：突变功能注释；s1881-16_type、s1882-16_type、s0454-16_type、s0416-16_type均表示对应样本的基因型；s1881-16_read、s1882-16_read、s0454-16_read、s0416-16_read均表示对应样本的reads数；omim_phenotype：omim：人类孟德尔疾病信息，omim表型；phenotype_omim_number：表型id；inheritance:遗传特点；phenotype_map_key：表型类型；avsnp147：该变异在dbsnp数数据库的id；cosmic70：癌症体细胞突变数据库，观察到的次数，以及观察到的癌组织，包括非编码突变。由此可知检测到的体细胞突变是否已被报导或观测到过，以及在哪些癌种中被报导、被报导的次数；all.2015_08：千人基因组数据库，参与人群包括afr(african)，amr(admixedamerican)，eas(eastasian)，eur(european)，sas(southasian)；esp6500siv2_all：通常采用0.01的标准进行过滤；exac*：exomeaggregationconsortium0.3版本，整合了60706个无亲缘关系个体的数据，这些个体来源于大量疾病研究和群体遗传学研究，能够用做严重疾病研究的参考数据库，目前数据库中包括65000exomeallelefrequencydataforall，afr(african)，amr(admixedamerican)，eas(eastasian)，fin(finnish)，nfe(non-finnisheuropean)，oth(other)，sas(southasian)，文献通常采用0.01的标准进行过滤；cl*：clindbn代表变异位点相关的疾病名称，clnacc代表变异在clinvar数据库中的accession号和版本号，clndsdb是疾病关联信息的数据库来源，clndsdbid是数据库中的编号；sift*：用来预测氨基酸改变是否影响功能，分值越小越可能“有害，如果得分小于0.05，被认为d(amaging)，否则被认为t(olerated)；polyphen2*：基于humandiv数据库和humanvar数据库用来预测非同义突变造成的氨基酸改变是否影响功能，数值越大越“有害”，包含hvar(孟德尔疾病诊断，d(probablydamaging，0.909-1)、p(possiblydamaging，0.447-0.908)和b(benign，0-0.446))和hdiv(在复杂表型可能的位点评估罕见的等位基因、d(probablydamaging，0.957-1)、p(possiblydamaging，0.453-0.956)和b(benign，0-0.452))；lrt*：预测snp导致蛋白结构和功能是否改变，包含d(deleterious)、n(neutral)和u(unknown)；mutationtaster*：分值越大，越可能是deleterious，包含a(disease_causing_automatic)、d(disease_causing)、n(polymorphism)和p(polymorphism_automatic)；mutationassessor*：预测是否有功能变异，越大越可能“有害”，包含功能性变异(h(high)和m(medium))和非功能性变异(l(low)和n(neutral))；fathmm*：根据疾病通路评估人类遗传性疾病突变，分值越小越可能“有害”，包含d(damaging，得分小于-1.5)和t(tolerated)；provean：预测snp或indel是否影响蛋白的功能，分值越小越可能“有害”，包含deleterious(得分小于-2.5)和neutral；vest3_score：采用机器学习方法来预测全基因组水平上非同义变异带来的功能显著性，得分越高越显著；cadd*：采用此方法进行预测，snp的有害性阈值为caddpred分值＞15，文章中常用10或15；dann_score：采用深度神经网络方法进行预测，大于0.96(0-1)表明是pathogenic，否则为benign；fathmm-mkl_coding*：类似于fathmm，采用mkl算法预测突变，得分越大越“有害”，大于0.5(0-1)表明是deleterious，否则为benign或neutral；metasvm*：采用svm方法预测突变，得分越大越“有害”，得分>0表明是deleterious，否则为tolerated；metalr*：类似于metasvm，利用逻辑回归(logisticregression)的预测方法，得分越高表明越“有害”，得分>0.5表明是deleterious，否则为tolerated；integrated_fitcons_score：fitcons预测分值，0-1，分值越大表示变异的重要性越强；integrated_confidence_value：fitcons预测分值可靠性评估，0表示非常显著(p<0.003)，1表示显著(p<0.05)，2表示有可能(p<0.25)，3表示预估(p≥0.25)；gerp++：通过定量替换情况识别多序列比对中的保守程度，数值越大越保守(-12.3-6.17)，分值大于2的位点认为是保守位点；phylop*:类似于gerp++，基于7个脊椎动物物种和20个哺乳动物物种的多序列比对得到位点的保守性分值，数值越大越保守；phastcons*：类似于phylop*，数值越大越保守(0-1)；siphy*：基于29种哺乳动物的多序列比对得到位点的保守性分值，数值越大越保守(0.0003-37.9718)；go_cc:go细胞组件，go_mf：go分子功能，go_bp：go生物进程，kegg_pathway：kegg通路。
步骤s2.2常染色体显性遗传模型(需提供家系模型),显性模式筛选是在突变位点过滤的基础之上，保留家系中患者常染色体共有的杂合突变(性染色体保留有突变的位点)。常染色体显性遗传模型筛选得到的结果如下表3所示：表3常染色体显性遗传模型筛选结果
步骤s2.3新生突变筛选，在没有患病史的家族中，由于自然突变等原因，后代往往会产生不来源于双亲的新生突变(denovomutation)。在这里，我们筛选获得正常双亲都没有，而患病后代特有的突变后，即为最终的新生突变。新生突变速率与疾病相关，编码区域的新生突变速率可以有效评估这些突变在患者群里中发生的偏好性。新生的snp(单核苷酸多态性)和indel(insertion-deletion，插入或缺失)的突变速率计算结果如表4所示：
表4新生snp突变速率
表5新生indel突变速率
表4和5中，proband表示患者编号，denovosnp表示新生snp数目，denovoindel表示新生indel数目，codingbase表示编码区域覆盖度4x以上的碱基数(单位bp)，mutationrate表示新生突变速率。
在获得新生突变后，利用annovar对突变进行注释，进一步地确定突变与蛋白功能之间的关系。
步骤s2.4共有突变基因筛选，对于散发样本，在过滤有害性样本的基础上，筛选样本间共有的突变基因，筛选比例为患者中至少20％共有且至少有2个患者共有，并且在90％以上的正常样本中不含有此基因的有害性突变。筛选得到的共有突变基因信息如表6所示：
表6共有突变基因信息
步骤s3基因功能筛选。疾病的发生往往与基因功能失调有关。基因功能筛选是基于步骤s2中显性遗传模式、隐形遗传模式、新生突变模式及共有突变模式进行突变基因筛选后，取上述各模式筛选到的突变基因并集进行go和kegg富集分析。具体包括：
步骤s3.1候选基因go富集分析：首先将上述各模式筛选到的突变基因取并集，将这些基因向go数据库(http://www.geneontology.org/)的各term映射，并计算每个term的基因数，得到具有某个go功能的基因列表及基因数目统计。然后应用超几何检验，找出与这个基因组背景相比，在基因中显著富集的go条目。候选基因的go富集结果以表格展示，点击表中的超链接可以查看详细结果，图2为候选基因分子功能(molecularfunction)的具体情况。候选基因的go分类柱状图如图3所示。
步骤s3.2候选基因kegg富集分析：以keggpathway为单位，应用超几何检验，找出与整个基因组背景相比，在基因中显著性富集的pathway。通过pathway显著性富集确定候选基因参与的最主要生化代谢途径和信号转导途径。图4为筛选得到的候选基因的kegg富集分析气泡图。
综上，本发明能够快速清晰地对人外显子测序数据进行有害突变位点的筛选；并将所检测到的snp以及indel等基因组变异和外部数据库进行注释分析，以确定变异的基因组位置、变异频率、蛋白有害性、基因型杂合性以及所在的功能通路等信息。此外，本发明利用常染色体隐形遗传模型筛选、常染色体显性遗传模型筛选、新生突变筛选及共有突变基因筛选，从而确定候选基因，并对候选基因进行go和kegg富集分析，为确定疾病相关的候选基因提供了强有力的证据。本发明检测数据生物信息分析步骤简洁清晰、分析内容全面，对遗传疾病的筛选和预防具有重要的指导意义。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。