用于确定DNA样本甲基化水平的方法、设备和存储介质与流程

本公开总体上涉及生物信息处理，并且具体地，涉及确定dna样本的甲基化水平的方法、计算设备和计算机存储介质。

背景技术：

dna甲基化是最早被发现、也是目前研究最深入的表观遗传调控机制之一。所谓的表观遗传是指在基因的dna序列不发生改变的情况下，基因的表达水平与功能发生改变，并可遗传给后代的现象。所谓的dna甲基化，一般是指在dna甲基化转移酶的作用下，在基因组cpg（即，胞嘧啶c-磷酸p-鸟嘌呤g的二核苷酸结构）二核苷酸的c碱基第五位碳原子（胞嘧啶c为一个六原子环状结构，环上含2个n原子和4个c原子，并按固定顺序标记为编号1-6，第五位碳原子即为编号为5的c原子）的位置以共价键的形式结合一个甲基基团的化学修饰过程。

人类基因组含有大约1%的已甲基化的胞嘧啶，因此其是最丰富、最广泛的dna修饰方式。基因区域内的cpg位点通过甲基化的方式影响基因转录活性，从而调控基因表达。在肿瘤细胞中，普遍存在着与正常细胞不同的dna甲基化水平的改变，主要特点是总体甲基化水平的降低与局部甲基化水平的升高。在肿瘤细胞中，原癌基因（即，细胞内与细胞增殖相关的基因）处于低甲基化水平而被激活，抑癌基因（即，肿瘤抑制基因，其是一类存在于正常细胞内可抑制细胞生长并具有潜在抑癌作用的基因）处于高甲基化水平而被抑制，从而导致肿瘤细胞的过度增殖。因此，准确测量和计算dna甲基化水平至关重要。

目前检测样本甲基化水平最常用的建库技术是重亚硫酸盐转化技术。即用重亚硫酸盐处理样本dna，所有未发生甲基化的胞嘧啶c被转化为尿嘧啶u，而已甲基化的胞嘧啶c则保持不变；然后经pcr扩增过后，尿嘧啶u会转换为胸腺嘧啶t，已甲基化的胞嘧啶c依然保持不变，从而将原始dna中未甲基化的胞嘧啶c和已甲基化的胞嘧啶c区分开来，进而统计各cpg位点的胞嘧啶c的甲基化水平。

传统的确定dna样本甲基化水平的方案例如是焦磷酸甲基化测序技术，焦磷酸测序技术是由4种酶（dna聚合酶、atp硫酸化酶、荧光素酶和三磷酸腺苷双磷酸酶）催化的同一反应体系中的酶级联化学发光反应，可准确测量每个位点的甲基化水平，同时可评估一段区域内的整体甲基化水平。焦磷酸甲基化测序技术能够快速地检测cpg位点甲基化的频率，对样品中的甲基化位点进行定性及定量检测。不过，焦磷酸甲基化测序及后期生信处理技术所存在的不足之处例如包括：测试周期长，正常情况下，焦磷酸甲基化确定技术例如需要针对比对结果中参考基因组所有位点逐一统计其碱基数量以用于甲基化水平计算，因此需要耗费较多时间用于数据统计与计算，例如从实验开始需要至少一周的时间才能拿到关于甲基化水平的计算结果。另外，焦磷酸甲基化测序及后期生信处理技术的稳定性不太理想，例如，同一样本、不同公司检测得出的结果会有所不同，有时甚至差异很大，而且同一公司、同一样本在不同时间检测得出的结果有时也会有所差异。此外，焦磷酸甲基化测序技术还存在长度限制和dna用量大的不足之处，例如，待检测序列的长度一般为几十个碱基，最多不超过一百个碱基。若需检测较长序列则需要多次实验或平行实验（即，同一样本同时进行多组实验，每一组的条件均有所差异；这里是指每一组检测的序列位置不同，从而达到同时检测较长序列的目的）。在dna用量方面，一般需至少1ug的dna片段才可用于实验。

综上，在传统的确定dna样本甲基化水平的方案中，存在确定周期长、结果的稳定性欠佳，长度限制和dna用量大的不足之处。

技术实现要素：

本公开提供一种确定dna样本的甲基化水平的方法、计算设备和计算机存储介质，能够快速准确计算得到dna样本的甲基化水平。

根据本公开的第一方面，提供了一种确定dna样本的甲基化水平的方法。该方法包括：过滤所接收的关于dna样本的测序数据，以便留下符合预定条件的测序数据；将所留下的测序数据分别与dna样本相对应的参考基因组和质控序列的基因组进行比对，以便生成比对到参考基因组的正链和负链的第一比对结果数据，以及比对到质控序列的基因组的正链和负链的第二比对结果数据；提取第一比对结果数据中预定标识符合第一预定阈值的预定位点的reads，以便生成第一提取信息；提取第一比对结果数据中预定标识符合第二预定阈值的预定位点的reads，以便生成第二提取信息；基于参考基因组的正链和负链的各预定位点，分别统计第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，预定位点属于预定的位点集合；基于第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，计算关于参考基因组的各预定位点的甲基化水平；以及基于第二比对结果数据，确定关于dna样本的甲基化的转化效率和错误率中的至少一个，以用于确定dna样本的甲基化水平。

根据本发明的第二方面，还提供了一种计算设备，该设备包括：存储器，被配置为存储一个或多个计算机程序；以及处理器，耦合至存储器并且被配置为执行一个或多个程序使装置执行本公开的第一方面的方法。

根据本公开的第三方面，还提供了一种非瞬态计算机可读存储介质。该非瞬态计算机可读存储介质上存储有机器可执行指令，该机器可执行指令在被执行时使机器执行本公开的第一方面的方法。

在一些实施例中，确定dna样本的甲基化水平的方法还包括：基于第二比对结果数据，统计质控序列的基因组的正链和负链的各预定位点的a、c、g、t四种碱基的数量，以便计算与质控序列的基因组相关联的各预定位点的甲基化水平。

在一些实施例中，统计质控序列的基因组的正链和负链的各预定位点的a、c、g、t四种碱基的数量包括：提取第二比对结果数据中预定标识符合第一预定阈值的预定位点的reads，以便生成第三提取信息；提取第二比对结果数据中预定标识符合第二预定阈值的预定位点的reads，以便生成第四提取信息；基于质控序列的基因组的正链和负链的各预定位点，分别统计第三提取信息中的a、c、g、t四种碱基的数量和第四提取信息中的a、c、g、t四种碱基的数量，以用于计算与质控序列的基因组相关联的各预定位点的甲基化水平。

在一些实施例中，分别统计第一提取信息中的碱基的数量和第二提取信息中的碱基的数量包括：统计第一提取信息中的a、c、g、t四种碱基的数量；以及统计第二提取信息中的a、c、g、t四种碱基的数量。

在一些实施例中，基于第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，计算关于参考基因组的各预定位点的甲基化水平包括：计算第二提取信息中的对应位点的c、t两种碱基的总数量；计算第一提取信息中的对应位点的g、a两种碱基的总数量；基于第二提取信息中的对应位点的c碱基的数量，第一提取信息中的对应位点的g碱基的数量，第二提取信息中的对应位点的c、t两种碱基的总数量，以及第一提取信息中的对应位点的g、a两种碱基的总数量，计算对应位点的甲基化水平，以便获得关于参考基因组的各预定位点的甲基化水平。

在一些实施例中，计算与质控序列的基因组相关联的各位点的甲基化水平包括：计算第四提取信息中的对应位点的c、t两种碱基的总数量；计算第三提取信息中的对应位点的g、a两种碱基的总数量；基于第四提取信息中的对应位点的c碱基的数量，第三提取信息中的对应位点的g碱基的数量，第四提取信息中的对应位点的c、t两种碱基的总数量，以及第三提取信息中的对应位点的g、a两种碱基的总数量，计算与质控序列的基因组相关联的各预定位点的甲基化水平。

在一些实施例中，第二比对结果数据确定关于dna样本的甲基化的转化效率和错误率中的至少一个包括：计算第四提取信息中所有未甲基化位点的c碱基的数量之和；计算第三提取信息中所有未甲基化位点的g碱基的数量之和；计算第四提取信息中所有未甲基化位点的c、t两种碱基的数量之和；计算第三提取信息中所有未甲基化位点的g、a两种碱基的数量之和；基于未甲基化位点的c碱基的数量之和，未甲基化位点的g碱基的数量之和，未甲基化位点的g、a两种碱基的数量之和，以及未甲基化位点的g、a两种碱基的数量之和，计算关于dna样本的甲基化的转化效率。

在一些实施例中，第二比对结果数据确定关于dna样本的甲基化的转化效率和错误率中的至少一个包括：计算第四提取信息中所有甲基化位点的c碱基的数量之和；计算第三提取信息中所有甲基化位点的g碱基的数量之和；计算第四提取信息中所有甲基化位点的c、t两种碱基的数量之和；计算第三提取信息中所有甲基化位点的g、a两种碱基的数量之和；基于甲基化位点的c碱基的数量之和，甲基化位点的g碱基的数量之和，甲基化位点的g、a两种碱基的数量之和，以及甲基化位点的g、a两种碱基的数量之和，计算关于dna样本的甲基化的错误率。

在一些实施例中，过滤测序数据以便留下符合预定条件的测序数据包括：在测序数据中将包含接头序列的测序数据移除；过滤掉测序质量值低于预定质量阈值的测序数据；以及过滤测序序列长度低于预定测序序列长度值的测序数据，以便将留下的测序数据作为符合预定条件的测序数据。

提供发明内容部分是为了以简化的形式来介绍对概念的选择，它们在下文的具体实施方式中将被进一步描述。发明内容部分无意标识本公开的关键特征或主要特征，也无意限制本公开的范围。

附图说明

图1示出了根据本公开的实施例的用于确定dna样本的甲基化水平的方法的系统的示意图。

图2示出了根据本公开的实施例的用于确定dna样本的甲基化水平的方法的流程图。

图3示出了根据本公开的实施例的用于计算与质控序列的基因组相关联的各预定位点的甲基化水平的方法的流程图。

图4示出了根据本公开的实施例的用于确定关于dna样本的甲基化的错误率的方法的流程图。

