CN115064211B - 一种基于全基因组甲基化测序的ctDNA预测方法及装置 - Google Patents

一种基于全基因组甲基化测序的ctDNA预测方法及装置 Download PDF

Info

Publication number
CN115064211B
CN115064211B CN202210977623.9A CN202210977623A CN115064211B CN 115064211 B CN115064211 B CN 115064211B CN 202210977623 A CN202210977623 A CN 202210977623A CN 115064211 B CN115064211 B CN 115064211B
Authority
CN
China
Prior art keywords
sample
methylation
block
principal component
training set
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210977623.9A
Other languages
English (en)
Other versions
CN115064211A (zh
Inventor
韩天澄
李苏星
李溪
李宇龙
洪媛媛
张琦
黄宇
陈维之
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhenyue Biotechnology Jiangsu Co ltd
Zhenhe Beijing Biotechnology Co ltd
Original Assignee
Zhenyue Biotechnology Jiangsu Co ltd
Zhenhe Beijing Biotechnology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhenyue Biotechnology Jiangsu Co ltd, Zhenhe Beijing Biotechnology Co ltd filed Critical Zhenyue Biotechnology Jiangsu Co ltd
Priority to CN202210977623.9A priority Critical patent/CN115064211B/zh
Publication of CN115064211A publication Critical patent/CN115064211A/zh
Application granted granted Critical
Publication of CN115064211B publication Critical patent/CN115064211B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/30—Detection of binding sites or motifs
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Biotechnology (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Theoretical Computer Science (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Molecular Biology (AREA)
  • Genetics & Genomics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明提供了一种基于全基因组甲基化测序的ctDNA预测方法及其应用。所述预测方法包括：根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，生成可计算甲基化水平的区块列表文件；计算训练集样本和待测样本在每个区块上的甲基化水平，对训练集样本和待测样本的甲基化水平进行降维，使用主成分分析法，根据训练集样本计算旋转矩阵，并利用旋转矩阵，生成训练集主成分矩阵和待测样本主成分矩阵，利用训练集主成分矩阵进行甲基化模型构建，用训练后的模型对待测样本中的ctDNA进行预测。所述方法克服了低深度甲基化测序中甲基化水平计算不准确、灵敏度低的缺陷，在肿瘤筛查或实时监控中具有重要应用价值。

Description

一种基于全基因组甲基化测序的ctDNA预测方法及装置
技术领域
本发明属于生物医学技术领域，具体涉及一种基于全基因组甲基化测序的ctDNA预测方法及其应用。
背景技术
循环肿瘤DNA（circulating tumor，ctDNA）由肿瘤细胞因分泌、凋亡或坏死而产生，是循环游离DNA（circulating cell-free DNA，cfDNA）中的一种。ctDNA在血液中的半衰期短，且携带有部分肿瘤细胞特有的特征，相比较于组织活检依从性好，可用于肿瘤患者的早期筛查或实时监控。除单核苷酸多态性（Single Nucleotide Polymorphisms，SNP）、插入缺失标记（insertion-deletion，InDel）与拷贝数变异（copy number variation，CNV）外，甲基化作为基因表达调控中的重要环节，同样能够影响基因组的稳定性。对于一些特定位点或区域的甲基化状况，肿瘤患者的ctDNA和健康人的cfDNA之间会存在明显差异。除特定位点或区域的甲基化状况外，也有文献报道肿瘤患者与健康人在全基因组范围内整体甲基化水平的差异。以此，通过从血浆中检测甲基化状况，能够在肿瘤发生的早期、或在肿瘤术后监测中识别出血浆中ctDNA的存在，为后续癌症的早期诊断或复发预测提供数据基础。
近年来甲基化测序的应用在一定程度上提高了ctDNA的检测灵敏度，但大多是依靠缩减测序的区域、提高测序深度达成的。在这种策略下，检测的甲基化位点数量与区域大小较为有限，容易遗漏一些具有检测意义的重要区域。全基因组甲基化测序虽然深度较低，但覆盖位点更全，更容易捕捉到全基因组范围内甲基化水平的变化。
由于各测序方式都需要在测序前对DNA进行重亚硫酸盐转化或酶转化，转化效率对于甲基化水平计算的准确性会产生较大的影响。由于患者血浆样本中ctDNA的占比通常很低，ctDNA的信号极易被淹没在噪音中，基于低深度全基因组甲基化测序的检测的灵敏度与特异性仍旧存在很大的优化空间。
此外，在一些研究中对锯齿末端现象进行了报道。该现象主要出现在cfDNA样本中，主要体现为DNA双链断裂位置的不一致。由于双链建库中，会将锯齿末端修剪或修复为平末端（使双链起始终止位置一致），而修复补齐时，使用的胞嘧啶均为无甲基基团修饰的胞嘧啶。若不进行额外的处理，这些胞嘧啶在测序中均会被判断为非甲基化，导致对应位置的甲基化水平被低估。
因此，提供一种能快速、简便和准确的ctDNA的预测方法，评估样本中存在ctDNA的概率，对肿瘤患者的早期筛查或实时监控具有重要意义。
发明内容
针对现有技术存在的不足，本发明的目的在于提供一种基于全基因组甲基化测序的ctDNA预测方法及其应用，所述ctDNA预测方法能快速简便地对全基因组范围内的甲基化水平进行准确计算，并评估样本中存在ctDNA的概率。所述预测方法可应用于血浆样本，克服了目前低深度甲基化测序中甲基化水平计算不准确、灵敏度不高的技术缺陷。
为达到此发明目的，本发明采用以下技术方案：
第一方面，本发明提供一种基于全基因组甲基化测序的ctDNA预测方法，所述预测方法包括如下步骤：
根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，生成可计算甲基化水平的区块列表文件；
计算训练集样本和待测样本在每个区块上的甲基化水平，对训练集样本和待测样本的甲基化水平进行降维，使用主成分分析法，根据训练集样本计算旋转矩阵，并利用旋转矩阵，生成训练集主成分矩阵和待测样本主成分矩阵；
利用训练集主成分矩阵进行甲基化模型构建，用训练后的模型对待测样本中的ctDNA进行预测。
优选地，所述将基因组划分为等长的区块采用如下步骤进行：
根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，根据全基因组范围内所有CpG位点的位置信息和多重比对区域的占比进行过滤，生成可计算甲基化水平的区块列表。
优选地，所述甲基化水平的计算方法包括：
分别获取训练集样本和待测样本的捕获测序FASTQ文件，进行序列比对，分别生成训练集样本和待测样本原始Bam文件，对原始Bam文件进行排序和标记重复处理，分别生成训练集样本和待测样本排序去重后Bam文件，记为样本j去重排序后的Bam文件；
将样本j去重排序后的Bam文件和所述可计算甲基化水平的区块列表作为输入文件，计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平，合并所有样本的校正后甲基化水平至一个I行J列的甲基化水平矩阵X。
优选地，所述训练集样本为正常人群血浆样本和患者的血浆样本；所述待测样本为被检测者的血浆样本；
优选地，计算甲基化水平前还包括将样本j在区块R i 上片段集合{S 1 ,⋯,S N }中的每一条片段S n ，切除每个片段XM标签上，序列2起始端的前A个碱基的信息。
本发明中，通过锯齿末端切除，以减少锯齿末端现象对区块R i 的平均甲基化水平beta i,j 、区块R i 的平均转化效率BS i,j 、片段S n 的甲基化水平beta i,j,n 、片段S n 的转化效率BS i,j,n 计算的影响。
优选地，计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平的方法包括：
方案一：根据非甲基化位点上胞嘧啶的状态，计算区块转化效率，并使用区块转化效率校正区块甲基化水平；
或方案二：根据每条原始基因片段非甲基化位点上胞嘧啶的状态，计算片段转化效率，过滤低转化率片段后再计算区块的甲基化水平。
优选地，所述方案一中，甲基化水平的计算公式如下所示：
Figure 372017DEST_PATH_IMAGE001
其中，beta_corrected i,j 为样本j在区块R i 的校正后平均甲基化水平；
beta i,j 为样本j在区块R i 的平均甲基化水平；
BS i,j 为样本j在区块R i 的平均转化效率。
优选地，所述方案二中，甲基化水平的计算公式如下所示：
Figure 329609DEST_PATH_IMAGE002
Figure 535462DEST_PATH_IMAGE003
其中，MFR i,j 为样本j在区块R i 的甲基化片段占比；
UFR i,j 为样本j在区块R i 的非甲基化片段占比；
N convertedi,j 为
Figure 527689DEST_PATH_IMAGE004
的片段的统计数量，BS i,j,n 为样本j在区块R i 上片段n的 转化效率；
N meth i,j 为甲基化片段的统计数量；
N unmeth i,j 为非甲基化片段的统计数量。
优选地，所述甲基化片段和非甲基化片段的判断标准为：
Figure 648878DEST_PATH_IMAGE005
的片段 判断为甲基化片段，
Figure 93766DEST_PATH_IMAGE006
的片段判断为非甲基化片段，beta i,j,n 为样本j在区块R i 上片段n的平均甲基化水平，
Figure 86999DEST_PATH_IMAGE007
。
优选地，所述训练集主成分矩阵采用如下方式获取得到：
将训练集样本的甲基化水平矩阵X作为输入文件，进行主成分分析，输出对应的旋 转矩阵V和训练集主成分矩阵Y；计算每个主成分Y p 对总方差的贡献λ p ，并计算前m个主成分 的累计贡献率，提取出累计贡献率超过t的前m个主成分到主成分列表中，
Figure 199311DEST_PATH_IMAGE008
；输 出旋转矩阵V、训练集主成分矩阵Y和主成分列表。
优选地，所述待测样本主成分矩阵采用如下方式获取得到：
将待测样本的甲基化水平矩阵X和所述旋转矩阵V作为输入文件，进行主成分分析，输出对应的待测样本主成分矩阵Y。
优选地，利用训练集主成分矩阵进行甲基化模型构建采用如下步骤进行：
首先根据所述主成分列表提取训练集主成分矩阵Y中所有样本的前m个主成分到新的矩阵Y'，将新矩阵Y'作为支持向量机模型的输入特征，采用装袋法进行建模，得到一个基于K个弱分类器建立的集成模型M，模型M中，每个弱分类器的权重相等，最终的预测值是K个弱分类器预测值的平均值。
优选地，对待测样本中的ctDNA进行预测采用如下步骤进行：
根据训练出的模型M中记录的模型变量，根据所述主成分列表提取待测样本主成分矩阵Y中所有待测样本的前m个主成分到新的矩阵Y'，并将Y'输入到训练出的模型M中，获得每个待测样本的预测值。
第二方面，本发明提供一种基于全基因组甲基化测序的ctDNA预测装置，所述预测装置包括：
甲基化区块划分模块：根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，生成可计算甲基化水平的区块列表文件；
区块甲基化水平计算与校正模块：用于计算训练集样本和待测样本在每个区块上的甲基化水平；
甲基化区块特征降维模块：对训练集样本和待测样本的甲基化水平进行降维，使用主成分分析法，根据训练集样本计算旋转矩阵，并利用旋转矩阵，生成训练集主成分矩阵和待测样本主成分矩阵；
甲基化模型构建模块：利用训练集主成分矩阵进行甲基化模型构建，对模型进行训练；
甲基化模型预测模块：用训练后的模型对待测样本中的ctDNA进行预测。
本发明中，所述甲基化区块划分模块的输入文件包括：人类参考基因组序列FASTA文件和参考基因组染色体BED文件。
所述甲基化区块划分模块将人类参考基因组按照输入参数L对应的区块长度进行划分，得到区块R i ，根据CpG位点数量和多重比对区域的占比进行过滤后输出。
本发明中，所述甲基化区块划分模块的输出文件包括：区块列表文件，所述区块列表文件为经过两步过滤后，保留下I个等长的区块{R 1 ,⋯,R I }。
优选地，所述甲基化区块划分模块中，所述将基因组划分为等长的区块采用如下步骤进行：
根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，根据全基因组范围内所有CpG位点的位置信息和多重比对区域的占比进行过滤，生成可计算甲基化水平的区块列表。
本发明中，所述区块甲基化水平计算与校正模块的输入文件包括：样本j去重排序后的Bam文件、由甲基化区块划分模块生成的区块列表文件，所述样本包括训练集样本和待测样本；
本发明中，所述区块甲基化水平计算与校正模块的输出文件包括：甲基化水平矩阵，所述甲基化水平矩阵X为一个I行J列的矩阵，其中，行代表区块，列代表样本。
优选地，所述区块甲基化水平计算与校正模块中，所述甲基化水平的计算方法包括：
分别获取训练集样本和待测样本的捕获测序FASTQ文件，进行序列比对，分别生成训练集样本和待测样本原始Bam文件，对原始Bam文件进行排序和标记重复处理，分别生成训练集样本和待测样本排序去重后Bam文件，记为样本j去重排序后的Bam文件；
将样本j去重排序后的Bam文件和所述可计算甲基化水平的区块列表作为输入文件，计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平，合并所有样本的校正后甲基化水平至一个I行J列的甲基化水平矩阵X；
优选地，所述训练集样本为正常人群血浆样本和患者的血浆样本；所述待测样本为被检测者的血浆样本；
优选地，计算甲基化水平前还包括将样本j在区块R i 上片段集合{S 1 ,⋯,S N }中的每一条片段S n ，切除每个片段XM标签上，序列2起始端的前A个碱基的信息。
优选地，计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平的方法包括：
方案一：根据非甲基化位点上胞嘧啶的状态，计算区块转化效率，并使用区块转化效率校正区块甲基化水平；
或方案二：根据每条原始基因片段非甲基化位点上胞嘧啶的状态计算片段转化效率，过滤低转化率片段后再计算区块的甲基化水平。
优选地，所述方案一中，甲基化水平的计算公式如下所示：
Figure 472161DEST_PATH_IMAGE010
其中，beta_corrected i,j 为样本j在区块R i 的校正后平均甲基化水平；
beta i,j 为样本j在区块R i 的平均甲基化水平；
BS i,j 为样本j在区块R i 的平均转化效率。
优选地，所述方案二中，甲基化水平的计算公式如下所示：
Figure 420657DEST_PATH_IMAGE011
Figure 968313DEST_PATH_IMAGE012
其中，MFR i,j 为样本j在区块R i 的甲基化片段占比；
UFR i,j 为样本j在区块R i 的非甲基化片段占比；
N convertedi,j 为
Figure 935131DEST_PATH_IMAGE004
的片段的统计数量，BS i,j,n 为样本j在区块R i 上片段n的 转化效率；
N meth i,j 为甲基化片段的统计数量；
N unmeth i,j 为非甲基化片段的统计数量。
优选地，所述甲基化片段和非甲基化片段的判断标准为：
Figure 175620DEST_PATH_IMAGE013
的片段 判断为甲基化片段，
Figure 844368DEST_PATH_IMAGE006
的片段判断为非甲基化片段，
Figure 461294DEST_PATH_IMAGE014
。
本发明中，所述甲基化区块特征降维模块的输入文件包括：训练集样本甲基化水平矩阵X；输出文件包括：训练集主成分矩阵Y、主成分列表和旋转矩阵V。
或者，所述甲基化区块特征降维模块的输入文件包括：待测样本甲基化水平矩阵X、旋转矩阵V；输出文件包括：待测样本主成分矩阵Y、主成分列表。
优选地，所述甲基化区块特征降维模块中，所述训练集主成分矩阵采用如下方式获取得到：
将训练集样本的甲基化水平矩阵X作为输入文件，进行主成分分析，输出对应的旋 转矩阵V和训练集主成分矩阵Y；计算每个主成分Y p 对总方差的贡献λ p ，并计算前m个主成分 的累计贡献率，提取出累计贡献率超过t的前m个主成分到主成分列表中，
Figure 282619DEST_PATH_IMAGE008
；输 出旋转矩阵V、训练集主成分矩阵Y和主成分列表。
累计贡献率的计算公式如公式8所示：
Figure 959588DEST_PATH_IMAGE015
优选地，所述甲基化区块特征降维模块中，所述待测样本主成分矩阵采用如下方式获取得到：
将待测样本的甲基化水平矩阵X和所述旋转矩阵V作为输入文件，进行主成分分析，输出对应的待测样本主成分矩阵Y。
本发明中，所述甲基化模型构建模块的输入文件包括：训练集主成分矩阵Y、主成分列表和样本列表，样本列表需要包含所有训练集样本的分组信息，包含样本名、样本分组两列，每一行记录一个样本。
本发明中，所述甲基化模型构建模块的输出文件包括：训练出的模型M、训练集的预测概率结果。
优选地，所述甲基化模型构建模块中，利用训练集主成分矩阵进行甲基化模型构建采用如下步骤进行：
首先根据所述主成分列表提取训练集主成分矩阵Y中所有样本的前m个主成分到新的矩阵Y'，将新矩阵Y'作为支持向量机模型的输入特征，采用装袋法进行建模，得到一个基于K个弱分类器建立的集成模型M，模型M中，每个弱分类器的权重相等，最终的预测值是K个弱分类器预测值的平均值。
本发明中，所述甲基化模型预测模块的输入文件包括：待测样本主成分矩阵Y、训练出的模型M和待测样本列表，待测样本列表文件作为输入文件，其中需要包含所有待测样本。
本发明中，所述甲基化模型预测模块的输出文件包括：待测样本的预测值，所述预测值是0-1之间的一个值，代表了该样本来源于患者血浆组的概率，即血浆中包含ctDNA的概率结果。
优选地，所述甲基化模型预测模块中，对待测样本中的ctDNA进行预测采用如下步骤进行：
根据训练出的模型M中记录的模型变量，根据所述主成分列表提取待测样本主成分矩阵Y中所有待测样本的前m个主成分到新的矩阵Y'，并将Y'输入到训练出的模型M中，获得每个待测样本的预测值。
本发明中，所述基于全基因组甲基化测序的ctDNA检测装置的关键点在于：
（1）针对低深度全基因组甲基化测序设计，捕捉全基因组范围内的甲基化水平波动
相较于高深度捕获测序，低深度全基因组甲基化测序的实验操作更简便，测序成本不高，覆盖的区域更全，能够避免因为捕获区间太少而遗漏信号的问题，更好地反映全基因组范围内的甲基化水平。
（2）以等长的区块为单位计算甲基化水平
本发明针对全基因组甲基化测序的数据特征进行了方案设计，考虑到全基因组甲基化测序的深度低，低深度下，对单个甲基化位点或较短的区间计算甲基化水平的准确度较低。本发明按照相同的长度对全基因组范围内的区域进行区块划分，每个区块的覆盖区域较广、甲基化位点数量较多，通过增加位点数量的方式，减少对深度的要求，弥补了传统分析方法中由于深度低而导致准确率低的缺陷。此外，相比于计算单个甲基化位点的甲基化水平，本发明计算区块甲基化水平的速度更快，且减少了特征的数量，大大节约了时间成本。
（3）去除包含多重比对区域的区块，降低比对错误的负面影响
在划分等长区块时，本发明还对每个区块上多重比对区域的长度进行了统计，去除了多重比对区域占比过高的区块，减少了比对错误带来的噪音，对甲基化水平计算和模型构建的影响，提高了模型准确性。
（4）针对双链建库，切除锯齿末端
本发明的区块甲基化水平计算与校正模块，在计算甲基化水平时考虑了锯齿末端的影响，预先对每条片段修复补齐端的碱基进行了切除，避免了修复步骤导致的对甲基化水平的低估，提升了甲基化水平计算的准确性，提高了模型的灵敏度。
（5）提供多种方案，根据转化效率校正甲基化水平
本发明的区块甲基化水平计算与校正模块，在传统甲基化水平计算方式的基础上，额外考虑了转化效率对甲基化水平产生的影响，提供了两种解决方案。方案一对区块平均甲基化水平进行校正，抵消了转化效率对甲基化水平计算的一部分影响；方案二通过对低转化率片段进行过滤，减少转化效率对判断片段甲基化状态的影响。两种方案均提高了甲基化水平计算的准确性与模型的灵敏度。
（6）主成分分析降维，减少模型构建时间，避免过拟合
相比常见的甲基化数据分析方法与模型，本发明不对甲基化区域进行差异筛选，不区分健康人血浆和患者血浆，在甲基化区块特征降维模块中使用主成分分析进行特征降维，提取在所有样本中对差异解释力最强的主成分进行建模，能够在一定程度上避免过拟合的问题，并大大缩短模型构建和预测的时间。
（7）使用血浆样本训练支持向量机，模型表现好，预测结果易于解读
本发明在甲基化模型构建模块中使用健康人血浆和患者血浆作为训练集，用主成分分析降维后获得的主成分作为变量，训练支持向量机集成模型，对待测样本与健康人或患者血浆的相似性进行量化。最终获得的模型同时具有较高的灵敏度和特异性，模型的泛化能力强。模型输出的预测值代表了待测样本含有ctDNA的概率，根据训练集中获得的阈值，可以对待测样本是否含有ctDNA进行判断。
第三方面，本发明提供第一方面所述的基于全基因组甲基化测序的ctDNA预测方法和/或第二方面所述的基于全基因组甲基化测序的ctDNA预测模型在制备ctDNA检测产品中的应用。
第四方面，本发明提供一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器运行所述计算机程序时实现第一方面所述的基于全基因组甲基化测序的ctDNA预测方法的步骤。
第五方面，本发明提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现第一方面所述的基于全基因组甲基化测序的ctDNA预测方法的步骤。
相对于现有技术，本发明具有以下有益效果：
本发明对于血浆样本低深度全基因组甲基化测序数据的甲基化水平计算进行了有针对性的优化：通过切除锯齿末端、使用校正后的区块甲基化水平，能够快速准确地计算出样本在全基因组范围内真实的甲基化水平；在优化甲基化水平计算的基础上，使用全基因组范围内的甲基化水平提取主成分进行降维，并根据较少数量的主成分构建支持向量机集成模型，能更快更准确地对待测血浆样本中存在ctDNA的概率进行预测，并避免过拟合。实例中使用回顾性队列的低深度血浆全基因组甲基化测序数据，对本发明进行了验证。测试结果显示，在低深度下，本发明可以在测试集中获得较高的灵敏度和特异性。
附图说明
图1是本发明技术方案的主流程图。
图2是测试集中肺癌的检测ROC曲线。
图3是测试集中结直肠癌的检测ROC曲线。
图4是测试集中胰腺癌的检测ROC曲线。
图5是测试集中胃癌的检测ROC曲线。
图6是测试集中肝癌的检测ROC曲线。
图7是测试集中乳腺癌的检测ROC曲线。
图8是测试集中食管癌的检测ROC曲线。
图9是测试集中卵巢癌的检测ROC曲线。
图10是测试集中全癌种的检测ROC曲线。
具体实施方式
下面通过具体实施方式来进一步说明本发明的技术方案。本领域技术人员应该明了，所述实施例仅仅是帮助理解本发明，不应视为对本发明的具体限制。
实施例中未注明具体技术或条件者，按照本领域内的文献所描述的技术或条件，或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可通过正规渠道商购获得的常规产品。
发明详述
本发明的技术原理：本发明针对全基因组甲基化测序数据的特点，将全基因组划分为等长的区块，在去除比对质量不佳的区块后，针对每个区块计算区块甲基化水平。为了克服锯齿末端对甲基化水平计算的影响，本发明专门针对双链建库、双端测序的样本进行了锯齿末端的切除。
为了克服转化效率对甲基化水平的影响，在计算区块甲基化水平时，本发明提供了两种方案：
方案一：根据非甲基化位点上胞嘧啶的状态，计算区块转化效率，并使用区块转化效率校正区块甲基化水平；
方案二：根据每条原始基因片段非甲基化位点上胞嘧啶的状态计算片段转化效率，过滤低转化率片段后再计算区块的甲基化水平。
获得校正后区块甲基化水平后，本发明使用一组健康人血浆样本和一组肿瘤患者血浆样本作为训练集，使用主成分分析的方法进行降维，并用降维后较少数量的特征作为支持向量机的输入变量，建立甲基化模型。
对每个待测样本，采用和训练集同样的方式计算校正后区块甲基化水平和降维，并使用甲基化模型对样本进行预测，获得待测样本中含有ctDNA的概率，提示样本中存在ctDNA的可能性。
本发明的实现过程：本发明的技术方案的主流程图如图1所示。
本发明的技术方案主要包括以下几个步骤：
（1）获取正常人群血浆样本（后续表述为基线样本）和被检测者的血浆样本（后续表述为待测样本）的捕获测序FASTQ文件，利用基因组比对工具Bismark进行序列比对，生成原始Bam文件；
（2）利用SAMtools和Bismark工具，对原始Bam文件进行排序和标记重复处理，生成排序去重后Bam文件；
（3）利用本发明的甲基化区块划分模块，根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，剔除CpG位点数量过少、比对质量较差的区块，生成可计算甲基化水平的区块列表；
（4）利用本发明的区块甲基化水平计算与校正模块，计算训练集样本和待测样本在每个区块上的校正后甲基化水平。计算过程中，本发明排除了末端锯齿影响，并依据转化效率对甲基化水平进行了校正。校正后的计算结果为两个矩阵：训练集样本校正后区块甲基化水平矩阵与待测样本校正后区块甲基化水平矩阵；
（5）利用本发明的甲基化区块特征降维模块，对训练集样本的校正后区块甲基化水平进行降维，使用主成分分析法，根据训练集样本计算旋转矩阵，并利用旋转矩阵，生成训练集主成分矩阵和待测样本主成分矩阵；
（6）利用本发明的甲基化模型构建模块，根据训练集样本的主成分矩阵对模型进行训练，构建甲基化模型；
（7）利用本发明的甲基化模型预测模块，使用甲基化模型对待测样本进行预测，最终给出待测样本中包含ctDNA的概率得分，所述概率得分能够反映出待测样本的甲基化水平是否存在异常，提示血浆中是否有ctDNA的存在。
本发明的技术方案主要包括以下几个模块：
（1）甲基化区块划分模块
该模块使用人类参考基因组序列FASTA文件和参考基因组染色体BED文件（记录基因组上每个染色体的起始和终止位置）作为输入文件，并要求提供区块长度作为输入参数L。将人类参考基因组按照一定的区块长度进行划分，根据一定的规则进行过滤后，输出记录区块染色体、起始和终止位置信息的区块列表文件。
首先，利用人类参考基因组序列FASTA文件，提取全基因组范围内所有CpG位点的位置信息。
将FASTA文件分别进行两次碱基转换：一次将胞嘧啶（cytosine，C）全部转换为胸腺嘧啶（thymine，T），其他碱基保持不变，记录为参考序列A；另一次将鸟嘌呤（Guanine，G）全部转换为腺嘌呤（Adenine，A），其他碱基保持不变，记录为参考序列B。以100个碱基为单位滑动，分别统计参考序列A（或参考序列B）中，每100个碱基序列在整个参考序列A（或参考序列B）中的出现次数，记录下存在多重比对现象（窗口序列出现超过2次）的区域位置。
将全基因组上的常染色体，按照输入的区块长度，划分为多个等长的区块。
对于每一个区块R i ，根据全基因组范围内所有CpG位点的位置信息，统计区块R i 上覆盖的CpG位点数量。若R i 上的CpG数量超过100个，则R i 被保留，否则R i 被剔除。
对于保留下的R i ，再根据存在多重比对现象的区域位置，统计R i 上多重比对区域的占比。若占比低于5%，则R i 被保留，否则R i 被剔除。
经过两步过滤后，保留下I个等长的区块{R 1 ,⋯,R I }，每个区块的染色体、起始和终止位置将按照“染色体：起始-终止”的格式，以列表形式被输出。
（2）区块甲基化水平计算与校正模块
该模块使用样本j去重排序后的Bam文件（比对软件须使用BisMark比对）和由本发明的甲基化区块划分模块生成的区块列表文件作为输入文件，根据特定的方法，计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平。
在计算样本j在每一个区块R i 的甲基化水平时，首先从样本j的Bam文件中提取出比对到区块R i 上的序列，并将来源于同一原始基因片段的两条序列（即序列1和序列2，由双端测序而产生，记录了同一条原始基因片段的序列）合并：根据其比对到参考基因组上的起始和终止位置，将序列1和序列2的XM标签进行合并为片段XM标签。
XM标签是BisMark在进行比对时同时输出到Bam文件中的标记，记录序列在每个胞嘧啶（cytosine，C）位置上的碱基状态。根据下游碱基的不同，C在基因组上的分布包含三种形式：CpG、CHH和CHG。其中，CpG代表C下游的一个碱基为鸟嘌呤（Guanine，G），CHH代表C下游的两个碱基都是H（除鸟嘌呤以外的其他碱基），CHG代表C下游的两个碱基是H和G。不同形式的C的标记字母不同，大写字母代表序列在对应位置上的碱基为C，小写字母代表对应位置的碱基为胸腺嘧啶（thymine，T）：CpG上的C由Z或z表示，CHH上的C由H或h表示，CHG上的C由X或x表示。
为了减少在双链建库方式下，锯齿末端现象对片段S n 的甲基化水平beta i,j,n 、片段S n 的转化效率BS i,j,n 、区块R i 的平均甲基化水平beta i,j 和区块R i 的平均转化效率BS i,j 计算产生的影响，对于每对序列1和序列2合并产生的片段XM标签，切除建库过程中可能发生补齐的一端的前A个碱基的信息。
对于区块上片段集合{S 1 ,⋯,S N }中的每一条片段S n ，根据去除碱基后的CpG标签，计算片段S n 的甲基化水平beta i,j,n ：
Figure 866364DEST_PATH_IMAGE016
根据去除碱基后的CHH和CHG标签，计算片段S n 的转化效率BS i,j,n ：
Figure 21402DEST_PATH_IMAGE017
其中Z n ，z n ，H n ，h n ，X n 和x n 分别代表片段S n 上Z，z，H，h，X和x的标签数量。
根据样本j在区块R i 上所有的{S 1 ,⋯,S N }，计算样本j在区块R i 的平均甲基化水平beta i,j ：
Figure 710616DEST_PATH_IMAGE018
和样本j在区块R i 的平均转化效率BS i,j ：
Figure 558486DEST_PATH_IMAGE019
样本j在区块R i 的校正后平均甲基化水平beta_corrected i,j 可以根据beta i,j 与BS i,j 计算：
Figure 952559DEST_PATH_IMAGE020
另外，区块甲基化水平的计算还可以通过计算样本j在区块R i 的甲基化片段占比MFR i,j 和非甲基化片段占比UFR i,j 获得。
在计算时，首先挑选出转化率高于阈值（以0.95阈值为例，即满足
Figure 911287DEST_PATH_IMAGE021
）的 片段，统计数量记录为N convertedi,j 。再根据N convertedi,j 条片段的beta i,j,n 是否高于阈值beta methy （
Figure 441626DEST_PATH_IMAGE022
，以0.9阈值为例，即满足
Figure 709665DEST_PATH_IMAGE023
）判断每条片段是否为甲基化 片段，或是否低于阈值beta unmethy （
Figure 591033DEST_PATH_IMAGE024
，以0.1阈值为例，即满足
Figure 353453DEST_PATH_IMAGE025
）判断 是否为非甲基化片段。甲基化片段的数量记录为N meth i,j ，非甲基化片段的数量记录为N unmeth i,j 。
最终，样本j在区块R i 的甲基化片段占比公式如下：
Figure 3877DEST_PATH_IMAGE026
样本j在区块R i 的非甲基化片段占比公式如下：
Figure 927971DEST_PATH_IMAGE027
对于每一样本j，最终输出的校正后区块甲基化水平为一个长度为I的向量。若使用平均甲基化水平，则输出：
Figure 47368DEST_PATH_IMAGE028
若使用甲基化片段占比，则输出：
Figure 347899DEST_PATH_IMAGE029
若使用非甲基化片段占比，则输出：
Figure 118409DEST_PATH_IMAGE030
合并所有样本的校正后区块甲基化水平至一个I行J列的矩阵X，其中，行代表区块，列代表样本，写入文件并输出至文件中。
（3）甲基化区块特征降维模块
在第一次使用甲基化区块特征降维模块时，需提供区块甲基化水平计算与校正模块输出的甲基化水平矩阵X文件作为输入文件。矩阵中需要包含所有后续将用于构建模型的训练集样本。
对I个区块的甲基化水平进行主成分分析，记录对应的旋转矩阵V和每个训练集样本在P个主成分上的值在主成分矩阵Y。
在获得旋转矩阵V和主成分矩阵Y后，计算每个主成分Y p 对总方差的贡献λ p ，并计算 前m个主成分的累计贡献率，累计贡献率的计算公式为公式8。记录提取出累计贡献率超过t
Figure 213404DEST_PATH_IMAGE031
的前m个主成分到主成分列表中。
Figure 69364DEST_PATH_IMAGE032
甲基化区块特征降维模块在第一次使用时输出旋转矩阵V、主成分矩阵Y和主成分列表三个文件。
当获得旋转矩阵V后，可在后续再次运行该模块时，在输入X矩阵文件的同时额外提供旋转矩阵V文件。此时，模块会直接使用V将甲基化水平矩阵X转换为主成分矩阵Y。此时，模块只输出主成分矩阵Y到输出文件中。
（4）甲基化模型构建模块
该模块使用甲基化区块特征降维模块输出的主成分矩阵Y、主成分列表和样本列表三个文件作为输入文件。矩阵中Y需要包含所有用于构建模型的训练集样本，包括健康人血浆和患者血浆。主成分列表包含了前m个主成分的名字，对使用的主成分数量进行指定。样本列表需要包含所有训练集样本的分组信息，包含样本名、样本分组两列，每一行记录一个样本。
在训练模型前，首先提取主成分矩阵Y中所有训练集样本的前m个主成分到新的矩阵Y'。将新矩阵Y'作为支持向量机模型的输入特征。
在建模时，采用装袋法，分别进行K次有放回的随机抽样，建立K个支持向量机弱分类器。对每个弱分类器单独进行参数寻优、训练和预测。最终的模型是一个基于K个弱分类器建立的集成模型M。模型M中，每个弱分类器的权重相等，最终的预测值是K个弱分类器预测值的平均值。
在训练模型的同时，甲基化模型构建模块会对每个训练集样本进行预测，并记录相应的预测值。预测值是0-1之间的一个值，代表了该样本来源于患者血浆组的概率，即血浆中包含ctDNA的概率。
甲基化模型构建模块最终会输出训练出的模型M、训练集的预测概率结果至文件中。其中，模型M采用特定的格式存储为二进制文件，训练集的预测概率是一个包含三列的表格，分别为样本名、样本所属分组和样本预测概率，每一行记录一个训练集样本的信息。
（5）甲基化模型预测模块
甲基化模型预测模块使用甲基化区块特征降维模块输出的主成分矩阵Y文件和甲基化模型构建模块输出的模型M、待测样本列表文件作为输入文件，矩阵中需要包含所有待测样本。
甲基化模型预测模块会根据模型M中记录的模型变量，提取主成分矩阵Y中所有待测样本的前m个主成分到新的矩阵Y'，并将Y'输入到模型M中，获得每个待测样本的预测值。预测值是0-1之间的一个值，代表了该样本来源于患者血浆组的概率，即血浆中包含ctDNA的概率。
甲基化模型预测模块最终会输出待测样本的预测概率结果至文件中。该文件包含两列的表格，分别为样本名、和样本预测概率，每一行记录一个待测样本的信息。
本发明的技术改进点：使用区块为单位计算甲基化水平，提高计算效率，使低深度样本中的甲基化水平计算更稳定；切除锯齿末端，减少末端修复对甲基化水平计算的影响；根据转化效率进行校正，提高血浆甲基化水平计算的准确性；在建模前先进行特征降维，避免过拟合，提高建模和预测速度；使用支持向量机集成模型，灵敏性、特异性较高。
为验证技术方案，回顾性选取了497例无癌症史的健康人血浆以及801例不同分期的多癌种癌症患者的血浆，并随机分组为训练集和验证集。患者的癌症种类包括了乳腺癌、结直肠癌、食管癌、胃癌、肝癌、肺癌、胰腺癌和卵巢癌。
训练集包括352例健康人及561例癌症患者（46例乳腺癌、105例结直肠癌、43例食管癌、80例胃癌、79例肝癌、110例肺癌、83例胰腺癌、15例卵巢癌）。
验证集包括145例健康人和240例癌症患者（20例乳腺癌、45例结直肠癌、18例食管癌、34例胃癌、34例肝癌、47例肺癌、36例胰腺癌、6例卵巢癌）。
实施例1验证技术方案的流程
实验流程如下所示：
1、血浆cfDNA提取
1.1每位受试者10mL全血存放在康为EDTA采血管中，在4℃以1600g转速离心10min使血浆、血细胞分层。将上层血浆转移至新离心管，再次以12000rpm转速4℃离心15min取上清以去除细胞碎屑。得到4mL血浆，-80℃冻存备用。
1.2血浆样本融化后，每1mL样本中加入15μL蛋白酶K（Proteinase K，20mg/mL，thermoscientific cat#EO0492）和50μL SDS（20%）。血浆量不足4mL，用PBS补足。
1.3翻转混匀，60℃孵育20min，然后冰浴5min。
1.4使用MagMAX Cell-Free DNA Isolation试剂盒（thermoscientific cat#A29319）提取cfDNA。
1.5使用Bioanalyzer 2100（Agilent Technologies）检测cfDNA的提取浓度和质量。
2、cfDNA文库构建
使用甲基化文库构建试剂盒NEBNext Enzymatic Methyl-seq Kit（NEB，cat#E7120），以5-30ng cfDNA起始量，通过TET2酶使5-甲基胞嘧啶（5-mC）转化为5-甲酰胞嘧啶（5-fC）和5-羧基胞嘧啶（5-caC），并且通过APOBEC酶，使非甲基化胞嘧啶（C）脱氨转化为尿嘧啶（U），然后进行扩增建库。
具体文库构建过程如下：
2.1内参准备
取50μL CpG全甲基化的pUC19 DNA和50μL CpG全非甲基化的lambda DNA混匀后加入100μL打断管中，使用M220打断仪（Covaris）打断。建库时，向待测cfDNA加入0.001ng的pUC19 DNA和0.02 ng的lambda DNA。
2.2cfDNA样本的准备
cfDNA样本起始量为5-30ng，不需要打断。
2.3末端修复
2.3.1在冰上混合，反应体系如表1所示；
表1
Figure 157275DEST_PATH_IMAGE033
表1中所述酶混合液为NEBNext Ultra II End Prep Enzyme Mix。
2.3.2反应体系置于PCR仪上，按表2所示的条件进行末端修复反应。
表2
Figure 782291DEST_PATH_IMAGE034
2.4连接Adaptor（接头）
2.4.1在冰上操作，按表3所示，将以下组分加入上步的60μL反应体系中；
表3
Figure 48188DEST_PATH_IMAGE035
2.4.220℃孵育15 min。
2.5连接后纯化
2.5.1上一步反应结束后，取出样本，加入110μL NEBNext Sample PurificationBeads，立即使用移液器吹打混匀。
2.5.2室温孵育5 min。
2.5.3离心管置于磁力架上5 min待液体澄清，弃去上清。
2.5.4加入200μL现配80%乙醇，孵育30 s后弃去。重复一次200μL 80%乙醇清洗步骤。
2.5.5用10μL移液器吸尽离心管底部的残留乙醇，室温干燥3-5min至乙醇完全挥发。
2.5.6从磁力架取下离心管，加入29 μL Elution Buffer，震荡混匀，室温孵育1min。
2.5.7短暂离心，离心管置于磁力架上3 min待液体澄清，取28 μL放进新的PCR管中。
2.6 5-甲基胞嘧啶和5-羟甲基胞嘧啶氧化反应
使用NEBNext Enzymatic Methyl-seq Kit（NEB，cat# E7120）进行以下反应操作。
2.6.1 TET2 Reaction Buffer Supplement干粉加入400μL TET2 ReactionBuffer，充分混合。
2.6.2在冰上按表4所示，将以下组分加入上述28 μL已连接adaptor的DNA：
表4
Figure 391444DEST_PATH_IMAGE036
2.6.3 将500 mM Fe（II）溶液按1:1250比例稀释。往上步混匀的产物中，加入已配好的Fe（II）溶液，各组分的用量如表5所示。
表5
Figure 33778DEST_PATH_IMAGE037
充分混合并在37℃孵育1h。
2.6.4 反应结束后移至冰上并加入1μL Stop Reagent（终止液），充分混合。各组分的用量如表6所示。
表6
Figure 778880DEST_PATH_IMAGE038
2.6.5将2.6.4中混合后的样品孵育，孵育条件如表7所示。
表7
Figure 981059DEST_PATH_IMAGE039
2.7氧化后纯化
2.7.1 上一步反应结束后，取出样本，加入90μL NEBNext Sample PurificationBeads，立即使用移液器吹打混匀。
2.7.2室温孵育5 min。
2.7.3 离心管置于磁力架上5 min待液体澄清，弃去上清。
2.7.4 加入200μL现配80%乙醇，孵育30 s后弃去。重复一次200μL 80%乙醇清洗步骤。
2.7.5用10μL移液器吸尽离心管底部的残留乙醇，室温干燥3-5min至乙醇完全挥发。
2.7.6从磁力架取下离心管，加入17 μL洗脱缓冲液（Elution Buffer），震荡混匀。室温孵育1 min。
2.7.7短暂离心，离心管置于磁力架上3 min待液体澄清，取16 μL放进新的PCR管中。
2.8 DNA变性
2.8.1 配制新鲜的0.1 MNaOH。
2.8.2 提前预热PCR仪到温度为50℃。
2.8.3加入4μL 0.1 M NaOH到上步16μL纯化产物中，充分混合。
2.8.4 50℃孵育10 min。
2.8.5 反应结束后立刻放入冰上。
2.9 胞嘧啶脱氨基
2.9.1在冰上将下列组分加入上步20μL变性DNA，充分混合。各组分的用量如表8所示。
表8
Figure 811611DEST_PATH_IMAGE040
2.9.2 在PCR仪上37℃孵育3h后转为4℃终止反应。
2.10脱氨后纯化
2.10.1 上一步反应结束后，取出样本，加入100μL磁珠（NEBNext SamplePurification Beads），立即使用移液器吹打混匀。
2.10.2室温孵育5 min。
2.10.3 离心管置于磁力架上5 min待液体澄清，弃去上清。
2.10.4 加入200μL现配80%乙醇，孵育30 s后弃去。重复一次200μL 80%乙醇清洗步骤。
2.10.5用10μL移液器吸尽离心管底部的残留乙醇，室温干燥3-5min至乙醇完全挥发。
2.10.6从磁力架取下离心管，加入21 μL洗脱缓冲液（Elution Buffer），震荡混匀，室温孵育1 min。
2.10.7短暂离心，离心管置于磁力架上3 min待液体澄清，取20 μL放进新的PCR管中。
2.11文库PCR扩增
2.11.1 在冰上将下列组分加入上步脱氨后的20μL DNA中，各组分的用量如表9所示。
表9
Figure 992057DEST_PATH_IMAGE041
2.11.2 充分混合后在PCR仪上进行PCR反应，PCR反应程序如表10所示。
表10
Figure 591666DEST_PATH_IMAGE042
2.12 PCR后纯化
2.12.1 上一步反应结束后，取出样本，加入45 μL磁珠（NEBNext SamplePurification Beads），立即使用移液器吹打混匀。
2.12.2室温孵育5 min。
2.12.3 离心管置于磁力架上5 min待液体澄清，弃去上清。
2.12.4 加入200μL现配80%乙醇，孵育30 s后弃去。重复一次200μL 80%乙醇清洗步骤。
2.12.5用10μL移液器吸尽离心管底部的残留乙醇，室温干燥3-5min至乙醇完全挥发。
2.12.6从磁力架取下离心管，加入21 μL Elution Buffer，震荡混匀。室温孵育1min。
2.12.7短暂离心，离心管置于磁力架上3 min待液体澄清，取20 μL放进新的PCR管中。
2.13 文库定量
使用Qubit高灵敏试剂（thermoscientific cat#Q32854）对所构建的文库进行定量，文库产量大于400ng进行后续上机测序。
3、文库测序
取100ng上述文库加入10% PhiX DNA（Illumina cat#FC-110-3001）混合成上机样品，在Novaseq 6000（Illumina）平台进行PE100测序。
实施例2生信分析流程
1、处理下机FASTQ数据为各模块可使用的Bam文件
1.1去接头
调用Trimmomatic-0.36将每一对FASTQ文件都作为paired reads比对到hg19人类参考基因组序列，除M参数与指定Reads Group的ID外，不使用其余参数选项，生成初始Bam文件。
1.2比对
调用Bismark-v0.19.0将去接头后的每一对FASTQ文件都作为paired reads比对到hg19人类参考基因组序列和Lambda DNA参考基因组序列，生成初始Bam文件。
1.3去重
调用Bismark-v0.19.0的deduplicate模块，对初始Bam文件进行去重复处理，生成去重后的Bam文件。
1.4排序标记
调用SAMtools-1.3的sort模块，对去重后的Bam文件进行排序，生成排序后的Bam文件。然后，调用Picard-2.1.0的AddOrReplaceReadGroups模块，对排序后的Bam文件进行标记分组。
1.5筛选
调用BamUtil-1.0.14的clipOverlap模块对标记分组后的Bam文件进行筛选，去除重叠的pair reads，生成Bam文件。并调用SAMtools-1.3 view对去除重叠的Bam文件的比对质量进行过滤，采用“-q 20”作为参数，生成过滤比对质量后Bam文件。
1.6 抽取配对序列
为保持各样本数据量一致，对每个样本，从过滤比对质量后Bam文件中随机抽取30,000,000对配对序列到新Bam文件中，作为最终Bam文件。
1.7建立索引
调用SAMtools-1.3的index模块对最终生成的Bam文件建立索引，生成与最终Bam文件配对的bai文件。
2、甲基化分析
2.1 划分甲基化区块
使用人类参考基因组序列FASTA文件和参考基因组染色体BED文件（记录基因组上每个染色体的起始和终止位置）作为甲基化区块划分模块的输入文件，并设定区块长度参数L=1,000,000。
首先，利用人类参考基因组序列FASTA文件，提取全基因组范围内所有CpG位点的位置信息，并计算出了全基因组范围内存在多重比对现象的所有区域的位置信息。
将人类参考基因组划分为2,734个区块，每个区块包含1,000,000个碱基。统计每个区块上CpG位点的数量和存在多重比对现象的区域长度，过滤CpG位点数量少于100个、或多重比对区域总长超过50000个碱基的区块，最终保留1,846个区块。这些区块将用于后续的分析，并保存在区块列表文件中。
2.2 分别生成训练集和测试集样本的校正后区块平均甲基化水平矩阵
将每个训练集内样本和测试集内样本的最终Bam文件和2.1中生成的区块列表文件作为区块甲基化水平计算与校正模块输入文件。
依次从每一个训练集或测试集样本j的Bam文件中，分别提取Bam中比对到每一个区块R i 的序列，并将来源于同一原始基因片段的两条配对序列的XM标签进行合并，作为片段XM标签。
XM标签是BisMark在进行比对时同时输出到Bam文件中的标记，记录序列在每个胞嘧啶（cytosine，C）位置上的碱基状态。根据下游碱基的不同，C在基因组上的分布包含三种形式：CpG、CHH和CHG。其中，CpG代表C下游的一个碱基为鸟嘌呤（Guanine，G），CHH代表C下游的两个碱基都是H（除鸟嘌呤以外的其他碱基），CHG代表C下游的两个碱基是H和G。不同形式的C的标记字母不同，大写字母代表序列在对应位置上的碱基为C，小写字母代表对应位置的碱基为胸腺嘧啶（thymine，T）：CpG上的C由Z或z表示，CHH上的C由H或h表示，CHG上的C由X或x表示。
对于样本j在区块R i 上片段集合{S 1 ,⋯,S N }中的每一条片段S n ，切除每个片段XM标签上，序列2起始端的前50个碱基的信息，以减少锯齿末端现象对区块R i 的平均甲基化水平beta i,j 和区块R i 的平均转化效率BS i,j 计算的影响。
根据样本j在区块R i 上片段集合{S 1 ,⋯,S N }，计算样本j在区块R i 的平均甲基化水平beta i,j ：
Figure 199365DEST_PATH_IMAGE043
和样本j在区块R i 的平均转化效率BS i,j ：
Figure 766481DEST_PATH_IMAGE044
其中Z n ，z n ，H n ，h n ，X n 和x n 分别代表片段S n 上Z，z，H，h，X和x的标签数量。
样本j在区块R i 的校正后平均甲基化水平beta_corrected i,j 可以根据beta i,j 与BS i,j 计算：
Figure 111137DEST_PATH_IMAGE045
对于每一样本j，最终输出的校正后区块平均甲基化水平为一个长度为I的向量：
Figure 752203DEST_PATH_IMAGE046
合并所有训练集样本的校正后区块平均甲基化水平至一个1,864行913列的矩阵，所有测试集样本的校正后区块平均甲基化水平至一个1,864行385列的矩阵，其中，行代表区块，列代表样本。两个矩阵分别输出为训练集校正后区块平均甲基化水平矩阵文件和测试集校正后区块平均甲基化水平矩阵文件。
2.3 分别生成训练集和测试集样本的校正后区块甲基化片段占比矩阵
将每个训练集内样本和测试集内样本的最终Bam文件和2.1中生成的区块列表文件作为区块甲基化水平计算与校正模块输入文件。
依次从每一个训练集或测试集样本j的Bam文件中，分别提取Bam中比对到每一个区块R i 的序列，并将来源于同一原始基因片段的两条配对序列的XM标签进行合并，作为片段XM标签。
对于样本j在区块R i 上片段集合{S 1 ,⋯,S N }中的每一条片段S n ，切除每个片段XM标签上，序列2起始端的前50个碱基的信息，以减少锯齿末端现象对片段S n 的甲基化水平beta i,j,n 、片段S n 的转化效率BS i,j 计算的影响。
根据去除碱基后片段XM标签上的CpG，计算片段S n 的甲基化水平beta i,j,n ：
Figure 530803DEST_PATH_IMAGE047
根据去除碱基后片段XM标签上的CHH和CHG标签，计算片段S n 的转化效率BS i,j,n ：
Figure 70369DEST_PATH_IMAGE048
其中Z n ，z n ，H n ，h n ，X n 和x n 分别代表片段S n 上Z，z，H，h，X和x的标签数量。
在计算样本j在区块R i 的甲基化片段占比MFR i,j 时，首先挑选出转化率高于阈值 （以0.95阈值为例，即满足
Figure 858196DEST_PATH_IMAGE004
）的片段，统计数量记录为N convertedi,j 。再根据N convertedi,j 条片段的beta i,j,n 是否高于阈值（以0.9阈值为例，即满足
Figure 166818DEST_PATH_IMAGE049
）判断每条 片段是否为甲基化片段。甲基化片段的数量记录为N meth i,j 。
最终，样本j在区块R i 的甲基化片段占比
Figure 864122DEST_PATH_IMAGE050
对于每一样本j，最终输出的校正后区块甲基化片段占比为一个长度为I的向量：
Figure 625404DEST_PATH_IMAGE051
合并所有训练集样本的校正后区块甲基化片段占比至一个1,864行913列的矩阵，所有测试集样本的校正后区块甲基化片段占比至一个1,864行385列的矩阵，其中行代表区块，列代表样本。两个矩阵分别输出为训练集校正后区块甲基化片段占比矩阵文件和测试集校正后甲基化片段占比矩阵文件。
2.4使用训练集对甲基化区块特征降维，获得旋转矩阵和训练集主成分矩阵
使用2.3中获得的训练集校正后区块甲基化片段占比矩阵文件，作为甲基化区块特征降维模块的输入文件。对训练集中913个样本1,864个区块的校正后甲基化片段占比进行主成分分析，共计产生913个主成分。该步骤同时生成了旋转矩阵V和一个913行913列的训练集主成分矩阵Y并生成训练集主成分文件。其中，行是训练集样本，列是主成分。
在获得旋转矩阵V和主成分矩阵Y后，计算每个主成分Y p 对总方差的贡献λ p ，并计算前m个主成分的累计贡献率，累计贡献率的计算公式为公式8。记录提取出累计贡献率超过95%的前66个主成分到主成分列表文件中。这66个主成分将用于甲基化模型的构建。
Figure 216923DEST_PATH_IMAGE052
2.5 获得降维后的测试集主成分矩阵
使用2.3中获得的测试集校正后区块甲基化片段占比矩阵文件和2.4中获得的旋转矩阵文件作为甲基化区块特征降维模块的输入文件。对测试集中385个样本1,864个区块的校正后甲基化片段占比进行转换，转换为913个主成分。该步骤同时生成了旋转矩阵V和一个385行913列的测试集主成分矩阵，并生成测试集主成分矩阵文件。其中，行是训练集样本，列是主成分。
2.6利用主成分进行甲基化模型构建
使用2.4中获得的训练集主成分矩阵文件和主成分列表文件作为甲基化模型构建模块的输入文件，同时提供训练集的分组信息列表。
首先提取主成分矩阵Y中所有训练集样本的前66个主成分到新的矩阵Y'中，将新矩阵Y'作为输入特征。训练集分组转化为二分类变量，健康人为0，患者为1，作为模型的应变量。
在建模时，采用装袋法，分别进行13次有放回的随机抽样，建立13个支持向量机弱分类器。对每个弱分类器单独进行参数寻优、训练和预测。最终的模型是一个基于13个弱分类器建立的集成模型M。模型M中，每个弱分类器的权重相等，最终的预测值是13个弱分类器预测值的平均值，代表样本中存在ctDNA的概率。
在训练模型的同时，对每个训练集样本进行预测，并记录和输出每个训练集样本的预测值到训练集预测结果文件中。
该步骤生成的支持向量机集成模型M会存储为一个二进制文件，用于后续对测试集样本的预测。
2.7用训练的模型对测试集样本进行预测
使用2.5中获得的测试集主成分矩阵文件和2.6中获得的集成模型M作为甲基化模型预测模块的输入文件。
首先提取训练集主成分矩阵中所有训练集样本的前66个主成分到新的矩阵中，并将新矩阵作为输入特征。支持向量机集成模型M中的13个弱分类器会分别依据66个主成分，对每个预测集样本进行预测，最终的预测值是13个弱分类器预测值的平均值，代表样本中存在ctDNA的概率。预测值被记录和输出到预测集预测结果文件中。
在训练集的特异性为95%的条件下，本发明对测试集中健康人的特异性达到97.24%，对单一癌种的检测灵敏度如表11所示，对全部8个癌种的检测灵敏度达到72.92%。
表11
Figure 380051DEST_PATH_IMAGE053
如图2-图10所示，本发明所述测试集中的单一癌种的检测ROC曲线面积（AUC）达到84.94%-100%，对全部8个癌种的检测AUC性能达到93.31%，显示出良好的癌症检测性能。
综上，本发明所述ctDNA检测装置能快速简便地对全基因组范围内的甲基化水平进行准确计算，并评估样本中存在ctDNA的概率。所述检测装置可应用于血浆样本，克服了目前低深度甲基化测序中甲基化水平计算不准确、灵敏度不高的技术缺陷，在肿瘤的早期筛查或实时监控中具有重要应用价值。
申请人声明，以上所述仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，所属技术领域的技术人员应该明了，任何属于本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，均落在本发明的保护范围和公开范围之内。

Claims (9)

1.一种基于全基因组甲基化测序的ctDNA预测方法，其特征在于，所述预测方法包括如下步骤：
根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，生成可计算甲基化水平的区块列表文件；
计算训练集样本和待测样本在每个区块上的甲基化水平，对训练集样本和待测样本甲基化水平进行降维，使用主成分分析法，根据训练集样本计算旋转矩阵，并利用旋转矩阵，生成训练集主成分矩阵和待测样本主成分矩阵；所述甲基化水平的计算方法包括：分别获取训练集样本和待测样本的捕获测序FASTQ文件，进行序列比对，分别生成训练集样本和待测样本原始Bam文件，对原始Bam文件进行排序和标记重复处理，分别生成训练集样本和待测样本排序去重后Bam文件，记为样本j去重排序后的Bam文件；将样本j去重排序后的Bam文件和所述可计算甲基化水平的区块列表作为输入文件，计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平，合并所有样本的校正后甲基化水平至一个I行J列的甲基化水平矩阵X；计算甲基化水平前还包括将样本j在区块R i 上片段集合{S 1 ,⋯,S N }中的每一条片段S n ，切除每个片段XM标签上，序列2起始端的前A个碱基的信息；
计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平的方法包括：方案一：根据非甲基化位点上胞嘧啶的状态，计算区块转化效率，并使用区块转化效率校正区块甲基化水平；或方案二：根据每条原始基因片段非甲基化位点上胞嘧啶的状态，计算片段转化效率，过滤低转化率片段后再计算区块的甲基化水平；
利用训练集主成分矩阵进行甲基化模型构建，用训练后的模型对待测样本中的ctDNA进行预测。
2.根据权利要求1所述的基于全基因组甲基化测序的ctDNA预测方法，其特征在于，所述将基因组划分为等长的区块采用如下步骤进行：
根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，根据全基因组范围内所有CpG位点的位置信息和多重比对区域的占比进行过滤，生成可计算甲基化水平的区块列表。
3.根据权利要求1所述的基于全基因组甲基化测序的ctDNA预测方法，其特征在于，所述训练集样本为正常人群血浆样本和患者的血浆样本；所述待测样本为被检测者的血浆样本。
4.根据权利要求1所述的基于全基因组甲基化测序的ctDNA预测方法，其特征在于，所述方案一中，甲基化水平的计算公式如下所示：
Figure DEST_PATH_IMAGE002
其中，beta_corrected i,j 为样本j在区块R i 的校正后平均甲基化水平；
beta i,j 为样本j在区块R i 的平均甲基化水平；
BS i,j 为样本j在区块R i 的平均转化效率；
所述方案二中，甲基化水平的计算公式如下所示：
Figure DEST_PATH_IMAGE004
Figure DEST_PATH_IMAGE006
其中，MFR i,j 为样本j在区块R i 的甲基化片段占比；
UFR i,j 为样本j在区块R i 的非甲基化片段占比；
N converted i,j 为
Figure DEST_PATH_IMAGE008
的片段的统计数量，BS i,j,n 为样本j在区块R i 上片段n的转化 效率；
N meth i,j 为甲基化片段的统计数量；
N unmeth i,j 为非甲基化片段的统计数量；
所述甲基化片段和非甲基化片段的判断标准为：
Figure DEST_PATH_IMAGE010
的片段判断为甲基化 片段，
Figure DEST_PATH_IMAGE012
的片段判断为非甲基化片段，beta i,j,n 为样本j在区块R i 上片段n的 平均甲基化水平，
Figure DEST_PATH_IMAGE014
。
5.根据权利要求1所述的基于全基因组甲基化测序的ctDNA预测方法，其特征在于，所述训练集主成分矩阵采用如下方式获取得到：
将训练集样本的甲基化水平矩阵X作为输入文件，进行主成分分析，输出对应的旋转矩 阵V和训练集主成分矩阵Y；计算每个主成分Y p 对总方差的贡献λ p ，并计算前m个主成分的累 计贡献率，提取出累计贡献率超过t的前m个主成分到主成分列表中，
Figure DEST_PATH_IMAGE016
；输出旋 转矩阵V、训练集主成分矩阵Y和主成分列表；
所述待测样本主成分矩阵采用如下方式获取得到：
将待测样本的甲基化水平矩阵X和所述旋转矩阵V作为输入文件，进行主成分分析，输出对应的待测样本主成分矩阵Y；
利用训练集主成分矩阵进行甲基化模型构建采用如下步骤进行：
首先根据所述主成分列表提取训练集主成分矩阵Y中所有样本的前m个主成分到新的矩阵Y'，将新矩阵Y'作为支持向量机模型的输入特征，采用装袋法进行建模，得到一个基于K个弱分类器建立的集成模型M，模型M中，每个弱分类器的权重相等，最终的预测值是K个弱分类器预测值的平均值；
对待测样本中的ctDNA进行预测采用如下步骤进行：
根据训练出的模型M中记录的模型变量，根据所述主成分列表提取待测样本主成分矩阵Y中所有待测样本的前m个主成分到新的矩阵Y'，并将Y'输入到训练出的模型M中，获得每个待测样本的预测值。
6.一种基于全基因组甲基化测序的ctDNA预测装置，其特征在于，所述预测装置包括：
甲基化区块划分模块：根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，生成可计算甲基化水平的区块列表文件；
区块甲基化水平计算与校正模块：用于计算训练集样本和待测样本在每个区块上的甲基化水平；所述甲基化水平的计算方法包括：分别获取训练集样本和待测样本的捕获测序FASTQ文件，进行序列比对，分别生成训练集样本和待测样本原始Bam文件，对原始Bam文件进行排序和标记重复处理，分别生成训练集样本和待测样本排序去重后Bam文件，记为样本j去重排序后的Bam文件；将样本j去重排序后的Bam文件和所述可计算甲基化水平的区块列表作为输入文件，计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平，合并所有样本的校正后甲基化水平至一个I行J列的甲基化水平矩阵X；计算甲基化水平前还包括将样本j在区块R i 上片段集合{S 1 ,⋯,S N }中的每一条片段S n ，切除每个片段XM标签上，序列2起始端的前A个碱基的信息；
计算样本j在列表文件中记录的I个等长的区块{R 1 ,⋯,R I }上的校正后甲基化水平的方法包括：方案一：根据非甲基化位点上胞嘧啶的状态，计算区块转化效率，并使用区块转化效率校正区块甲基化水平；或方案二：根据每条原始基因片段非甲基化位点上胞嘧啶的状态，计算片段转化效率，过滤低转化率片段后再计算区块的甲基化水平；
甲基化区块特征降维模块：对训练集样本和待测样本的甲基化水平进行降维，使用主成分分析法，根据训练集样本计算旋转矩阵，并利用旋转矩阵，生成训练集主成分矩阵和待测样本主成分矩阵；
甲基化模型构建模块：利用训练集主成分矩阵进行甲基化模型构建，对模型进行训练；
甲基化模型预测模块：用训练后的模型对待测样本中的ctDNA进行预测。
7.根据权利要求6所述的基于全基因组甲基化测序的ctDNA预测装置，其特征在于，所述甲基化区块划分模块中，所述将基因组划分为等长的区块采用如下步骤进行：
根据人类参考基因组FASTA文件和BED文件，将基因组划分为等长的区块，根据全基因组范围内所有CpG位点的位置信息和多重比对区域的占比进行过滤，生成可计算甲基化水平的区块列表；
所述甲基化区块特征降维模块中，所述训练集主成分矩阵采用如下方式获取得到：
将训练集样本的甲基化水平矩阵X作为输入文件，进行主成分分析，输出对应的旋转矩 阵V和训练集主成分矩阵Y；计算每个主成分Y p 对总方差的贡献λ p ，并计算前m个主成分的累 计贡献率，提取出累计贡献率超过t的前m个主成分到主成分列表中，
Figure 761379DEST_PATH_IMAGE016
；输出旋 转矩阵V、训练集主成分矩阵Y和主成分列表；
所述甲基化区块特征降维模块中，所述待测样本主成分矩阵采用如下方式获取得到：
将待测样本的甲基化水平矩阵X和所述旋转矩阵V作为输入文件，进行主成分分析，输出对应的待测样本主成分矩阵Y。
8.一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器运行所述计算机程序时实现权利要求1-5所述的基于全基因组甲基化测序的ctDNA预测方法的步骤。
9.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1-5所述的基于全基因组甲基化测序的ctDNA预测方法的步骤。
CN202210977623.9A 2022-08-15 2022-08-15 一种基于全基因组甲基化测序的ctDNA预测方法及装置 Active CN115064211B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210977623.9A CN115064211B (zh) 2022-08-15 2022-08-15 一种基于全基因组甲基化测序的ctDNA预测方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210977623.9A CN115064211B (zh) 2022-08-15 2022-08-15 一种基于全基因组甲基化测序的ctDNA预测方法及装置

Publications (2)

Publication Number Publication Date
CN115064211A CN115064211A (zh) 2022-09-16
CN115064211B true CN115064211B (zh) 2023-01-24

Family

ID=83208023

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210977623.9A Active CN115064211B (zh) 2022-08-15 2022-08-15 一种基于全基因组甲基化测序的ctDNA预测方法及装置

Country Status (1)

Country Link
CN (1) CN115064211B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115323058A (zh) * 2022-10-17 2022-11-11 深圳市睿法生物科技有限公司 一种基于ctDNA甲基化模式的癌种定位方法
CN115376616B (zh) * 2022-10-24 2023-04-28 臻和(北京)生物科技有限公司 一种基于cfDNA多组学的多分类方法及装置
CN115678964B (zh) * 2022-11-08 2023-07-14 广州女娲生命科技有限公司 基于胚胎培养液的植入前胚胎的无创筛选方法
CN116153418B (zh) * 2023-04-18 2023-07-18 臻和(北京)生物科技有限公司 校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质
CN117423388B (zh) * 2023-12-19 2024-03-22 北京求臻医疗器械有限公司 一种基于甲基化水平的多癌种检测系统及电子设备

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104899474A (zh) * 2015-06-09 2015-09-09 大连三生科技发展有限公司 基于岭回归矫正MB-seq甲基化水平的方法及系统
CN109852672A (zh) * 2017-11-30 2019-06-07 深圳豪石生物科技有限公司 一种筛选急性髓系白血病dna甲基化预后标志物的方法
CN110322928A (zh) * 2019-08-16 2019-10-11 河海大学常州校区 Dna甲基化谱检测方法
CN112397150A (zh) * 2021-01-20 2021-02-23 臻和(北京)生物科技有限公司 基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法
CN112397151A (zh) * 2021-01-21 2021-02-23 臻和(北京)生物科技有限公司 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置
CN112951418A (zh) * 2021-05-17 2021-06-11 臻和(北京)生物科技有限公司 基于液体活检的连锁区域甲基化评估方法和装置、终端设备及存储介质
CN113257350A (zh) * 2021-06-10 2021-08-13 臻和(北京)生物科技有限公司 基于液体活检的ctDNA突变程度分析方法和装置、ctDNA性能分析装置
CN113539355A (zh) * 2021-07-15 2021-10-22 云康信息科技(上海)有限公司 预测cfDNA的组织特异性来源及相关疾病概率评估系统及应用
CN113903401A (zh) * 2021-12-10 2022-01-07 臻和(北京)生物科技有限公司 基于ctDNA长度的分析方法和系统
CN114045345A (zh) * 2022-01-07 2022-02-15 臻和(北京)生物科技有限公司 基于游离dna的基因组癌变信息检测系统和检测方法
CN114703284A (zh) * 2022-04-15 2022-07-05 北京莱盟君泰国际医疗技术开发有限公司 一种血液游离dna甲基化定量检测方法及其应用
CN114868191A (zh) * 2019-10-11 2022-08-05 格瑞尔有限责任公司 利用起源组织阈值的癌症分类

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010085343A1 (en) * 2009-01-23 2010-07-29 Cold Spring Harbor Laboratory Methods and arrays for profiling dna methylation
NZ717423A (en) * 2012-09-20 2017-08-25 Univ Hong Kong Chinese Non-invasive determination of methylome of fetus or tumor from plasma
US10174375B2 (en) * 2013-09-20 2019-01-08 The Chinese University Of Hong Kong Sequencing analysis of circulating DNA to detect and monitor autoimmune diseases
EP3067432A1 (en) * 2015-03-11 2016-09-14 Deutsches Krebsforschungszentrum Stiftung des Öffentlichen Rechts DNA-methylation based method for classifying tumor species of the brain
CA3168463A1 (en) * 2015-12-17 2017-06-22 Illumina, Inc. Distinguishing methylation levels in complex biological samples
EP3494234A4 (en) * 2016-08-05 2020-03-04 The Broad Institute, Inc. METHOD FOR GENOM CHARACTERIZATION
AU2019253112A1 (en) * 2018-04-13 2020-10-29 Grail, Llc Multi-assay prediction model for cancer detection
GB201818159D0 (en) * 2018-11-07 2018-12-19 Cancer Research Tech Ltd Enhanced detection of target dna by fragment size analysis
EP4127231A1 (en) * 2020-03-31 2023-02-08 Grail, LLC Cancer classification with genomic region modeling
CN111755072B (zh) * 2020-08-04 2021-02-02 深圳吉因加医学检验实验室 一种同时检测甲基化水平、基因组变异和插入片段的方法及装置
EP4268231A1 (en) * 2021-01-20 2023-11-01 Genecast Biotechnology Co., Ltd Dna methylation sequencing analysis methods

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104899474A (zh) * 2015-06-09 2015-09-09 大连三生科技发展有限公司 基于岭回归矫正MB-seq甲基化水平的方法及系统
CN109852672A (zh) * 2017-11-30 2019-06-07 深圳豪石生物科技有限公司 一种筛选急性髓系白血病dna甲基化预后标志物的方法
CN110322928A (zh) * 2019-08-16 2019-10-11 河海大学常州校区 Dna甲基化谱检测方法
CN114868191A (zh) * 2019-10-11 2022-08-05 格瑞尔有限责任公司 利用起源组织阈值的癌症分类
CN112397150A (zh) * 2021-01-20 2021-02-23 臻和(北京)生物科技有限公司 基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法
CN112397151A (zh) * 2021-01-21 2021-02-23 臻和(北京)生物科技有限公司 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置
CN112951418A (zh) * 2021-05-17 2021-06-11 臻和(北京)生物科技有限公司 基于液体活检的连锁区域甲基化评估方法和装置、终端设备及存储介质
CN113257350A (zh) * 2021-06-10 2021-08-13 臻和(北京)生物科技有限公司 基于液体活检的ctDNA突变程度分析方法和装置、ctDNA性能分析装置
CN113539355A (zh) * 2021-07-15 2021-10-22 云康信息科技(上海)有限公司 预测cfDNA的组织特异性来源及相关疾病概率评估系统及应用
CN113903401A (zh) * 2021-12-10 2022-01-07 臻和(北京)生物科技有限公司 基于ctDNA长度的分析方法和系统
CN114045345A (zh) * 2022-01-07 2022-02-15 臻和(北京)生物科技有限公司 基于游离dna的基因组癌变信息检测系统和检测方法
CN114703284A (zh) * 2022-04-15 2022-07-05 北京莱盟君泰国际医疗技术开发有限公司 一种血液游离dna甲基化定量检测方法及其应用

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
人类基因组DNA甲基化数据分析的研究现状;凡时财等;《中国科学:生命科学》;20150520;第45卷(第05期);全文 *

Also Published As

Publication number Publication date
CN115064211A (zh) 2022-09-16

Similar Documents

Publication Publication Date Title
CN115064211B (zh) 一种基于全基因组甲基化测序的ctDNA预测方法及装置
US11371074B2 (en) Method and system for determining copy number variation
JP2020537487A (ja) メチローム解析を用いる癌の検出及び分類
CN110846411B (zh) 一种基于二代测序的单独肿瘤样本区分基因突变类型的方法
CN114045345B (zh) 基于游离dna的基因组癌变信息检测系统和检测方法
CN112397151B (zh) 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置
CN112218957A (zh) 用于确定在无细胞核酸中的肿瘤分数的系统及方法
CN112397150B (zh) 基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法
WO2020237184A1 (en) Systems and methods for determining whether a subject has a cancer condition using transfer learning
CN112941180A (zh) 一组肺癌dna甲基化分子标志物及其在制备用于肺癌早期诊断试剂盒中的应用
WO2022253288A1 (zh) 一种甲基化测序方法和装置
CN115803447A (zh) 染色体邻近实验中的结构变异检测
CN117524301B (zh) 一种拷贝数变异的检测方法、装置以及计算机可读介质
CN109712671B (zh) 基于ctDNA的基因检测装置、存储介质及计算机系统
CN110970091A (zh) 标签质控的方法及装置
CN109461473B (zh) 胎儿游离dna浓度获取方法和装置
CN110373458B (zh) 一种地中海贫血检测的试剂盒及分析系统
KR102347463B1 (ko) 핵산 서열 분석에서 위양성 변이를 검출하는 방법 및 장치
CN112837748A (zh) 一种区分不同解剖学起源肿瘤的系统及其方法
CN115831355A (zh) 多癌种wgs的肿瘤早期筛查方法
CN110305945A (zh) 一种基于二代测序技术的游离线粒体dna突变检测技术
CN109979534B (zh) 一种c位点提取方法及装置
US20200216888A1 (en) Method for increasing accuracy of analysis by removing primer sequence in amplicon-based next-generation sequencing
CN110684830A (zh) 一种石蜡切片组织rna分析方法
CN117637020B (zh) 一种基于深度学习的四倍体牡蛎全基因组snp分型方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant