阅读说明：本技术 一种利用基因组数据探究疾病亚型亲缘性的方法 (Method for exploring disease subtype affinity by using genome data ) 是由 侯群星 袁卫兰 高军晖 林灵 吴昊天 蒋丽莎 李无霜 王瑶瑶 吴守信 许骋 于 2020-05-28 设计创作，主要内容包括：本发明提供了一种利用基因组数据探究疾病亚型亲缘性的方法,所述方法包括采用超几何分布检验计算基因非沉默突变富集值的步骤；所述基因非沉默突变富集值的计算公式为：<Image he="120" wi="248" file="DDA0002513874290000011.GIF" imgContent="drawing" imgFormat="GIF" orientation="portrait" inline="no"></Image>其中,n<Sub>f</Sub>为基因在疾病亚型中发生非沉默突变的样本数,N<Sub>f</Sub>为基因在所有样本中发生非沉默突变的样本数,n为疾病亚型样本数,N为样本总数；所述疾病亚型的分组数量不小于3。本发明在进行疾病亚型亲缘性分析之前,首先利用超几何分布检验计算基因在各肿瘤亚型中的非沉默突变富集分数,再利用此富集分数进行亲缘性分析,减小了样本总数、疾病亚型样本数等背景因素对分析结果造成的影响,提高了方法的准确性。(The invention provides a method for exploring disease subtype affinity by using genome data, which comprises the steps of calculating a gene non-silent mutation enrichment value by adopting a hyper-geometric distribution test; the calculation formula of the gene non-silent mutation enrichment value is as follows: wherein n is f Number of samples of genes with non-silent mutations in disease subtypes, N f The number of samples of the gene which has non-silent mutation in all samples, N is the number of samples of disease subtype, and N is the total number of samples; the number of subgroups of disease subtypes is not less than 3. The invention is about to treat diseasesBefore disease subtype affinity analysis, the non-silent mutation enrichment fraction of the gene in each tumor subtype is calculated by using hyper-geometric distribution test, and affinity analysis is performed by using the enrichment fraction, so that the influence of background factors such as total number of samples, disease subtype sample number and the like on an analysis result is reduced, and the accuracy of the method is improved.)

一种利用基因组数据探究疾病亚型亲缘性的方法

技术领域

本发明属于生物信息分析技术领域，涉及一种利用基因组数据探究疾病亚型亲缘性的方法。

背景技术

癌症是一类由于细胞分类和调往机制失常而导致的疾病，通常表现为恶性肿瘤。由于癌症的早期诊断准确性差，复发率和死亡率高，已经成为严重威胁人类健康的疾病之一。近年来人们已经认识到肿瘤的发生和转移是多基因、多步骤相互作用连续发展的结果，对肿瘤进行整体、全面、动态的研究是防治肿瘤的根本途径。不同的肿瘤存在不同的亚型，由于肿瘤的临床异质性，在临床上对于不同的肿瘤亚型需要采用不同的治疗策略。尽管如此，探究肿瘤亚型亲缘性对于肿瘤的临床治疗和预后同样具有重要的意义。

目前，研究肿瘤亚型亲缘性的方法的主要步骤为：1)根据病人样本获取WES测序数据；2)根据测序数据分析病人的体细胞突变情况；3)对突变位点进行注释；4)对注释结果进行格式转换并对突变类型进行注释；5)突变类型筛选；6)统计各基因在各样本中是否存在非沉默突变；7)根据样本的非沉默突变统计结果，使用层次聚类方法计算各疾病亚型之间的亲缘性。

但是现有技术只统计了样本中非沉默突变情况的有或无，并直接根据此统计结果聚类分析疾病亚型的亲缘性，并未考虑样本数量等背景因素对结果造成的影响。因此，现有方法得出的结果可能存在准确性差的问题。

因此，提供一种更加精准的疾病亚型亲缘性的分析方法，在肿瘤的临床治疗和预后监测领域具有重要意义。

发明内容

针对现有技术的不足和实际需求，本发明提供了一种利用基因组数据探究疾病亚型亲缘性的方法，所述方法在进行疾病亚型亲缘性分析之前，首先利用超几何分布检验计算基因在各肿瘤亚型中的非沉默突变富集分数，再利用此富集分数进行亲缘性分析，减小了样本总数、疾病亚型样本数等因素对分析结果造成的影响。

为达此目的，本发明采用以下技术方案：

第一方面，本发明提供了一种疾病亚型亲缘性的分析方法，所述方法包括采用超几何分布检验计算基因非沉默突变富集值的步骤；

所述基因非沉默突变富集值的计算公式为：

其中，n_f为基因在疾病亚型中发生非沉默突变的样本数，N_f为基因在所有样本中发生非沉默突变的样本数，n为疾病亚型样本数，N为样本总数；

所述疾病亚型的分组数量不小于3。