图5示出了根据本公开实施例的用于呈现指示各预定位点甲基化水平分布情况的图像的示意图。

图6示意性示出了适于用来实现本公开实施例的电子设备的框图。

在各个附图中，相同或对应的标号表示相同或对应的部分。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施例。虽然附图中显示了本公开的优选实施例，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

在本文中使用的术语“包括”及其变形表示开放性包括，即“包括但不限于”。除非特别申明，术语“或”表示“和/或”。术语“基于”表示“至少部分地基于”。术语“一个示例实施例”和“一个实施例”表示“至少一个示例实施例”。术语“另一实施例”表示“至少一个另外的实施例”。术语“第一”、“第二”等等可以指代不同的或相同的对象。

在本文中使用的术语“甲基化水平”表示甲基化的胞嘧啶c占所有胞嘧啶c的比例。术语“转化效率”表示未甲基化的胞嘧啶c经重亚硫酸盐处理后转化为尿嘧啶u的碱基比例。术语“错误率”表示甲基化的胞嘧啶c经重亚硫酸盐处理后转化为尿嘧啶u的碱基比例，正常情况不应该出现这种情况的。术语“质控序列”表示在甲基化建库过程中添加到dna文库中的一段外源dna片段，这段dna序列包含多个已知的完全甲基化的cpg位点和完全未甲基化的胞嘧啶c位点。根据未甲基化的胞嘧啶c的转化情况，可评估样本整体的转化效率。根据甲基化的胞嘧啶c的转化情况，可评估样本整体的错误率。常用的质控序列有lambdadna序列和puc19c序列等。术语“reads”代表测序数据的每一条记录称为一条read，是测序结果的通用称呼。另外，鸟嘌呤、胞嘧啶、腺嘌呤、胸腺嘧啶均为单核甘酸类型。

如前文提及，在传统的确定dna样本甲基化水平的方案中，例如需要针对比对结果中参考基因组所有位点逐一统计其碱基数量以用于甲基化水平计算，因此需要耗费较多数据统计与计算时间。因此传统的确定dna样本甲基化水平的方案存在确定周期长、结果的稳定性欠佳等不足之处。

为了至少部分地解决上述问题以及其他潜在问题中的一个或者多个，本公开的示例实施例提出了一种用于确定dna样本的甲基化水平的方案。该方案包括：过滤所接收的关于dna样本的测序数据，以便留下符合预定条件的测序数据；将所留下的测序数据分别与dna样本相对应的参考基因组和质控序列的基因组进行比对，以便生成比对到参考基因组的正链和负链的第一比对结果数据，以及比对到质控序列的基因组的正链和负链的第二比对结果数据；提取第一比对结果数据中预定标识符合第一预定阈值的预定位点的reads，以便生成第一提取信息；提取第一比对结果数据中预定标识符合第二预定阈值的预定位点的reads，以便生成第二提取信息；基于参考基因组的正链和负链的各预定位点，分别统计第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，预定位点属于预定的位点集合；基于第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，计算关于参考基因组的各预定位点的甲基化水平；以及基于第二比对结果数据，确定关于dna样本的甲基化的转化效率和错误率中的至少一个，以用于确定dna样本的甲基化水平。

在上述方案中，通过将经过滤符合预定条件的测序数据分别比对到dna样本对应的参考基因组和质控序列基因组，以便生成第一比对结果数据和第二比对结果数据，以及基于所提取的、“第一比对结果数据中预定标识符合第一预定阈值的预定位点的reads”的碱基统计数量，以及所提取的、“第二比对结果数据中预定标识符合第二预定阈值的预定位点的reads”的碱基统计数量来计算关于质控序列的转化效率和错误率，以用于确定dna样本的甲基化水平；本公开无需针对参考基因组中的所有位点逐一进行碱基统计和甲基化水平计算，仅需针对预定的位点集合中的预定位点进行reads提取、碱基统计和甲基化水平计算，因此，能够有效缩短确定dna样本的甲基化水平的周期。另外，通过基于第二比对结果数据的碱基统计所计算的转化效率或错误率来确认dna样本的甲基化计算结果，能够提高dna样本的甲基化计算结果的稳定性和准确性。再者，本公开通过提取预定标识符合第一预定阈值、第二预定阈值的预定位点的reads，能够提取正链和负链的对比结果均为准确状态下的预定位点的reads，有利于提高甲基化水平计算的准确性。因此，本公开能够快速准确计算得到dna样本中各位点的甲基化水平。

图1示出了根据本公开的实施例的用于确定dna样本的甲基化水平的方法的系统100的示意图。如图1所示，系统100例如包括计算设备110、生信服务器150、测序单元160和网络170。计算设备110可以通过网络170以有线或者无线的方式与生信服务器150和测序单元160进行数据交互。

