CN112767999A - 一种全基因组测序数据的分析方法及装置 - Google Patents

Abstract

本申请公开了一种全基因组测序数据的分析方法及装置，方法包括：从预设数据库收集来自同一受检对象的全基因组测序数据；计算受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性；根据每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型；通过BIC算法获取最优P值；获取最优P值对应的各个片段的基因型、非正常组织样本纯度γ、非正常组织细胞倍性κ、各个片段的拷贝数以及等位基因拷贝数。采用本申请所提供的方案，从全基因组测序数据获得更全面的信息，利用这些参数，能够使得拷贝数计算结果更加精确。

Description

一种全基因组测序数据的分析方法及装置

技术领域

本申请涉及生物科学领域，特别涉及一种全基因组测序数据的分析方法及装置。

背景技术

非正常组织的研究是生命医学中的重要研究领域，对于肿瘤组织这类非正常组织的研究，对人类健康生活有重大影响。恶性肿瘤又称为癌症，是一类细胞恶性增殖的疾病，因其病理十分复杂，人类尚无法攻克这类疾病。精准医疗的提出，为攻克这一疾病带来了希望。精准医疗可以针对每个病人的遗传信息和癌症发展状况有目标地准确治疗病人。二代测序(next generation sequencing)为快速检测病人遗传信息提供了可能。目前，有很多与癌症相关的基因被证实与拷贝数变异相关。例如，与拷贝数增加相关的原癌基因有CDK4，CKD6，PIK3CA，EGFR等，与拷贝数减少相关的原癌基因有RB1，NF1，SMAD4等。对基因拷贝数的分析提供了这个基因的总拷贝数，但对其等位基因拷贝数的分析能提供关于这个基因的更多细节。

等位基因拷贝数的分析包含了基因拷贝数分析，同时它多了以下几个优点。首先，它能检测出更多类型的拷贝数变异事件。例如，某个基因的基因型从AB变为了AA，此时它的拷贝数没有变，仍然为2，但是基因型却从1∶1变为了2∶0，这种拷贝数没有变化，因此杂合性丢失的变异事件不能通过拷贝数分析检测到，但是可以通过等位基因拷贝数分析检测到。其次，拷贝数分析只能得到某个基因拷贝数的变化，而分析等位基因拷贝数能得到某个基因特定等位基因拷贝数的变化，这提供了这个基因的更多信息。然而测序需要从病人组织中提取样本，但通常癌症组织并不是只单纯地包含癌症细胞，它还有非常丰富的微环境。另外，样本纯度对于基因拷贝数的影响也较大，因此，得到肿瘤样本纯度之后，通过肿瘤样本的纯度可以对基因的拷贝数进行校正从而得到准确的癌细胞基因拷贝数，其次，拷贝数的准确率还受非正常组织细胞倍性的影响，可见，等位基因、基因型、非正常组织样本纯度、非正常组织细胞倍性等多个参数都会影响拷贝数计算结果的准确度，可见，在分析拷贝数时，需要从测序数据中获得更为全面的数据，才能使拷贝数计算结果更加准确。

因此，亟需提供一种全基因组测序数据的分析方法，用以从全基因组测序数据获得更全面的信息。

发明内容

本申请实施例的目的在于提供一种全基因组测序数据的分析方法及装置，用以从全基因组测序数据获得更全面的信息。

为了解决上述技术问题，本申请的实施例采用了如下技术方案：一种全基因组测序数据的分析方法，包括：

从预设数据库收集来自同一受检对象的全基因组测序数据，所述全基因组测序数据包括正常组织样本数据和非正常组织样本数据；

计算所述受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性；其中，所述P值为相邻拷贝数对应的基因组片段的TRE的差值；

根据所述每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型；

通过BIC算法获取最优P值；

获取所述最优P值对应的各个片段的基因型、非正常组织样本纯度γ、非正常组织细胞倍性κ、各个片段的拷贝数以及等位基因拷贝数。

本申请的有益效果在于：在计算拷贝数时，还得到了等位基因、基因型、非正常组织样本纯度、非正常组织细胞倍性等多个对拷贝数产生影响的参数，从而从全基因组测序数据获得更全面的信息，利用这些参数，能够使得拷贝数计算结果更加精确。

在一个实施例中，所述P值通过如下方式计算：

统计各基因组片段内的TRE的均值方差和窗口数量；

基于基因组片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值P。

在一个实施例中，所述基因组片段通过如下方式形成：

根据所述测序序列将全基因组测序数据划分为特定碱基数的窗口；

计算每个窗口的TRE值；

利用TRE值将相邻窗口进行合并，以形成新的基因组片段。

在一个实施例中，所述计算所述受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性，包括：

根据所述差值P推算出拷贝数为2的片段对应的TRE位点Q；

根据所述TRE位点Q计算出所有P取值下非正常组织样本纯度和非正常细胞倍性，并计算所有P取值下各个片段的拷贝数。

在一个实施例中，还包括：

根据如下公式计算各个P取值下的非正常组织样本纯度和非正常组织细胞倍性：

其中σ_i表示拷贝数为i的所有片段中，其TRE的标准差。S_s为片段S的TRE观测值，Sⁱ为在某一P值下，理论上拷贝数为i的片段的TRE值。p_i表示片段S的拷贝数为i的权重，取值均为1。

在一个实施例中，根据所述每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型，包括：

根据所述各个P取值下的非正常组织样本纯度和非正常组织细胞倍性、所述各个片段的拷贝数以及变异位点的LAR数据计算每个P值对应的各个片段的基因型估计值。

在一个实施例中，所述变异位点的LAR数据通过以下方式获取：

将所述全基因组测序数据与参考基因组进行比对，以获取比对文件；

从比对文件中获取每一个测序序列的位置信息、变异位点的LAR数据。

本实施例还提供一种全基因组测序数据的分析装置，包括：

收集模块，用于从预设数据库收集来自同一受检对象的全基因组测序数据，所述全基因组测序数据包括正常组织样本数据和非正常组织样本数据；

第一计算模块，用于计算所述受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性；其中，所述P值为相邻拷贝数对应的基因组片段的TRE的差值；

第二计算模块，用于根据所述每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型；

第一获取模块，用于通过BIC算法获取最优P值；

第二获取模块，用于获取所述最优P值对应的各个片段的基因型、非正常组织样本纯度γ、非正常组织细胞倍性κ、各个片段的拷贝数以及等位基因拷贝数。

在一个实施例中，所述P值通过如下方式计算：

统计各基因组片段内的TRE的均值方差和窗口数量；

基于基因组片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值P。

在一个实施例中，所述基因组片段通过如下方式形成：

根据所述测序序列将全基因组测序数据划分为特定碱基数的窗口；

计算每个窗口的TRE值；

利用TRE值将相邻窗口进行合并，以形成新的基因组片段。

在一个实施例中，所述第一计算模块，包括：

推算子模块，用于根据所述差值P推算出拷贝数为2的片段对应的TRE位点Q；

第一计算子模块，用于根据所述TRE位点Q计算出所有P取值下非正常组织样本纯度和非正常细胞倍性，并计算所有P取值下各个片段的拷贝数。

在一个实施例中，还包括：

第三计算模块，用于根据如下公式计算各个P取值下的非正常组织样本纯度和非正常组织细胞倍性：

其中σ_i表示拷贝数为i的所有片段中，其TRE的标准差。S_s为片段S的TRE观测值，Sⁱ为在某一P值下，理论上拷贝数为i的片段的TRE值。p_i表示片段S的拷贝数为i的权重，取值均为1。

在一个实施例中，所述第二计算模块，包括：

第二计算子模块，用于根据所述各个P取值下的非正常组织样本纯度和非正常组织细胞倍性、所述各个片段的拷贝数以及变异位点的LAR数据计算每个P值对应的各个片段的基因型估计值。

在一个实施例中，所述变异位点的LAR数据通过以下方式获取：

将所述全基因组测序数据与参考基因组进行比对，以获取比对文件；

从比对文件中获取每一个测序序列的位置信息、变异位点的LAR数据。

附图说明

图1为本申请一实施例中一种全基因组测序数据的分析方法的流程图；

图2为本申请一实施例中癌症组织样本和正常组织样本中细胞组分的示例图；

图3为本申请一实施例中TRE和LAR的示意图；

图4为本申请一实施例中拷贝数的计算过程示意图；

图5为本申请一实施例中拷贝数输出结果的示意图，输出结果中包含肿瘤组织样本的纯度和倍性，以及完全拷贝数图谱和主要等位基因拷贝数图谱；

图6为本申请另一实施例中一种全基因组测序数据的分析方法的流程图；

图7为本申请一总的实施例中一种全基因组测序数据的分析方法的流程图；

图8为本申请一实施例中一种全基因组测序数据的分析装置的框图；

图9为本申请另一实施例中一种全基因组测序数据的分析装置的框图。

具体实施方式

此处参考附图描述本申请的各种方案以及特征。

应理解的是，可以对此处申请的实施例做出各种修改。因此，上述说明书不应该视为限制，而仅是作为实施例的范例。本领域的技术人员将想到在本申请的范围和精神内的其他修改。

包含在说明书中并构成说明书的一部分的附图示出了本申请的实施例，并且与上面给出的对本申请的大致描述以及下面给出的对实施例的详细描述一起用于解释本申请的原理。

通过下面参照附图对给定为非限制性实例的实施例的优选形式的描述，本申请的这些和其它特性将会变得显而易见。

还应当理解，尽管已经参照一些具体实例对本申请进行了描述，但本领域技术人员能够确定地实现本申请的很多其它等效形式。

当结合附图时，鉴于以下详细说明，本申请的上述和其他方面、特征和优势将变得更为显而易见。

此后参照附图描述本申请的具体实施例；然而，应当理解，所申请的实施例仅仅是本申请的实例，其可采用多种方式实施。熟知和/或重复的功能和结构并未详细描述以避免不必要或多余的细节使得本申请模糊不清。因此，本文所申请的具体的结构性和功能性细节并非意在限定，而是仅仅作为权利要求的基础和代表性基础用于教导本领域技术人员以实质上任意合适的详细结构多样地使用本申请。

本说明书可使用词组“在一种实施例中”、“在另一个实施例中”、“在又一实施例中”或“在其他实施例中”，其均可指代根据本申请的相同或不同实施例中的一个或多个。

图1为本申请实施例的一种全基因组测序数据的分析方法的流程图，该方法包括以下步骤S11-S15：

在步骤S11中，从预设数据库收集来自同一受检对象的全基因组测序数据，全基因组测序数据包括正常组织样本数据和非正常组织样本数据；

在步骤S12中，计算受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性；其中，P值为相邻拷贝数对应的基因组片段的TRE的差值；

在步骤S13中，根据每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型；

