本发明涉及基因组学高通量测序数据分析领域，具体包括对靶向基因文库测序数据进行质量控制、扩增子效率评估和过滤、基因组比对、突变识别及注释，进而结合已记录突变及病例资料完成统计分析，为肿瘤患者提供一整套非定制肿瘤基因突变检测分析解决方案，为肿瘤个体化医疗提供技术支持。
背景技术：
肿瘤在本质上是基因病。各种环境的和遗传的致癌因素以协同的方式引起dna损害，从而激活原癌基因和(或)灭活肿瘤抑制基因，加之凋亡调节基因和(或)dna修复基因的改变，继而引起表达水平的异常，使靶细胞逐步向癌细胞发生转化。被转化的细胞先呈现多克隆性的增生，经过一个漫长的多阶段的演进过程，其中一个克隆相对无限制的扩增，通过累积附加突变，选择性地形成具有不同特点的亚克隆(异质化)，从而获得浸润和转移的能力(恶性转化)，形成恶性肿瘤。
肿瘤基因检测是提取人体细胞内的遗传物质，通过测序检测人体内的肿瘤基因或易感基因，用于肿瘤的预防、诊断、预后预测、靶向用药、术后监控等。
靶向测序是对特定长度的pcr产物或者捕获的片段进行测序，分析序列中的变异。靶向测序能够根据不同需求，对目标区域进行高覆盖度的测序，还可以检测到低频突变。随着测序成本的减少以及人类功能基因组学研究的深入，靶向测序已经从研究机构走向临床，用于遗传筛查、疾病风险评估、肿瘤诊治以及精准用药等多个领域。
靶向测序数据分析的问题：第一，尚无靶向测序数据工具能结合已知突变数据库或病例资料进行统计分析，不能给出与疾病显著相关的突变分析结果。第二，一般分析软件只满足固定panel文库测序数据的分析，如illumina公司的trueseqamplicon，无法满足不同需求的用户的靶向文库分析。第三，现有方法中，没有评估扩增子的扩增效率并对引物序列进行修剪的操作，如果不做处理两者均会导致后续的突变分析得到的snv结果具有较高的假阳性，从而影响分析结论。
技术实现要素：
针对现有分析方法存在的缺陷，本发明为自主设计的靶向基因测序数据分析提供了一整套解决方案。
本发明目的是通过下述技术方案实现的：
一种肿瘤靶向基因测序数据解析方法，其特征在于，包括以下步骤：
步骤一：获取含有突变信息的读段序列，即高通量的测序数据；
步骤二：测序数据的质量控制，通过fastqc软件对步骤一获取的所有测序数据进行质量分析，获得测序数据质量控制报告，并过滤掉被报告为低质量的数据；
步骤三：统计不同扩增区域的扩增效率，删除扩增异常的数据；
步骤四：删除测序数据读段序列中的引物序列，即获得读段序列中真正的目标区域dna序列；
步骤五：将步骤四中得到的所有序列数据比对到靶向区域上，获得比对结果数据；
步骤六：应用突变识别工具从比对结果数据中检测所有突变；
步骤七：统计扩增区域内所有覆盖碱基的测序深度，根据突变位点的测序深度筛选可靠性高的突变；
步骤八：结合已注释在癌症相关数据库中的突变筛选差异显著的突变；
步骤九：结合病例的临床资料信息，对各种突变进行统计分析，识别与性状显著相关的胚系突变(germlinemutations)及体细胞突变(somaticmutations)；
步骤十：以图表的形式生成数据分析报告。
所述的测序数据来自包括illuminamiseq/hiseq在内的高通量测序平台的靶向扩增文库，靶向区域可定制，即在首次分析时提供所有靶位点基因组定位信息。
步骤二中所述的低质量数据是指单碱基平均测序质量得分低于20的测序数据。
步骤三中扩增区域异常数据获取及删除过程如下：
(1)将测序得到的读段数据比对到靶向参考基因组上；
(2)判断读段的两端序列所对应的扩增子引物序列是否来自同一个引物对，即允许前引物和后引物与测序片段的5’和3’端各有2个错配，去除不符合条件的读段序列；
(3)统计覆盖到每对扩增子对应的靶向扩增区域的读段序列，并应用扩增数来衡量和比较它们的扩增效率；
(4)当扩增数小于所有扩增子对应的扩增数均值的1/3时，则判断该扩增子所对应的扩增区域存在异常，并将其扩增出来的所有读段序列在分析中删除。
所述步骤五中的比对过程需要根据首次提供的靶向区域基因组位置信息，提取这些靶向区域在基因组中的核酸序列信息，并生成索引。
步骤七中对于所有突变位点，只将在该位点测序深度超过100×的病例纳入统计分析。
步骤八中差异显著的突变的筛选是基于目前已知的癌症相关数据库中存储的数据来进行的。
步骤九中的临床资料信息的关联分析，具体应用r统计软件中的卡方检验函数对患者的临床资料，包括年龄、性别、癌症tnm分期、肿瘤体积、肿块大小、是否有淋巴结侵袭、ki67恶性程度、组织亚型，进行关联分析，找出与特定基因突变发生相关的风险因素，即是否具有某一临床特征的病人倾向于发生哪些特异的突变。
所述参考序列是来源于ucsc公共数据库的参考肿瘤基因或人类参考基因组序列。
本发明的有益效果：本发明一种肿瘤靶向基因测序数据解析方法包括⑴获取含有突变信息的读段序列，即测序数据；⑵测序质量控制；⑶靶向区域扩增效率质量控制；⑷删除读段序列中的引物序列；⑸将读段序列与参考靶位点序列进行比对，获得比对上的读段序列；⑹识别读段序列中的突变；⑺根据突变位点的测序深度进行样本筛选；⑻结合已注释在癌症相关数据库中的突变筛选差异显著的突变；⑼结合病例的临床资料信息进行关联分析；⑽以图表的形式生成数据分析报告。满足非定制靶向文库的分析需求，结合已知突变数据库或病例资料给出与疾病显著相关的突变分析结果。该解析方法中将结合已知的参考突变数据库和患者的临床资料信息，应用不同的统计检验方法筛选出具有显著性差异的突变。突变数据库包括胚系突变数据库和体细胞突变数据库两类。其中，常用的胚系突变数据库包括千人基因组数据库(http://www.1000genomes.org/)和6万人exac人类外显子组整合数据库(http://exac.broadinstitute.org/)等。常用的体细胞突变数据库包括美国肿瘤基因组图谱tcga数据库(http://cancergenome.nih.gov/)和国际癌症基因组联盟icgc数据库(https://dcc.icgc.org/)等。通常需要用到四种对象，第一是突变人数比例，即携带突变的病人数量；第二，群体中的突变比例及群体等位基因频率；第三，纯合突变人数比例；第四，杂合突变人数比例。上述四种对象的数据拿到之后，我们就可以使用fisher精确检验统计的方法筛选差异显著的突变(即与肿瘤发生相关的基因突变)。本方法应用r统计软件中的卡方检验函数对患者的临床资料信息(包括年龄、性别、癌症tnm分期、肿瘤体积、肿块大小、是否有淋巴结侵袭、ki67恶性程度、组织亚型等信息)进行关联分析，找出在肺癌中与特定基因突变发生相关的风险因素，即是否具有某一临床特征的病人倾向于发生哪些特异的突变。一方面，本发明一种肿瘤靶向基因测序数据解析方法相较目前存在的方法更具有通用性。另一方面，本方法中特别针对用户定制的扩增子靶向文库的扩增效率进行全面评估，并对引物序列进行修剪，尽可能的保证不同扩增子的扩增效率保持在一个大致相同的水平，来规避由于不同扩增子的扩增效率导致的snv结果的假阳性问题。总而言之，本方法不仅可以对不同用户定制的不同扩增子文库进行通用的分析处理并获得结合临床信息的与疾病显著相关的突变分析结果，还自主增加了对于靶向文库扩增效率的全面评估和引物序列修剪步骤，使得整个分析方法在兼具新颖性的同时又提高了分析结果的可靠性。
附图说明
图1为本发明方法实现流程图。
具体实施方式
现有的高通量测序平台有多种，包括illuminanextseq、miseq和hiseq等。本发明中的实施例以illuminahiseq/miseq测序平台作说明。
本发明提供的方法适用于靶向dna或rna中突变检测，因此将分别以实施例作阐述。实施例中样本dna/rna提取、构建文库、高通量测序等均可利用现有技术进行。
实施例中未注明具体条件的，按照常规条件或制造商建议的条件进行；所用试剂或仪器未注明生产厂商的，均可以通过市面购买获得的常规产品。
实施例一：10例肺癌患者胸水样本靶向基因测序数据解析：
本实施例中的文库为10例肺癌患者胸水样本游离dna构建的靶向扩增子文库。文库构建的具体步骤如下：
(1)靶向基因的选择：选择肿瘤热点突变基因、原癌基因、抑癌基因、与靶向药物作用的基因，具体是abl1、egfr、gnas、mlh1、ret、akt1、erbb2、hnf1a、alk、erbb4、hras、notch1、smarcb1、apc、fbxw7、idh1、npm1、smo、atm、fgfr1、jak2、nras、src、braf、fgfr2、jak3、pdgfra、stk11、cdh1、fgfr3、kdr、pik3ca、tp53、cdkn2a、flt3、kit、pten、vhl、csf1r、gna11、kras、ptpn11、ezh2、tnnb1、gnaq、met、rb1、idh2这48个基因，作为我们研究的靶基因。
(2)游离dna的提取及定量：对于肺癌患者的胸水样本，我们先进行低速离心(3,000rpm)5分钟取上清，在进行高速离心(14,000rpm)10分钟取上清，得到了胸水样本中的游离dna(平均长度约为166bp)；并应用qbuit2.0(invitrogen公司)仪器进行定量。
(3)扩增子设计：通过在线引物设计工具designstudio，针对48个靶基因进行引物设计。最终，我们得到了覆盖48个靶基因全部外显子区域的2,158对扩增子，每对扩增子片段的大小约为150bp。由于不同靶基因的序列长度不同，我们的每对扩增子的片段大小又是几乎固定的，因此每个靶基因都对应着不同数目的扩增子引物对。靶基因与扩增子引物对数目的对应列表，见表1。
表1
(4)多重pcr扩增靶向基因的外显子：扩增子引物设计完成后，根据设计报告提供的引物序列，合成引物核酸，并以多重pcr的形式扩增靶向基因的全部外显子序列。
(5)illumina测序接头的连接和文库pcr扩增：对于上述扩增产物，我们连接illumina测序仪的测序接头。测序接头序列如下：
上游序列：5'p-nnn……nnngatcggaagagcacacgtctgaa-3’
下游序列：5'-acactctttccctacacgacgctcttccgatcnnn……nnnt-3’接头连接完成后，我们会根据模板起始量的不同，使用kapahifihotstartpcrkit对文库进行6-15个循环数的pcr扩增。
(6)文库质检和q-pcr定量：通过琼脂糖凝胶电泳来检测文库质量，使用2％的琼脂糖凝胶，120v，30分钟，凝胶成像，目标条带为270bp。通过agilent公司的2100bioanalyzer对文库片段大小精确定量，并通过q-pcr对文库浓度精确定量。
(7)miseq测序仪上机测序：在illuminamiseq测序平台下获得了读段序列长度为75bp的双端测序数据。
请参考图1，本实施例的具体步骤包括：
s101：通过扩增子文库的构建和上机测序，我们可以获取覆盖48个靶基因全部外显子区域核酸序列信息的读段序列(即双端75bp的测序数据)。
s102：使用fastqc软件对测序数据中的所有读段序列进行质量控制，对于单个碱基平均测序质量得分低于20的读段序列数据定为测序质量差，并在分析中予以删除。
s103:使用扩增子引物序列对测序数据进行过滤，即提取配对的两个读段序列中的扩增子引物序列对，去除扩增子引物序列对不是来源于同一个配对引物的读段序列，进而统计覆盖到每对扩增子引物对应的靶向扩增区域的读段序列数目，比较不同扩增子的扩增效率，删除扩增异常的扩增子所对应的读段序列数据。
本实施例中，如果读段的两端序列所对应的扩增子引物序列是来自我们设计好的同一个引物对，则认为该读段是来源于该扩增子引物的扩增，对所有符合上述条件的扩增子对应的读段序列进行计数，比较所有扩增子的扩增效率(这里用扩增数来衡量扩增效率)，当扩增数小于所有扩增子对应的扩增数均值的1/3时，则判断该扩增子存在异常，并将其扩增出来的所有读段序列在分析中删除。
s104：删除读段序列中的扩增子引物序列，提高突变识别的准确性及比对效率。
本实例中，我们撰写python程序将配对的读段序列与已知的扩增子引物对的序列进行比对，并将读段序列中与引物序列匹配的部分从读段序列中删除，从而得到真正的靶基因dna序列。这里，我们对读段序列进行截取从而剔除引物序列部分的目的在于，该做法一方面能够避免引物合成错误的碱基被当作突变识别出来。另一方面，精简后的序列长度能缩减后续的比对时间。
s105：首先，我们从ucscgenomebrowser数据库中
(http://genome.ucsc.edu/cgi-bin/hgtracks？db＝hg19)下载人类参考基因组序列hg19。其次，我们编写程序将48个靶基因的参考序列从整个人类基因组序列
(hg19)中提取出来。再次，我们将上步中得到的所有读段序列比对到靶基因的参考序列上，从而得到记录着比对结果的bam文件。
本实施例中，利用读段序列与靶基因参考序列进行比对而不是与整个人类基因组序列进行比对，从而提高比对准确性及比对效率。真核生物基因由外显子及内含子拼接而成，直接与靶基因的参考序列进行比对可以较为直接、准确。比对过程应用bwa比对工具，在其他实施案例中，亦可以使用其他的比对软件，比如bowtie、soap2等。
s106：对于上述比对后的bam文件，我们应用突变识别工具探测位于靶基因区域的单核苷酸突变。
在本实施例中，突变识别过程应用varscan2和mutect两种突变识别工具，将分别得到的突变列表取交集，作为用于后续分析的结果数据。
s107：筛选测序深度超过100×的突变位点。
在本实施例中，我们首先应用samtools软件中的depth子程序获得每个突变位点的测序深度。然后，对于那些测序深度低于一定阈值的突变位点予以剔除。本领域人员知晓，目前利用高通量测序进行某区域的snv检测，一般需要该区域30×的测序数据，测序深度越高，获得的等位基因频率越可靠，根据测序深度设置阈值，阈值越大，留下的snv准确程度越高，越可靠，但后续可用数据减少；阈值越小，后续数据量越大，但数据可靠性低。利用这些混有可靠性低的snv进行统计分析得到的假阳性snv多。这里，我们筛选突变位点的阈值设定为100×，在其他实施案例中，亦可根据具体情况修改该阈值。
s108：根据突变频率，参考已有的突变数据库与数据库记录的snv作显著性差异检验，筛选差异显著的突变(p<0.05)。最后使用annovar工具对差异显著的突变进行功能注释，将snv注释到gene上以及各种突变数据库中，从而阐明这些突变的所属类型(同义突变、非同义突变、无义突变等)是否会影响该基因所编码的蛋白质的功能，来进一步揭示这些突变在肺癌形成和发展中的作用。
在本实施例中，突变频率来源于snv识别工具的输出结果。参考的突变数据库包括胚系突变数据库和体细胞突变数据库两类。其中，常用的胚系突变数据库包括千人基因组数据库(http://www.1000genomes.org/)和6万人exac人类外显子组整合数据库(http://exac.broadinstitute.org/)等。常用的体细胞突变数据库包括美国肿瘤基因组图谱tcga数据库(http://cancergenome.nih.gov/)和国际癌症基因组联盟icgc数据库(https://dcc.icgc.org/)等。使用fisher精确检验统计的方法筛选差异显著的突变(即与肿瘤发生相关的基因突变)。共包括四种对象，第一是突变人数比例，即携带突变的病人数量；第二，群体中的突变比例及群体等位基因频率；第三，纯合突变人数比例；第四，杂合突变人数比例。
s109：对本实施案例中的10例肺癌患者的临床资料信息(包括年龄、性别、癌症tnm分期、吸烟史、肿瘤体积、肿块大小、是否有淋巴结侵袭、ki67恶性程度、组织亚型等信息)进行关联分析，找出在肺癌中与特定基因突变发生相关的风险因素，即是否具有某一临床特征的病人倾向于发生哪些特异的突变。
这里，我们应用r统计软件中的卡方检验函数进行关联分析。
实施例一的snv统计结果，见表2
表2
实施例二：11例肺癌患者血清样本肿瘤靶向基因测序数据解析。
本实施例中的文库为11例肺癌患者血清游离dna构建的靶向扩增子文库。文库构建的具体步骤如下：
(1)靶向基因的选择：选择肿瘤热点突变基因、原癌基因、抑癌基因、一些靶向药物作用的基因，具体abl1、egfr、gnas、mlh1、ret、akt1、erbb2、hnf1a、alk、erbb4、hras、notch1、smarcb1、apc、fbxw7、idh1、npm1、smo、atm、fgfr1、jak2、nras、src、braf、fgfr2、jak3、pdgfra、stk11、cdh1、fgfr3、kdr、pik3ca、tp53、cdkn2a、flt3、kit、pten、vhl、csf1r、gna11、kras、ptpn11、ezh2、tnnb1、gnaq、met、rb1、idh2这48个基因，作为我们研究的靶基因。
(2)游离dna的提取及定量：对于肺癌患者的胸水样本，我们先进行低速离心(3,000rpm)5分钟取上清，在进行高速离心(14,000rpm)10分钟取上清，得到了胸水样本中的游离dna(平均长度约为166bp)；并应用qbuit2.0(invitrogen公司)仪器进行定量。
(3)扩增子设计：通过在线引物设计工具designstudio，针对48个靶基因进行引物设计。最终，我们得到了覆盖48个靶基因全部外显子区域的2,158对扩增子，每对扩增子片段的大小约为150bp。由于不同靶基因的序列长度不同，我们的每对扩增子的片段大小又是几乎固定的，因此每个靶基因都对应着不同数目的扩增子引物对。靶基因与扩增子引物对数目的对应列表同表1。
(4)多重pcr扩增靶向基因的外显子：扩增子引物设计完成后，根据设计报告提供的引物序列，合成引物核酸，并以多重pcr的形式扩增靶向基因的全部外显子序列。
(5)illumina测序接头的连接和文库pcr扩增：对于上述扩增产物，我们连接illumina测序仪的测序接头。测序接头序列如下：
上游序列：5'p-nnn……nnngatcggaagagcacacgtctgaa-3’
下游序列：5'-acactctttccctacacgacgctcttccgatcnnn……nnnt-3’接头连接完成后，我们会根据模板起始量的不同，使用kapahifihotstartpcrkit对文库进行6-15个循环数的pcr扩增。
(6)文库质检和q-pcr定量：通过琼脂糖凝胶电泳来检测文库质量，使用2％的琼脂糖凝胶，120v，30分钟，凝胶成像，目标条带为270bp。通过agilent2100bioanalyzer对文库片段大小精确定量，并通过q-pcr对文库浓度精确定量。
(7)hiseq4000测序仪上机测序：在illuminahiseq4000测序平台下获得了读段序列长度为125bp的双端测序数据。
请参考图1，本实施例的具体步骤包括：
s101：通过扩增子文库的构建和上机测序，我们可以获取覆盖48个靶基因全部外显子区域核酸序列信息的读段序列(即双端125bp的测序数据)。
s102：使用fastqc软件对测序数据中的所有读段序列进行质量控制，对于单个碱基平均测序质量低于20的读段序列数据定为测序质量差，并在分析中予以删除。
s103:使用扩增子引物序列对测序数据进行过滤，即提取配对的两个读段序列中的扩增子引物序列对，去除扩增子引物序列对不是来源于同一个配对引物的读段序列，进而统计覆盖到每对扩增子引物对应的靶向扩增区域的读段序列数目，比较不同扩增子的扩增效率，删除扩增异常的扩增子所对应的读段序列数据。
本实施例中，判断该扩增子存在异常的根据和条件同实施例一。
s104：删除读段序列中的扩增子引物序列，提高突变识别的准确性及比对效率。
本实例中，我们撰写python程序将配对的读段序列与已知的扩增子引物对的序列进行比对，并将读段序列中与引物序列匹配的部分从读段序列中删除，从而得到真正的靶基因dna序列。这里，我们对读段序列进行截取从而剔除引物序列部分的目的在于，该做法一方面能够避免引物合成错误的碱基被当作突变识别出来。另一方面，精简后的序列长度能缩减后续的比对时间。
s105：首先，我们从ucscgenomebrowser数据库中下载人类参考基因组序列hg19。其次，我们编写程序将48个靶基因的参考序列从整个人类基因组序列(hg19)中提取出来。再次，我们将上步中得到的所有读段序列比对到靶基因的参考序列上，从而得到记录着比对结果的bam文件。
本实施例中，比对过程应用bowtie作为比对工具。
s106：对于上述比对后的bam文件，我们应用突变识别工具探测出现在靶向基因区域的单核苷酸突变。
在本实施例中，突变识别过程应用varscan2和freebayes两种snv识别工具，将分别得到的突变列表取交集，作为用于后续分析的结果数据。
s107：筛选测序深度超过100×的突变位点，具体操作同实施例一。
s108：统计突变频率，参考已有的突变数据库与数据库记录的snv作显著性差异检验，筛选差异显著的突变(p<0.05)，并使用annovar工具对筛选出的突变进行功能注释。
在本实施例中，突变频率来源于snv识别工具的输出结果。参考的突变数据库以及筛选差异显著的突变的方法同实施例一。
s109：对本实施案例中的11例肺癌患者的临床资料信息(包括年龄、性别、癌症tnm分期、吸烟史、肿瘤体积、ki67恶性程度、组织亚型等信息)进行关联分析，找出在肺癌中与特定基因突变发生相关的风险因素，即是否具有某一临床特征的病人倾向于发生哪些特异的突变。
实施例二的snv统计结果，见表3.
表3
本领域技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。
以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换。