计算设备110用于计算关于参考基因组的各预定位点的甲基化水平，以及确定关于dna样本的甲基化的转化效率和错误率中的至少一个，以用于确定dna样本的甲基化水平。在一些实施例中，计算设备110可以具有一个或多个处理单元，包括诸如gpu、fpga和asic等的专用处理单元以及诸如cpu的通用处理单元。计算设备110例如是将基于linux操作系统，其上配置有shell编程语言。在一些实施例中，计算设备110例如而不限于配置有bismark软件、bsmap软件、methylextract软件。另外，在每个计算设备上也可以运行着一个或多个虚拟机。计算设备110例如包括数据接收单元112、过滤单元114、参考基因组单元120和质控序列基因组单元130。上述数据接收单元112、过滤单元114、参考基因组单元120和质控序列基因组单元130可以配置在一个或者多个计算设备110上。

关于参考基因组单元120，其用于评估关于参考基因组的各预定位点的甲基化水平。参考基因组单元120例如包括：参考基因组比对单元116、第一提取信息生成单元118、第二提取信息生成单元122、参考基因组碱基统计单元124和参考基因组甲基化水平计算单元126。

关于质控序列基因组单元130，其用于评估关于dna样本的甲基化的转化效率和错误率。质控序列基因组单元130例如包括：质控序列基因组比对单元132、第三提取信息生成单元134、第四提取信息生成单元136、质控序列基因组碱基统计单元138、质控序列基因组甲基化水平计算单元140和转化效率和错误率确定单元142。

关于数据接收单元112，其用于接收来自测序单元160（例如测序仪）或者生信服务器150的关于dna样本的测序数据。

关于过滤单元114，其用于过滤数据接收单元112所接收的关于dna样本的测序数据，以便留下符合预定条件的测序数据。

关于参考基因组比对单元116，其用于将经由过滤单元114过滤而留下的测序数据与dna样本相对应的参考基因组进行比对，以便生成比对到参考基因组的正链和负链的第一比对结果数据。

关于第一提取信息生成单元118，其用于提取第一比对结果数据中预定标识符合第一预定阈值的预定位点的reads，以便生成第一提取信息。

关于第二提取信息生成单元122，其用于提取第一比对结果数据中预定标识符合第二预定阈值的预定位点的reads，以便生成第二提取信息。

关于参考基因组碱基统计单元124，其用于基于参考基因组的正链和负链的各预定位点，分别统计第一提取信息中的碱基的数量和第二提取信息中的碱基的数量。

关于参考基因组甲基化水平计算单元126，其用于基于第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，计算关于参考基因组的各预定位点的甲基化水平。

关于质控序列基因组比对单元132，其用于将所留下的测序数据与质控序列的基因组进行比对，以便生成比对到质控序列的基因组的正链和负链的第二比对结果数据。

关于第三提取信息生成单元134，其用于提取第二比对结果数据中预定标识符合第一预定阈值的预定位点的reads，以便生成第三提取信息。

关于第四提取信息生成单元136，其用于提取第二比对结果数据中预定标识符合第二预定阈值的预定位点的reads，以便生成第四提取信息。

关于质控序列基因组碱基统计单元138，其用于基于质控序列的基因组的正链和负链的各预定位点，分别统计第三提取信息中的a、c、g、t四种碱基的数量和第四提取信息中的a、c、g、t四种碱基的数量。

关于质控序列基因组甲基化水平计算单元140，其用于计算与质控序列的基因组相关联的各预定位点的甲基化水平。例如，质控序列基因组甲基化水平计算单元140计算第四提取信息中的对应位点的c、t两种碱基的总数量；计算第三提取信息中的对应位点的g、a两种碱基的总数量；以及基于第四提取信息中的对应位点的c碱基的数量，第三提取信息中的对应位点的g碱基的数量，第四提取信息中的对应位点的c、t两种碱基的总数量，以及第三提取信息中的对应位点的g、a两种碱基的总数量，计算与质控序列的基因组相关联的各预定位点的甲基化水平。

关于转化效率和错误率确定单元142，其用于基于第二比对结果数据，确定关于dna样本的甲基化的转化效率和错误率中的至少一个，以用于确定dna样本的甲基化水平。例如，转化效率和错误率确定单元142计算第四提取信息中所有未甲基化位点的c碱基的数量之和；计算第三提取信息中所有未甲基化位点的g碱基的数量之和；计算第四提取信息中所有未甲基化位点的c、t两种碱基的数量之和；计算第三提取信息中所有未甲基化位点的g、a两种碱基的数量之和；以及基于未甲基化位点的c碱基的数量之和，未甲基化位点的g碱基的数量之和，未甲基化位点的g、a两种碱基的数量之和，以及未甲基化位点的g、a两种碱基的数量之和，计算关于dna样本的甲基化的转化效率。

再例如，转化效率和错误率确定单元142计算第四提取信息中所有甲基化位点的c碱基的数量之和；计算第三提取信息中所有甲基化位点的g碱基的数量之和；计算第四提取信息中所有甲基化位点的c、t两种碱基的数量之和；计算第三提取信息中所有甲基化位点的g、a两种碱基的数量之和；以及基于甲基化位点的c碱基的数量之和，甲基化位点的g碱基的数量之和，甲基化位点的g、a两种碱基的数量之和，以及甲基化位点的g、a两种碱基的数量之和，计算关于dna样本的甲基化的错误率。