在步骤S14中，通过BIC算法获取最优P值；

在步骤S15中，获取最优P值对应的各个片段的基因型、非正常组织样本纯度γ、非正常组织细胞倍性κ、各个片段的拷贝数以及等位基因拷贝数。

为了更好地理解本发明，下面提供相关的解释和说明：

whole genome sequencing(WGS)：全基因组测序；

Read：基因组测序得到的测序片段读长，bam文件的基本输出单位；

Window(窗口)将基因组按照一定长度分割成的片段单位；

Tumor Read Enrichment(TRE)：肿瘤片段富集程度，指基因组某一段序列在癌症组织样本中的片段富集程度与该段序列在正常样本中的片段富集程度的比值。

定义公式如下：

公式中

和

表示在癌症样本中覆盖片段S的read数量和相匹配的正常样本中覆盖片段S的read数量，N_t和N_n分别表示癌症样本与匹配的正常样本的测序read总量。

Heterozygous germline single nucleotide variants(HGSNV)：杂合生殖系细胞单碱基变异，由于人类染色体属于二倍体，理论上由胚胎细胞发育而来HGSNV位点只有两种碱基类型A和B，其中一种来源于父本，另一种来源于母本。

Log ratio of Allelic coverage Ratios of HGSNVs(LAR)：HGSNV位点上片段覆盖数的对数比例。公式定义如下：

公式中

和

分别代表正常样本中reference base(参考碱基)和alternativebase(备选碱基)的碱基数目。

种

分别代表癌症组织中reference base和alternativebase的碱基数目。

Major allele copy number：主要等位基因的拷贝数，指在拷贝数为i的片段中，主要等位基因的拷贝数的取值，它的值应该是大于等于

的整数。

Peak：基因组中的各个片段具有不同的TRE，具有相似TRE片段聚集为一个簇，每个簇定义为peak，TRE最小的簇的认为是first peak

本实施例中，从预设数据库收集来自同一受检对象的全基因组测序数据，全基因组测序数据包括正常组织样本数据和非正常组织样本数据；

该非正常组织样本数据可以是指恶性肿瘤样本数据，也可以称之为癌症样本数据。具体的，从预设数据库中获取来自同一癌症病人的全基因组测序数据，该全基因组测序数据包括癌症样本数据和正常组织样本数据。该预设数据库可以是如TCGA(The CancerGenome Atlas，癌症基因信息)数据库这类存放有大量癌症信息的数据库。

在得到全基因组测序数据之后，计算受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性；其中，P值为相邻拷贝数对应的基因组片段的TRE的差值；

具体的，将测序数据比对到人类基因组参考序列。从比对到参考基因组的BAM格式的测序数据中提取每一个测序序列的位置信息，提取各HGSNV(Heterozygous germlinesingle nucleotide variants，杂合生殖系细胞单碱基变异)位点上的LAR(Log ratioofAllelic coverage Ratios of HGSNV，HGSNV位点上片段覆盖数的对数比例)信息和覆盖该位点的序列数量。

图2上述BAM格式的测序数据示意图，其用于表示癌症组织样本和正常组织样本中细胞组分的示例图，其中癌症组织样本中除了含有N标注的正常细胞外，还很有C1和C2这两种不同的癌细胞，它们各自含有的突变类型进行了标注(扩增和缺失)，通过对这两种细胞进行全基因组测序(WGS)会得到Accucopy的输入数据(BAM格式)。

统计样本有多少个序列，由于序列总数过大，所以需要划分窗口，划分窗口，以500个碱基为单位划分窗口，统计窗口的序列数量。

计算每一个窗口的TRE(Tumor Read Enrichment，肿瘤片段富集程度)值，然后运用TRE，对基因组进行片段化，即将TRE值相似的窗口拼接成新的基因组片段。

以拼接后的基因组片段为单位，统计片段内的TRE的均值方差和窗口数量。基于片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值，即P值。

其中，在建立自回归模型时，使用上述片段内的TRE的均值和方差，对基因组每个片段的窗口数量进行平滑处理，使TRE的分布更平滑。平滑处理后的TRE分布图如图3所示。

图3为TRE和LAR的示意图，图上写了3个与该图有关的步骤，TRE分布图的x轴表示TRE的值，y轴表示该值上的丰度(即有多少个TRE值与x轴坐标值相等)，LAR图的x轴表示LAR的值，同样地，y轴表示该值上的丰度。TRE图上每一个峰上的片段都会有一个LAR图与之对应。

对于自回归模型计算所得P值，选择[P-0.2，P+0.2]作为P值的取值范围。遍历所有P(例如，以千分之一为间隔进行遍历)，寻找每个P对应的first peak(第一个峰)，即P对应的第一个整数拷贝数所对应的TRE值，也即每个P点在图4中TRE柱状图(平滑处理后的TRE分布图)的第一个柱形对应的值。P值和first peak位点，推断出在first peak之前可能存在的peak(峰、即柱状图中的各个柱形)的个数。根据TRE分布规律，虚拟的peak与实际peak的可能排布情况，推算出拷贝数为2的片段对应的TRE位点，即为Q。

在得到每个P值对应的TRE位点，即Q值之后，可以通过公式

计算出每个P取值下的癌症样本纯度γ；可以通过公式

计算出每个P取值下的癌细胞倍性κ；

计算其他peak与Q距离的周期数p，如图4所示，如果这个peak在Q的右边，则为正值，在左边则为负值，将p加2则得到这个peak代表的拷贝数，落在这个peak上的片段均为这个拷贝数。如果拷贝数为浮点数，即当拷贝数不是整数时，则说明这个peak上的片段为亚克隆片段。

