CN108021788B

CN108021788B - 基于细胞游离dna的深度测序数据提取生物标记物的方法和装置

Info

Publication number: CN108021788B
Application number: CN201711276870.1A
Authority: CN
Inventors: 万季; 宋麒; 潘有东; 崔好; 韩艳; 高峰
Original assignee: Beijing Xinhe Ruien Biomedical Technology Co ltd
Current assignee: Beijing Xinhe Ruien Biomedical Technology Co ltd
Priority date: 2017-12-06
Filing date: 2017-12-06
Publication date: 2022-08-05
Anticipated expiration: 2037-12-06
Also published as: CN108021788A

Abstract

提供一种基于细胞游离DNA的深度测序数据提取生物标记物的方法和装置。方法包括由处理器执行的如下步骤：获取不同类别的样本的细胞游离DNA的深度测序数据；计算每个样本数据在每个碱基位置的各变异的频率值；基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布；确定频率分布在类别间具有充分差异的变异；以及基于所确定的变异提取所述生物标记物。该方法可以极大限度地挖掘cfDNA数据中的信息。

Description

基于细胞游离DNA的深度测序数据提取生物标记物的方法和装置

技术领域

本发明涉及分子生物学技术领域，具体涉及基于细胞游离DNA的深度测序数据提取生物标记物的方法以及装置。

背景技术

人类的疾病直接或间接与基因有关(文中又称此种类型的疾病为基因疾病)，如所有的遗传病，都是由于基因结构异常或基因表达异常所引起的；有些疾病则是环境因素与遗传因素综合作用的结果，如吸烟引起肺癌。基因异常所致疾病(即，基因疾病)大致分为以下三类：1.单基因病：在一个基因位点上存在缺陷，如镰刀红细胞贫血。2.多基因病：涉及一个以上的基因以及基因与环境因素的相互作用，如肿瘤、心血管病、代谢性疾病、神经和精神类疾病、免疫性疾病等。3.获得性基因病：这类疾病由病原微生物感染引起，不符合孟德尔遗传规律，但可发生病原微生物基因组与人类基因组的相互作用，这类疾病多数涉及人类基因组结构及表达功能的改变，如爱滋病、鼻咽癌等。

cfDNA、ctDNA和CTC的检测作为当前比较热门的液体活检的手段，可应用于疾病的诊断。其中，cfDNA(cell-free DNA，细胞游离DNA)，即血液或其他体液中游离的DNA片段。在基因病，例如肿瘤的发展过程中，肿瘤细胞会通过各种机制将其DNA释放到循环系统中(如血液、尿液等)。肿瘤DNA上面携带了大量与肿瘤发生发展相关的突变信息。这些来自肿瘤细胞的DNA片段和血液(或者其他体液)里来自各种正常细胞的DNA片段混合即cfDNA，而来自肿瘤的DNA片段被称为ctDNA(circulating tumor DNA, 循环肿瘤DNA)片段。

ctDNA由于其非侵入式的特点，可用于检测多种类型肿瘤的疾病状况。当前ctDNA的检测主要分为数字PCR(digital polymerase chain reaction，数字聚合酶链式反应)和NGS(next-generation sequencing，下一代基因测序) 两大主要分支。数字PCR可以从突变比例极低的样本中检测到指定的突变以及突变的绝对含量信息，具有很高的灵敏度和可靠性，是当前ctDNA检测的金标准。但是这一方法需要预先指定待检测的突变位点，同时检测的位点数目受仪器的承载容量限制。与数字PCR不同，NGS可以对基因组的指定区域进行扫描，检测范围不再局限于指定的具体突变，这样有助于覆盖更多的检测区域和突变。但是，现有的以肿瘤组织突变数据为参照从血液(或者其他体液)中寻找肿瘤DNA生物标记物的方案并不完善。首先，即便是同一肿瘤类型，相关的突变在病人间的重现率并不高。其次，同一病人肿瘤组织与ctDNA突变信息的一致性关系仍有待大规模的实验验证。由于生物上的和检测技术上的诸多因素，在病人肿瘤组织样本中检出的突变无法在(血液或其他体液)ctDNA样本中检出。肿瘤突变的低重现率和ctDNA突变缺失的两大问题，限制了ctDNA在临床，尤其是早期诊断方面的推广。

发明内容

针对ctDNA临床应用中存在的上述问题，本发明综合考虑了血液中游离 DNA的来源和深度测序的数据特点，提出了一套提取cfDNA(而非ctDNA) 特有的生物标记物的方法。

根据本发明的一个方案，提供了一种由计算机实现的基于不同类别的细胞游离DNA的深度测序数据提取生物标记物的方法，包括由处理器执行的如下步骤：获取不同类别的样本的细胞游离DNA的深度测序数据；计算每个样本数据在每个碱基位置的各变异的频率值；基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布；确定频率分布在类别间具有充分差异的变异；以及基于所确定的变异提取所述生物标记物。

优选地，所述样本为体液样本。

优选地，基于所确定的变异提取所述生物标记物的步骤可以包括：提取所确定的变异所在碱基位置的所有变异的频率分布，作为所述生物标记物。

作为替代，基于所确定的变异提取所述生物标记物的步骤可以包括：提取所确定的变异的频率分布，作为所述生物标记物。

优选地，确定频率分布在类别间具有充分差异的变异的步骤可以包括：计算变异的频率分布在不同类别之间的差异程度；在所述变异的差异程度满足预定条件的情况下，认为所述变异的频率分布在类别间具有充分差异。

优选地，所述差异程度可以包括频率分布的距离，所述预定条件可以包括：所述频率分布的距离大于第一预定阈值。

作为替代，所述差异程度可以包括P值，所述预定条件可以包括：所述P值小于第二预定阈值。

优选地，所述碱基位置可以为靶向区域中的碱基位置。

优选地，所述变异包括如下类型中的一种或更多种：单核苷酸变异类型、插入变异类型、删除变异类型。

优选地，所述类别为肿瘤类别、关于预后的类别、与关于治疗反应的类别中的任何一个类别。