关于生信服务器150，其用于管理生物信息数据库，提供生物信息。计算设备110可以经由网络170和生信服务器150获得关于dna样本的测序数据。在一些实施例中，生信服务器150可以具有一个或多个处理单元，包括诸如gpu、fpga和asic等的专用处理单元以及诸如cpu的通用处理单元。另外，在每个计算设备上也可以运行着一个或多个虚拟机。

关于测序单元160，其用于生成关于dna样本的测序数据。测序单元160例如是基于ngs的微生物测序技术，例如，在进行dna的提取后，利用高通量测序仪建库测序，即获得原始reads数据，存储格式通常为fastq文件，该fastq文件例如包含碱基质量信息。测序单元160所检测的关于dna样本的测序数据可以经由网络170发送至计算设备110，以便计算设备110基于关于dna样本的测序数据计算关于参考基因组的各预定位点的甲基化水平，以及确定关于dna样本的甲基化的转化效率和错误率。

以下将结合图2描述根据本公开的实施例的用于确定dna样本的甲基化水平的方法。图2示出了根据本公开的实施例的用于确定dna样本的甲基化水平的方法200的流程图。应当理解，方法200例如可以在图6所描述的电子设备600处执行。也可以在图1所描述的计算设备110处执行。应当理解，方法200还可以包括未示出的附加动作和/或可以省略所示出的动作，本公开的范围在此方面不受限制。

在步骤202处，计算设备110可以过滤所接收的关于dna样本的测序数据，以便留下符合预定条件的测序数据。例如，计算设备110获得原始的dna样本的测序数据后，使用adapterremoval软件进行数据质控，以过滤掉不符合要求及低质量的测序数据。

所接收的关于dna样本的测序数据例如是经由测序仪所获得的针对dna样本的测序结果数据。

在一些实施例中，计算设备110可以在测序数据中将包含接头序列的测序数据移除；过滤掉测序质量值低于预定质量阈值的测序数据；以及过滤测序序列长度低于预定测序序列长度值的测序数据，以便将留下的测序数据作为符合预定条件的测序数据。

关于移除包含接头序列的测序数据的方式例如包括：计算设备110可以基于相关参数（例如：--adapter1agatcggaagagc--adapter2agatcggaagagc），在测序数据中将包含接头序列的测序数据移除。

关于测序质量的过滤方式例如包括：计算设备110可以使用“--trimqualities–minquality30”，移除测序质量值低于30的测序数据。

关于测序序列长度的过滤方式例如包括：计算设备110首先确定预定测序序列长度值为150；然后计算设备110可以使用“--minlength150”，移除测序序列长度低于150的测序数据。在一些实施例中，可以将预定测序序列长度设定为软件的默认值。

在步骤204处，计算设备110将所留下的测序数据分别与dna样本相对应的参考基因组和质控序列的基因组进行比对，以便生成比对到参考基因组的正链和负链的第一比对结果数据，以及比对到质控序列的基因组的正链和负链的第二比对结果数据。例如，第一比对结果数据例如是比对结果bam文件。

关于质控序列，其是在甲基化建库过程中添加到dna文库中的一段外源dna片段，这段dna序列包含多个已知的完全甲基化的cpg位点和完全未甲基化的胞嘧啶c位点。根据未甲基化的胞嘧啶c的转化情况，可评估样本整体的转化效率；以及根据甲基化的胞嘧啶c的转化情况，可评估样本整体的错误率。在甲基化实验时，质控序列就被加入至样本中，可以通过该预先合成的质控序列中未甲基化的胞嘧啶c的转化情况。来检测前序的甲基化实验质量及其测序数据是否符合预定要求。

关于比对的方式，例如，计算设备110使用软件bismark的数据比对模块，进行测序数据与dna样本相对应的参考基因组的比对，以及测序数据与和质控序列的基因组的比对。比对的参数设置方式例如包括：关于比对算法设置为：“--path_to_bowtie2”（即，采用bowtie2的比对算法）、以及关于测序方向设置为：“--non_directional”。关于测序方向的上述设置，主要是因为所采用的pe150(即双端测序paired-end150测序，序列长度为150个碱基)的测序策略是双端测序，无方向性，因此需要针对测序方向进行如上参数设置。其它的软件参数例如可以均设置为默认值。

在步骤206处，计算设备110提取第一比对结果数据中预定标识符合第一预定阈值的预定位点的reads，以便生成第一提取信息。

关于预定位点，其例如是预先确定的药物相关的甲基化位点，即具有临床意义或者药物使用具有意义的一些位点。可以通过公开信息或者生信数据库收集关于甲基化具有临床意义或者药物使用意义的所有位点，并将这些位点组成预定的位点集合。该预定的位点集合中的各位点即各预定位点。本公开通过提取第一比对结果数据中预定标识符合第一预定阈值的预定位点的reads，而非所有参考基因组的所有位点的reads，可以在保证甲基化数据的计算准确性的前提下大大提高数据处理的效率。

关于第一预定阈值，其例如为83或163。该第一预定阈值用于指示与参考基因组的正链和负链的中的一个链的对比结果的准确状态。

例如，计算设备110提取数据比对结果bam文件中比对到与dna样本相对应的参考基因组的正链和负链的测序数据；然后将比对结果bam文件中flag标识（flag标识代表reads比对情况的数字标识，flag标识通常位于比对结果文件bam文件的第二列）为83或163的reads提取到第一提取信息中，该第一提取信息例如为fc.bam文件。