通过图4所示，通过求解P，再找到Q，就可以推测出其它峰所代表的拷贝数。LAR图中峰的特点是关于x＝0对称，有些图只有一个峰是因为两个峰在x＝0处重合了，表现为只有一个峰了。通过对LAR的数据进行处理后，会得到每个柱形代表的拷贝数，其中较大的那个拷贝数我们称为主要等位基因拷贝数。

得到所有片段的拷贝数后，根据如下公式计算各个P取值下的非正常组织样本纯度和非正常组织细胞倍性：

其中σ_i表示拷贝数为i的所有片段中，其TRE的标准差。S_s为片段S的TRE观测值，Sⁱ为在某一P值下，理论上拷贝数为i的片段的TRE值。p_i表示片段S的拷贝数为i的权重，取值均为1，L(S_s；γ，κ)为各个P取值下的非正常组织样本纯度和非正常组织细胞倍性。

根据每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型；

具体的，根据各个P取值下的非正常组织样本纯度和非正常组织细胞倍性、各个片段的拷贝数以及变异位点的LAR数据计算每个P值对应的各个片段的基因型估计值。

通过BIC算法获取最优P值；获取最优P值对应的各个片段的基因型、非正常组织样本纯度γ、非正常组织细胞倍性κ、各个片段的拷贝数以及等位基因拷贝数。

具体的，在得到某个P取值下的非正常组织样本纯度γ和非正常组织细胞倍性κ，片段总拷贝数，和每个片段上HGSNVs位点上的LAR数据，使用EM算法得出在该P取值下每个片段的基因型估计值

和最大的LAR似然值L(r；k，l)。

得到所有可能P取值下的L(S_s；γ，κ)和L(r；k，l)后，通过BIC算法获取最优P值，以及该最优P值下的癌症样本纯度γ和癌细胞倍性κ，各个片段的拷贝数和等位基因拷贝数。

在具体应用过程中，上述步骤S11-S15可通过软件的形式来实现，具体设计的软件包括PICARD软件、GADA软件、strelka2软件和回归模型。具体通过以下步骤A1-A10来实现：

在步骤A1中，在TCGA中下载全基因组测序数据；

在步骤A2中，比对到参考基因组，并将步骤一中癌症和正常样本的bam文件，通过PICARD软件(一个用来处理高通量测序数据的程序)的SamToFastq方法，转换为fastq文件。将得到的fastq文件使用bwa mem方法比对到参考基因组上，得到新的bam文件。

在步骤A3中，使用strelka2(一种快速、准确的小变异体调用程序，用于分析小型世代的种系变异和肿瘤/正常样本对的体细胞变异)对癌症和正常样本的bam文件做双样本的人类基因组变异分析，同时输入参考HGSNVs位点文件。strelka2输出vcf格式的人类基因组变异信息文件，再使用rust语言的rust-htslib包处理该文件，提取HGSNVs位点上的LAR信息。

在步骤A4中，使用rust语言的rust-htslib包处理癌症样本和正常样本的bam文件，以500bp为一个窗口，计算出每个窗口的read数量，过滤掉MAPQ不高于30和重复的序列。

在步骤A5中，对每一个窗口，使用read数量计算TRE。并依据TRE，通过GADA对基因组进行片段化。

在步骤A6中，对步骤A5得到的每个片段计算片段上所有TRE的均值和方差，以及该片段内包含的窗口数量。使用正态分布的方法，以每个片段的TRE均值和方差作为正态分布的均值和方差，将片段中的窗口按照正态分布进行平滑化。汇总所有片段平滑化后的TRE以及对应窗口数的信息。

在步骤A7中，使用自回归模型与GADA相结合的方法，找出P的2个候选值：420和262。

在步骤A8中，对步骤A7中的每一个P候选值，计算出Q值，根据P和Q，计算出癌症样本纯度γ和癌细胞倍性κ，并计算出所有片段的拷贝数

在步骤A9中，根据步骤A8计算出的片段拷贝数，使用EM算法对LAR数据进行计算，挑选出每个片段最佳的基因型。

在步骤A10中，根据步骤A8和A9的结果，计算出步骤A7找到的2个P候选值对应的BIC值，选择BIC值最小的作为最终结果，选择的P值为420，并输出该P值下的片段拷贝数及主要等位基因拷贝数。具体输出结果如图5所示。

本申请的有益效果在于：在计算拷贝数时，还得到了等位基因、基因型、非正常组织样本纯度、非正常组织细胞倍性等多个对拷贝数产生影响的参数，从而从全基因组测序数据获得更全面的信息，利用这些参数，能够使得拷贝数计算结果更加精确。

在一个实施例中，P值通过如下步骤B1-B2计算：

在步骤B1中，统计各基因组片段内的TRE的均值方差和窗口数量；

在步骤B2中，基于基因组片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值P。

本实施例中，以拼接后的基因组片段为单位，统计片段内的TRE的均值方差和窗口数量。基于片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值，即P值。

其中，在建立自回归模型时，使用上述片段内的TRE的均值和方差，对基因组每个片段的窗口数量进行平滑处理，使TRE的分布更平滑。平滑处理后的TRE分布图如图3所示。

在一个实施例中，基因组片段通过如下步骤C1-C3形成：

在步骤C1中，根据测序序列将全基因组测序数据划分为特定碱基数的窗口；