根据本发明的另一方案，提供一种由计算机实现的基于不同类别的细胞游离DNA的深度测序数据提取生物标记物的装置，包括：样本数据获取模块，其配置为获取不同类别的样本的细胞游离DNA的深度测序数据；频率值计算模块，其配置为计算每个样本数据在每个碱基位置的各变异的频率值；频率值分布提取模块，其配置为基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布；分布差异程度计算模块，其配置为确定频率分布在类别间具有充分差异的变异；以及生物标记物提取模块，其配置为基于所确定的变异提取所述生物标记物。

优选地，所述生物标记物提取模块进一步配置为：提取所确定的变异所在碱基位置的所有变异的频率分布，作为所述生物标记物；或者提取所确定的变异的频率分布，作为所述生物标记物。

本发明再一方面还提供一种计算设备，包括处理器，所述处理器配置为运行计算机指令以实现以下步骤：获取不同类别的样本的细胞游离DNA的深度测序数据；计算每个样本数据在每个碱基位置的各变异的频率值；基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布；确定频率分布在类别间具有充分差异的变异；基于所确定的变异提取所述生物标记物。

优选地，所述处理器进一步配置为执行计算机指令以实现以下步骤：利用所提取的生物标记物作为特征，来训练分类模型。

优选地，所述处理器进一步配置为执行计算机指令以实现以下步骤：基于待分类的体液样本的细胞游离DNA深度测序数据，利用训练好的所述分类模型，来确定分类结果。

本发明还一方面提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序使计算机执行以下步骤：获取不同类别的样本的细胞游离 DNA的深度测序数据；计算每个样本数据在每个碱基位置的各变异的频率值；基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布；确定频率分布在类别间具有充分差异的变异；基于所确定的变异提取所述生物标记物。

优选地，所述计算机程序使计算机进一步执行以下步骤：将所提取的生物标记物作为特征训练分类模型。

本发明的方案充分考虑了cfDNA等位基因频率(allele frequency)的重要生物意义，并将等位基因频率的考察对象范围拓宽至基因区域上面每个碱基位以及各种变异类型，因此可以极大限度地挖掘cfDNA数据中的信息。

与现有技术相比，本发明的方案具有如下优势：

1、从来源上讲，通过本发明的方案发现的变异既可以来自肿瘤释放的 DNA，也可以来自其他各类游离DNA片段(如各正常组织释放的DNA、肿瘤所在组织的正常部分释放的DNA、以及其他外源的DNA片段)的变异，因此可以用于发现多种基因疾病的变异，而不局限于肿瘤的变异。当前常用的cfDNA分析方法主要沿袭肿瘤组织DNA的处理流程，通过体细胞变异识别(somatic variant call)找出与正常组织序列信息不同的变异位点。这实质上是将分析对象从cfDNA片段限制在了携带体细胞突变的ctDNA片段。

2、本发明的方案不需要体细胞变异识别步骤，因此免去了用作对照的正常组织样本(一般为血液白细胞)，降低了成本和实施的复杂度。

3、本发明对cfDNA库制备过程中引入的碱基置换错误不敏感。当前的游离DNA库制备步骤均会以一定的模式引入微量的碱基置换(与碱基位置和参考碱基类型相关)。若分析程序没有考虑这种情况，结果可能会产生变异识别的假阳性，如果分析程序考虑了背景的碱基错误(且不论实现的复杂度和所需的训练数据量)，可能会导致变异识别的假阴性。本发明直接讨论不同肿瘤类型间的变异分布的迁移，其变化差异受碱基错误影响小，这是因为：来自不同肿瘤的样本受同样模式的碱基置换错误影响，计算差异时这些错误可以相互抵消。

附图说明

图1是本发明实施例的提取生物标记物的方法的流程图。

图2是本发明另一实施例的提取生物标记物(例如，有关肿瘤类别的生物标记物)的方法的流程示意图。

图3(a)和图3(b)是示例性示出评估图2所示实施例中提取的生物标记物对肿瘤状况的区分效果的图。

图4是评估本发明实施例的方法的技术效果的工作特征曲线(ROC)指标的图。

图5是本发明再一实施例的提取生物标记物的装置的结构框图。

图6示出本发明实施例的提取生物标记物的方法中频率值的计算的示意图。

具体实施方式

为使本领域技术人员更好地理解本发明，下面参照附图对本发明的实施例进行详细说明。本文具体实施例中所用的“首先”、“其次”、“再次”、“然后”、“接下来”等表达并不旨在对先后顺序进行限定。

图1是本发明实施例的由计算机实现的基于不同类别的细胞游离DNA 的深度测序数据提取生物标记物的方法的流程图。如图所示，该方法包括由处理器执行的如下步骤：

S101，获取不同类别的样本的细胞游离DNA的深度测序数据。

具体地，获取来自不同类别的多个样本(例如，体液样本)的细胞游离 DNA深度测序数据作为样本数据，样本数据中包括多个重叠或部分重叠的短读序列；

其中，“不同类别”的具体内容在不同情况下可以是不同的，并且例如，“不同类别”可以是指根据所期望提取的生物标记物的特性而预先确定的多个类别。例如，当期望提取的生物标记物具有能够表达有无肿瘤存在的特性时，该多个类别可以包括：无肿瘤和肿瘤这两个类别；又例如，当期望提取的生物标记物能够更细致地表达肿瘤存在并且细分肿瘤的类型时，该多个类别可以包括：无肿瘤；肿瘤1；肿瘤2；…肿瘤N(N为自然数)；此外，例如，当期望提取的生物标记物具有能够表达预后(具体地，例如为生存时间)的特性时，该多个类别可以包括：生存时间1；生存时间2；…；生存时间N；另外，例如，当期望提取的生物标记物具有能够表达治疗反应(具体地，例如为治疗效果)的特性时，该多个类别可以包括：治疗效果良好；治疗效果一般；治疗效果差；等等。

利用提取患者的细胞游离DNA并对其深度测序的方法，可以获取体液样本的细胞游离DNA深度测序数据，所获取的样本数据中包括多个重叠或部分重叠的短读序列，重叠程度的不同与测序的深度有关，由于篇幅所限，在此不一一赘述。

S102，计算每个样本数据在每个碱基位置的各变异的频率值。

具体地，针对每一体液样本的样本数据，计算在多个碱基位置中的每一碱基位置处，各变异在当前样本数据中的包含当前碱基位置的短读序列中出现的次数所占包含当前碱基位置的短读序列总数的频率值；