本发明中，在进行疾病亚型亲缘性分析之前，首先利用超几何分布检验计算基因在各肿瘤亚型中的非沉默突变富集分数，再利用此富集分数进行亲缘性分析，减小了样本总数、疾病亚型样本数等因素对分析结果造成的影响，提高了方法的准确性。

优选地，所述方法包括以下步骤：

(1)对肿瘤和正常样本进行测序，获取全外显子测序数据；

(2)根据测序数据分析样本的体细胞突变情况；

(3)对突变位点进行注释；

(4)将注释结果进行格式转换，并对突变类型进行注释；

(5)根据注释后的突变类型，筛选出非沉默突变类型；

(6)计算基因非沉默突变富集值；

(7)根据基因的非沉默突变富集值，使用层次聚类方法计算疾病亚型间的亲缘性。

优选地，步骤(2)所述分析的步骤为：

1)对获得的全外显子测序数据进行过滤，筛选Q20≥90％、Q30≥80％的测序数据；

2)对参考基因组构建比对索引，并将步骤1)筛选后的测序数据比对到参考基因组，得到比对后的数据；

3)统计测序数据中比对到参考基因组的测序序列的比例；

4)计算步骤2)获得的比对后的数据的深度、平均比对质量和覆盖度；

5)统计步骤2)获得的比对后的数据在参考基因组的目标区域的比例、平均覆盖深度和覆盖度；

6)标记步骤2)获得的比对后的数据中的PCR重复序列，并对PCR重复序列进行去重；

7)对步骤6)获得的去重后的数据进行位点矫正；

8)对步骤7)矫正后的数据进行分组；

9)对步骤8)分组后的数据进行过滤，得到变异频率大于5％的初始体细胞突变。

优选地，步骤(3)包括：对步骤(2)获得的变异频率大于5％的初始体细胞突变进行筛选，并对筛选后的体细胞突变进行注释的步骤。

优选地，步骤(5)所述非沉默突变类型包括移码缺失突变(Frame_Shift_Del)、移码***突变(Frame_Shift_Ins)、框内缺失(In_Frame_Del)、框内***(In_Frame_Ins)、错义突变(Missense_Mutation)、无义突变(Nonsense_Mutation)、终止密码子突变(Nonstop_Mutation)或剪接位点(Splice_Site)中的任意一种或至少两种的组合。

优选地，在步骤(6)之前还包括：根据筛选出的非沉默突变类型数据，统计基因在样本中的非沉默突变情况，得到基因非沉默突变情况与样本的矩阵数据。

第二方面，本发明提供了一种疾病亚型亲缘性的分析装置，所述装置包括基因非沉默突变富集值计算模块，用于统计基因在样本中的非沉默突变情况，并根据基因非沉默突变富集值的计算公式

计算基因在样本中的基因非沉默突变富集值；

其中，n_f为基因在疾病亚型中发生非沉默突变的样本数，N_f为基因在所有样本中发生非沉默突变的样本数，n为疾病亚型样本数，N为样本总数；

所述疾病亚型的分组数量不小于3。

优选地，所述装置还包括：

测序模块，用于获取肿瘤和正常样本的全外显子测序数据；

样本体细胞突变情况分析模块；

突变位点注释模块，用于对初始体细胞突变结果进行筛选，并对筛选后的体细胞突变位点进行注释；

格式转换和突变类型注释模块；

非沉默突变类型筛选模块，用于筛选突变类型为移码缺失突变(Frame_Shift_Del)、移码***突变(Frame_Shift_Ins)、框内缺失(In_Frame_Del)、框内***(In_Frame_Ins)、错义突变(Missense_Mutation)、无义突变(Nonsense_Mutation)、终止密码子突变(Nonstop_Mutation)或剪接位点(Splice_Site)中的任意一种或至少两种的组合；

疾病亚型间的亲缘性计算模块，用于根据基因非沉默突变富集值，使用层次聚类方法计算疾病亚型间的亲缘性。

优选地，所述样本的基因突变情况分析模块包括：

测序数据质控单元，用于对获得的全外显子组测序数据进行过滤，并筛选Q20≥90％、Q30≥80％的测序数据；

序列比对单元，用于构建参考基因组比对索引，将质控后的数据比对到参考基因组，并获得比对后的数据；

比对数据分析单元，用于统计测序数据中比对到参考基因组的测序序列的比例，计算比对后的数据的深度、平均比对质量和覆盖度，统计比对后的数据在参考基因组的目标区域的比例、平均覆盖深度和覆盖度；

比对数据处理单元，用于获得比对后的数据中的PCR重复序列，并对PCR重复序列进行去重，对去重后的数据进行位点矫正，对矫正后的数据进行分组；

初始体细胞突变位点获取单元，用于对分组后的数据进行过滤，得到变异频率大于5％的初始体细胞突变。