在步骤C2中，计算每个窗口的TRE值；

在步骤C3中，利用TRE值将相邻窗口进行合并，以形成新的基因组片段。

本实施例中，由于序列总数过大，所以需要划分窗口，以500个碱基为单位划分窗口，统计窗口的序列数量。计算每一个窗口的TRE(Tumor Read Enrichment，肿瘤片段富集程度)值，然后运用TRE，对基因组进行片段化，即将TRE值相似的窗口拼接成新的基因组片段。

在一个实施例中，如图6所示，上述步骤S12可被实施为如下步骤S61-S62：

在步骤S61中，根据差值P推算出拷贝数为2的片段对应的TRE位点Q；

在步骤S62中，根据TRE位点Q计算出所有P取值下非正常组织样本纯度和非正常细胞倍性，并计算所有P取值下各个片段的拷贝数。

本实施例中，根据差值P推算出拷贝数为2的片段对应的TRE位点Q；

对于自回归模型计算所得P值，选择[P-0.2，P+0.2]作为P值的取值范围。遍历所有P(例如，以千分之一为间隔进行遍历)，寻找每个P对应的first peak(第一个峰)，即P对应的第一个整数拷贝数所对应的TRE值，也即每个P点在图4中TRE柱状图(平滑处理后的TRE分布图)的第一个柱形对应的值。P值和first peak位点，推断出在first peak之前可能存在的peak(峰、即柱状图中的各个柱形)的个数。根据TRE分布规律，虚拟的peak与实际peak的可能排布情况，推算出拷贝数为2的片段对应的TRE位点，即为Q。

在计算出TRE位点Q之后，根据TRE位点Q计算出所有P取值下非正常组织样本纯度和非正常细胞倍性，并计算所有P取值下各个片段的拷贝数；具体的，在得到每个P值对应的TRE位点，即Q值之后，可以通过公式

计算出每个P取值下的癌症样本纯度γ；可以通过公式

计算出每个P取值下的癌细胞倍性κ。

在一个实施例中，方法还可被实施为如下步骤：

根据如下公式计算各个P取值下的非正常组织样本纯度和非正常组织细胞倍性：

其中σ_i表示拷贝数为i的所有片段中，其TRE的标准差。S_s为片段S的TRE观测值，Sⁱ为在某一P值下，理论上拷贝数为i的片段的TRE值。p_i表示片段S的拷贝数为i的权重，取值均为1。

本实施例中，计算其他peak与Q距离的周期数p，如图4所示，如果这个peak在Q的右边，则为正值，在左边则为负值，将p加2则得到这个peak代表的拷贝数，落在这个peak上的片段均为这个拷贝数；得到所有片段的拷贝数后，根据如下公式计算各个P取值下的非正常组织样本纯度和非正常组织细胞倍性：

其中σ_i表示拷贝数为i的所有片段中，其TRE的标准差。S_s为片段S的TRE观测值，Sⁱ为在某一P值下，理论上拷贝数为i的片段的TRE值。p_i表示片段S的拷贝数为i的权重，取值均为1。

在一个实施例中，上述步骤S13可被实施为如下步骤：

根据各个P取值下的非正常组织样本纯度和非正常组织细胞倍性、各个片段的拷贝数以及变异位点的LAR数据计算每个P值对应的各个片段的基因型估计值。

在一个实施例中，变异位点的LAR数据通过以下步骤D1-D2获取：

在步骤D1中，将全基因组测序数据与参考基因组进行比对，以获取比对文件；

在步骤D2中，从比对文件中获取每一个测序序列的位置信息、变异位点的LAR数据。

图7为本申请一总的实施例中，一种全基因组测序数据的分析方法的流程图，该方法包括以下步骤S701-S710：

在步骤S701中，获取来自同一癌症病人的癌症样本和正常组织样本的WGS数据，并将测序数据比对到参考基因组hs37d5或hg38。

在步骤S702中，从比对到参考基因组的BAM文件中提取每一个测序序列的位置信息，HGSNVs位点上的LAR数据和覆盖该位点的序列数量。

在步骤S703中，将全基因组划分为500碱基大小的窗口；

在步骤S704中，计算每一个窗口的TRE值，然后运用TRE将相邻窗口进行合并，形成新的基因组片段；

在步骤S705中，以新的基因组片段为单位，统计新的基因组片段内的TRE的均值方差和窗口数量；

在步骤S706中，基于片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值即P值；其中，建立自回归模型之前，使用新的基因组片段内的TRE的均值方差，对基因组每个片段的窗口数量进行平滑处理，使TRE的分布更平滑；

在步骤S707中，对于自回归模型计算所得P值，选择[P-0.2，P+0.2]作为P值的取值范围，遍历所有P，寻找每个P对应的first peak；

在步骤S708中，通过P值和firstpeak位点，推断出在firstpeak之前可能存在的peak的个数；

在步骤S709中，通过某个P取值下的γ和κ，片段总拷贝数，和每个片段上HGSNVs位点上的LAR数据，使用EM算法对公式(22)进行求解，得出在各P取值下每个片段的基因型估计值

和最大的LAR似然值L(r；k，l)。

在步骤S710中，通过I步骤得到所有可能P取值下的L(S_s；γ，κ)和L(r；k，l)后，通过对公式(24)求解，得出最优的P值，以及该最优P值下的癌症样本纯度γ和癌细胞倍性k，各个片段的拷贝数和等位基因拷贝数。

其中，上述S708包括以下步骤：

H1、根据TRE分布规律，虚拟的peak与实际peak的可能排布情况，推算出拷贝数为2的片段对应的TRE位点，即为Q。计算出各P取值下的癌症样本纯度γ和癌细胞倍数k。

H2、计算其他peak与Q距离的周期数p，如果这个peak在Q的右边，则为正值，在左边则为负值，将p加2则得到这个peak代表的拷贝数，落在这个peak上的片段均为这个拷贝数。如果拷贝数为浮点数，则说明这个peak上的片段为亚克隆片段。

H3、得到所有片段的拷贝数后，计算出该P值下的L(S_s；γ，κ)

在上述步骤S701中，采用了1000基因组计划第三期(phase 3)项目使用的参考基因组hs37d5。比对软件使用Burrows-Wheeler Aligner(BWA)，比对方法使用BWA中的mem，最终获得癌症和正常样本的比对结果bam格式文件。

上述步骤S702中，采用了rust语言的rust-htslib包提取测序read的位置信息。在提取序列的位置信息过程中，我们过滤掉序列比对质量(MAPQ)不高于30的序列，过滤掉重复序列，并仅仅提取正向的序列。在提取HGSNVs信息时，使用strelka2软件输出包含SNV信息的vcf文件，并输入提前准备一个包含SNP位点信息的bed格式文件。随后使用rust语言的rust-htslib包对vcf文件进行处理，提取出HGSNVs信息。本发明方法提前收集了1000基因组计划中根据大量样本统计出来的杂合等位基因位点，并且过滤掉B-allele frequence至小于0.05的位点，然后做成bed文件。该方法为后续步骤提供了充足的HGSNV位点，并大大加快了HGSNV位点的提取速度，提高了软件运行效率。

在上述步骤S703中，将全基因组划分为以500bp长度为单位的window(窗口)。合适大小的窗口可以使read数量的统计更精确。

上述步骤S704中，计算每个window的TRE取值。然后运用TRE的值，使用GADA软件对全基因组进行片段化。

上述步骤S706中，定义TRE的取值范围为[0，3]，以TRE的取值为X，以某个TRE对应的window的数量为Y。计算L＝Y_x×Y_p+x，从结果可以知道L的值随着P的改变程现出多峰分布。第二高峰对应的P值就是最合理的相邻拷贝数片段之间的TRE差值。自回归模型在找到P值得同时，会依据P值以及L的分布，给可能的P值划定范围，一般的范围为[P-0.2，P+0.2]，该范围内的TRE差值，都有可能成为实际的P值。

而在建立回归模型时，会对步骤S705中产生的片段的TRE进行平滑处理，针对每一个基因组片段，以TRE的均值作为正态分布的均值μ，以TRE的方差作为正态分布的方差σ，计算出TRE在[μ-2σ，μ+2σ]范围内的window分布。将所有片段的window根据TRE值smooth了之后，可以得到一个平滑的TRE多峰分布图。

上述步骤S707中，首先对[P-0.2，P+0.2]内的每一个TRE差值进行遍历。通过P值在[10，1250]的范围内，检索TRE分布图中的first peak，具体通过检索公式

进行检索，计算所得

即为TRE分布的first peak，在TRE的分布图中每一个peak对应一种拷贝数。其中，该检索公式中，P₀为从[P-0.2，P+0.2]中选取的P值，i指第几个峰，已知第一个峰的TRE值为

周期为P，那么第i个峰的TRE值为S_i＝l+P₀×i。将S_i带入下文中所提及的公式(11)中，可以构建TRE高斯分布模型。

上述步骤S708中，根据公式n＝floor(i/P)可计算出在first peak之前存在的虚拟的峰的个数。其中i表示first peak位点，P表示上一步计算出的TRE差值，floor表示对小数向下取整。例如当I＝980，P＝360时，n＝2，表示在first peak之前，TRE＝980-360＝620和TRE＝980-360*2＝260的位点可能存在两个虚拟的峰。

上述推断癌症样本拷贝数的方法步骤H1中，根据步骤S708中假设的情况，可推断拷贝数为0的片段可以对应TRE为260，620，980的peak，那么相应的拷贝数为2的片段对应的TRE的位点也只有三种可能。同时，在拷贝数为2的片段对应的TRE位点一定的情况下，TRE分布图上每一个峰都有唯一的拷贝数对应。而在步骤H2中，通过P和Q，以及其它peak距离Q的周期数，可以非常快速准确地计算出每个peak上片段的总拷贝数，以及识别片段是否为亚克隆片段。

在步骤S709中，对用于求解的EM算法进行了优化，μ₁与μ₂存在互为相反数的关系，对于每一个高斯组分，方差σ是相同的，μ_m的可能取值是有限个，最多为

个，在实际应用中，μ_m的取值不会超过6个，因此我们极大的缩减了EM算法的求解空间，提高了EM算法的求解速度，降低了EM算法不收敛及陷入局部最优解的概率，使得最终的结果更加可靠。另外，我们还对μ_m进行了校正，使得在低测序深度低肿瘤样本纯度的条件下，该算法依然能够正常运行。

上述步骤S710中，使用BIC来避免模型出现过度拟合的问题，提高了结果准确性。

下面，结合具体公式，对于本申请的具体实施方式进行详细描述如下：

对于基因组某一片断S，TRE的期望值为以下公式计算值：

公式(3)中，γ为癌症样本纯度(purity)，C_s为片段S的拷贝数，κ为癌症样本基因组的倍性(ploidy)。