在步骤208处，计算设备110提取第一比对结果数据中预定标识符合第二预定阈值的预定位点的reads，以便生成第二提取信息。

关于第二预定阈值，其例如为99或147。该99或147用于指示与参考基因组的正链和负链的中的另一个链的对比结果的准确状态。本公开通过提取预定标识符合第一预定阈值、第二预定阈值的预定位点的reads，能够提取参考基因组正链和负链的对比结果均为准确状态下的预定位点的reads，以便提高后续甲基化水平计算的准确性。

例如，计算设备110提取数据比对结果bam文件中比对到与dna样本相对应的参考基因组的正链和负链的测序数据；然后将比对结果bam文件中flag标识为99或147的reads提取到第二提取信息中，该第二提取信息例如为fw.bam文件。上述第二预定阈值例如为99或147。

在步骤210处，计算设备110基于参考基因组的正链和负链的各预定位点，分别统计第一提取信息中的碱基的数量和第二提取信息中的碱基的数量。在一些实施例中，计算设备110可以针对预定位点中的所有位点分别统计第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，也可以针对预定位点中的部分位点分别统计第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，以便进一步提高数据处理的效率。

关于统计第一提取信息中的碱基的数量和第二提取信息中的碱基的数量的方法例如包括：统计第一提取信息中的a、c、g、t四种碱基的数量；以及统计第二提取信息中的a、c、g、t四种碱基的数量。例如，计算设备110对步骤208所提取出的两个文件，例如是fc.bam文件和fw.bam文件，使用igvtools软件分别统计各位点的a、c、g、t四种碱基的数量。igvtools软件的相关参数为：count，使用计数统计功能；-w1，按照单个位点逐一进行统计。igvtools软件的其它参数可以均设定为默认值。

在步骤212处，计算设备110基于第一提取信息中的碱基的数量和第二提取信息中的碱基的数量，计算关于参考基因组的各预定位点的甲基化水平。

计算关于参考基因组的各预定位点的甲基化水平的方法例如包括：计算设备110计算第二提取信息中的对应位点的c、t两种碱基的总数量；计算第一提取信息中的对应位点的g、a两种碱基的总数量；以及基于第二提取信息中的对应位点的c碱基的数量，第一提取信息中的对应位点的g碱基的数量，第二提取信息中的对应位点的c、t两种碱基的总数量，以及第一提取信息中的对应位点的g、a两种碱基的总数量，计算对应位点的甲基化水平，以便获得关于参考基因组的各预定位点的甲基化水平。以下结合公式（1）来说明计算关于参考基因组的各预定位点的甲基化水平的方式。

（1）

在上述公式（1）中，代表第二提取信息（例如是fw.bam文件）中的对应位点的c碱基的数量。代表第一提取信息（例如是fc.bam文件）中的对应位点的g碱基的数量。代表第二提取信息（例如是fw.bam文件）中的对应位点的c、t两种碱基的总数量。代表第一提取信息（例如是fc.bam文件）中的对应位点的g、a两种碱基的总数量。代表关于参考基因组的各预定位点的甲基化水平。上述对应位点为预定位点中的各个位点。

在步骤214处，计算设备110基于第二比对结果数据，确定关于dna样本的甲基化的转化效率和错误率中的至少一个，以用于确定dna样本的甲基化水平。例如，计算设备110如果确定关于dna样本的甲基化的转化效率达到99%，则认为甲基化实验的质量符合要求的，所计算的dna样本的甲基化水平为准确的。

确定关于dna样本的甲基化的转化效率的方法例如包括：计算设备110计算第四提取信息中所有未甲基化位点的c碱基的数量之和；计算第三提取信息中所有未甲基化位点的g碱基的数量之和；计算第四提取信息中所有未甲基化位点的c、t两种碱基的数量之和；计算第三提取信息中所有未甲基化位点的g、a两种碱基的数量之和；基于未甲基化位点的c碱基的数量之和，未甲基化位点的g碱基的数量之和，未甲基化位点的g、a两种碱基的数量之和，以及未甲基化位点的g、a两种碱基的数量之和，计算关于dna样本的甲基化的转化效率。以下结合公式（2）来说明计算关于dna样本的甲基化的转化效率的方式。

（2）

在上述公式（2）中，代表第四提取信息（关于质控序列fw.bam文件）中所有（例如n个）未甲基化位点的c碱基的数量之和。代表第三提取信息（关于质控序列fc.bam文件）中所有（例如n个）未甲基化位点的g碱基的数量之和。代表第四提取信息（关于质控序列fw.bam文件）中所有（例如n个）未甲基化位点的c、t两种碱基的数量之和。代表第三提取信息（关于质控序列fc.bam文件）中所有（例如n个）未甲基化位点的g、a两种碱基的数量之和。n代表质控序列所有未甲基化位点的数量。代表关于dna样本的甲基化的转化效率，即未甲基化的胞嘧啶c经重亚硫酸盐处理后转化为尿嘧啶u的碱基比例。