具体而言，假设对于当前碱基位置P1计算该频率值f_A。图6示例性地示出本发明实施例的提取生物标记物的方法中频率值的计算的示意图。

首先，对包含当前碱基位置P1的短读序列的数目进行计数，得到包含当前碱基位置P1的短读序列总数m1，假设图中m1＝10(该数值仅为示例)。

其次，对于当前碱基位置P1的具体的变异，比如在前碱基位置P1，如图所示，当碱基为A(这里，A与参考碱基相同，这种非变异类型在计算时也被当成是一种变异类型进行处理)时，对在当前碱基位置P1处变异A在这m1个短读序列中出现的次数c_1A计数，图中可以看到由于在该碱基位置 P1处，在这m1个短读序列中变异为A出现了两次，因此c_1A＝2。

然后，计算c_1A/m1(即，＝0.2)作为当前碱基位置P1处变异为A时的频率值。

接下来，以同样的方法，分别计算当前碱基位置P1处，变异为C、G、 T时的频率值。此外，还计算在当前碱基位置P1处的m1短读序列中，如果有插入类型的变异(图中未示出)以及删除类型的变异(图中未示出)，它们各自的频率值。在一优选实施例中，对于插入类型的变异，不同个数的插入碱基作为不同类型的变异来计算频率值，例如，插入一个碱基和插入两个碱基是两种类型的插入变异，分别计算频率值。可以理解，对于每一碱基位置，所计算出的各频率值之和为1。

对于所要计算的多个碱基位置的每一碱基位置，进行如上的同样的计算。

接下来，在S103，基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布。

具体地，对于不同类别中的每一类别，对属于该类别的每一体液样本的样本数据在每一碱基位置上各变异的频率值进行复合计算，得到每一类别在每一碱基位置上各变异的频率值分布。

例如，假设对于第一类别有1000个体液样本的样本数据，则对于碱基位置P1，将这1000个体液样本的样本数据中对应于碱基位置P1的上面计算的变异A的频率值均提取出来，进行诸如求平均，或者用这些样本数据中变异 A在碱基位置P1上的总出现次数除以该碱基位置P1上的短读序列总数，再除以样本数，来得到第一类别下碱基位置P1的变异为A的频率值分布。同样地，将这1000个体液样本的样本数据中对应于碱基位置P1的上面计算的变异C的频率值均提取出来作为第一类别下碱基位置P1的变异为C的频率值分布；等等。

接下来，在S104，确定频率分布在类别间具有充分差异的变异。

具体地，基于频率值分布，通过计算在不同类别之间的每一对应碱基位置处各变异的分布差异程度，来确定频率分布在类别间具有充分差异的变异。

更具体地，对于变异C，在各个类别下计算变异C所有碱基位置上的频率分布，从而得到变异C在各个类别下的分布，计算变异C在第一类别下的分布与第二类别下的分布之间的差异程度。计算两个分布之间的差异程度可使用计算两个分布之间的距离、或计算两个分布之间的P值等等来进行。更具体地，例如可采用对两分布距离的计算可以采用Kolmogorov-Smirnov test (KS test)进行或者计算巴氏距离、欧式距离等等，本发明对此不做限制。

当存在多个类别时，分别计算多个类别中任意两类别之间，各变异在每一对应碱基位置处的分布的差异程度。

S105，基于所确定的变异提取所述生物标记物。

具体地，根据分布差异程度满足预定条件的变异提取生物标记物。

更具体地，如果某一变异在所有碱基位置上的频率分布的差异程度满足预定条件，则认为该变异能够很好地区分两种类型，提取该变异所在碱基位置的所有变异的频率分布作为生物标记物。

或者，在一可选实施例中，如果某一变异在所有碱基位置上的频率分布的差异程度满足预定条件，提取该变异的该频率分布作为生物标记物。

如上，衡量分布差异程度的参数有距离、P值等。当参数为距离时，预定条件例如为：分布的距离大于第一预定阈值(可根据经验设定第一阈值的具体数值，在后面的实施例中会描述)，或者距离的P值小于第二预定阈值 (可根据经验设定第二阈值的具体数值，在后面的实施例中会描述)。

通过以上的实施例，能够提取得到对于区分不同类别具有显著意义的生物标记物(特征)。

实施例使用了细胞游离DNA深度测序数据作为样本数据，可用于发现与基因疾病的具体类型对应的生物标记物，也可用于发现与基因疾病的其他方面(如预后、治疗反应、疾病有无)的类型相关的生物标记物。

图2是本发明另一实施例的提取诊断性生物标记物(例如，有关肿瘤类别的生物标记物)的方法的流程示意图。

注意的是，本实施例的意图不在于限制所提取的生物标记物为有关肿瘤类别的生物标记物，显而易见的是，本发明的方案也可适用于提取其他相关的生物标记物，如预后生物标记物(prognostic biomarker)或者预测治疗反应的生物标记物(predictivebiomarker)。

下面以从cfDNA中寻找肿瘤类别相关的生物标记物为例，详细地介绍本发明的一方法实施例，但是本发明不限于肿瘤类别。

S201、准备一批不同肿瘤类型患者的cfDNA原始测序数据，将短读序列定位到基因组，并进行数据清理，如图2中左上部分所示。

不同肿瘤类型的cfDNA原始数据指的是从患有不同肿瘤的患者的血液或者组织液样本中提取游离DNA并进行测序而得到的短读序列文件。该序列文件的常见格式为FASTQ，但是也可以为bam或者其他格式。

测序时定位的基因组靶向区域(target region)需要预先指定，以便引物的设计。通常可以选择常见的高频突变肿瘤相关基因的外显子及启动子区域作为靶向区域。常规的实验流程包括库制备(以Capp-Seq为金标准)和深度测序(例如，>5000x)。本发明实施例的方法优选地要求采集这些数据的原始实验条件完全相同，即，使用同样的样本提取、库制备流程、试剂和测序仪器设定，以确保后续提取的生物标记物的精度。