当片段S的拷贝数为i时和i+1时：

通过公式(4)和公式(5)，我们可以发现对于相邻的拷贝数对应的片段，他们的TRE的差值如下：

此外，对于I＝2的片段，相应的TRE应该为：

通过上述(6)和(7)，可以解得癌症样本的纯度(γ)和倍性(κ)分别为：

可见，通过P和Q可以直接解出癌症样本纯度γ和倍性κ。本发明通过一种自回归模型与GADA相结合的方法，找出P的有限个候选取值。考虑到具有相同拷贝数的基因组片段的TRE值服从正态分布，并且HGSNV的等位基因频率也近似服从高斯分布。本方法通过对片段TRE的高斯分布模型，以及HGSNV的高斯模型如下所示，求混合模型的极大似然值，对P与Q进行极大似然估计，可以对TRE差值P进行较好的校正。TRE的高斯分布模型如公式11所示：

其中σ_i表示拷贝数为i的所有片段中，其TRE的标准差。S_s为片段S的TRE观测值，Sⁱ为在某一P值下，理论上拷贝数为i的片段的TRE值。p_i表示片段S的拷贝数为i的权重，取值均为1。在确定P和Q后，我们通过计算每一个片段S距离Q的周期数+2来得出片段S的总拷贝数。

另一方面，假设在某个HGSNV位点上癌细胞和正常样本中正常细胞的基因型分别为(k，l)和(1，1)，k和1分布代表主要等位基因拷贝数和次要等位基因拷贝数，将癌症样本纯度γ考虑进去，则这个HGSNV位点上癌细胞基因型观测值为(γk+(1-γ)×1，γl+(1-γ)×1)，而正常样本中正常细胞的基因型观测值仍然为(1，1)，由于不能从测序数据中识别出主要等位基因是否是参考基因(reference allele)，因此LAR的朴素期望值(

expectations)有两个，定义为：

由于公式(2)中没有考虑到样本中HGSNVs存在reads读数为0的情况，导致公式(12)的LAR朴素估计值会有较大的偏差，因此我们需要对公式(12)进行校正。假设基因型的测序深度服从泊松分布，λ_k+l＝λ_k+λ_l代表总的片段测序深度，λ_k和λ_l分别代表主要等位基因和次要等位基因的测序深度，则：

我们以癌症样本中HGSNVs位点上测序深度的中位数值来估计λ_k+l。用d^k和d^l分别代表主要基因和次要基因的观测读数，并且均服从泊松分布，则：

d^k～P_o(λ_k)，d^k～Po(λ_l) 公式(16)

排除读数为0的数据后，得到校正后的LAR期望值μ₁和μ₂分别定义为：

在实际数据中，我们无法得知观测到的LAR的期望值是μ₁还是μ₂，因为我们无法得知这个基因是主要等位基因还是次要等位基因。因此我们假设LAR的观测值服从两组分的高斯混合分布，并且这两个组分有相同的方差，均值互为相反数：

其中α_m代表组分m所占的比例且满足α₁+α₂＝1，

代表这两个组分在片段s上的方差。

我们使用EM算法对公式(18)中的α_m和

进行求解，为此，我们对一个LAR数据i，引入一个隐变量Δ_i，代表这个LAR数据属于哪一个组分。已知片段s总拷贝数为C_s，在这个片段上有N^s个LAR数据，对于每一个可能的等位基因型组合

μ_m使用公式(17)进行计算，然后计算LAR数据属于哪一个组分的E步为：

用于更新各个参数的M步为：

E步和M步反复迭代直至参数收敛，随后我们计算等位基因型(k，l)的似然值：

EM算法用于求解最优的α_m和

并得出每个可能等位基因型组合(k，l)的似然值。我们挑选似然值最大的等位基因型组合作为片段s的等位基因型

接下来，我们使用BIC将所有片段最大LAR似然值和TRE的似然值进行组合，求出最可能的癌症样本纯度，癌症细胞倍性，基因片段的总拷贝数和等位基因拷贝数：

其中，I代表可能的TCNs的总数，M代表观测到的TREs的总数，J代表可能的ASCNs的总数，N代表HGSNVs的总数。本算法通过使用自回归求得P和Q的可能取值后，通过公式(9)(10)得到对应的可能的癌症样本纯度γ和癌细胞倍性κ，进一步推出可能的片段拷贝数，得到L(S_s；γ，κ)，同时对于可能的片段拷贝数，使用EM求得可能的等位基因型及L(r；k，l)，最后使用最小化公式(23)的参数最为最终的估计结果：

图8为本申请实施例的一种全基因组测序数据的分析装置的框图，如图8所示，该装置包括以下模块：

收集模块81，用于从预设数据库收集来自同一受检对象的全基因组测序数据，全基因组测序数据包括正常组织样本数据和非正常组织样本数据；

第一计算模块82，用于计算受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性；其中，P值为相邻拷贝数对应的基因组片段的TRE的差值；

第二计算模块83，用于根据每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型；

第一获取模块84，用于通过BIC算法获取最优P值；

第二获取模块85，用于获取最优P值对应的各个片段的基因型、非正常组织样本纯度γ、非正常组织细胞倍性κ、各个片段的拷贝数以及等位基因拷贝数。

在一个实施例中，P值通过如下方式计算：

统计各基因组片段内的TRE的均值方差和窗口数量；

基于基因组片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值P。

在一个实施例中，基因组片段通过如下方式形成：