关于确定关于dna样本的甲基化的错误率的方法，下文将结合图4加以说明，在此，不再赘述。

在上述方案中，通过将经过滤符合预定条件的测序数据分别比对到dna样本对应的参考基因组和质控序列基因组，以便生成第一比对结果数据和第二比对结果数据，以及基于所提取的、“第一比对结果数据中预定标识符合第一预定阈值的预定位点的reads”的碱基统计数量，以及所提取的、“第二比对结果数据中预定标识符合第二预定阈值的预定位点的reads”的碱基统计数量来计算关于质控序列的转化效率和错误率，以用于确定dna样本的甲基化水平；本公开无需针对参考基因组中的所有位点逐一进行碱基统计和甲基化水平计算，仅需针对预定的位点集合中的预定位点进行reads提取、碱基统计和甲基化水平计算，因此，能够有效缩短确定dna样本的甲基化水平的周期。另外，通过基于第二比对结果数据的碱基统计所计算的转化效率或错误率来确认dna样本的甲基化计算结果，能够提高dna样本的甲基化计算结果的稳定性和准确性。再者，本公开通过提取预定标识符合第一预定阈值、第二预定阈值的预定位点的reads，能够提取正链和负链的对比结果均为准确状态下的预定位点的reads，有利于提高甲基化水平计算的准确性。因此，本公开能够快速准确计算得到dna样本中各位点的甲基化水平。

另外，实验数据表明：从获取样本组织到最终得到dna样本的甲基化计算结果，平均三天即可完成。而且本公开的确定dna样本的甲基化水平的方法稳定性好、可重复性高，在多个样本的生物学重复实验中，各位点的甲基化水平基本一致，差异不超过1%。在同一样本不同时间的重复试验中，得到的结果的一致性也较高。例外，本公开的方法不存在长度限制，最长可对整个基因组进行测序分析，也可根据实际需求针对某一段序列进行检测，灵活性高。再者，本公开对dna使用量的要求较小，50ng的dna片段即可满足实验要求，因此，对组织取样困难或样本十分珍贵的情况极为适用。另外，本公开的方法通过加入质控序列的方式评估转化效率和错误率，进一步保证了dna样本的甲基化计算结果的准确性。

在一些实施例中，方法200还包括用于呈现指示各预定位点甲基化水平分布情况的图像的方法。例如，计算设备110基于各预定位点在基因组上的位置信息和所计算的各预定位点的甲基化水平，绘制用于指示各位点甲基化水平分布情况的图像，以用于呈现图像。例如，图5示出了根据本公开实施例的用于呈现指示各预定位点甲基化水平分布情况的图像500。在图5中，每一个峰代表一个预定位点。图5的横坐标代表该预定位点在基因组上的位置，纵坐标代表该位点所计算的甲基化水平，峰顶端的数值为该预定位点对应的纵坐标具体数值。

在一些实施例中，方法200还包括用于计算与质控序列的基因组相关联的各位点的甲基化水平的方法300。例如，计算设备110基于第二比对结果数据，统计质控序列的基因组的正链和负链的各预定位点的a、c、g、t四种碱基的数量，以便计算与质控序列的基因组相关联的各位点的甲基化水平。下文将结合图3具体说明用于计算与质控序列的基因组相关联的各位点的甲基化水平的方法300。

图3示出了根据本公开的实施例的用于计算与质控序列的基因组相关联的各预定位点的甲基化水平的方法300的流程图。应当理解，方法300例如可以在图6所描述的电子设备600处执行。也可以在图1所描述的计算设备110处执行。应当理解，方法300还可以包括未示出的附加动作和/或可以省略所示出的动作，本公开的范围在此方面不受限制。

在步骤302处，计算设备110提取第二比对结果数据中预定标识符合第一预定阈值的预定位点的reads，以便生成第三提取信息。关于第一预定阈值，其例如为83或163。该第一预定阈值用于指示与质控序列的基因组的正链和负链的中的一个链的对比结果的准确状态。在步骤304处，计算设备110提取第二比对结果数据中预定标识符合第二预定阈值的预定位点的reads，以便生成第四提取信息。关于第二预定阈值，其例如为99或147。该99或147用于指示与质控序列的基因组的正链和负链的中的另一个链的对比结果的准确状态。

在步骤306处，计算设备110基于质控序列的基因组的正链和负链的各预定位点，分别统计第三提取信息中的a、c、g、t四种碱基的数量和第四提取信息中的a、c、g、t四种碱基的数量。

上述第三提取信息中的a、c、g、t四种碱基的数量和第四提取信息中的a、c、g、t四种碱基的数量被用于计算与质控序列的基因组相关联的各预定位点的甲基化水平。具体的计算方式例如包括以下的步骤308至步骤312。

在步骤308处，计算设备110计算第四提取信息中的对应位点的c、t两种碱基的总数量。

在步骤310处，计算设备110计算第三提取信息中的对应位点的g、a两种碱基的总数量。

在步骤312处，计算设备110基于第四提取信息中的对应位点的c碱基的数量，第三提取信息中的对应位点的g碱基的数量，第四提取信息中的对应位点的c、t两种碱基的总数量，以及第三提取信息中的对应位点的g、a两种碱基的总数量，计算与质控序列的基因组相关联的各预定位点的甲基化水平。