本发明实施例的方法的肿瘤类型既可以包含肿瘤大类，比如肺癌、结肠癌、肾癌等，也可以包含肿瘤的亚类，比如非小细胞肺癌中的肺腺癌(LUAD) 与肺鳞状细胞癌(LUSC)。本发明实施例的方法将“正常(无肿瘤)”也视作一种肿瘤类型。本发明实施例的方法的不同肿瘤的种类数目没有具体限定，在有正常人cfDNA数据的情况下，可以通过比较某类肿瘤患者与正常人的 cfDNA数据，找出的生物标记物适用于区分该类肿瘤病症与正常状况；如果没有正常人cfDNA数据，则对肿瘤类型按两两成对分析，所找出的生物标记物适用于区分这两种肿瘤，在有N种肿瘤的场合下，可以得到N(N-1)/2组生物标记物。

不同肿瘤来源的肿瘤ctDNA原始数据的数量没有具体限定，但是理论上来说，每种肿瘤的样本数量越多，结果的可靠性和可重复性越能够得到保证。

随后可以对上述肿瘤ctDNA原始数据与参考基因组数据进行序列比对(alignment)。参考基因组数据是指人类各染色体上的碱基序列数据，通常为 FASTA格式。参考基因组数据既可以通过UCSC，Broad Institute等机构下载 (常见版本有hg19/GRCh37,hg38/GRCh38)，也可以自己添加。

进行序列比对是为了定位短读序列在参考基因组上面的位置，例如可以采用比对软件进行序列比对，例如可以使用bwa、bowtie，类似的软件还包括NovoAlign等。

例如，对bwa，其示例命令为：

bwa mem-M-t 4ref.fa R1.fastq.gz R2.fastq.gz|samtools view-Shb-osample.aln.bam– (R1.fastq.gz和R2.fastq.gz为来自同一样本的配对测序(paired-endsequencing) 的原始文件)

或者

bwa mem-M-t 4ref.fa single_read.fastq.gz|samtools view-Shb-osample.aln.bam–(single_read.fastq.gz为单短读测序(single read sequencing)的原始文件)

其中，

mem命令表示bwa mem算法，这是现在推荐使用的主流算法；

-M表示生成格式会兼容下游PICARD工具；

-t表示调用的线程数目；

ref.fa指的是参考基因组数据文件；

samtools view表示samtools工具的view命令，此处用于对bwa mem结果做进一步的筛选；

-Shb表示输出文件的注释调整，对结果分析没有直接影响；

–o表示输出文件名称；

aln.bam为进行序列比对和基因组定位生成的bam格式文件。

上述命令的输入文件为测序原始文件(FASTQ或其压缩格式)和参考基因组ref.fa(FASTA格式)。原始文件的格式与实验设置相关，对分析流程没有直接影响。最后生成的文件为bam格式文件。

为了用于清除或修复数据中的错误，例如可以对序列比对和基因组定位所得文件进行数据清理。常规的肿瘤组织DNA数据分析数据清理流程包括：去重复(Deduplicates)、插入删除重比对(Indel realignment)和碱基重校正 (base recalibration)。由于深度测序数据从统计上会倾向形成重复的短读序列，因此一些学者不建议对cfDNA数据使用去重复这一步骤。本发明方法也遵循了这一观点，仅使用插入删除重比对和碱基重校正。以下以GATKv3.7 工具集为例，来说明各步骤的参数设置。

插入删除重比对用于处理序列比对过程中因碱基插入/删除突变造成的短读序列定位混乱的问题。此步骤一般使用GATK的IndelRealignment或者 ABRA。以IndelRealignment为例，示例代码如下：

java-jar gatk.jar-T RealignerTargetCreator

-R ref.fa

-I sample.aln.bam

-o realigner.intervals

-dcov 15000

--knownSites 1000G_phase1.indels.hg19.sites.vcf

--knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf

java-jar gatk.jar-T IndelRealigner

-I sample.aln.bam

--targetIntervals realigner.intervals

-o sample.realign.bam

--maxReadsForRealignment 60000

--knownSites 1000G_phase1.indels.hg19.sites.vcf

--knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf

其中gatk.jar为GATK工具可执行文件，RealignerTargetCreator一步的-R 指明参考基因组FASTA文件，-I指明输入文件为序列对比步骤得到的sample.aln.bam文件，-o指明该步骤生成的插入删除重对比的检查区域，-dcov 设定了downsampling执行的短读序列覆盖数目，对cfDNA深度测序而言，该值有必要设定为一个较高的值，避免downsample造成的信息丢失， --knownSites指明公共的含插入和删除变异的数据库。IndelRealigner一步的-I指明输入文件序列对比步骤生成的sample.aln.bam文件，--targetIntervals 指明检查的区域为RealignerTargetCreator步骤生成的realigner.intervals，-o 指明该步骤输出为sample.realign.bam，--maxReadsForRealignment指明用于重比对的最大的短读序列数目，对于深度测序而言，该值有必要设定为一较高的值。--knownSites指明公共的含插入和删除变异的数据库。

此外，GATK的IndelRealigner对较长(>30bp)的插入删除识别率不高，使用ABRA来替代GATK IndelResligner可以提供更多的对插入删除变异的识别。

碱基重校正用于调整短读序列的碱基质量分数，使之能更准确反映测序过程中的各种错误。这一步骤可以使用GATK的BaseRecalibration工具实现。示例代码如下：

java-jar gatk.jar-T BaseRecalibrator

-I input.realign.bam

-L target_region.bed

-R ref.fa

-o recal.table

--knownSites 1000G_phase1.indels.hg19.sites.vcf

--knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf

--knownSites dbsnp_138.hg19.vcf

--indels_context_size 6

--mismatches_context_size 6

-rf BadCigar

-rf MappingQualityZero

-rf MappingQualityUnavailable

-rf NotPrimaryAlignment

Java-jar gatk.jar-T PrintReads

-R ref.fa

-I sample.realign.bam

-BQSR recal.table

-o sample.bqsr.bam

其中gatk.jar为GATK程序。在BaseRecalibrator一步中，-I表明输入文件为插入删除重校正步骤得到的bam文件sample.realign.bam，-L是靶向区域文件(一般为bed或者interval格式),-R指明参考基因组的FASTA文件， -o指明输出的重校正表格记录，该记录仅用作中间临时文件，--knownSites 指明需要使用的公共变异数据库，--indels_context_size指明评估碱基插入和删除所使用的k-mer的长度，而--mismatches_context_size指明评估碱基不匹配所使用的k-mer的长度。在PrintReads一步中，-R是参考基因组的FASTA 文件，-I是插入删除重校正步骤得到的bam文件sample.realign.bam，-BQSR 指明BaseRecalibrator一步得到的表格记录，-o指明本步生成的bam文件 sample.bqsr.bam