根据测序序列将全基因组测序数据划分为特定碱基数的窗口；

计算每个窗口的TRE值；

利用TRE值将相邻窗口进行合并，以形成新的基因组片段。

在一个实施例中，如图9所示，第一计算模块82，包括：

推算子模块91，用于根据差值P推算出拷贝数为2的片段对应的TRE位点Q；

第一计算子模块92，用于根据TRE位点Q计算出所有P取值下非正常组织样本纯度和非正常细胞倍性，并计算所有P取值下各个片段的拷贝数。

在一个实施例中，还包括：

第三计算模块，用于根据如下公式计算各个P取值下的非正常组织样本纯度和非正常组织细胞倍性：

其中σ_i表示拷贝数为i的所有片段中，其TRE的标准差。S_s为片段S的TRE观测值，Sⁱ为在某一P值下，理论上拷贝数为i的片段的TRE值。p_i表示片段S的拷贝数为i的权重，取值均为1。

在一个实施例中，第二计算模块，包括：

第二计算子模块，用于根据各个P取值下的非正常组织样本纯度和非正常组织细胞倍性、各个片段的拷贝数以及变异位点的LAR数据计算每个P值对应的各个片段的基因型估计值。

在一个实施例中，变异位点的LAR数据通过以下方式获取：

将全基因组测序数据与参考基因组进行比对，以获取比对文件；

从比对文件中获取每一个测序序列的位置信息、变异位点的LAR数据。

以上实施例仅为本申请的示例性实施例，不用于限制本申请，本申请的保护范围由权利要求书限定。本领域技术人员可以在本申请的实质和保护范围内，对本申请做出各种修改或等同替换，这种修改或等同替换也应视为落在本申请的保护范围内。

Claims (10)

1.一种全基因组测序数据的分析方法，其特征在于，包括：

从预设数据库收集来自同一受检对象的全基因组测序数据，所述全基因组测序数据包括正常组织样本数据和非正常组织样本数据；

计算所述受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性；其中，所述P值为相邻拷贝数对应的基因组片段的TRE的差值；

根据所述每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型；

通过BIC算法获取最优P值；

获取所述最优P值对应的各个片段的基因型、非正常组织样本纯度γ、非正常组织细胞倍性κ、各个片段的拷贝数以及等位基因拷贝数。

2.如权利要求1所述的方法，其特征在于，所述P值通过如下方式计算：

统计各基因组片段内的TRE的均值方差和窗口数量；

基于基因组片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值P。

3.如权利要求2所述的方法，其特征在于，所述基因组片段通过如下方式形成：

根据所述测序序列将全基因组测序数据划分为特定碱基数的窗口；

计算每个窗口的TRE值；

利用TRE值将相邻窗口进行合并，以形成新的基因组片段。

4.如权利要求3所述的方法，其特征在于，所述计算所述受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性，包括：

根据所述差值P推算出拷贝数为2的片段对应的TRE位点Q；

根据所述TRE位点Q计算出所有P取值下非正常组织样本纯度和非正常细胞倍性，并计算所有P取值下各个片段的拷贝数。

5.如权利要求1所述的方法，其特征在于，还包括：

根据如下公式计算各个P取值下的非正常组织样本纯度和非正常组织细胞倍性：

其中σ_i表示拷贝数为i的所有片段中，其TRE的标准差。S_s为片段S的TRE观测值，Sⁱ为在某一P值下，理论上拷贝数为i的片段的TRE值。p_i表示片段S的拷贝数为i的权重，取值均为1。

6.如权利要求5所述的方法，其特征在于，根据所述每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型，包括：

根据所述各个P取值下的非正常组织样本纯度和非正常组织细胞倍性、所述各个片段的拷贝数以及变异位点的LAR数据计算每个P值对应的各个片段的基因型估计值。

7.如权利要求6所述的方法，其特征在于，所述变异位点的LAR数据通过以下方式获取：

将所述全基因组测序数据与参考基因组进行比对，以获取比对文件；

从比对文件中获取每一个测序序列的位置信息、变异位点的LAR数据。

8.一种全基因组测序数据的分析装置，其特征在于，包括：

收集模块，用于从预设数据库收集来自同一受检对象的全基因组测序数据，所述全基因组测序数据包括正常组织样本数据和非正常组织样本数据；

第一计算模块，用于计算所述受检对象全基因组测序数据中每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性；其中，所述P值为相邻拷贝数对应的基因组片段的TRE的差值；

第二计算模块，用于根据所述每个P值下各个片段的拷贝数、非正常组织样本纯度和非正常组织细胞倍性计算每个P值对应的各个片段的基因型；

第一获取模块，用于通过BIC算法获取最优P值；

第二获取模块，用于获取所述最优P值对应的各个片段的基因型、非正常组织样本纯度γ、非正常组织细胞倍性κ、各个片段的拷贝数以及等位基因拷贝数。

9.如权利要求8所述的装置，其特征在于，所述P值通过如下方式计算：

统计各基因组片段内的TRE的均值方差和窗口数量；

基于基因组片段内的TRE的均值方差和窗口数量建立自回归模型，计算相邻拷贝数对应的基因组片段的TRE的差值P。

10.如权利要求9所述的装置，其特征在于，所述基因组片段通过如下方式形成：

根据所述测序序列将全基因组测序数据划分为特定碱基数的窗口；

计算每个窗口的TRE值；

利用TRE值将相邻窗口进行合并，以形成新的基因组片段。

Images