第三方面，本发明提供了一种如第二方面所述的装置在分析疾病亚型亲缘性中的应用。

与现有技术相比，本发明具有如下有益效果：

本发明在进行疾病亚型亲缘性分析之前，首先利用超几何分布检验计算基因在各肿瘤亚型中的非沉默突变富集分数，再利用此富集分数进行亲缘性分析，减小了样本总数、疾病亚型样本数等因素对分析结果造成的影响，提高了方法的准确性。

附图说明

图1为一种疾病亚型亲缘性的分析方法的流程图；

图2为一种疾病亚型亲缘性的分析装置的结构图；

图3A为仅根据统计的各基因在疾病亚型中的样本突变个数进行聚类分析的结果，图3B为根据本发明所述方法计算的各基因突变富集分数进行聚类分析的结果。

具体实施方式

为进一步阐述本发明所采取的技术手段及其效果，以下结合实施例和附图对本发明作进一步地说明。可以理解的是，此处所描述的具体实施方式仅仅用于解释本发明，而非对本发明的限定。

实施例中未注明具体技术或条件者，按照本领域内的文献所描述的技术或条件，或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可通过正规渠道商购获得的常规产品。

实施例1

图1为一种疾病亚型亲缘性的分析方法的流程图，本实施例的方法可以由疾病亚型亲缘性的分析装置来执行，该装置可通过硬件和/或软件的方式实现，并一般可集成于用于实现疾病亚型亲缘性分析功能的服务器中，本实施例的方法具体包括：

110、对样本进行测序，获取全外显子测序数据(WES)

将Illumina测序仪测序得到的原始图像数据进行碱基识别(Base Calling)，转化为原始序列数据，即一条条的读长(reads)，称为raw data，结果以fastq文件格式储存(文件名：*.fq)，其中每条reads的分子标签以fastq文件格式储存；

120、根据测序数据分析样本的体细胞突变情况

利用fastqc软件统计raw data的碱基质量得分(Q20和Q30)、序列长度、N字符数量、数据量、GC含量和PCR重复率，筛选Q20≥90％、Q30≥80％的测序数据；

利用bwasw算法对参考基因组构建比对索引，使用bwa软件的BWA-MEM算法将测序数据中的目标序列比对到参考基因组(hg19/hg38)，得到比对后的数据(格式为bam文件)；

利用samtools软件的flagstat模块统计测序数据中比对到参考基因组的测序序列的比例；

利用qulimap软件计算bam文件(即比对后的数据)的深度、平均比对质量和覆盖度；

利用GATK CollectHsmetrics模块统计bam文件在参考基因组的目标区域的数据量比例、平均覆盖深度和1×、20×、50×和100×覆盖度百分比；

使用GATK4 MarkDuplicates模块标记bam文件中的PCR重复序列，并对PCR重复序列进行去重；

使用GATK4 BaseRecalibrator模块、外部公共位点数据库和GATK4ApplyBQSR模块对去重后的数据进行位点矫正，得到BQSR矫正后数据；

使用GATK4 AddOrReplaceReadGroups模块对比对后的数据进行分组，即增加bam文件的分组信息；

使用GATK4 mutect2模块和GATK4 FilterMutectCalls模块检测过滤，得到变异频率大于5％的初始体细胞突变；

130、对突变位点进行注释

筛选初始体细胞突变结果文件(vcf格式)中FILTER为“PASS”的结果；

利用ANNOVAR软件以及相应的数据库(refGene、avsnp138、clinvar_20170130、cosmic70、1000g2015aug_all、ljb26_all)对筛选后的体细胞突变结果进行注释；

140、将注释结果进行格式转换，并对突变类型进行注释

利用R语言(v3.5.2)及相关程序将所有肿瘤样本注释后的体细胞突变结果的vcf文件转换为maf文件，并合并为一个文件；

注释体细胞变异类型，类型注释标准如下：

1)如果在Func.refGene列中匹配到表1指定的几种类型，将ExonicFunc.refGene列的内容更改为相对应的类型；若Func.refGene列中的类型未在表1中指定，则ExonicFunc.refGene列原有类型保持不变；

2)根据表2ExonicFunc.refGene列的类型对相应的变异信息进行注释；

表1

表2

150、根据注释后的突变类型，筛选出非沉默突变类型

根据注释后的突变类型结果，筛选出非沉默突变类型，包括移码缺失突变(Frame_Shift_Del)、移码***突变(Frame_Shift_Ins)、框内缺失(In_Frame_Del)、框内***(In_Frame_Ins)、错义突变(Missense_Mutation)、无义突变(Nonsense_Mutation)、终止密码子突变(Nonstop_Mutation)和剪接位点(Splice_Site)中的任意一种或至少两种的组合；

160、计算基因非沉默突变富集值