经数据清理后的所得文件sample.bqsr.bam用于后续步骤。

ST202、计算每个样本数据在每个碱基位置的各变异的频率值。

在靶向区域上，检查前述步骤清理后的bam文件在每个碱基位置上的各变异类型的计数。

靶向区域指的是实验引物设计阶段优化的基因组区域，一般是肿瘤突变高发的区域。

本发明实施例所讨论的变异类型包括：、单核苷酸变异类型(SNV)、插入/删除(indel)变异类型。对于结构变异(structural variation,SV)，拷贝数变异(copy numbervariation,CNV)等，本发明实施例的方法不做讨论。

计数(在每一碱基位置处，各变异在当前样本数据中的包含当前碱基位置的短读序列中出现的次数)例如可以采用多种工具实现，例如，可以直接使用samtools从bam文件生成mpileup文件，然后用单独的脚本对mpileup 结果进行计数，也可以使用其他基于samtools的第三方工具，如iDES软件里的ides-bam2freq.pl(https:// cappseq.stanford.edu/ides/download.php)脚本， R库PlasmaMutationDetector等。上述所使用的工具的主要差别在于应用的数据类型(如PlasmaMutationDetector主要用于iontorrent产生的数据)和程序界面，对结果本身影响不大。

以ides-bam2freq.pl为例，程序的界面为：

perl ides-bam2freq.pl[options]input.bam(s)ref.fa targets.bed

ides-bam2freq.pl为统计变异类型计数的perl脚本，可从https:// cappseq.stanford.edu/ides/download.php下载；

[options]表示可用于筛选符合一定质量要求的碱基和短读序列的参数；

input.bam(s)表示数据清理步骤之后得到的bam文件；

ref.fa为参考基因组FASTA格式的数据文件；

targets.bed提供需要检查的基因组区域。

程序的输入为bam文件或者bam文件存放的文件夹，ides-bam2freq.pl 最后生成文本文件，包括如下信息：染色体编号，基因组位置，该位置的短读序列数目，参考碱基，支持参考碱基的来自+/-链的短读序列数目，支持其他每类变异的来自+/-链的短读序列数目。

使用下列公式将每个碱基位置上支持变异(其中，在计算时,参考碱基(即，非变异类型)也被视为变异的一种)的短读序列的数目值转化为频率值：

频率值＝支持当前变异的短读序列数目(不区分+/-链)/该位置的短读序列总数目。

具体而言，假设当前碱基位置，参考碱基为A，单核苷酸变异类型为C、 G或T，则分别计算A、C、G、T的频率值；此外，计算插入变异(优选地，按照插入碱基的数目不同，分别计算各种插入变异，如插入一个碱基的插入变异、插入两个碱基的插入变异)以及删除变异的频率值，如图2中右上部分所示。

ST203、对上述结果文件进行整理，得到每种肿瘤类型在每个碱基位置上的各变异的频率值分布。

对前一步得到的结果文件，提取每个碱基位置上的各变异的频率值。每个样本数据在每个碱基位置上的每种变异(变异包括非变异类型、单核苷酸变异类型、插入变异类型、删除变异类型)对应一个频率值。将频率值按照样本所属的肿瘤类型、碱基位置和变异类型(这里，非变异也视为一种变异类型)进行归类，这样可以得到每种肿瘤类型在每个碱基位置上的各变异的频率值的分布D(i,p,t),其中i代表肿瘤类型，p代表碱基位置，t代表变异。

ST204、寻找肿瘤类型间分布差异显著的碱基位置，并将涉及的肿瘤类型在该碱基位置上各变异的频率值分布作为提取的生物标记物，如图2中的特征提取部分所示。

前述步骤提供了肿瘤类型i与j在每个碱基位置和每种变异上的分布D(i, p,t)与D(j,p,t)。计算这两个分布的距离ΔD(ij,p,t)和距离的P值P(ij,p,t)，其中p代表碱基位置，t代表变异的类型。

对两分布距离的计算可以采用Kolmogorov-Smirnov检验(KS test)来进行。例如，在R语言里面，KS test的用法如下：

ks.test(D(i,p,t),D(j,p,t))

该函数返回两分布间的距离ΔD(ij,p,t)和距离的P值P(ij,p,t)，其中i,j, p,t的注释同上。

ΔD越大，代表两个分布差异越大；P值越小，代表两个分布的差异越显著。优选选取距离大而P值小的变异。如果某一变异的分布的差异达到一定阈值(例如P<0.05，D>0.95)，则将该变异所在碱基位置的各变异的频率值分布作为生物标记物提取出来。

本步骤的改动调整会影响最后生物标记物的发现。本步骤一些可能的变动方案包括：

1.提取不同肿瘤类型间差异最大的变异所在碱基位的所有变异的频率值分布(当前方案)。

2.提取不同肿瘤类型间差异最大的变异的频率值分布。

第二个方案最为直观，但是从公有多样本数据集中提取的变异的频率值生成的热图来讲，方案2不如方案1效果好。

步骤S204主要考虑两种肿瘤类型间生物标记物的提取。实际应用时，其中一种肿瘤类别如果为“正常(无肿瘤)”，这样得出的生物标记物可以直接反映肿瘤病症的发生发展。若两种肿瘤类别均对应病理上的两种肿瘤，得出的生物标记物可以反应两种肿瘤的差异，这种情况可以适用相似肿瘤类别的区分，如非小细胞肺癌中的腺癌和鳞状细胞癌的区分。

若存在多种肿瘤类别，根据本发明实施例的方法步骤S201-203，可以得到一张数据的二维表，其中每一行为某个碱基位置的某个变异，每一列为一个样本。这里的碱基位置涵盖了靶向区域上面的所有碱基位置。该问题为常见的多类别分类问题。可以使用例如多元逻辑回归(multi-nomial logistic regression)模型配合特征的贡献值来筛选合适的碱基(变异)类型，例如在 R语言里，使用caret库的示例代码如下：

上述代码示范了从不同肿瘤样本数据中筛选生物标记物。trainControl函数设定了模型训练的基本框架，本例使用10倍交叉验证，重复5次。train() 函数用于设定模型的信息，本例中使用glmnet模型用于多元逻辑回归。该模型也可以替换为xgboost,glm等各种分类模型的名称。varImp()提供了每个变异在模型中的贡献。根据实际数据情况，可以保留贡献前20或者40的变异，并将其频率分布用作生物标记物。

本发明实施例的方法与当前ctDNA分析方法的不同体现在以下几点：1) 本发明实施例的方法直接从cfDNA中发现与肿瘤相关的变异信息用作生物标记物，而不是寻找在肿瘤组织中已经被发现的各种变异；2)本发明实施例的方法在流程上不需要体细胞变异识别步骤；3)其他使用cfDNA来追踪肿瘤来源的方案(如利用HPV DNA来鉴定肺癌，利用ctDNA序列特征来追踪肿瘤来源)流程复杂，往往需要大量肿瘤组织样本的验证工作。本发明实施例的方法流程简单，可以直接提取单点突变，对不同肿瘤类型的适用性较强，同时并不需要肿瘤组织DNA序列的验证工作(本发明实施例的方法并不假设发现的变异直接来自肿瘤)。4)当前cfDNA分析流程沿袭肿瘤组织DNA 分析流程，而肿瘤组织DNA数据中的突变等位基因频率(mutant allele frequency,MAF)并不与肿瘤负担关联(受肿瘤取样偏差以及肿瘤异质性影响)。因此，传统肿瘤组织DNA分析中，往往并不把MAF纳入生物标记物的考虑。最近有研究报道特定的肿瘤突变MAF与治疗之间的关联(ref),但是并没有涉及新的突变MAF的发现，也没有涉及对肿瘤类型的区分。