在一些实施例中，计算设备110可以基于上述公式（1）来计算与质控序列的基因组相关联的各位点的甲基化水平。

图4示出了根据本公开的实施例的用于确定关于dna样本的甲基化的错误率的方法400的流程图。应当理解，方法400例如可以在图6所描述的电子设备600处执行。也可以在图1所描述的计算设备110处执行。应当理解，方法400还可以包括未示出的附加动作和/或可以省略所示出的动作，本公开的范围在此方面不受限制。

在步骤402处，计算设备110计算第四提取信息中所有甲基化位点的c碱基的数量之和。

在步骤404处，计算设备110计算第三提取信息中所有甲基化位点的g碱基的数量之和。

在步骤406处，计算设备110计算第四提取信息中所有甲基化位点的c、t两种碱基的数量之和。

在步骤408处，计算设备110计算第三提取信息中所有甲基化位点的g、a两种碱基的数量之和。

在步骤410处，计算设备110基于甲基化位点的c碱基的数量之和，甲基化位点的g碱基的数量之和，甲基化位点的g、a两种碱基的数量之和，以及甲基化位点的g、a两种碱基的数量之和，计算关于dna样本的甲基化的错误率。以下结合公式（3）来说明计算关于dna样本的甲基化的错误率的方式。

（3）

在上述公式（3）中，代表第四提取信息（例如关于质控序列的fw.bam文件）中所有（例如m个）甲基化位点的c碱基的数量之和。代表第三提取信息（例如关于质控序列的fc.bam文件）中所有（例如m个）甲基化位点的g碱基的数量之和。代表第四提取信息（例如关于质控序列的fw.bam文件）中所有（例如m个）甲基化位点的c、t两种碱基的数量之和。代表第三提取信息例如关于质控序列的fc.bam文件）中所有（例如m个）甲基化位点的g、a两种碱基的数量之和。m代表质控序列所有甲基化位点的数量。代表关于dna样本的甲基化的错误率，即甲基化的胞嘧啶c经重亚硫酸盐处理后转化为尿嘧啶u的碱基比例。

图6示意性示出了适于用来实现本公开实施例的电子设备600的框图。设备600可以是用于实现执行图2至4所示的方法200至400的设备。如图6所示，设备600包括中央处理单元（cpu）601，其可以根据存储在只读存储器（rom）602中的计算机程序指令或者从存储单元608加载到随机存取存储器（ram）603中的计算机程序指令，来执行各种适当的动作和处理。在ram中，还可存储设备600操作所需的各种程序和数据。cpu、rom以及ram通过总线604彼此相连。输入/输出（i/o）接口605也连接至总线604。

设备600中的多个部件连接至i/o接口605，包括：输入单元606、输出单元607、存储单元608，中央处理单元601执行上文所描述的各个方法和处理，例如执行方法200至400。例如，在一些实施例中，方法200至400可被实现为计算机软件程序，其被存储于机器可读介质，例如存储单元608。在一些实施例中，计算机程序的部分或者全部可以经由rom和/或通信单元609而被载入和/或安装到设备600上。当计算机程序加载到ram并由cpu执行时，可以执行上文描述的方法200至400的一个或多个操作。备选地，在其他实施例中，cpu可以通过其他任何适当的方式（例如，借助于固件）而被配置为执行方法200至400的一个或多个动作。

需要进一步说明的是，本公开可以是方法、装置、系统和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于执行本公开的各个方面的计算机可读程序指令。

计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子（非穷举的列表）包括：便携式计算机盘、硬盘、随机存取存储器（ram）、只读存储器（rom）、可擦式可编程只读存储器（eprom或闪存）、静态随机存取存储器（sram）、便携式压缩盘只读存储器（cd-rom）、数字多功能盘（dvd）、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波（例如，通过光纤电缆的光脉冲）、或者通过电线传输的电信号。

这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令，并转发该计算机可读程序指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。

用于执行本公开操作的计算机程序指令可以是汇编指令、指令集架构（isa）指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，该编程语言包括面向对象的编程语言—诸如smalltalk、c++等，以及常规的过程式编程语言—诸如“c”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络—包括局域网(lan)或广域网(wan)—连接到用户计算机，或者，可以连接到外部计算机（例如利用因特网服务提供商来通过因特网连接）。在一些实施例中，通过利用计算机可读程序指令的状态信息来个性化定制电子电路，例如可编程逻辑电路、现场可编程门阵列（fpga）或可编程逻辑阵列（pla），该电子电路可以执行计算机可读程序指令，从而实现本公开的各个方面。

这里参照根据本公开实施例的方法、设备（系统）、和计算机程序产品的流程图和/或框图描述了本公开的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。

这些计算机可读程序指令可以提供给语音交互装置中的处理器、通用计算机、专用计算机或其它可编程数据处理装置的处理单元，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理单元执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。

也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。

附图中的流程图和框图显示了根据本公开的多个实施例的设备、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，该模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

以上已经描述了本公开的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。

以上该仅为本公开的可选实施例，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等效替换、改进等，均应包含在本公开的保护范围之内。