根据筛选出的非沉默突变类型数据maf文件，统计基因在样本中的非沉默突变情况，得到基因非沉默突变情况与样本的矩阵数据；

根据样本的疾病亚型分组(组数不小于3组)、样本数量等信息，使用超几何分布检验，计算基因非沉默突变富集值，以此将基因与样本非沉默突变情况的矩阵数据转化为基因与疾病亚型的非沉默突变富集矩阵数据；

基因突变富集值的计算公式如下：

其中，n_f为基因在疾病亚型中发生非沉默突变的样本数，N_f为基因在所有样本中发生非沉默突变的样本数，n为疾病亚型样本数，N为样本总数；

所述疾病亚型的分组数量不小于3。

170、疾病亚型间的亲缘性分析

根据基因在各疾病亚型中的基因突变富集值，使用层次聚类的方法对样本的亲缘性进行计算，得到疾病亚型间的亲缘关系，利用R语言的pheatmap数据包进行可视化展示。

实施例2

图2为一种疾病亚型亲缘性的分析装置的结构图，所述装置包括：

测序模块210、样本的体细胞突变情况分析模块220、突变位点注释模块230、格式转换和突变类型注释模块240、非沉默突变类型筛选模块250、基因非沉默突变富集值计算模块260和疾病亚型间的亲缘性计算模块270；

测序模块210，用于获取样本的全外显子测序数据；

样本的体细胞突变情况分析模块220，包括：

测序数据质控单元，用于对获得的测序数据进行过滤，并筛选Q20≥90％、Q30≥80％的测序数据；

序列比对单元，用于构建参考基因组比对索引，将质控后的数据比对到参考基因组，并获得比对后的数据；

比对数据分析单元，用于统计测序数据中比对到参考基因组的测序序列的比例，计算比对后的数据的深度、平均比对质量和覆盖度，统计比对后的数据在参考基因组的目标区域的比例、平均覆盖深度和覆盖度；

比对数据处理单元，用于获得比对后的数据中的PCR重复序列，并对PCR重复序列进行去重，对去重后的数据进行位点矫正，对矫正后的数据进行分组；

初始体细胞突变位点获取单元，用于对分组后的数据进行过滤，得到变异频率大于5％的初始体细胞突变；

突变位点注释模块230，用于对初始体细胞突变结果进行筛选，并对筛选后的体细胞突变位点进行注释；

格式转换和突变类型注释模块240，用于对突变类型进行注释；

非沉默突变类型筛选模块250，用于筛选突变类型为移码缺失突变、移码***突变、框内缺失、框内***、错义突变、无义突变、终止密码子突变或剪接位点中的任意一种或至少两种的组合；

基因突变富集值计算模块260，用于统计基因在样本中的非沉默突变情况，并根据基因非沉默突变富集值的计算公式

计算基因在样本中的基因非沉默突变富集值；

其中，n_f为基因在疾病亚型中发生非沉默突变的样本数，N_f为基因在所有样本中发生非沉默突变的样本数，n为疾病亚型样本数，N为样本总数；所述疾病亚型的分组数量不小于3。

疾病亚型间的亲缘性计算模块270，用于根据基因非沉默突变富集值，使用层次聚类方法计算疾病亚型间的亲缘性。

实施例3

本实施例采用临床数据验证方法的可行性，下载TCGA数据库乳腺癌基因层面数据及临床数据，并对其进行整理；

根据临床对乳腺癌亚型的划分标准，即Luminal A型(ER+/PR+,HER-2-)、LuminalB型(ER+/PR+,HER-2+)、HER-2+型(ER-,PR-,HER-2+)和Basal-like型(ER-,PR-,HER-2-)，对整理好的样本数据进行亚型分组；

根据亚型分组情况统计各基因在疾病亚型中的样本突变个数，并仅根据统计的突变样本数进行聚类分析，结果如图3A所示；

根据亚型分组情况使用本文所述方法统计各基因在各亚型中的样本突变富集分数，并根据此突变富集分数进行聚类分析，结果如图3B所示。

在临床上，由于Luminal A型与Luminal B型关系更接近，通过对比图3A与图3B可以看出，利用富集分数进行亲缘性分析，明显提高了方法的准确性。

综上所述，本发明在进行疾病亚型亲缘性分析之前，首先利用超几何分布检验计算基因在各肿瘤亚型中的非沉默突变富集分数，再利用此富集分数进行亲缘性分析，减小了样本总数、疾病亚型样本数等因素对分析结果造成的影响，提高了方法的准确性。

申请人声明，本发明通过上述实施例来说明本发明的详细方法，但本发明并不局限于上述详细方法，即不意味着本发明必须依赖上述详细方法才能实施。所属技术领域的技术人员应该明了，对本发明的任何改进，对本发明产品各原料的等效替换及辅助成分的添加、具体方式的选择等，均落在本发明的保护范围和公开范围之内。