图2示例性地显示了对胰腺癌和肺癌患者样本采用基于肿瘤患者游离 DNA深度测序数据的生物标记物发现方法寻找生物标记物的基本流程。

--搜集胰腺癌和肺癌患者的ctDNA深度测序原始数据；

--将生成的fastq文件使用bwa进行序列对比，然后使用PICARD和GATK 相关工具进行数据清理(碱基重校正和插入删除的重比对)；

--使用samtools mpileup统计每个碱基位的各变异类型的频率；

--统计胰腺癌与肺癌两组样本在每个碱基位的各变异的分布，选出差异显著的分布，并提取出对应的碱基位下的各碱基(变异)类型的频率作为生物标记物(特征)；

对于筛选出来的生物标记物(特征)，可以使用热图或者协同聚类初步展示特征的区分效果。如图2所示，例如，可以利用线性分类器对生物标记物 x进行分类，分类函数y＝f(x^T)的值为0.1时，分类结果是胰腺癌，分类函数y 的值为0.9时，分类结果是肺癌。

图3(a)和图3(b)是示例性示出评估图2所示实施例中提取的诊断性生物标记物对肿瘤状况的区分效果的图。可以使用图3(a)中的热图或者图 3(b)中的协同聚类图初步展示特征的区分效果。热图中，每一行代表筛选出来的某个碱基位置上的某个变异，每一列代表一个样本。热图中的每列用不同颜色深度区分该列对应的肿瘤类型(浅灰色为肺癌，深灰色为胰腺癌)。热图中的某行某列处的栅格对应值(例如灰度值)为某样本在某碱基位置上特定变异类型的频率值。频率值接近的样本间通过层次聚类会被安排到相近的位置。理想情况下，同一形状标记的色块应当各自连成一片。最坏情况下，不同形状的标记混杂在一起。

图3(b)的协同聚类图中每一行和每一列都代表一个样本。图中位于(i, j)的小方块的颜色越深，代表样本i与j越相似。理想情况下，协同聚类图中深色小方块会连在一起形成数个深色的大方块。每个大方块代表一种肿瘤类型或者亚型。最坏情况下，深浅不同的小方块混杂在一起。列上方的相同颜色色块的聚集程度就反映了生物标记物的肿瘤类型区分效果。由图3(b) 表明，可以很好地区分两种肿瘤类型，并且提示可能具有肿瘤亚型。

此外，该实施例中进一步可以包括对所提取的生物标记物的后续使用的步骤，包括：

--训练步骤，即，利用训练数据集对分类模型进行训练，训练数据由从已知肿瘤类型的样本的细胞游离DNA的深度测序数据所提取的特征(也就是生物标记物)和所述肿瘤类型构成；

--(可选地)测试步骤，基于测试数据集，对训练好的分类模型进行性能评估。

--预测步骤，对于未做分类的(待预测的)样本数据，可以重复上述生物标记物的提取步骤，并将提取出的生物标记物作为特征导入分类模型，从而得到分类模型对肿瘤类型的预测值。例如，可以采用各种分类模型，包括线性分类器和非线性分类器。再例如，分类模型可以采用支持向量机、逻辑回归算法、人工神经网络中的任何一种或多种来实现。例如，分类模型的计算结果接近0，则预测样本的类别为胰腺癌，接近1，则预测样本的类别为肺癌。

图4是评估本发明实施例的方法所建立的模型的技术效果的工作特征曲线(ROC)指标的示意图。其中，横坐标为假阳性率，纵坐标为真阳性率。对角线代表随机猜测结果的模型。

评估本发明实施例的方法所建模型的效果的方法如下：

--如果该肿瘤信息为有序变量(比如治疗反应，肿瘤发展期)，则例如计算出这批样本的均方根误差(RMSE)值。RMSE能够很好地反映测量的精度，其计算公式如下：

其中n为样本个数，

为肿瘤状况预测值，y_i为实际值。RMSE越低表明模型预测效果越好。由于RMSE与预测值在同一量级上。举例，对于预测治疗反应的生物标记物，如果预测值为1.8(考虑这样的场合，0代表肿瘤恶化，1代表肿瘤稳定，2代表肿瘤对治疗部分反应，3代表完全反应)，如果 RMSE＝0.1，则实际结果很大可能是2[1.8–0.1,1.8+0.1]。

--如果该肿瘤信息为二值分类变量(比如两种肿瘤类别)，则可以使用工作特征曲线(ROC)图的曲线下面面积(AUC)值来评估分类好坏。AUC 值越接近1，代表模型效果越可靠，而AUC值越接近0.5，代表模型效果越接近随机选择(模型不起作用)。如果该肿瘤信息为多值分类变量(比如多种肿瘤类别)，则可以考虑精准率(precision)，召回率(recall)和F1分值。其公式如下：

其中TP代表真阳性个数(模型预测为真且实际为真)，FP代表假阳性个数(模型预测为真但实际为假)，FN代表假阴性个数(模型预测为假而实际为真)。这些分值越接近1代表模型的性能越好。如图4中虚线所示，可以看到模型的分类效果较好。

下面，描述本发明的装置实施例。如图5所示，提供一种基于不同类别的细胞游离DNA的深度测序数据提取生物标记物的装置500，包括：样本数据获取模块501，其获取不同类别的样本的细胞游离DNA的深度测序数据；频率值计算模块502，其计算每个样本数据在每个碱基位置的各变异的频率值；频率值分布提取模块503，其基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布；分布差异程度计算模块504，其确定频率分布在类别间具有充分差异的变异；以及生物标记物提取模块505，其基于所确定的变异提取所述生物标记物。

具体地，在生物标记物提取模块505中，提取所确定的变异所在碱基位置的所有变异的频率分布，作为生物标记物；或者，可替换地，提取所确定的变异的频率分布，作为生物标记物。

在优选实施例中，分布差异程度可以包括分布的距离，预定条件可以包括：分布的距离大于第一预定阈值。分布差异程度还可以包括P值，预定条件还包括：P值小于第二预定阈值。此外，优选地，碱基位置为靶向区域中的碱基位置。此外，还优选的是，变异包括如下类型：单核苷酸变异类型、插入变异类型、删除变异类型。另外，优选的是，类别为肿瘤类别、关于预后的类别或与关于治疗反应的类别。

上述各个模块501-505实现为装置500上的程序模块。为了便于描述，文中已经定义了系统400的功能构建模块的配置和边界。只要适当地执行指定的功能和关系，就可以定义其他可替代的边界，从而将以上各个程序模块整合或分离以得到新的程序模块。基于文中所包含的教导，替代方案(包括文中所述的那些方案的等同形式、扩展形式、变化形式、偏离形式等在内) 对于相关领域的技术人员来说将是显而易见的。这样的替代方案落入所公开的实施例的范围之内。上述装置500可以包括处理器，所述处理器可通信地连接到存储有可执行计算机指令的存储器，当所述处理器执行存储器上的对应于上述各个程序模块的计算机可执行指令时，则执行各个程序模块的相应功能。

本文中的处理器可以是包括诸如微处理器、中央处理单元(CPU)、图形处理单元(GPU)等的一个或更多通用处理装置的处理器电路。更具体地，处理器可以是复杂指令集计算(CISC)微处理器、精简指令集计算(RISC) 微处理器、超长指令字(VLIW)微处理器、实现其他指令集的处理器、或实现指令集的组合的处理器。

处理器还可以是一个或更多专用处理装置，诸如专用集成电路(ASIC)、现场可编程门阵列(FPGA)、数字信号处理器(DSP)、片上系统(SoC)等。如本领域技术人员将领会的，在一些实施例中，处理器装置可以是专用处理器，而不是通用处理器。处理器可以包括一个或更多已知的处理装置，诸如由Intel^TM制造的Pentium^TM或Xeon^TM系列的微处理器、由AMD^TM制造的 Turion^TM系列的微处理器、或由其他供应商诸如Oracle^TM制造的任何各种处理器(例如，SPARC^TM架构处理器)。处理器还可以包括由Nvidia^TM制造的图形处理单元。所公开的实施例不限于与所公开的实施例相符的以其他方式被配置为满足识别、分析、维护、生成和/或提供大量样本数据或任何其他类型的数据的计算需求的任何类型的(一个或数个)处理器。

本文中的存储器可以包括被配置为存储由处理器使用的计算机可执行指令的一个或更多存储设备，以执行与所公开的实施例相关的功能。例如，存储器可以存储用于操作系统、样本数据获取模块501、频率值计算模块502、频率值分布提取模块503、分布差异程度计算模块504、生物标记物提取模块 505的任何一个的计算机可执行软件指令。

在下面一个实施例中，提供了一种计算设备，包括处理器，处理器配置为运行计算机指令以实现以下步骤：取不同类别的样本的细胞游离DNA的深度测序数据；计算每个样本数据在每个碱基位置的各变异的频率值；基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布；确定频率分布在类别间具有充分差异的变异；基于所确定的变异提取生物标记物。

其中，处理器进一步配置为执行计算机指令以实现以下步骤：利用所提取的生物标记物作为特征，对分类模型进行训练。

如此得到的分类模型可以用于对新的样本数据的类别进行预测。

在计算设备中，处理器可以进一步配置为执行计算机指令以实现以下步骤：将待分类的体液样本的细胞游离DNA深度测序数据输入分类模型中，输出分类结果。

例如，当提取的生物标记物是预后生物标记物时，可以对新的样本数据进行预后方面的预测。具体的预测过程在图2所示实施例中已进行了说明和例示，此处由于篇幅所限，不展开描述。

注意的是，本发明装置实施例的未详尽描述之处，可以参考本发明上述方法实施例的相关说明。

还在另一实施例中，还提供了一种计算机可读存储介质，如光介质、磁介质等，其上存储有计算机程序，计算机程序被处理器执行时，实现以下步骤：获取不同类别的样本的细胞游离DNA的深度测序数据；计算每个样本数据在每个碱基位置的各变异的频率值；基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布；确定频率分布在类别间具有充分差异的变异；基于所确定的变异提取所述生物标记物。

其中计算机程序被处理器执行时，可以进一步执行以下步骤：利用所提取的生物标记物作为特征，对分类模型进行训练。所得到的分类模型可以一并存储在该计算机可读存储介质中，或者也可以存储在计算机的存储单元中 (例如，ROM)，或者也可以以任何其他的形式存储。

此外，所述计算机可读存储介质可以可通信地连接到存储有分类模型的其他装置，例如其他的服务器、数据库、缓存装置等，由此可以将所提取的生物标记物传输到分类模型所在的其他装置，以便在其他装置处对分类模型进行训练。

所得到的分类模型可供对新的样本数据在疾病类型、治疗反应、预后等方面的预测。

本文描述的上述方法示例至少部分可以是机器或计算机实现的。一些示例可以包括编码有指令的计算机可读介质或机器可读介质，该指令可经操作以配置电子装置来执行如上述示例中所描述的方法。这种方法的实现可以包括诸如微代码、汇编语言代码、较高级语言代码等的代码。这种代码可以包括用于执行各种方法的计算机可读指令。这种代码可以形成计算机程序产品的一部分。此外，在一个示例中，例如在运行期间或在其它时间，代码可以有形地存储在一个或更多易失性、非暂时性或非易失性有形计算机可读介质上。这些有形的计算机可读介质的示例可以包括但不限于硬盘、可移动磁盘、可移动光盘(例如，压缩光盘和数字视频盘)、磁带、存储卡或存储棒、随机存取存储器(RAM)、只读存储器(ROM)等。

在实施例中，虽然计算机可读存储介质可以是单一介质，但是术语“计算机可读存储介质”应当被理解为包括单一介质或多介质(例如，集中式或分布式数据库，和/或相关的缓存和服务器)，其存储一个或更多的计算机可执行指令集或数据集。术语“计算机可读存储介质”还应当被理解为包括能够对由计算机执行并且使得计算机执行本公开的方法中的任何一个或更多的指令集进行存储或编码的任何介质。

本发明不局限于上述特定实施例，虽然上述特定实施例以一定次序对各个步骤和构件进行了说明，但这仅仅作为示例，说明的各个步骤和构件对于本发明并不是必需的，次序也并非必要的，只要能够实现所附权利要求书限定的技术方案即可。在不背离本发明精神及其实质情况下，熟悉本领域的技术人员可根据本发明作出各种相应改变和变形，例如对不同实施例中的步骤和构件进行合并、调序、分离和组合，这些相应改变和变形都应属于本发明所附权利要求的保护范围之内。

Claims

1.一种基于细胞游离DNA的深度测序数据提取生物标记物的方法，所述方法由计算机实现，包括由处理器执行的如下步骤：

获取不同肿瘤类别的样本的细胞游离DNA的深度测序数据，其中，所述样本为体液样本；

计算每个样本数据在每个碱基位置的各变异的频率值；

基于所计算的每个碱基位置的各变异的频率值，得到每个类别在碱基位置上的变异的频率分布，具体包括：对于不同类别中的每一类别，对属于该类别的每一体液样本的样本数据在每一碱基位置上各变异的频率值进行复合计算，得到每一类别在每一碱基位置上各变异的频率值分布；

确定频率分布在类别间具有充分差异的变异；以及

基于所确定的变异提取所述生物标记物而无需肿瘤组织DNA序列的验证。

2.根据权利要求1所述的方法，其特征在于，基于所确定的变异提取所述生物标记物的步骤包括：提取所确定的变异所在碱基位置的所有变异的频率分布，作为所述生物标记物。

3.根据权利要求1所述的方法，其特征在于，基于所确定的变异提取所述生物标记物的步骤包括：提取所确定的变异的频率分布，作为所述生物标记物。

4.根据权利要求1所述的方法，其特征在于，确定频率分布在类别间具有充分差异的变异的步骤包括：

计算变异的频率分布在不同类别之间的差异程度；

在所述变异的差异程度满足预定条件的情况下，认为所述变异的频率分布在类别间具有充分差异。

5.根据权利要求4所述的方法，其中：

所述差异程度包括频率分布的距离，所述预定条件包括：所述频率分布的距离大于第一预定阈值。

6.根据权利要求4所述的方法，其中：

所述差异程度包括P值，所述预定条件包括：所述P值小于第二预定阈值。

7.根据权利要求1所述的方法，其中：

所述碱基位置为靶向区域中的碱基位置。

8.根据权利要求1～7中任一项所述的方法，其中：

所述变异包括如下类型中的一种或更多种：单核苷酸变异类型、插入变异类型、删除变异类型。

9.根据权利要求1～7中任一项所述的方法，其中：

所述类别为肿瘤类别、关于预后的类别、与关于治疗反应的类别中的任何一个类别。

10.一种基于细胞游离DNA的深度测序数据提取生物标记物的装置，包括：

样本数据获取模块，其配置为获取不同肿瘤类别的样本的细胞游离DNA的深度测序数据，其中，所述样本为体液样本；

频率值计算模块，其配置为计算每个样本数据在每个碱基位置的各变异的频率值；

频率值分布提取模块，其配置为基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布，具体包括：对于不同类别中的每一类别，对属于该类别的每一体液样本的样本数据在每一碱基位置上各变异的频率值进行复合计算，得到每一类别在每一碱基位置上各变异的频率值分布；

分布差异程度计算模块，其配置为确定频率分布在类别间具有充分差异的变异；以及

生物标记物提取模块，其配置为基于所确定的变异提取所述生物标记物而无需肿瘤组织DNA序列的验证。

11.根据权利要求10所述的装置，其中所述生物标记物提取模块进一步配置为：

提取所确定的变异所在碱基位置的所有变异的频率分布，作为所述生物标记物；或者

提取所确定的变异的频率分布，作为所述生物标记物。

12.一种计算设备，包括处理器，所述处理器配置为运行计算机指令以实现以下步骤：

计算每个样本数据在每个碱基位置的各变异的频率值；

基于所计算的每个碱基位置的各变异的频率值，得到每个类别在所有碱基位置上的变异的频率分布，具体包括：对于不同类别中的每一类别，对属于该类别的每一体液样本的样本数据在每一碱基位置上各变异的频率值进行复合计算，得到每一类别在每一碱基位置上各变异的频率值分布；

确定频率分布在类别间具有充分差异的变异；以及

基于所确定的变异提取生物标记物而无需肿瘤组织DNA序列的验证。

13.根据权利要求12所述的计算设备，所述处理器进一步配置为执行计算机指令以实现以下步骤：

利用所提取的生物标记物作为特征，来训练分类模型。

14.根据权利要求12所述的计算设备，所述处理器进一步配置为执行计算机指令以实现以下步骤：

基于待分类的体液样本的细胞游离DNA深度测序数据，利用训练好的所述分类模型，确定分类结果。

15.一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序使计算机执行以下步骤：

计算每个样本数据在每个碱基位置的各变异的频率值；

确定频率分布在类别间具有充分差异的变异；以及

16.根据权利要求15所述的计算机可读存储介质，其中所述计算机程序使计算机进一步执行以下步骤：

利用所提取的生物标记物作为特征，来训练分类模型。