CN112397151B - 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置 - Google Patents

基于靶向捕获测序的甲基化标志物筛选与评价方法及装置 Download PDF

Info

Publication number
CN112397151B
CN112397151B CN202110078570.2A CN202110078570A CN112397151B CN 112397151 B CN112397151 B CN 112397151B CN 202110078570 A CN202110078570 A CN 202110078570A CN 112397151 B CN112397151 B CN 112397151B
Authority
CN
China
Prior art keywords
methylation
matrix
sample
level
site
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
CN202110078570.2A
Other languages
English (en)
Other versions
CN112397151A (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.)
Wuxi Zhenhe Biotechnology Co.,Ltd.
Zhenhe (Beijing) Biotechnology Co.,Ltd.
Original Assignee
Wuxi Zhenhe Biotechnology 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 Wuxi Zhenhe Biotechnology Co ltd, Zhenhe Beijing Biotechnology Co ltd filed Critical Wuxi Zhenhe Biotechnology Co ltd
Priority to CN202110078570.2A priority Critical patent/CN112397151B/zh
Publication of CN112397151A publication Critical patent/CN112397151A/zh
Application granted granted Critical
Publication of CN112397151B publication Critical patent/CN112397151B/zh
Priority to EP21920475.7A priority patent/EP4268231A1/en
Priority to PCT/CN2021/091761 priority patent/WO2022156089A1/en
Priority to US17/490,549 priority patent/US20220228209A1/en
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/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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
    • 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
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding

Abstract

本发明提供了一种基于靶向捕获测序的甲基化标志物筛选与评价方法及装置，方法包括：分别获取N个待测样本捕获测序的FASTQ文件，并生成Bam文件；计算每个甲基化位点的甲基化水平和覆盖深度，合并得到甲基化水平矩阵和位点深度矩阵；针对每个甲基化位点，计算其与下一甲基化位点之间的距离及线性相关系数，并根据结果合并得到甲基化连锁区域；计算连锁区域甲基化水平均值矩阵和位点深度均值矩阵，筛选出与正常人群组存在设定差异的特异连锁区域；根据得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据甲基化得分对甲基化标志物进行评价。经过本发明筛选与评价的标志物能有效发现血浆中的ctDNA甲基化信号，获得更高的灵敏度。

Description

基于靶向捕获测序的甲基化标志物筛选与评价方法及装置
技术领域
本发明涉及生物医学技术领域，尤其涉及一种甲基化标志物筛选与评价方法及装置。
背景技术
循环肿瘤DNA（circμLating tumor ，ctDNA）由肿瘤细胞因分泌、凋亡或坏死而产生，是循环游离DNA（circμLating cell-free DNA，cfDNA）中的一种。ctDNA在血液中的半衰期短，且携带有部分肿瘤细胞特有的特征，可用于肿瘤患者的早期筛查或实时监控。除单核苷酸多态性（Single Nucleotide Polymorphisms，SNP）、插入缺失标记（insertion-deletion，InDel）与拷贝数变异（copy number variation，CNV）外，甲基化作为基因表达调控中的重要环节，同样能够影响基因组的稳定性。对于一些特定位点或区域的甲基化状况，肿瘤患者的ctDNA和健康人的cfDNA之间会存在明显差异，以此，通过从血浆中检测这些特定位点或区域的甲基化状况，能够在肿瘤发生的早期识别出血浆中ctDNA的存在，为后续癌症的早期诊断或复发预测提供数据基础。
近年来甲基化测序的应用一定程度上提高了ctDNA的检测灵敏度，但是这些技术大多将检测样本限定在粪便、痰液等与肿瘤的发生器官高度相关的样本中，且只能检测某一特定类型的ctDNA进行检测。目前，广泛使用的甲基化测序技术大多为聚合酶链式反应（PCR）、全基因组甲基化测序（WGBS）或靶向捕获测序。其中，PCR虽然成本较低，但能检测的甲基化位点数量较为有限，进而影响检测的灵敏度和特异性。WGBS虽然覆盖位点更全，但成本高、深度低，不利于从血浆样本中发现ctDNA的甲基化信号。由于各测序方式都需要在测序前对DNA进行重亚硫酸盐转化，甲基化位点的甲基化水平计算的准确性会受到转化效率的影响，不利于甲基化标志物的筛选，同时也影响对样本进行标志物检验时的准确性。
发明内容
针对上述问题，本发明提供了一种基于靶向捕获测序的甲基化标志物筛选与评价方法及装置，有效解决现有甲基化测序中依从性差、适用范围窄、灵敏度不高的技术问题。
本发明提供的技术方案如下：
一方面，本发明提供了一种基于靶向捕获测序的甲基化标志物筛选与评价方法，包括：
分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件，所述待测样本为血浆样本；
依次计算各待测样本Bam文件目标区域内
Figure 100002_DEST_PATH_IMAGE001
个甲基化位点上的甲基化水平和覆盖 深度，并合并得到甲基化水平矩阵和位点深度矩阵；
基于所述甲基化水平矩阵和位点深度矩阵，针对待测样本中的每一个甲基化位 点，分别计算其与下一甲基化位点之间的距离及甲基化水平的线性相关系数，并根据计算 结果依次对相邻甲基化位点进行合并得到甲基化连锁区域，将包含有预设数量甲基化位点 的
Figure 100002_DEST_PATH_IMAGE002
个甲基化连锁区域作为甲基化标志物候选区域输出；
依次计算各待测样本中
Figure 920299DEST_PATH_IMAGE002
个甲基化连锁区域内的甲基化水平均值和位点深度均 值，生成连锁区域甲基化水平均值矩阵和位点深度均值矩阵；
根据所述连锁区域甲基化水平均值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组存在设定差异的特异连锁区域，得到甲基化标志物；
根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据所述甲基化得分对甲基化标志物进行评价。
进一步优选地，在所述分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件之后，还包括：
根据预先设定的C-T转化率对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件；和/或，
根据目标区域Bed文件及预先设定的各reads中覆盖甲基化位点的数量对Bam文件进行过滤，得到过滤后的Bam文件。
进一步优选地，在所述分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件之后，还包括：
根据预先设定的C-T转化率对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件；和/或，
根据目标区域Bed文件及预先设定的各reads中覆盖甲基化位点的数量对Bam文件进行过滤，得到过滤后的Bam文件。
进一步优选地，在所述依次计算各待测样本Bam文件目标区域内
Figure 100002_DEST_PATH_IMAGE003
个甲基化位点上 的甲基化水平和覆盖深度，并合并到甲基化水平矩阵和位点深度矩阵中，包括：
依次提取各待测样本Bam文件中各个甲基化位点的正链信息和负链信息；
依次计算各待测样本中各个甲基化位点的甲基化水平和覆盖深度，其中，样本n在 甲基化位点
Figure 100002_DEST_PATH_IMAGE004
上的甲基化水平
Figure 100002_DEST_PATH_IMAGE005
和覆盖深度
Figure 100002_DEST_PATH_IMAGE006
分别为：
Figure 100002_DEST_PATH_IMAGE007
其中，
Figure 100002_DEST_PATH_IMAGE008
表示样本n在甲基化位点
Figure 192622DEST_PATH_IMAGE004
上的正链甲基化水平，
Figure 100002_DEST_PATH_IMAGE009
表示样本n在甲基化位点
Figure 447628DEST_PATH_IMAGE004
上的负链甲基化水平，
Figure 100002_DEST_PATH_IMAGE010
表示样本n 在甲基化位点
Figure 802255DEST_PATH_IMAGE004
上的正链覆盖深度，
Figure 100002_DEST_PATH_IMAGE011
表示样本n在甲基化位点
Figure 472795DEST_PATH_IMAGE004
上的负 链覆盖深度，
Figure 100002_DEST_PATH_IMAGE012
，
Figure 100002_DEST_PATH_IMAGE013
；
将计算得到的各甲基化位点的甲基化水平和覆盖深度进行合并得到甲基化水平 矩阵
Figure 100002_DEST_PATH_IMAGE014
和位点深度矩阵
Figure 100002_DEST_PATH_IMAGE015
，所述甲基化水平矩阵
Figure 100002_DEST_PATH_IMAGE016
和位点深度矩 阵
Figure 100002_DEST_PATH_IMAGE017
均为I行N列矩阵，其中，行对应甲基化位点，列对应待测样本。
进一步优选地，在所述基于甲基化水平矩阵和位点深度矩阵，针对待测样本中的 每一个甲基化位点，分别计算其与下一甲基化位点之间的距离及甲基化水平的线性相关系 数，并根据计算结果依次对相邻甲基化位点进行合并得到甲基化连锁区域，将包含有预设 数量甲基化位点的
Figure 100002_DEST_PATH_IMAGE018
个甲基化连锁区域作为甲基化标志物候选区域输出中，针对样本n甲基 化位点
Figure 216103DEST_PATH_IMAGE004
的甲基化连锁区域合并步骤包括：
判断甲基化位点
Figure 447102DEST_PATH_IMAGE004
及其下一甲基化位点
Figure 100002_DEST_PATH_IMAGE019
的覆盖深度是否均在预设最低深度 要求
Figure 100002_DEST_PATH_IMAGE020
之上；
若是，计算甲基化位点
Figure 886261DEST_PATH_IMAGE004
和甲基化位点
Figure 911724DEST_PATH_IMAGE019
之间的距离
Figure 100002_DEST_PATH_IMAGE021
；
计算甲基化位点
Figure 721872DEST_PATH_IMAGE004
上各待测样本的甲基化水平
Figure 100002_DEST_PATH_IMAGE022
与甲基化位点
Figure 634857DEST_PATH_IMAGE019
上各待测样本的甲基化水平
Figure 100002_DEST_PATH_IMAGE023
之间的线性相关系数
Figure 100002_DEST_PATH_IMAGE024
；其 中，
Figure 100002_DEST_PATH_IMAGE025
表示样本1在甲基化位点
Figure 554622DEST_PATH_IMAGE004
上的甲基化水平，
Figure 100002_DEST_PATH_IMAGE026
表示样本N在甲基化位点
Figure 891931DEST_PATH_IMAGE004
上的甲基化水平，
Figure 100002_DEST_PATH_IMAGE027
表示样本1在甲基化位点
Figure 173613DEST_PATH_IMAGE019
上的甲基化水平，
Figure 100002_DEST_PATH_IMAGE028
表 示样本N在甲基化位点
Figure 437760DEST_PATH_IMAGE019
上的甲基化水平；
判断是否同时满足条件
Figure 100002_DEST_PATH_IMAGE029
和
Figure 100002_DEST_PATH_IMAGE030
，其中，
Figure 100002_DEST_PATH_IMAGE031
为预设最大 位点间间距，
Figure 100002_DEST_PATH_IMAGE032
为预设最小相关系数；
若是，将甲基化位点
Figure 100002_DEST_PATH_IMAGE033
并入甲基化位点
Figure 475947DEST_PATH_IMAGE004
当前所在的甲基化连锁区域
Figure 100002_DEST_PATH_IMAGE034
形成新 的甲基化连锁区域
Figure 46475DEST_PATH_IMAGE034
，否则断开甲基化位点
Figure 429832DEST_PATH_IMAGE004
当前所在的甲基化连锁区域
Figure 407323DEST_PATH_IMAGE034
，
Figure 100002_DEST_PATH_IMAGE035
。
进一步优选地，所述依次计算各待测样本中
Figure 494490DEST_PATH_IMAGE018
个甲基化连锁区域内的甲基化水平 均值和位点深度均值，生成连锁区域甲基化水平均值矩阵和位点深度均值矩阵中：
样本
Figure 100002_DEST_PATH_IMAGE036
于甲基化连锁区域
Figure 100002_DEST_PATH_IMAGE037
中的甲基化水平均值
Figure 100002_DEST_PATH_IMAGE038
为：
Figure 100002_DEST_PATH_IMAGE039
样本
Figure 522005DEST_PATH_IMAGE036
于甲基化连锁区域
Figure 183538DEST_PATH_IMAGE037
中的位点深度均值
Figure 100002_DEST_PATH_IMAGE040
为：
Figure 100002_DEST_PATH_IMAGE041
在根据计算的甲基化水平均值
Figure 100002_DEST_PATH_IMAGE042
和位点深度均值
Figure 100002_DEST_PATH_IMAGE043
形成的连锁区域甲 基化水平均值矩阵
Figure 100002_DEST_PATH_IMAGE044
和位点深度均值矩阵
Figure 100002_DEST_PATH_IMAGE045
中，行对应甲基化连 锁区域，列对应样本。
进一步优选地，所述根据所述连锁区域甲基化水平均值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组存在设定差异的特异连锁区域中，包括：
将所述甲基化水平均值矩阵
Figure 100002_DEST_PATH_IMAGE046
进行
Figure 100002_DEST_PATH_IMAGE047
转换得到矩阵
Figure 100002_DEST_PATH_IMAGE048
；
根据预先设定的样本分组信息表将矩阵
Figure 100002_DEST_PATH_IMAGE049
、甲基化水平均值矩阵
Figure 100002_DEST_PATH_IMAGE050
和位点深度均值矩阵
Figure 100002_DEST_PATH_IMAGE051
分别按照control组和case组进行分割， 得到矩阵
Figure 100002_DEST_PATH_IMAGE052
、矩阵
Figure 100002_DEST_PATH_IMAGE053
、矩阵
Figure 100002_DEST_PATH_IMAGE054
、 矩阵
Figure 100002_DEST_PATH_IMAGE055
、矩阵
Figure 100002_DEST_PATH_IMAGE056
和矩阵
Figure 100002_DEST_PATH_IMAGE057
6个独立矩 阵，其中，control组中矩阵
Figure 100002_DEST_PATH_IMAGE058
、矩阵
Figure 100002_DEST_PATH_IMAGE059
和矩阵
Figure 100002_DEST_PATH_IMAGE060
包含正常人群组的样本数据，case组中矩阵
Figure 100002_DEST_PATH_IMAGE061
、矩 阵
Figure 100002_DEST_PATH_IMAGE062
和矩阵
Figure 22046DEST_PATH_IMAGE057
包含异常人群组的样本数据；
遍历所有甲基化连锁区域，筛选出与正常人群组存在设定差异的特异连锁区域， 其中，针对甲基化连锁区域
Figure 291616DEST_PATH_IMAGE037
的筛选过程为：
依次计算control组与case组的错误发现率
Figure 100002_DEST_PATH_IMAGE063
、甲基化水平差值
Figure 100002_DEST_PATH_IMAGE064
、case组 差异样本占比
Figure 100002_DEST_PATH_IMAGE065
及control组中的低甲基化基线样本占比
Figure 100002_DEST_PATH_IMAGE066
，其中，所述case组差异样本占 比
Figure 100002_DEST_PATH_IMAGE067
表示case组
Figure 100002_DEST_PATH_IMAGE068
转换后甲基化水平在control组一倍标准差范围外的样本数量占 case组样本总数量的比值，所述control组中的低甲基化基线样本占比
Figure 139966DEST_PATH_IMAGE066
表示control组中 甲基化水平低于预设低甲基化阈值
Figure 100002_DEST_PATH_IMAGE069
的样本数量占control组样本总数量的比值，
Figure 100002_DEST_PATH_IMAGE070
；
判断是否同时满足条件
Figure 100002_DEST_PATH_IMAGE071
、
Figure 100002_DEST_PATH_IMAGE072
、
Figure 100002_DEST_PATH_IMAGE073
及
Figure 100002_DEST_PATH_IMAGE074
，其 中，
Figure 100002_DEST_PATH_IMAGE075
为预设最大错误发现率，
Figure 100002_DEST_PATH_IMAGE076
为预设最小甲基化水平差值，
Figure 100002_DEST_PATH_IMAGE077
为预设差异 样本占比阈值，
Figure 100002_DEST_PATH_IMAGE078
为预设低甲基化基线样本占比阈值；
若是，判断甲基化连锁区域
Figure 100002_DEST_PATH_IMAGE079
与正常人存在设定差异。
进一步优选地，所述control组与case组的错误发现率
Figure 100002_DEST_PATH_IMAGE080
由control组所有待测 样本
Figure 100002_DEST_PATH_IMAGE081
转换后的甲基化水平
Figure 100002_DEST_PATH_IMAGE082
与case组所有待测样 本
Figure 709418DEST_PATH_IMAGE081
转换后的甲基化水平
Figure 100002_DEST_PATH_IMAGE083
经过校正t检验得到的差异检验 值
Figure 100002_DEST_PATH_IMAGE084
进一步经过Benjaminiand Hochberg法校正后得到，其中，
Figure 100002_DEST_PATH_IMAGE085
表示control组 样本1在甲基化连锁区域
Figure 100002_DEST_PATH_IMAGE086
中
Figure 100002_DEST_PATH_IMAGE087
转换后的甲基化水平，
Figure 100002_DEST_PATH_IMAGE088
表示control 组样本N在甲基化连锁区域
Figure 49218DEST_PATH_IMAGE086
中
Figure 100002_DEST_PATH_IMAGE089
转换后的甲基化水平；
Figure 100002_DEST_PATH_IMAGE090
表示case组样本1在 甲基化连锁区域
Figure 100002_DEST_PATH_IMAGE091
中
Figure 235611DEST_PATH_IMAGE081
转换后的甲基化水平，
Figure 100002_DEST_PATH_IMAGE092
表示case组样本N在甲基 化连锁区域
Figure 100002_DEST_PATH_IMAGE093
中
Figure 100002_DEST_PATH_IMAGE094
转换后的甲基化水平；
和/或，所述甲基化水平差值
Figure 100002_DEST_PATH_IMAGE095
表示control组平均甲基化水平与case组平均 甲基化水平之差：
Figure 100002_DEST_PATH_IMAGE096
其中，
Figure 100002_DEST_PATH_IMAGE097
表示case组甲基化连锁区域
Figure 189616DEST_PATH_IMAGE093
中样本1的平均甲基化水平，
Figure 100002_DEST_PATH_IMAGE098
表示case组甲基化连锁区域
Figure 982997DEST_PATH_IMAGE086
中样本N的平均甲基化水平，
Figure 100002_DEST_PATH_IMAGE099
表示control组甲基化连锁区域
Figure 242334DEST_PATH_IMAGE086
中样本1的平均甲基化水平，
Figure 100002_DEST_PATH_IMAGE100
表示 control组甲基化连锁区域
Figure 713592DEST_PATH_IMAGE086
中样本N的平均甲基化水平；
和/或，所述case组差异样本占比
Figure 100002_DEST_PATH_IMAGE101
为：
Figure 100002_DEST_PATH_IMAGE102
其中，
Figure 100002_DEST_PATH_IMAGE103
表示case组甲基化连锁区域
Figure 282239DEST_PATH_IMAGE086
中
Figure 100002_DEST_PATH_IMAGE104
转换后样本
Figure 100002_DEST_PATH_IMAGE105
的甲基化水平，
Figure 100002_DEST_PATH_IMAGE106
表示control组甲基化连锁区域
Figure 52266DEST_PATH_IMAGE086
中
Figure 100002_DEST_PATH_IMAGE107
转 换后的平均甲基化水平，
Figure 100002_DEST_PATH_IMAGE108
表示control组甲基化连锁区域
Figure 336530DEST_PATH_IMAGE086
中
Figure 497516DEST_PATH_IMAGE107
转换后 的标准差，
Figure 100002_DEST_PATH_IMAGE109
表示case组样本总数；
和/或，control组中的低甲基化基线样本占比
Figure 100002_DEST_PATH_IMAGE110
为：
Figure 100002_DEST_PATH_IMAGE111
其中，
Figure 100002_DEST_PATH_IMAGE112
表示control组甲基化连锁区域
Figure 736865DEST_PATH_IMAGE086
中样本
Figure 100002_DEST_PATH_IMAGE113
的甲基化水平，
Figure 100002_DEST_PATH_IMAGE114
表示control组样本总数，
Figure 100002_DEST_PATH_IMAGE115
表示预设 的甲基化水平背景噪音最大值。
进一步优选地，在所述根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据所述甲基化得分对甲基化标志物进行评价中，包括：
将所述甲基化水平均值矩阵
Figure 100002_DEST_PATH_IMAGE116
进行
Figure 100002_DEST_PATH_IMAGE117
转换得到矩阵
Figure 100002_DEST_PATH_IMAGE118
；
从所述矩阵
Figure 100002_DEST_PATH_IMAGE119
和位点深度均值矩阵
Figure 100002_DEST_PATH_IMAGE120
中提取筛选得到的
Figure 100002_DEST_PATH_IMAGE121
个特异连锁区域的数据，并根据预先设定的样本分组信息表将其分割为矩阵
Figure 100002_DEST_PATH_IMAGE122
、矩阵
Figure 100002_DEST_PATH_IMAGE123
、矩阵
Figure 100002_DEST_PATH_IMAGE124
和矩阵
Figure 100002_DEST_PATH_IMAGE125
，其中，矩阵
Figure 100002_DEST_PATH_IMAGE126
和矩阵
Figure 100002_DEST_PATH_IMAGE127
包 含正常人群组样本数据，矩阵
Figure 100002_DEST_PATH_IMAGE128
和矩阵
Figure 100002_DEST_PATH_IMAGE129
包含待测样本 数据；
根据分割得到的矩阵分别计算各待测样本的甲基化得分，并判断是否存在甲基化得分大于预设分数阈值的待测样本；若是，判定该待测样本中包含筛选获得的甲基化标志物；
其中，待测样本
Figure 100002_DEST_PATH_IMAGE130
甲基化得分
Figure 100002_DEST_PATH_IMAGE131
为：
Figure 100002_DEST_PATH_IMAGE132
其中，
Figure DEST_PATH_IMAGE133
表示待测样本
Figure DEST_PATH_IMAGE134
在特异连锁区域
Figure DEST_PATH_IMAGE135
的平均深度，
Figure DEST_PATH_IMAGE136
；
Figure DEST_PATH_IMAGE137
表示待测样本
Figure DEST_PATH_IMAGE138
在特异连锁区域
Figure DEST_PATH_IMAGE139
的p-value值，所述p- value值为待测样本
Figure 446718DEST_PATH_IMAGE134
在特异连锁区域
Figure DEST_PATH_IMAGE140
上
Figure DEST_PATH_IMAGE141
转换后甲基化水平
Figure DEST_PATH_IMAGE142
的 Z-score值
Figure DEST_PATH_IMAGE143
转换为标准正态分布的分位数：
Figure DEST_PATH_IMAGE144
其中，
Figure DEST_PATH_IMAGE145
为正常人群组
Figure 203451DEST_PATH_IMAGE141
转换后的甲基化水平均值，
Figure DEST_PATH_IMAGE146
为正常人群组
Figure 910288DEST_PATH_IMAGE107
转换后的方差。
进一步优选地，在所述根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据所述甲基化得分对甲基化标志物进行评价中，还包括：根据每个待测样本的已知分组情况计算检出的灵敏度和特异性，或根据待测样本的已知ctDNA浓度，计算甲基化得分与ctDNA浓度的线性相关系数，进而根据灵敏度、特异性与线性相关系数的高低对筛选的甲基化标志物进行评价的步骤。
另一方面，本发明还提供了一种基于靶向捕获测序的甲基化标志物筛选与评价装置，应用于上述甲基化标志物筛选与评价方法，所述甲基化标志物筛选与评价装置中包括：
Bam文件生成模块，用于分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件，所述待测样本为血浆样本；
位点甲基化水平提取模块，用于依次计算各待测样本Bam文件目标区域内
Figure DEST_PATH_IMAGE147
个甲基 化位点上的甲基化水平和覆盖深度，并合并得到甲基化水平矩阵和位点深度矩阵；
甲基化连锁区域合并模块，用于基于所述甲基化水平矩阵和位点深度矩阵，针对 待测样本中的每一个甲基化位点，分别计算其与下一甲基化位点之间的距离及甲基化水平 的线性相关系数，并根据计算结果依次对相邻甲基化位点进行合并得到甲基化连锁区域， 将包含有预设数量甲基化位点的
Figure 559751DEST_PATH_IMAGE002
个甲基化连锁区域作为甲基化标志物候选区域输出；
区域甲基化平均水平提取模块，用于依次计算各待测样本中
Figure 598377DEST_PATH_IMAGE002
个甲基化连锁区域 内的甲基化水平均值和位点深度均值，生成连锁区域甲基化水平均值矩阵和位点深度均值 矩阵；
差异区域筛选模块，用于根据所述连锁区域甲基化水平均值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组存在设定差异的特异连锁区域，得到甲基化标志物；
基线构建与得分计算模块，用于根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据所述甲基化得分对甲基化标志物进行评价。
另一方面，本发明提供了一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器运行所述计算机程序时实现上述基于靶向捕获测序的甲基化标志物筛选与评价方法的步骤。
另一方面，本发明提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现上述任一项所述基于靶向捕获测序的甲基化标志物筛选与评价方法的步骤。
本发明提供的基于靶向捕获测序的甲基化标志物筛选与评价装置及方法中，至少能够带来以下有益效果：
1.相较于甲基化测序中常用的聚合酶链式反应（PCR）或全基因组甲基化测序（WGBS）来说，本发明针对的高深度捕获测序，能够更好地平衡成本、深度和覆盖度的问题，更有利于发现血浆中的ctDNA甲基化信号，获得更高的灵敏度，且检测结果可靠易于解读。
2.基于普遍认为的相邻甲基化位点之间存在高相关性的连锁关系，本发明根据得到的甲基化水平矩阵和位点深度矩阵，对相邻甲基化位点进行合并得到甲基化连锁区域，以此用连锁区域的甲基化水平替代单个位点的甲基化水平，在计算连锁区域甲基化水平时综合考虑多个甲基化位点的甲基化情况，从而提升甲基化水平计算的稳定性。尤其对于待测血浆样本中的低甲基化水平区域，实验或测序中产生的误差对这些区域的甲基化水平计算影响很大，甲基化水平计算的稳定性对于甲基化标志物筛选的准确性至关重要，具体表现在两方面：一，优化低甲基化区域的甲基化水平计算，减少实验和测序过程中引入的误差；二，使用连锁区域特征作为输入进行建模时，将多个存在高相关性的甲基化位点特征合并为一个区域特征，从而减少模型的变量数，去除冗余变量，优化模型的表现。
3.由于血浆中cfDNA的释放来源十分复杂，肿瘤患者血浆中的片段，有一大部分的来源可能与健康人相类似。扣除这部分的甲基化背景信号，对于特定甲基化标志物的检测至关重要。以此，本发明在对甲基化连锁区域进行合并之后，进一步对差异性大的区域进行筛选，保留高甲基化区域，且尽量保证筛选出的甲基化标志物在健康人中的信号尽可能地少，从而保证在后续分析时，与病灶相关的ctDNA信号不会被来源于健康组织的cfDNA信号淹没。
4.在对甲基化得分进行计算中，建立了正常人群甲基化水平基线数据，并对待测样本与基线之间的差异进行量化分析，而不是直接使用待测样本的甲基化水平进行分析，以此待测样本与参考人群组基线水平的偏离程度能够被更加准确地衡量，用于构建基线的样本越多，优势越明显，准确率越高，使用该得分进行甲基化标志物评价时也越准确。
5.在构建的用于评价标志物优劣的得分计算模型中，训练阶段仅使用健康人血浆样本即可，无需患者样本。且建立的模型不同于常见的逻辑回归、随机森林、支持向量机等机器学习方法，仅对待测样本与基线样本在不同甲基化连锁区域上的甲基化水平差异进行简单的加权平均，权重为待测样本在区域的平均测序深度。该模型使用的计算方法更简便，能够克服过拟合问题，降低对训练模型所需的患者样本数量的要求，且同时支持使用的甲基化标志物数量也比常见的机器学习方法更多。
附图说明
下面将以明确易懂的方式，结合附图说明优选实施方式，对上述特性、技术特征、优点及其实现方式予以进一步说明。
图1为本发明中基于靶向捕获测序的甲基化标志物筛选与评价方法流程示意图；
图2为本发明中基于靶向捕获测序的甲基化标志物筛选与评价装置结构示意图；
图3为本发明一实例中甲基化连锁区域划分流程图；
图4为本发明一实例中对甲基化连锁区域进行差异性筛选流程图；
图5为本发明一实例中梯度稀释血浆样本得分计算流程图；
图6为本发明一实例中每个原始样本在不同稀释梯度对应的甲基化得分变化图；
图7为本发明一实例中12个梯度稀释样本理论ctDNA占比和甲基化得分的散点图；
图8为本发明中终端设备结构示意图。
附图标记：
100-甲基化标志物筛选与评价装置，110-Bam文件生成模块，120-位点甲基化水平提取模块，130-甲基化连锁区域合并模块，140-区域甲基化平均水平提取模块，150-差异区域筛选模块，160-基线构建与得分计算模块。
具体实施方式
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对照附图说明本发明的具体实施方式。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，并获得其他的实施方式。
如图1所示为本发明提供的基于靶向捕获测序的甲基化标志物筛选与评价方法流程示意图，从图中可以看出，该甲基化标志物筛选与评价方法中包括：
S10 分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件，待测样本为血浆样本；
S20 依次计算各待测样本Bam文件目标区域内
Figure 46282DEST_PATH_IMAGE147
个甲基化位点上的甲基化水平和 覆盖深度，并合并得到甲基化水平矩阵和位点深度矩阵；
S30 基于甲基化水平矩阵和位点深度矩阵，针对待测样本中的每一个甲基化位 点，分别计算其与下一甲基化位点之间的距离及甲基化水平的线性相关系数，并根据计算 结果依次对相邻甲基化位点进行合并得到甲基化连锁区域，将包含有预设数量甲基化位点 的
Figure 856893DEST_PATH_IMAGE002
个甲基化连锁区域作为甲基化标志物候选区域输出；
S40 依次计算各待测样本中
Figure 872997DEST_PATH_IMAGE002
个甲基化连锁区域内的甲基化水平均值和位点深度 均值，生成连锁区域甲基化水平均值矩阵和位点深度均值矩阵；
S50 根据连锁区域甲基化水平均值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组存在设定差异（显著差异）的特异连锁区域，得到甲基化标志物；
S60 根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据甲基化得分对甲基化标志物进行评价。
从甲基化水平的计算角度考虑，各类测序技术都需要在测序前对DNA进行重亚硫酸盐转化，转化不完全会导致位点甲基化状态的误判。常用的甲基化水平衡量指标为甲基化位点的beta值，即位点上甲基化reads占总reads的比例。考虑到早期肿瘤患者的血浆中，ctDNA占比通常较低，特异的甲基化标志物特征信号较为微弱，转化效率会极大地影响beta值的计算。是以，在一实施例中，步骤S10，分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件之后，还包括：根据预先设定的C-T转化率（原始序列非CpG位点的C碱基转化为T碱基的比例）对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件步骤；和/或，根据目标区域Bed文件及预先设定的各reads中覆盖甲基化位点的数量对Bam文件进行过滤，得到过滤后的Bam文件，以提升后续的筛选准确率。
在根据C-T转化率进行过滤中，包括：获取待测样本捕获测序的的FASTQ文件之后，利用去接头软件Trimmomatic去除接头和低质量reads，得到过滤后的FASTQ文件，并利用FASTQC（一种高通量测序数据的质控软件，用于数据质量的评估）软件对待测样本的数据量、碱基质量分布、碱基含量比例进行统计分析。之后，利用基因组比对工具Bismark（一种比对方法软件，用于查找测序序列在基因参考序列中的位置，可输出Bam格式结果文件）将FASTQ文件中的基因序列分别与人类参考基因组和内参lambda DNA参考基因组进行比对并去重，生成人类参考基因组的Bam文件、去重前的比对报告和去重后的比对报告，及内参lambda DNA参考基因组Bam文件、去重前的比对报告和去重后的比对报告；并利用SAMtools和Picard工具对比对后的Bam文件进行排序和标记重复处理。接着，逐行读取Bam文件中的reads，根据Bam文件每条read中原始序列为C碱基的位点的实际碱基情况，对其non-CpGcontext模式下甲基化和非甲基化的碱基数量进行统计；并根据发生甲基化的non-CpGcontext的碱基数及non-CpG context碱基数总和（甲基化和非甲基化的碱基数量之和）对每条reads的C-T转化率进行计算；最后，将Bam文件中C-T转化率（如90%、95%等）小于预先设定的C-T转化率的reads滤除，以此过滤掉不符合non-CpG context C-T转化率最低要求的reads，输出过滤后的Bam格式文件、过滤后待测样本C-T转化率及过滤后待测样本的reads数据量。
更具体来说，FASTQ为一种常见的高通量测序文件类型。reads为测序读长，测序仪测到的基因组或转录组序列片段。根据甲基化的C碱基所处的上下文环境，分为三种类型CpG、CHG和CHH，H代表除了G碱基之外的其他碱基，即A、C、T中的任意一种；CpG为甲基化的C的下游是1个G碱基，CHG代表甲基化的C下游的2个碱基是H和G，CHH表示甲基化的C下游的两个碱基都是H，CHG和CHH可以合称为non CpG context。Bam文件用来存储测序序列回贴到参考基因组的结果。
在对reads中覆盖甲基化位点数量进行过滤中，包括：利用BisSNP软件（一种分析甲基化数据的软件，可用于鉴定甲基化位点和预测甲基化水平）根据目标区域Bed文件对dbSNP数据库中已知的SNP位点及因特定变异原因（如，结构变异、染色体拷贝数变异等）产生的SNP位点进行过滤得到待测样本的CpG位点之后，根据过滤得到的CpG位点及预先设定的各reads中覆盖CpG位点的数量（即上述每条reads上覆盖CpG位点的最低要求）对Bam文件进行过滤，将不满足覆盖CpG位点数量最低要求的reads滤除。
由于基因组上部分区域的甲基化位点，倾向于同时发生甲基化或不发生甲基化，可称为甲基化位点间的连锁关系，存在连锁关系的甲基化位点所在的区域，可被定义为甲基化连锁区域。基于甲基化的该种特性，本实施例根据待测样本相邻甲基化位点的甲基化水平值的相关性及距离，将目标区域上存在高相关性的相邻甲基化位点进行合并得到甲基化连锁区域。
在这一过程中，步骤S20 依次计算各待测样本Bam文件目标区域内I个甲基化位点上的甲基化水平和覆盖深度，并合并到甲基化水平矩阵和位点深度矩阵中，包括：
S21 依次提取各待测样本Bam文件中各个甲基化位点的正链信息和负链信息。这里，N个待测样本对应生成有N个Bam文件，每个待测样本中目标区域中包括I个甲基化位点，以此这里使用BisSNP结果文件（记录位点对应的甲基化水平）列表作为输入文件，分别从样本1到样本N的BisSNP结果文件中，提取出样本在I个甲基化位点上的正链信息和负链信息。
S22 依次计算各待测样本中各个甲基化位点的甲基化水平（覆盖位点的测序序列 中，发生甲基化的位点所占比例）和覆盖深度，其中，样本n在甲基化位点
Figure 398919DEST_PATH_IMAGE004
上的甲基化水平
Figure 355959DEST_PATH_IMAGE005
和覆盖深度
Figure DEST_PATH_IMAGE148
分别如式(1)和式(2)：
Figure DEST_PATH_IMAGE149
(1)
Figure DEST_PATH_IMAGE150
(2)
其中，
Figure DEST_PATH_IMAGE151
表示样本n在甲基化位点
Figure 893513DEST_PATH_IMAGE004
上的正链甲基化水平，
Figure 199293DEST_PATH_IMAGE009
表示样本n在甲基化位点
Figure 147264DEST_PATH_IMAGE004
上的负链甲基化水平，
Figure 2219DEST_PATH_IMAGE010
表示样本n 在甲基化位点
Figure 849739DEST_PATH_IMAGE004
上的正链覆盖深度，
Figure DEST_PATH_IMAGE152
表示样本n在甲基化位点
Figure 896061DEST_PATH_IMAGE004
上的负 链覆盖深度，
Figure 145644DEST_PATH_IMAGE012
，
Figure 817672DEST_PATH_IMAGE013
；
S23 将计算得到的各甲基化位点的甲基化水平和覆盖深度进行合并得到甲基化 水平矩阵
Figure DEST_PATH_IMAGE153
和位点深度矩阵
Figure DEST_PATH_IMAGE154
，甲基化水平矩阵
Figure DEST_PATH_IMAGE155
和位点深度矩 阵
Figure DEST_PATH_IMAGE156
均为I行N列矩阵，其中，行对应甲基化位点，列对应待测样本。
得到甲基化水平矩阵
Figure DEST_PATH_IMAGE157
和位点深度矩阵
Figure DEST_PATH_IMAGE158
之后，进入步骤S30对 甲基化连锁区域进行合并的步骤。具体，在S30 在基于甲基化水平矩阵和位点深度矩阵，针 对待测样本中的每一个甲基化位点，分别计算其与下一甲基化位点之间的距离及甲基化水 平的线性相关系数，并根据计算结果依次对相邻甲基化位点进行合并得到甲基化连锁区 域，将包含有预设数量甲基化位点的
Figure 901122DEST_PATH_IMAGE018
个甲基化连锁区域作为甲基化标志物候选区域输出 中，针对样本n甲基化位点
Figure 560423DEST_PATH_IMAGE004
的甲基化连锁区域合并步骤包括：
S31 判断甲基化位点
Figure 279724DEST_PATH_IMAGE004
及其下一甲基化位点
Figure DEST_PATH_IMAGE159
的覆盖深度（
Figure DEST_PATH_IMAGE160
和
Figure DEST_PATH_IMAGE161
）是否均在预设最低深度要求
Figure DEST_PATH_IMAGE162
之上；
S32 若是，计算甲基化位点
Figure 805644DEST_PATH_IMAGE004
和甲基化位点
Figure 562510DEST_PATH_IMAGE159
之间的距离
Figure 318677DEST_PATH_IMAGE021
；若覆盖深度
Figure DEST_PATH_IMAGE163
或
Figure DEST_PATH_IMAGE164
低于预设最低深度要求
Figure 905035DEST_PATH_IMAGE162
，则不将样本n计入线性相 关系数的计算，以确保数据的可靠性。
S33 计算甲基化位点
Figure 840236DEST_PATH_IMAGE004
上各待测样本的甲基化水平
Figure 513925DEST_PATH_IMAGE022
与甲基化 位点
Figure 526748DEST_PATH_IMAGE159
上各待测样本的甲基化水平
Figure DEST_PATH_IMAGE165
之间的线性相关系数
Figure DEST_PATH_IMAGE166
； 其中，
Figure 538085DEST_PATH_IMAGE025
表示样本1在甲基化位点
Figure 214660DEST_PATH_IMAGE004
上的甲基化水平，
Figure 609433DEST_PATH_IMAGE026
表示样本N在甲基化位 点
Figure 650945DEST_PATH_IMAGE004
上的甲基化水平，
Figure 569485DEST_PATH_IMAGE027
表示样本1在甲基化位点
Figure 117928DEST_PATH_IMAGE159
上的甲基化水平，
Figure 497700DEST_PATH_IMAGE028
表示样本N在甲基化位点
Figure 837677DEST_PATH_IMAGE159
上的甲基化水平；
S34 判断是否同时满足条件
Figure DEST_PATH_IMAGE167
和
Figure DEST_PATH_IMAGE168
，其中，
Figure 225538DEST_PATH_IMAGE031
为预设 最大位点间间距，
Figure 312092DEST_PATH_IMAGE032
为预设最小相关系数；
S35 若是，将甲基化位点
Figure DEST_PATH_IMAGE169
并入甲基化位点
Figure 172470DEST_PATH_IMAGE004
当前所在的甲基化连锁区域
Figure 753187DEST_PATH_IMAGE034
形 成新的甲基化连锁区域
Figure 705706DEST_PATH_IMAGE034
，否则断开甲基化位点
Figure 205739DEST_PATH_IMAGE004
当前所在的甲基化连锁区域
Figure 563033DEST_PATH_IMAGE034
，
Figure 445145DEST_PATH_IMAGE035
。
应当注意，在上述甲基化连锁区域合并过程中，若后续甲基化位点
Figure DEST_PATH_IMAGE170
和甲基化 位点
Figure DEST_PATH_IMAGE171
之间的距离
Figure DEST_PATH_IMAGE172
和线性相关系数
Figure DEST_PATH_IMAGE173
满足条件
Figure DEST_PATH_IMAGE174
和
Figure DEST_PATH_IMAGE175
，则进一步将甲基化位点
Figure 764522DEST_PATH_IMAGE170
并入甲基化位点
Figure 384423DEST_PATH_IMAGE171
当前所在的甲基化连锁区 域
Figure DEST_PATH_IMAGE176
形成新的甲基化连锁区域
Figure 973055DEST_PATH_IMAGE176
，否则断开甲基化连锁区域
Figure 760489DEST_PATH_IMAGE176
。后续满足条件的甲基化位 点合并入甲基化连锁区域
Figure DEST_PATH_IMAGE177
。完成所有甲基化位点的计算合并过程之后，输出包含有预 设数量（如包含3个和3个以上）甲基化位点的
Figure 268960DEST_PATH_IMAGE018
个甲基化连锁区域，并形成包含目标区间的 甲基化位点位置、相邻位点间的距离和线性相关系数及合并后位点所属连锁区域信息的特 异连锁区域列表文件。
甲基化连锁区域合并完成并输出后，在步骤S40 依次计算各待测样本中
Figure DEST_PATH_IMAGE178
个甲基 化连锁区域内的甲基化水平均值和位点深度均值，生成连锁区域甲基化水平均值矩阵和位 点深度均值矩阵，其中，
样本
Figure 124177DEST_PATH_IMAGE036
于甲基化连锁区域
Figure 518380DEST_PATH_IMAGE037
中的甲基化水平均值
Figure 804612DEST_PATH_IMAGE038
如式(3)：
Figure 956370DEST_PATH_IMAGE039
(3)
样本
Figure 918073DEST_PATH_IMAGE036
于甲基化连锁区域
Figure 101536DEST_PATH_IMAGE037
中的位点深度均值
Figure 296019DEST_PATH_IMAGE040
如式(4)：
Figure 764434DEST_PATH_IMAGE041
(4)
之后根据计算的甲基化水平均值
Figure 913785DEST_PATH_IMAGE042
和位点深度均值
Figure DEST_PATH_IMAGE179
形成的连锁区域 甲基化水平均值矩阵
Figure DEST_PATH_IMAGE180
和位点深度均值矩阵
Figure 571994DEST_PATH_IMAGE045
中，行对应甲基化 连锁区域，列对应样本，
Figure DEST_PATH_IMAGE181
。
之后对差异性模块进行筛选，在步骤S50 根据连锁区域甲基化水平均值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组存在设定差异的特异连锁区域中，包括：
S51 将所述甲基化水平均值矩阵
Figure 560548DEST_PATH_IMAGE046
进行
Figure 627511DEST_PATH_IMAGE047
转换得到矩阵
Figure DEST_PATH_IMAGE182
；
S52 根据预先设定的样本分组信息表将矩阵
Figure DEST_PATH_IMAGE183
、甲基化水平均值矩阵
Figure 203172DEST_PATH_IMAGE050
和位点深度均值矩阵
Figure DEST_PATH_IMAGE184
分别按照control组和case组进行分割， 得到矩阵
Figure 597114DEST_PATH_IMAGE052
、矩阵
Figure 152641DEST_PATH_IMAGE053
、矩阵
Figure DEST_PATH_IMAGE185
、 矩阵
Figure 318350DEST_PATH_IMAGE055
、矩阵
Figure 228144DEST_PATH_IMAGE056
和矩阵
Figure 181318DEST_PATH_IMAGE057
6个独立矩 阵，其中，control组中矩阵
Figure 297960DEST_PATH_IMAGE058
、矩阵
Figure 989448DEST_PATH_IMAGE059
和矩阵
Figure DEST_PATH_IMAGE186
包含正常人群组的样本数据，case组中矩阵
Figure 612189DEST_PATH_IMAGE061
、矩 阵
Figure 416940DEST_PATH_IMAGE062
和矩阵
Figure 560608DEST_PATH_IMAGE057
包含异常人群组的样本数据；
S53 遍历所有甲基化连锁区域，筛选出与正常人群组存在设定差异的特异连锁区 域，其中，针对甲基化连锁区域
Figure 89459DEST_PATH_IMAGE037
的筛选过程为：
依次计算control组与case组的错误发现率
Figure 154105DEST_PATH_IMAGE063
、甲基化水平差值
Figure 869820DEST_PATH_IMAGE064
、case组 差异样本占比
Figure 932192DEST_PATH_IMAGE067
及control组中的低甲基化基线样本占比
Figure DEST_PATH_IMAGE187
，其中，所述case组差异样本占 比
Figure 433492DEST_PATH_IMAGE067
表示case组
Figure DEST_PATH_IMAGE188
转换后甲基化水平在control组一倍标准差范围外的样本数量占 case组样本总数量的比值，所述control组中的低甲基化基线样本占比
Figure 114878DEST_PATH_IMAGE187
表示control组中 甲基化水平低于预设低甲基化阈值
Figure DEST_PATH_IMAGE189
的样本数量占control组样本总数量的比值，
Figure 659603DEST_PATH_IMAGE070
。具体：
control组与case组的错误发现率
Figure 145073DEST_PATH_IMAGE063
由R软件中的limma包生成。计算中，先对甲 基化连锁区域
Figure 445254DEST_PATH_IMAGE037
上control组所有待测样本
Figure DEST_PATH_IMAGE190
转换后的甲基化水平
Figure DEST_PATH_IMAGE191
与case组所有待测样本
Figure 282062DEST_PATH_IMAGE188
转换后的甲基化水平
Figure DEST_PATH_IMAGE192
进行校正t检验，得到甲基化连锁区域
Figure 482623DEST_PATH_IMAGE037
的差异检验值
Figure 762164DEST_PATH_IMAGE084
，其 中，
Figure DEST_PATH_IMAGE193
表示control组样本1在甲基化连锁区域
Figure DEST_PATH_IMAGE194
中
Figure 889256DEST_PATH_IMAGE117
转换后的甲基化水平，
Figure DEST_PATH_IMAGE195
表示control组样本N在甲基化连锁区域
Figure 853057DEST_PATH_IMAGE086
中
Figure 601134DEST_PATH_IMAGE089
转换后的甲基化水 平；
Figure 366090DEST_PATH_IMAGE090
表示case组样本1在甲基化连锁区域
Figure 369424DEST_PATH_IMAGE091
中
Figure 389114DEST_PATH_IMAGE081
转换后的甲基化水平，
Figure 997556DEST_PATH_IMAGE092
表示case组样本N在甲基化连锁区域
Figure 948062DEST_PATH_IMAGE093
中
Figure DEST_PATH_IMAGE196
转换后的甲基化水平。计算 出所有的
Figure 2475DEST_PATH_IMAGE018
个甲基化连锁区域的差异检验值
Figure DEST_PATH_IMAGE197
后，进一步对
Figure DEST_PATH_IMAGE198
使用 Benjamini-Hochberg法计算错误发现率，得到错误发现率
Figure 183445DEST_PATH_IMAGE063
。
甲基化水平差值
Figure DEST_PATH_IMAGE199
表示control组平均甲基化水平与case组平均甲基化水平 之差如式(5)：
Figure DEST_PATH_IMAGE200
(5)
其中，
Figure 765168DEST_PATH_IMAGE097
表示case组甲基化连锁区域
Figure DEST_PATH_IMAGE201
中样本1的平均甲基化水平，
Figure DEST_PATH_IMAGE202
表示case组甲基化连锁区域
Figure 138772DEST_PATH_IMAGE086
中样本N的平均甲基化水平，
Figure DEST_PATH_IMAGE203
表示control组甲基化连锁区域
Figure 883743DEST_PATH_IMAGE086
中样本1的平均甲基化水平，
Figure DEST_PATH_IMAGE204
表示control组甲基化连锁区域
Figure 789776DEST_PATH_IMAGE086
中样本N的平均甲基化水平；
case组差异样本占比
Figure 438057DEST_PATH_IMAGE101
如式(6)：
Figure DEST_PATH_IMAGE205
(6)
其中，
Figure DEST_PATH_IMAGE206
表示case组甲基化连锁区域
Figure 162387DEST_PATH_IMAGE086
中
Figure 630540DEST_PATH_IMAGE104
转换后样本
Figure 648918DEST_PATH_IMAGE105
的甲基化水平，
Figure DEST_PATH_IMAGE207
表示control组甲基化连锁区域
Figure 854770DEST_PATH_IMAGE086
中
Figure 238085DEST_PATH_IMAGE107
转 换后的平均甲基化水平，
Figure DEST_PATH_IMAGE208
表示control组甲基化连锁区域
Figure 632681DEST_PATH_IMAGE086
中
Figure 448891DEST_PATH_IMAGE107
转换后 的标准差，
Figure DEST_PATH_IMAGE209
表示case组样本总数；
control组中的低甲基化基线样本占比
Figure 120700DEST_PATH_IMAGE110
如式(7)：
Figure DEST_PATH_IMAGE210
(7)
其中，其中，
Figure DEST_PATH_IMAGE211
表示control组甲基化连锁区域
Figure 66395DEST_PATH_IMAGE086
中样本
Figure 446823DEST_PATH_IMAGE113
的甲基化水平，
Figure 134900DEST_PATH_IMAGE114
表示control组样本总数，
Figure 349630DEST_PATH_IMAGE115
表示预设 的甲基化水平背景噪音最大值（低于该值的样本被认为在该区域的甲基化水平低）。
S54 判断是否同时满足条件
Figure 74747DEST_PATH_IMAGE071
、
Figure DEST_PATH_IMAGE212
、
Figure 559781DEST_PATH_IMAGE073
及
Figure 661336DEST_PATH_IMAGE074
， 其中，
Figure 730573DEST_PATH_IMAGE075
为预设最大错误发现率，
Figure 957417DEST_PATH_IMAGE076
为预设最小甲基化水平差值，
Figure 106245DEST_PATH_IMAGE077
为预设差 异样本占比阈值，
Figure 779252DEST_PATH_IMAGE078
为预设低甲基化基线样本占比阈值；
S55 若是，判断甲基化连锁区域
Figure 27962DEST_PATH_IMAGE079
与正常人存在设定差异。
另外，在步骤S53 遍历所有甲基化连锁区域，筛选出与正常人群组存在设定差异 的特异连锁区域的过程中，为进一步保证指标计算的可靠性，若某一样本
Figure 485095DEST_PATH_IMAGE036
在甲基化连锁 区域
Figure 64762DEST_PATH_IMAGE079
的甲基化水平均值
Figure DEST_PATH_IMAGE213
低于预先设定的深度阈值
Figure DEST_PATH_IMAGE214
，则样本
Figure 499330DEST_PATH_IMAGE036
不会被用 于甲基化连锁区域
Figure 730110DEST_PATH_IMAGE079
各项指标的计算。
最后进入基线构建和甲基化得分的步骤，在步骤S60 根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据甲基化得分对甲基化标志物进行评价中，包括：
S61 将所述甲基化水平均值矩阵
Figure DEST_PATH_IMAGE215
进行
Figure 874257DEST_PATH_IMAGE117
转换得到矩阵
Figure DEST_PATH_IMAGE216
；
S62 从所述矩阵
Figure 567319DEST_PATH_IMAGE119
和位点深度均值矩阵
Figure DEST_PATH_IMAGE217
中提取筛选得到的
Figure DEST_PATH_IMAGE218
个特异连锁区域的数据，并根据预先设定的样本分组信息表将其分割为矩阵
Figure 182409DEST_PATH_IMAGE122
、矩阵
Figure 74885DEST_PATH_IMAGE123
、矩阵
Figure DEST_PATH_IMAGE219
和矩阵
Figure DEST_PATH_IMAGE220
，其中，矩阵
Figure DEST_PATH_IMAGE221
和矩阵
Figure 458110DEST_PATH_IMAGE127
包 含正常人群组样本数据，矩阵
Figure 68827DEST_PATH_IMAGE128
和矩阵
Figure DEST_PATH_IMAGE222
包含待测样本 数据。
S63 根据分割得到的矩阵分别计算各待测样本的甲基化得分，并判断是否存在甲 基化得分大于预设分数阈值的待测样本；若是，判定该待测样本中包含筛选获得的甲基化 标志物，其中，待测样本
Figure 610579DEST_PATH_IMAGE130
（为了对公式8和9表述清晰，这里用下标
Figure 623141DEST_PATH_IMAGE130
表示待测样 本n）甲基化得分
Figure 321539DEST_PATH_IMAGE131
如式(8)：
Figure DEST_PATH_IMAGE223
(8)
其中，
Figure 373546DEST_PATH_IMAGE133
表示待测样本
Figure 89960DEST_PATH_IMAGE134
在特异连锁区域
Figure DEST_PATH_IMAGE224
的平均深度，
Figure 523740DEST_PATH_IMAGE136
；
Figure 557423DEST_PATH_IMAGE137
表示待测样本
Figure 814835DEST_PATH_IMAGE138
在特异连锁区域
Figure DEST_PATH_IMAGE225
的p-value值，所述p- value值为待测样本
Figure 698390DEST_PATH_IMAGE134
在特异连锁区域
Figure DEST_PATH_IMAGE226
上
Figure 373959DEST_PATH_IMAGE141
转换后甲基化水平
Figure 555105DEST_PATH_IMAGE142
的 Z-score值
Figure 362131DEST_PATH_IMAGE143
转换为标准正态分布的分位数如式(9)：
Figure DEST_PATH_IMAGE227
(9)
其中，
Figure DEST_PATH_IMAGE228
为正常人群组
Figure 734862DEST_PATH_IMAGE141
转换后的甲基化水平均值，
Figure DEST_PATH_IMAGE229
为正常人群组
Figure 282516DEST_PATH_IMAGE107
转换后的方差。
预设分数阈值可以为基线样本得分
Figure DEST_PATH_IMAGE230
的最大值
Figure DEST_PATH_IMAGE231
或95%分位数
Figure DEST_PATH_IMAGE232
，基线样本得分
Figure DEST_PATH_IMAGE233
由正常健康人群计算得到， 计算步骤与待测样本相同。得到所有待测样本
Figure DEST_PATH_IMAGE234
的甲基化得分
Figure DEST_PATH_IMAGE235
之后，可以根据每个待测样本的已知分组情况计算检出的灵敏度和 特异性，或根据待测样本的已知ctDNA浓度，计算甲基化得分与ctDNA浓度的线性相关系数， 进而根据灵敏度、特异性与线性相关系数的高低对筛选的甲基化标志物进行评价。
相对应地，本发明还提供了一种基于靶向捕获测序的甲基化标志物筛选与评价装 置，应用于上述甲基化标志物筛选与评价方法，如图2所示，该甲基化标志物筛选与评价装 置100中包括：Bam文件生成模块110，用于分别获取N个待测样本捕获测序的FASTQ文件，并 分别与参考基因组进行比对生成Bam文件，待测样本为血浆样本；位点甲基化水平提取模块 120，用于依次计算各待测样本Bam文件目标区域内
Figure 713192DEST_PATH_IMAGE147
个甲基化位点上的甲基化水平和覆盖 深度，并合并得到甲基化水平矩阵和位点深度矩阵；甲基化连锁区域合并模块130，基于甲 基化水平矩阵和位点深度矩阵，针对待测样本中的每一个甲基化位点，分别计算其与下一 甲基化位点之间的距离及甲基化水平的线性相关系数，并根据计算结果依次对相邻甲基化 位点进行合并得到甲基化连锁区域，将包含有预设数量甲基化位点的
Figure 13373DEST_PATH_IMAGE002
个甲基化连锁区域 作为甲基化标志物候选区域输出；区域甲基化平均水平提取模块140，用于依次计算各待测 样本中
Figure 547385DEST_PATH_IMAGE002
个甲基化连锁区域内的甲基化水平均值和位点深度均值，生成连锁区域甲基化水 平均值矩阵和位点深度均值矩阵；差异区域筛选模块150，用于根据连锁区域甲基化水平均 值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组 存在设定差异的特异连锁区域，得到甲基化标志物；基线构建与得分计算模块160，用于根 据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据甲基化得分对甲基 化标志物进行评价。
在一实施例中，Bam文件生成模块110中还用于：根据预先设定的C-T转化率对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件步骤；和/或，根据目标区域Bed文件及预先设定的各reads中覆盖甲基化位点的数量对Bam文件进行过滤，得到过滤后的Bam文件，以提升后续的筛选准确率。
由于基因组上部分区域的甲基化位点，倾向于同时发生甲基化或不发生甲基化，可称为甲基化位点间的连锁关系，存在连锁关系的甲基化位点所在的区域，可被定义为甲基化连锁区域。基于甲基化的该种特性，本实施例根据待测样本相邻甲基化位点的甲基化水平值的相关性及距离，将目标区域上存在高相关性的相邻甲基化位点进行合并得到甲基化连锁区域。
在位点甲基化水平提取模块120中，将N个待测样本甲基化水平Bed文件中的数据 进行合并生成甲基化水平矩阵和位点深度矩阵。该模块中，使用BisSNP结果文件（记录位点 对应的甲基化水平）列表作为输入文件，分别从样本1到样本N的BisSNP结果文件中，提取出 样本在I个甲基化位点上的正链信息和负链信息，并根据式(1)和式(2)计算各样本中所有 甲基化位点的甲基化水平矩阵和位点深度矩阵合并得到甲基化水平矩阵
Figure DEST_PATH_IMAGE236
和位点 深度矩阵
Figure DEST_PATH_IMAGE237
。
在甲基化连锁区域合并模块130中，输入文件包括位点甲基化水平提取模块120输 出的甲基化水平矩阵
Figure DEST_PATH_IMAGE238
、位点深度矩阵
Figure DEST_PATH_IMAGE239
及用于合并甲基化连锁区域的样 本信息表，同时要求输入预先配置的最低深度要求
Figure DEST_PATH_IMAGE240
、最大位点间距
Figure DEST_PATH_IMAGE241
和最小相 关性
Figure DEST_PATH_IMAGE242
三个参数。对于待测样本n中的甲基化位点
Figure 629173DEST_PATH_IMAGE004
，采用步骤S31~S35对其是否合并入 同一甲基化连锁区域进行判断，根据该步骤完成待测样本中所有甲基化位点的判断之后， 输出包含有预设数量（如包含3个和3个以上）甲基化位点的
Figure DEST_PATH_IMAGE243
个甲基化连锁区域，并形成包 含目标区间的甲基化位点位置、相邻位点间的距离和线性相关系数及合并后位点所属连锁 区域信息的甲基化连锁区域列表文件。
在待测样本为肿瘤组织样本的实例中，结合位点甲基化水平提取模块120和甲基化连锁区域合并模块130，甲基化连锁区域划分过程如图3所示。甲基化连锁区域划分开始之后，位点甲基化水平提取模块120根据肿瘤组织样本BisSNP输出的bed文件生成肿瘤组织位点beta值矩阵和肿瘤组织位点深度矩阵；之后，甲基化连锁区域合并模块130根据样本信息列表对甲基化连锁区域进行合并得到甲基化连锁区域列表，完成甲基化连锁区域的划分。
在区域甲基化平均水平提取模块140中，输入文件为甲基化连锁区域合并模块130 输出的甲基化连锁区域列表文件、位点甲基化水平提取模块120输出的甲基化水平矩阵
Figure DEST_PATH_IMAGE244
和位点深度矩阵
Figure DEST_PATH_IMAGE245
。根据式(3)和式(4)分别对各待测样本中每个甲基化 连锁区域计算甲基化水平均值和位点深度均值后，形成连锁区域甲基化水平均值矩阵
Figure DEST_PATH_IMAGE246
和位点深度均值矩阵
Figure DEST_PATH_IMAGE247
。
在差异区域筛选模块150中，输入文件包括记录连锁区域甲基化水平均值矩阵
Figure DEST_PATH_IMAGE248
的文件、记录位点深度均值矩阵
Figure DEST_PATH_IMAGE249
的文件和用于筛选差异甲基化区 域的样本分组信息表，同时输入预先配置的最低平均深度要求
Figure DEST_PATH_IMAGE250
、最大错误发现率
Figure DEST_PATH_IMAGE251
、最小甲基化水平差值
Figure DEST_PATH_IMAGE252
、差异样本占比阈值
Figure 53855DEST_PATH_IMAGE077
、低甲基化阈值
Figure DEST_PATH_IMAGE253
与低 甲基化基线样本占比阈值
Figure DEST_PATH_IMAGE254
，并根据步骤S51~S55对与正常人存在设定差异的甲基化连 锁区域
Figure DEST_PATH_IMAGE255
进行筛选，并将筛选后得到的特异连锁区域输出。
在待测样本为肿瘤组织样本的实例中，结合位点甲基化水平提取模块120和差异区域筛选模块150，对甲基化连锁区域进行差异性筛选的过程如图4所示。肿瘤差异连锁区域筛选开始之后，位点甲基化水平提取模块120根据肿瘤组织样本BisSNP输出bed文件和基线血浆样本BisSNP输出bed文件生成组织与基线血浆位点beta值矩阵及组织与基线血浆位点深度矩阵，并根据甲基化连锁区域列表进一步生成组织与基线血浆连锁区域平均beta值矩阵及组织与基线血浆连锁区域平均深度矩阵；之后，差异区域筛选模块150根据基线血浆与肿瘤组织样本分组信息表对甲基化连锁区域进行筛选得到肿瘤特异连锁区域，并形成肿瘤特异连锁区域列表，完成肿瘤特异连锁区域的筛选。
在基线构建与得分计算模块160中，输入文件包括记录连锁区域甲基化水平均值 矩阵
Figure DEST_PATH_IMAGE256
的文件、记录位点深度均值矩阵
Figure DEST_PATH_IMAGE257
的文件、差异区域筛选模块 150筛选出来的甲基化区域列表及用于进行得分计算的样本分组信息表作为输入文件，并 根据步骤S61~S63判断该待测样本中是否包含ctDNA的甲基化信号。且得到所有待测样本
Figure DEST_PATH_IMAGE258
的甲基化得分
Figure DEST_PATH_IMAGE259
之后，可以根据每个待测样本的已知分组 情况计算检出的灵敏度和特异性，或根据待测样本的已知ctDNA浓度，计算甲基化得分与 ctDNA浓度的线性相关系数，进而根据灵敏度、特异性与线性相关系数的高低对筛选的甲基 化标志物进行评价。
在一实例中，选取健康人血浆50例，泛癌种肿瘤组织FFPE 166例（其中肺癌肿瘤11例）以及灵敏度测试血浆4例分别进行如下操作：
1.血浆样本的准备
1.1 血浆样本融化后，每1mL样本中加入15μL蛋白酶K（Proteinase K）(20mg/mL)和50 μL十二烷基硫酸钠（SDS）溶液(20%)。当血浆量不足4mL时，用磷酸缓冲盐（PBS）溶液补足。
1.2 翻转混匀，60℃孵育20min，然后冰浴5min。
1.3 向深孔板内加入如表1示的试剂。
表1：深孔板中加入的试剂列表
Figure DEST_PATH_IMAGE260
1.4 运行KingFisher FLEX磁珠提取仪。
程序运行前需将干净磁头套放入检测程序指定位置，运行程序，检测磁头套是否会掉落。深孔板加好后，点击自动提取仪上SATRT键，按照显示屏要求依次放入磁头套和对应深孔板。再次点击SATRT键，自动提取仪开始运行。程序时长为49min左右。
1.5 吸出DNA样品：
自动提取仪运行结束后，先取出7号深孔板，然后点击STOP键。用移液器将DNA样本吸出至对应的贴标标签的离心管中。
石蜡包埋组织（FFPE）样本的准备
与血浆样本准备步骤不同的是，FFPE样本准备中基因组DNA样本需要单管打断，打断后使用Qubit荧光定量仪测浓度和质检，其余步骤一致，这里不做赘述。
梯度稀释样本的准备
血浆样本的DNA提取步骤与1中一致。
将完成提取后的4例血浆样本cfDNA样本分别按照1/27、1/81、1/243的比例稀释于健康人的血浆样本中。
内参准备
取Lamdba DNA加入50μL打断管中，使用M220打断仪打断，将打断的内参DNA稀释，建库时加入样本中。Lamdba为参考品，用于确定样本的转化情况。
文库制备
5.1使用EZ DNA Methylation-LightningTM试剂盒（Zymo Research公司生产）对DNA进行转化
样品起始体积为20μL，不足20μL时，用水补足。取130μL试剂盒中的LightningConversion Reagent加入DNA样本中，震荡混匀，短暂离心，置于PCR仪上，并按表2进行PCR反应。
表2：PCR反应的条件
Figure DEST_PATH_IMAGE261
向试剂盒中的Zymo-Spin™ IC Column中加入600μL试剂盒中的M-BindingBuffer，将上步骤反应后的产物加入含有M-Binding Buffer的Zymo-Spin™ IC Column中，用枪吹打混匀，静置2min。12000rpm离心1min。
将收集管中的液体重新加回吸附柱中，静置2min，12000rpm离心1min，弃废液。
加入100μL试剂盒中的M-Wash Buffer，12000rpm离心1min，弃废液。
加入200μL试剂盒中的L-DesμLphonation Buffer室温(20-30°C)孵育15-20min，孵育完成后，12000rpm离心1min，弃废液。
加入200μL试剂盒中的M-Wash Buffer，12000rpm离心1min，弃废液，重复两次。
将吸附柱放回收集管中，12,000 rpm离心2 min，倒掉废液。将吸附柱开盖置于室温放置2-5 min，以彻底晾干吸附材料中残余的漂洗液。
将吸附柱转入一个干净的离心管中，向吸附膜的中间部位悬空滴加20μL洗脱缓冲液TE洗脱，室温放置2-5min，12000 rpm离心1 min。
将收集管中的液体重新加回吸附柱中，室温放置2-5min，12000 rpm离心1 min，将收集有转化后DNA的离心管-20℃保存。
5.2 DNA预处理
PCR仪提前95℃预热，热盖温度105℃。
取转化后的片段化DNA放入0.2ml的PCR管中，加入低浓度乙二胺四乙酸TE缓冲液（Low EDTA TE）稀释总体积到15μL。
将PCR管放入PCR仪中，进行95℃孵育2min后，立即放置到冰上，静置2min。
5.3加T7接头
PCR仪提前37℃预热，热盖温度105℃。
按照表3配置反应体系，表格中的试剂均来自ACCEL-NGS® METHYL-SEQ DNALIBRARY KIT试剂盒（Swift Biosciences公司生产）。
表3：反应试剂列表
Figure DEST_PATH_IMAGE262
加25μL试剂到冰上放置的预处理DNA样本PCR管中，使用移液器进行吹打混匀，瞬时离心。
将PCR管置于PCR仪中，进行反应，条件如表4所示。
表4：反应条件
Figure DEST_PATH_IMAGE263
5.4二链合成反应（Second strand synthesis reaction）
PCR仪提前98℃预热，热盖温度105℃。
按照表5配置反应试剂，表格中的试剂均来自剂来自ACCEL-NGS® METHYL-SEQDNA LIBRARY KIT试剂盒（Swift Biosciences公司生产）。
表5：反应试剂列表
Figure DEST_PATH_IMAGE264
加44μL表5试剂到上一步反应体系中，使用移液器进行吹打混匀，瞬时离心。
将PCR管置于PCR仪中，进行二链合成反应，反应条件如表6所示。
表6：二链合成反应条件
Figure DEST_PATH_IMAGE265
提前将纯化磁珠从4℃取出，室温平衡半小时。
待上一步反应结束后，在产物中加入101μL磁珠，吹打混匀。
室温静置5min，置于磁力架上至液体澄清，弃去上清。
加入200μL 80％乙醇孵育30sec后弃去。注意：80%乙醇现用现配。重复一次200μL80％乙醇清洗步骤。
用10μL枪头弃去离心管底部的残留乙醇，室温干燥至乙醇完全挥发。
从磁力架取下离心管，加入16μL超纯水，振荡混匀。室温孵育2min。
短暂离心，置于磁力架上至液体澄清，将15μL样本转入新的离心管中。
5.5加T5接头
按照表7配置反应试剂，表格中的试剂均来自ACCEL-NGS® METHYL-SEQ DNALIBRARY KIT试剂盒（Swift Biosciences公司生产）。加15μL反应体系到上一步的样本中，使用移液器进行吹打混匀，瞬时离心。
表7：反应试剂列表
Figure DEST_PATH_IMAGE266
将PCR管置于PCR仪中，按表8的条件进行PCR反应。
表8：PCR反应的条件
Figure DEST_PATH_IMAGE267
提前将纯化磁珠从4℃取出，室温平衡半小时。
连接反应结束后，加入36μL磁珠，吹打混匀。
室温静置5min，置于磁力架上至液体澄清，弃去上清。
加入200μL 80％乙醇孵育30sec后弃去。注意：80%乙醇现用现配。重复一次200μL80％乙醇清洗步骤。
用10μL枪头弃去离心管底部的残留乙醇，室温干燥至乙醇完全挥发。
从磁力架取下离心管，加入20μL超纯水，振荡混匀。室温孵育2min。
短暂离心，置于磁力架上至液体澄清，将20μL样本转入新的离心管中。
5.6 扩增
按照表9配置反应试剂，加30μL反应体系到上一步的样本中，使用移液器进行吹打混匀，瞬时离心。表格中的试剂来自ACCEL-NGS® METHYL-SEQ DNA LIBRARY KIT试剂盒（Swift Biosciences公司生产）。
表9：反应试剂列表
Figure DEST_PATH_IMAGE268
将PCR管置于PCR仪中，按表10的条件进行PCR反应。
表10：PCR反应的条件
Figure DEST_PATH_IMAGE269
提前将纯化磁珠从4℃取出，室温平衡半小时。
连接反应结束后，加入60μL磁珠，吹打混匀。
室温静置5min，置于磁力架上至液体澄清，弃去上清。
加入200μL 80％乙醇孵育30sec后弃去。注意：80%乙醇现用现配。重复一次200μL80％乙醇清洗步骤。
用10μL枪头弃去离心管底部的残留乙醇，室温干燥至乙醇完全挥发。
从磁力架取下离心管，加入50μL超纯水，振荡混匀。室温孵育2min。
短暂离心，置于磁力架上至液体澄清，将50μL样本转入新的离心管中。
文库捕获
6.1 混合文库：
按每个捕获总量1ug捕获。向上述体系中加入杂交试剂，震荡混匀，短暂离心。
用封口膜封住EP管，放入真空离心浓缩仪中蒸干（60℃，约20min-1hr）。注意随时查看是否已蒸干。
6.2 DNA变性：
样本完全蒸干后，每个捕获中加入7.5μL 2×Hybridization Buffer (vial5)和3μL Hybridization Component A (vial 6)，震荡混匀，短暂离心。置于95℃变性10min。该步骤中的两种试剂都来自SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）。
6.3 文库与探针杂交：
取出探针，短暂离心。
将变性的DNA（始终保持在95℃）快速转移至含有探针的PCR管中，震荡混匀，短暂离心。
置于PCR仪中，47℃杂交。
6.4 配制纯化试剂：
一个捕获所需纯化试剂的配制方法如表11所示，根据捕获的个数按下表配制缓冲液。表格中试剂均来自SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）。
表11：捕获所需纯化试剂的配制试剂列表
Figure DEST_PATH_IMAGE270
孵育捕获磁珠（Capture Beads）和清洗缓冲液（Wash Buffer）工作液。其中，Capture Beads使用前须室温平衡30min，Wash Buffer使用前须47℃孵育2hr。
6.5 杂交后纯化：
每个捕获分装100μL捕获磁珠，将100μL捕获磁珠置于磁力架上至液体澄清，弃去上清。
加入200μL 1×Bead Wash Buffer (vial 7)，震荡混匀。置于磁力架上至液体澄清，弃去上清。1×Bead Wash Buffer (vial 7)来自SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）。
再次加入200μL 1×Bead Wash Buffer (vial 7)，震荡混匀。置于磁力架上至液体澄清，弃去上清。
再次加入100μL 1×Bead Wash Buffer(vial 7)，震荡混匀。置于磁力架上至液体澄清，彻底弃去上清。此时磁珠预处理完成，立即进行下一步试验。
将捕获过夜的杂交液体转入清洗好的磁珠中，移液器吹打十次。置于PCR仪中47℃孵育45min（PCR热盖温度设为57℃），每隔15min震荡一次保证磁珠悬浮。
6.6 清洗
该步骤所使用的试剂均来自SeqCap® Hyb and Wash Kit试剂盒（Roche公司生产）。
孵育完成后，每管加入100μL 47℃预热的1×Wash Buffer I(vial 1)，震荡混匀。置于磁力架上至液体澄清，弃去上清。
加入200μL 47℃预热的1×Stringent Wash Buffer (vial 4)，移液器吹打十次混匀。47℃孵育5min，置于磁力架上至液体澄清，弃去上清。
加入200μL 47℃预热的1×Stringent Wash Buffer (vial 4)，移液器吹打十次混匀。47℃孵育5min，置于磁力架上至液体澄清，弃去上清。
加入200μL室温放置的1×Wash Buffer I (vial 1)，振荡2min，短暂离心，置于磁力架上至液体澄清，弃去上清。
加入200μL室温放置的1×Wash Buffer II (vial 2)，震荡1min，短暂离心，放置磁力架上至液体澄清，弃去上清。
加入200μL室温放置的1×Wash Buffer III (vial 3)，震荡30sec，短暂离心，放置磁力架上至液体澄清，弃去上清。
向离心管中加入36μL超纯水洗脱，震荡混匀，进行下一步扩增试验。
6.7PCR反应：
根据捕获个数，按照表12配制混合液，震荡混匀。表中试剂均来自SeqCap® Hyband Wash Kit试剂盒（Roche公司生产）。
表12：混合液的配置试剂列表
Figure DEST_PATH_IMAGE271
短暂离心，将混合液分装至PCR管中，每管30μL。每个捕获样本分为两管进行PCR扩增，每管样本20μL。上述样本转入PCR反应中，震荡混匀，短暂离心。
置于PCR仪上，按表13的条件进行PCR反应。
表13：PCR反应的条件
Figure DEST_PATH_IMAGE272
6.8 扩增后纯化
取出纯化磁珠，室温平衡30min备用。
取180μL纯化磁珠于1.5mL离心管中，加入100μL扩增后的捕获DNA文库，振荡混匀，室温孵育15min。
置于磁力架上至液体澄清，弃去上清。
加入200μL 80％乙醇孵育30sec后弃去。注意：80%乙醇现用现配。重复一次200μL80％乙醇清洗步骤。
用10μL枪头弃去离心管底部的残留乙醇，室温干燥至乙醇完全挥发。
从磁力架取下离心管，加入120μL超纯水，振荡混匀。室温孵育2min。
短暂离心，置于磁力架上至液体澄清，将捕获样本转入新的离心管中。
文库混库和测序
将上述每个捕获按照数据量比例计算混库质量，按照数据量比例将不同捕获混合成一个样本。加入Phix文库混合成上机样本，进行测序。Phix为一种噬菌体，能够改善碱基不平衡，其作为参考品也可以对测序质量进行评估。
将下机FASTQ文件处理为各模块和软件可使用的输入文件
数据下机后，首先将下机数据从FASTQ文件处理成Bam文件，具体使用的软件和步骤如下：
8.1 去接头
调用Trimmomatic-0.36将每一对FASTQ文件都作为配对序列（paired reads）进行去除接头，之后切去剩余部分开头和结尾处碱基质量低于20的碱基，从reads的5’端开始，以大小为5的窗口进行划窗计算平均质量，如果窗口内平均碱基质量低于20，则切除该窗口，并要求切除后剩余碱基数量超过75，生成去接头后的FASTQ文件。
8.2 比对
调用Bismark-v0.19.0将去接头后的每一对FASTQ文件都作为配对序列比对到hg19人类参考基因组序列，生成初始bam文件；
8.3 去重
调用Bismark- v0.19.0的deduplicate模块，对初始Bam文件进行去重复处理，生成去重后的Bam文件；
8.4排序标记
调用SAMtools-1.3的sort模块，对去重后的Bam文件进行排序，生成排序后的Bam文件；调用Picard-2.1.0的AddOrReplaceReadGroups模块，对排序后的Bam文件进行标记分组；
8.5筛选
调用BamUtil-1.0.14的clipOverlap模块对标记分组后的Bam文件进行筛选，去除配对序列中的重叠部分，并调用SAMtools-1.3 view 对去除重合部分的序列的Bam文件的比对质量进行过滤，去除比对质量低于20的序列。完成该步骤后，统计每条序列中非甲基化位点中C碱基的转化率，剔除转化率低于95%的序列，并输出到最终的bam文件中；
8.6 建立索引
调用SAMtools-1.3的index模块对最终生成的Bam文件建立索引，生成与最终Bam文件配对的bai文件；
8.7 计算样本的位点甲基化水平
调用BisSNP-0.82.2对最终生成的Bam文件进行处理。首先调用BisμLfiteCountCovariates和BisμLfiteTableRecalibration模块，进行碱基质量校正，生成校正后的csv文件和Bam文件；再利用BisμLfiteGenotyper模块，鉴定待测样本的SNP位点和CpG位点，生成SNP和CpG的原始VCF文件；根据生成的VCF文件，调用VCFpostprocess模块对CpG位点进行过滤，得到最终的CpG位点及其甲基化水平，输出各样本的甲基化水平结果文件。该文件为本发明所需的输入文件之一。
用本发明方法划分甲基化连锁区域
本实例中一共使用了166个肿瘤样本，以将距离较近、甲基化水平较相关的甲基化位点合并为甲基化连锁区域，流程如图3所示。
9.1 生成肿瘤组织样本的位点矩阵文件
甲基化连锁区域划分开始之后，位点甲基化水平提取模块根据输入的肿瘤组织样本BisSNP输出的bed文件生成2个txt格式的肿瘤组织位点beta值矩阵文件和肿瘤组织位点深度矩阵文件。
9.2 划分甲基化连锁区域
获得两个矩阵文件后，将矩阵文件和肿瘤组织样本列表文件同时输入甲基化连锁 区域合并模块，其中，肿瘤组织样本列表文件仅包含一列，无标题行，其中记录有166个肿瘤 样本的名称。模块的最低深度
Figure DEST_PATH_IMAGE273
设置为100，最大位点间距
Figure DEST_PATH_IMAGE274
设置为100，最小相 关性
Figure DEST_PATH_IMAGE275
设置为0.95。该模块最终合并出了6042个甲基化连锁区域，记录在输出的甲基化 连锁区域列表文件中，完成对甲基化连锁区域的划分。输出列表中包含了所有甲基化连锁 区域的详细信息：在目标区间内的甲基化位点位置、位点与相邻位点的距离和相关系数、是 否满足甲基化连锁区域的标准、以及满足标准的甲基化位点所属的甲基化连锁区域起始、 终止位置。
用本发明方法筛选肿瘤特异连锁区域
本实例中共使用11个肺癌肿瘤样本和50个基线健康人血浆样本，用于筛选可用于得分计算的肿瘤差异连锁区域，流程如图4所示。
10.1生成肿瘤组织样本与基线健康人血浆样本的位点矩阵文件
肿瘤差异连锁区域筛选开始之后，位点甲基化水平提取模块根据输入的1个肺癌肿瘤样本、50个基线健康人血浆样本的BisSNP输出bed文件和基线血浆样本BisSNP输出txt格式的bed文件生成组织与基线血浆位点beta值矩阵文件及组织与基线血浆位点深度矩阵文件。
10.2生成肿瘤组织样本与基线健康人血浆样本的连锁区域矩阵文件
获得两个矩阵文件后，将矩阵文件和9.2中获得的连锁区域列表文件输入区域甲基化平均水平提取模块，对每个样本在6042个甲基化连锁区域上的平均beta值与平均深度依次进行计算，生成txt格式的组织与基线血浆连锁区域平均beta值矩阵文件及组织与基线血浆连锁区域平均深度矩阵文件。
10.3 筛选肿瘤特异连锁区域
将10.2中生成的连锁区域平均beta值矩阵文件、平均深度矩阵文件及肿瘤组织样 本分组信息表作为差异区域筛选模块输入文件，其中，样本分组信息表包括两列，分别记录 11例肺癌肿瘤样本和50例基线血浆样本的名称和所属分组，且肿瘤和血浆样本的组别分别 设置为case组和control组。参数中，最低深度要求
Figure DEST_PATH_IMAGE276
为100，最大调整后错误发现率
Figure DEST_PATH_IMAGE277
为0.05，最小甲基化水平差值
Figure DEST_PATH_IMAGE278
为0.1，差异样本占比阈值
Figure DEST_PATH_IMAGE279
为0.8，低甲基 化阈值
Figure DEST_PATH_IMAGE280
为0.02，低甲基化基线样本占比阈值
Figure DEST_PATH_IMAGE281
为0.8。筛选过后输出的特异连 锁区域列表文件中，包含了所有特异连锁区域的平均beta值均值、方差、组间均值差、错误 发现率等信息。共筛选出208个满足筛选条件的特异连锁区域。
用本发明方法对梯度稀释血浆样本计算甲基化得分并进行评价
实例中共使用了50个健康人血浆样本作为基线，计算了这50个样本和12个梯度稀释血浆样本的得分，并对梯度稀释血浆样本进行了评价，实施流程图如图5所示。
11.1生成基线健康人血浆样本与梯度稀释血浆样本的连锁区域矩阵文件
将基线血浆BisSNP输出bed文件和梯度稀释血浆BisSNP输出bed文件输入位点甲基化水平提取模块；之后将生成的基线与梯度稀释血浆位点beta值矩阵、基线与梯度稀释血浆位点深度矩阵及步骤9.2中获得的交际花连锁区域列表文件同时输入区域甲基化平均水平提取模块，生成基线与梯度稀释血浆区域平均beta值矩阵文件和平均深度矩阵文件。
11.2 计算基线健康人血浆样本与梯度稀释血浆样本的甲基化得分
将11.1中生成的两个矩阵、10.3中筛选获得的特异连锁区域列表和基线与梯度稀释血浆样本信息表作为基线构建与得分计算模块的输入文件，其中，基线与梯度稀释血浆样本信息表中需要包含两列信息：50个基线血浆样本与12个梯度稀释血浆样本的样本名与所属分组，前者在表中的组别为正常人群组，后者为待测样本组。
本步骤中使用50个健康人血浆作为基线，对每一特异连锁区域Logit转化后的平均beta值建立分布，并计算每个样本与分布的差异值p-value。计算完p-value后，使用208个区域的平均测序深度作为权重，对处理后的p-value计算加权平均数作为样本的得分，用于量化样本与健康人血浆的差别。最终会输出健康人血浆和梯度稀释血浆样本的甲基化得分列表。
11.3 对梯度稀释血浆样本进行预测并评估其准确性
根据11.2中输出的得分列表，最终计算出的基线血浆样本的甲基化得分在1.25-3.60之间，Sample A~ Sample D梯度稀释后得到的12个样本的甲基化得分见表14。
表14：梯度稀释样本的甲基化得分见表
Figure DEST_PATH_IMAGE282
其中，理论ctDNA占比由原始血浆样本的ctDNA占比乘稀释梯度计算获得。
通过计算线性相关系数，评估本发明计算出的甲基化得分与理论ctDNA占比的相关性。本实例中计算得到的线性相关系数为0.82，相关性明显，说明筛选出的甲基化标志物能够提示血浆中ctDNA的存在。每个原始样本在不同稀释梯度对应的甲基化得分变化如图6所示，其中，横坐标为稀释梯度，纵坐标为甲基化得分；12个梯度稀释样本理论ctDNA占比和甲基化得分的散点图如图7所示，其中，横坐标为理论ctDNA占比，纵坐标为甲基化得分，黑色虚线为健康人中最高分（3.60）。当使用虚线值作为阈值时，理论ctDNA占比高于0.5%的样本得分均高于该值，特异性与灵敏度均为100%，说明采用本发明筛选出的甲基化标志物检出率高，拥有高特异性、高灵敏度、低检测下限等优势。
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各程序模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的程序模块完成，即将装置的内部结构划分成不同的程序单元或模块，以完成以上描述的全部或者部分功能。实施例中的各程序模块可以集成在一个处理单元中，也可是各个单元单独物理存在，也可以两个或两个以上单元集成在一个处理单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件程序单元的形式实现。另外，各程序模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。
图8是本发明一个实施例中提供的终端设备的结构示意图，如所示，该终端设备200包括：处理器220、存储器210以及存储在存储器210中并可在处理器220上运行的计算机程序211，例如：基于靶向捕获测序的甲基化标志物筛选与评价关联程序。处理器220执行计算机程序211时实现上述各个基于靶向捕获测序的甲基化标志物筛选与评价方法实施例中的步骤，或者，处理器220执行计算机程序211时实现上述基于靶向捕获测序的甲基化标志物筛选与评价装置实施例中各模块的功能。
终端设备200可以为笔记本、掌上电脑、平板型计算机、手机等设备。终端设备200可包括，但不仅限于处理器220、存储器210。本领域技术人员可以理解，图8仅仅是终端设备200的示例，并不构成对终端设备200的限定，可以包括比图示更多或更少的部件，或者组合某些部件，或者不同的部件，例如：终端设备200还可以包括输入输出设备、显示设备、网络接入设备、总线等。
处理器220可以是中央处理单元(Central Processing Unit，CPU)，还可以是其他通用处理器、数字信号处理器 (Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器220可以是微处理器或者该处理器也可以是任何常规的处理器等。
存储器210可以是终端设备200的内部存储单元，例如：终端设备200的硬盘或内存。存储器210也可以是终端设备200的外部存储设备，例如：终端设备200上配备的插接式硬盘，智能TF存储卡(Smart Media Card，SMC)，安全数字(Secure Digital，SD)卡，闪存卡(Flash Card)等。进一步地，存储器210还可以既包括终端设备200的内部存储单元也包括外部存储设备。存储器210用于存储计算机程序211以及终端设备200所需要的其他程序和数据。存储器210还可以用于暂时地存储已经输出或者将要输出的数据。
在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详细描述或记载的部分，可以参见其他实施例的相关描述。
本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
在本申请所提供的实施例中，应该理解到，所揭露的装置/终端设备和方法，可以通过其他的方式实现。例如，以上所描述的装置/终端设备实施例仅仅是示意性的，例如，模块或单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通讯连接可以是通过一些接口，装置或单元的间接耦合或通讯连接，可以是电性、机械或其他的形式。
作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本申请各个实施例中的各功能单元可能集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读存储介质中。基于这样的理解，本发明实现上述实施例方法中的全部或部分流程，也可以通过计算机程序211发送指令给相关的硬件完成，的计算机程序211可存储于一计算机可读存储介质中，该计算机程序211在被处理器220执行时，可实现上述各个方法实施例的步骤。其中，计算机程序211包括：计算机程序代码，计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。计算机可读存储介质可以包括：能够携带计算机程序211代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器 (ROM，Read-Only Memory)、随机存取存储器(RAM，RandomAccess Memory)、电载波信号、电信信号以及软件分发介质等。需要说明的是，计算机可读存储介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如：在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括电载波信号和电信信号。
应当说明的是，上述实施例均可根据需要自由组合。以上仅是本发明的优选实施方式，应当指出，对于本技术领域的普通相关人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (12)

1.一种基于靶向捕获测序的甲基化标志物筛选与评价方法，其特征在于，包括：
分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件，所述待测样本为血浆样本；
依次计算各待测样本Bam文件目标区域内
Figure DEST_PATH_IMAGE001
个甲基化位点上的甲基化水平和覆盖深度， 并合并得到甲基化水平矩阵和位点深度矩阵；
基于所述甲基化水平矩阵和位点深度矩阵，针对待测样本中的每一个甲基化位点，分 别计算其与下一甲基化位点之间的距离及甲基化水平的线性相关系数，并根据计算结果依 次对相邻甲基化位点进行合并得到甲基化连锁区域，将包含有预设数量甲基化位点的
Figure DEST_PATH_IMAGE002
个 甲基化连锁区域作为甲基化标志物候选区域输出；
依次计算各待测样本中
Figure 199505DEST_PATH_IMAGE002
个甲基化连锁区域内的甲基化水平均值和位点深度均值，生 成连锁区域甲基化水平均值矩阵和位点深度均值矩阵；
根据所述连锁区域甲基化水平均值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组存在设定差异的特异连锁区域，得到甲基化标志物；
根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据所述甲基化得分对甲基化标志物进行评价。
2.如权利要求1所述的甲基化标志物筛选与评价方法，其特征在于，在所述分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件之后，还包括：
根据预先设定的C-T转化率对生成的Bam文件中的reads进行逐条过滤，得到过滤后的Bam文件；和/或，
根据目标区域Bed文件及预先设定的各reads中覆盖甲基化位点的数量对Bam文件进行过滤，得到过滤后的Bam文件。
3.如权利要求1或2所述的甲基化标志物筛选与评价方法，其特征在于，在所述依次计 算各待测样本Bam文件目标区域内
Figure 802131DEST_PATH_IMAGE001
个甲基化位点上的甲基化水平和覆盖深度，并合并到甲 基化水平矩阵和位点深度矩阵中，包括：
依次提取各待测样本Bam文件中各个甲基化位点的正链信息和负链信息；
依次计算各待测样本中各个甲基化位点的甲基化水平和覆盖深度，其中，样本n在甲基 化位点
Figure DEST_PATH_IMAGE003
上的甲基化水平
Figure DEST_PATH_IMAGE004
和覆盖深度
Figure DEST_PATH_IMAGE005
分别为：
Figure DEST_PATH_IMAGE006
其中，
Figure DEST_PATH_IMAGE007
表示样本n在甲基化位点
Figure 412496DEST_PATH_IMAGE003
上的正链甲基化水平，
Figure DEST_PATH_IMAGE008
表 示样本n在甲基化位点
Figure 110501DEST_PATH_IMAGE003
上的负链甲基化水平，
Figure DEST_PATH_IMAGE009
表示样本n在甲基化位点
Figure 180000DEST_PATH_IMAGE003
上 的正链覆盖深度，
Figure DEST_PATH_IMAGE010
表示样本n在甲基化位点
Figure 204061DEST_PATH_IMAGE003
上的负链覆盖深度，
Figure DEST_PATH_IMAGE011
，
Figure DEST_PATH_IMAGE012
；
将计算得到的各甲基化位点的甲基化水平和覆盖深度进行合并得到甲基化水平矩阵
Figure DEST_PATH_IMAGE013
和位点深度矩阵
Figure DEST_PATH_IMAGE014
，所述甲基化水平矩阵
Figure DEST_PATH_IMAGE015
和位点深度矩阵
Figure DEST_PATH_IMAGE016
均为I行N列矩阵，其中，行对应甲基化位点，列对应待测样本。
4.如权利要求1或2所述的甲基化标志物筛选与评价方法，其特征在于，在所述基于甲 基化水平矩阵和位点深度矩阵，针对待测样本中的每一个甲基化位点，分别计算其与下一 甲基化位点之间的距离及甲基化水平的线性相关系数，并根据计算结果依次对相邻甲基化 位点进行合并得到甲基化连锁区域，将包含有预设数量甲基化位点的
Figure DEST_PATH_IMAGE017
个甲基化连锁区域 作为甲基化标志物候选区域输出中，针对样本n甲基化位点
Figure DEST_PATH_IMAGE018
的甲基化连锁区域合并步骤 包括：
判断甲基化位点
Figure 468165DEST_PATH_IMAGE018
及其下一甲基化位点
Figure DEST_PATH_IMAGE019
的覆盖深度是否均在预设最低深度要求
Figure DEST_PATH_IMAGE020
之上；
若是，计算甲基化位点
Figure 976986DEST_PATH_IMAGE018
和甲基化位点
Figure 150260DEST_PATH_IMAGE019
之间的距离
Figure DEST_PATH_IMAGE021
；
计算甲基化位点
Figure 455208DEST_PATH_IMAGE018
上各待测样本的甲基化水平
Figure DEST_PATH_IMAGE022
与甲基化位点
Figure 577579DEST_PATH_IMAGE019
上 各待测样本的甲基化水平
Figure DEST_PATH_IMAGE023
之间的线性相关系数
Figure DEST_PATH_IMAGE024
；其中，
Figure DEST_PATH_IMAGE025
表示样本1在甲基化位点
Figure 878373DEST_PATH_IMAGE018
上的甲基化水平，
Figure DEST_PATH_IMAGE026
表示样本N在甲基化位点
Figure 793238DEST_PATH_IMAGE018
上的甲基化 水平，
Figure DEST_PATH_IMAGE027
表示样本1在甲基化位点
Figure 778075DEST_PATH_IMAGE019
上的甲基化水平，
Figure DEST_PATH_IMAGE028
表示样本N在甲基 化位点
Figure 271897DEST_PATH_IMAGE019
上的甲基化水平；
判断是否同时满足条件
Figure DEST_PATH_IMAGE029
和
Figure DEST_PATH_IMAGE030
，其中，
Figure DEST_PATH_IMAGE031
为预设最大位点间间 距，
Figure DEST_PATH_IMAGE032
为预设最小相关系数；
若是，将甲基化位点
Figure 729379DEST_PATH_IMAGE019
并入甲基化位点
Figure 996019DEST_PATH_IMAGE018
当前所在的甲基化连锁区域
Figure DEST_PATH_IMAGE033
形成新的甲 基化连锁区域
Figure 586312DEST_PATH_IMAGE033
，否则断开甲基化位点
Figure 719DEST_PATH_IMAGE018
当前所在的甲基化连锁区域
Figure 595098DEST_PATH_IMAGE033
，
Figure DEST_PATH_IMAGE034
。
5.如权利要求1或2所述的甲基化标志物筛选与评价方法，其特征在于，所述依次计算 各待测样本中
Figure 404660DEST_PATH_IMAGE017
个甲基化连锁区域内的甲基化水平均值和位点深度均值，生成连锁区域甲 基化水平均值矩阵和位点深度均值矩阵中：
样本
Figure DEST_PATH_IMAGE035
于甲基化连锁区域
Figure 556067DEST_PATH_IMAGE033
中的甲基化水平均值
Figure DEST_PATH_IMAGE036
为：
Figure DEST_PATH_IMAGE037
样本
Figure 435600DEST_PATH_IMAGE035
于甲基化连锁区域
Figure 941274DEST_PATH_IMAGE033
中的位点深度均值
Figure DEST_PATH_IMAGE038
为：
Figure DEST_PATH_IMAGE039
在根据计算的甲基化水平均值
Figure 109737DEST_PATH_IMAGE036
和位点深度均值
Figure 755349DEST_PATH_IMAGE038
形成的连锁区域甲基化 水平均值矩阵
Figure DEST_PATH_IMAGE040
和位点深度均值矩阵
Figure DEST_PATH_IMAGE041
中，行对应甲基化连锁区 域，列对应样本。
6.如权利要求1或2所述的甲基化标志物筛选与评价方法，其特征在于，所述根据所述连锁区域甲基化水平均值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组存在设定差异的特异连锁区域中，包括：
将所述甲基化水平均值矩阵
Figure DEST_PATH_IMAGE042
进行
Figure DEST_PATH_IMAGE043
转换得到矩阵
Figure DEST_PATH_IMAGE044
；
根据预先设定的样本分组信息表将矩阵
Figure DEST_PATH_IMAGE045
、甲基化水平均值矩阵
Figure DEST_PATH_IMAGE046
和 位点深度均值矩阵
Figure DEST_PATH_IMAGE047
分别按照control组和case组进行分割，得到矩阵
Figure DEST_PATH_IMAGE048
、矩阵
Figure DEST_PATH_IMAGE049
、矩阵
Figure DEST_PATH_IMAGE050
、矩阵
Figure DEST_PATH_IMAGE051
、 矩阵
Figure DEST_PATH_IMAGE052
和矩阵
Figure DEST_PATH_IMAGE053
6个独立矩阵，其中，control组中矩阵
Figure 565522DEST_PATH_IMAGE048
、矩阵
Figure DEST_PATH_IMAGE054
和矩阵
Figure DEST_PATH_IMAGE055
包含正常人群组的样 本数据，case组中矩阵
Figure DEST_PATH_IMAGE056
、矩阵
Figure DEST_PATH_IMAGE057
和矩阵
Figure DEST_PATH_IMAGE058
包含异常人 群组的样本数据；
遍历所有甲基化连锁区域，筛选出与正常人群组存在设定差异的特异连锁区域，其中， 针对甲基化连锁区域
Figure 42597DEST_PATH_IMAGE033
的筛选过程为：
依次计算control组与case组的错误发现率
Figure DEST_PATH_IMAGE059
、甲基化水平差值
Figure DEST_PATH_IMAGE060
、case组差异 样本占比
Figure DEST_PATH_IMAGE061
及control组中的低甲基化基线样本占比
Figure DEST_PATH_IMAGE062
，其中，所述case组差异样本占比
Figure 595061DEST_PATH_IMAGE061
表示case组
Figure DEST_PATH_IMAGE063
转换后甲基化水平在control组一倍标准差范围外的样本数量占case组 样本总数量的比值，所述control组中的低甲基化基线样本占比
Figure 76553DEST_PATH_IMAGE062
表示control组中甲基化 水平低于预设低甲基化阈值
Figure DEST_PATH_IMAGE064
的样本数量占control组样本总数量的比值，
Figure DEST_PATH_IMAGE065
；
判断是否同时满足条件
Figure DEST_PATH_IMAGE066
、
Figure DEST_PATH_IMAGE067
、
Figure DEST_PATH_IMAGE068
及
Figure DEST_PATH_IMAGE069
，其 中，
Figure DEST_PATH_IMAGE070
为预设最大错误发现率，
Figure DEST_PATH_IMAGE071
为预设最小甲基化水平差值，
Figure DEST_PATH_IMAGE072
为预设差异 样本占比阈值，
Figure DEST_PATH_IMAGE073
为预设低甲基化基线样本占比阈值；
若是，判断甲基化连锁区域
Figure DEST_PATH_IMAGE074
与正常人存在设定差异。
7.如权利要求6所述的甲基化标志物筛选与评价方法，其特征在于，
所述control组与case组的错误发现率
Figure DEST_PATH_IMAGE075
由control组所有待测样本
Figure DEST_PATH_IMAGE076
转换后 的甲基化水平
Figure DEST_PATH_IMAGE077
与case组所有待测样本
Figure 563902DEST_PATH_IMAGE076
转换后的 甲基化水平
Figure DEST_PATH_IMAGE078
经校正t检验得到的差异检验值
Figure DEST_PATH_IMAGE079
进一步经过 Benjaminiand Hochberg法校正后得到，其中，
Figure DEST_PATH_IMAGE080
表示control组样本1在甲基化 连锁区域
Figure DEST_PATH_IMAGE081
中
Figure 239077DEST_PATH_IMAGE076
转换后的甲基化水平，
Figure DEST_PATH_IMAGE082
表示control组样本N在甲基 化连锁区域
Figure 428137DEST_PATH_IMAGE081
中
Figure 329797DEST_PATH_IMAGE076
转换后的甲基化水平；
Figure DEST_PATH_IMAGE083
表示case组样本1在甲基化连锁区 域
Figure 991592DEST_PATH_IMAGE081
中
Figure 710936DEST_PATH_IMAGE076
转换后的甲基化水平，
Figure DEST_PATH_IMAGE084
表示case组样本N在甲基化连锁区域
Figure 204103DEST_PATH_IMAGE081
中
Figure 96840DEST_PATH_IMAGE076
转换后的甲基化水平；
和/或，所述甲基化水平差值
Figure DEST_PATH_IMAGE085
表示control组平均甲基化水平与case组平均甲基化 水平之差：
Figure DEST_PATH_IMAGE086
其中，
Figure DEST_PATH_IMAGE087
表示case组甲基化连锁区域
Figure 750325DEST_PATH_IMAGE081
中样本1的平均甲基化水平，
Figure DEST_PATH_IMAGE088
表示case组甲基化连锁区域
Figure 556914DEST_PATH_IMAGE081
中样本N的平均甲基化水平，
Figure DEST_PATH_IMAGE089
表示control组甲基化连锁区域
Figure 48463DEST_PATH_IMAGE081
中样本1的平均甲基化水平，
Figure DEST_PATH_IMAGE090
表示control组甲基化连锁区域
Figure 720621DEST_PATH_IMAGE081
中样本N的平均甲基化水平；
和/或，所述case组差异样本占比
Figure DEST_PATH_IMAGE091
为：
Figure DEST_PATH_IMAGE092
其中，
Figure DEST_PATH_IMAGE093
表示case组甲基化连锁区域
Figure 610472DEST_PATH_IMAGE081
中
Figure 668689DEST_PATH_IMAGE076
转换后样本
Figure DEST_PATH_IMAGE094
的甲基化水平，
Figure DEST_PATH_IMAGE095
表示control组甲基化连锁区域
Figure 612812DEST_PATH_IMAGE081
中
Figure 95353DEST_PATH_IMAGE076
转换后的平均甲基化水平，
Figure DEST_PATH_IMAGE096
表示control组甲基化连锁区域
Figure 707118DEST_PATH_IMAGE081
中
Figure 506709DEST_PATH_IMAGE076
转换 后的标准差，
Figure DEST_PATH_IMAGE097
表示case组样本总数；
和/或，control组中的低甲基化基线样本占比
Figure DEST_PATH_IMAGE098
为：
Figure DEST_PATH_IMAGE099
其中，
Figure DEST_PATH_IMAGE100
表示control组甲基化连锁区域
Figure 159053DEST_PATH_IMAGE081
中样本
Figure DEST_PATH_IMAGE101
的甲基化水平，
Figure DEST_PATH_IMAGE102
表示control组样本总数，
Figure DEST_PATH_IMAGE103
表示预设 的甲基化水平背景噪音最大值。
8.如权利要求6所述的甲基化标志物筛选与评价方法，其特征在于，在所述根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据所述甲基化得分对甲基化标志物进行评价中，包括：
将所述甲基化水平均值矩阵
Figure DEST_PATH_IMAGE104
进行
Figure 497981DEST_PATH_IMAGE076
转换得到矩阵
Figure DEST_PATH_IMAGE105
；
从所述矩阵
Figure 733395DEST_PATH_IMAGE105
和位点深度均值矩阵
Figure DEST_PATH_IMAGE106
中提取筛选得到的
Figure DEST_PATH_IMAGE107
个特 异连锁区域的数据，并根据预先设定的样本分组信息表将其分割为矩阵
Figure DEST_PATH_IMAGE108
、矩阵
Figure DEST_PATH_IMAGE109
、矩阵
Figure DEST_PATH_IMAGE110
和矩阵
Figure DEST_PATH_IMAGE111
，其中，矩阵
Figure DEST_PATH_IMAGE112
和矩阵
Figure DEST_PATH_IMAGE113
包含 正常人群组样本数据，矩阵
Figure DEST_PATH_IMAGE114
和矩阵
Figure DEST_PATH_IMAGE115
包含待测样本数据；
根据分割得到的矩阵分别计算各待测样本的甲基化得分，并判断是否存在甲基化得分大于预设分数阈值的待测样本；若是，判定该待测样本中包含筛选获得的甲基化标志物；
其中，待测样本
Figure DEST_PATH_IMAGE116
甲基化得分
Figure DEST_PATH_IMAGE117
为：
Figure DEST_PATH_IMAGE118
其中，
Figure DEST_PATH_IMAGE119
表示待测样本
Figure DEST_PATH_IMAGE120
在特异连锁区域
Figure DEST_PATH_IMAGE121
的平均深度，
Figure DEST_PATH_IMAGE122
；
Figure DEST_PATH_IMAGE123
表示待测样本
Figure DEST_PATH_IMAGE124
在特异连锁区域
Figure DEST_PATH_IMAGE125
的p-value值，所述p- value值为待测样本
Figure DEST_PATH_IMAGE126
在特异连锁区域
Figure DEST_PATH_IMAGE127
上
Figure 327491DEST_PATH_IMAGE076
转换后甲基化水平
Figure DEST_PATH_IMAGE128
的 Z-score值
Figure DEST_PATH_IMAGE129
转换为标准正态分布的分位数：
Figure DEST_PATH_IMAGE130
其中，
Figure DEST_PATH_IMAGE131
为正常人群组
Figure 320242DEST_PATH_IMAGE076
转换后的甲基化水平均值，
Figure DEST_PATH_IMAGE132
为正常人群组
Figure 637261DEST_PATH_IMAGE076
转换后的方差。
9.如权利要求8所述的甲基化标志物筛选与评价方法，其特征在于，在所述根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据所述甲基化得分对甲基化标志物进行评价中，还包括：根据每个待测样本的已知分组情况计算检出的灵敏度和特异性，或根据待测样本的已知ctDNA浓度，计算甲基化得分与ctDNA浓度的线性相关系数，进而根据灵敏度、特异性与线性相关系数的高低对筛选的甲基化标志物进行评价的步骤。
10.一种基于靶向捕获测序的甲基化标志物筛选与评价装置，其特征在于，应用于如权利要求1-9任意一项所述的甲基化标志物筛选与评价方法，所述甲基化标志物筛选与评价装置中包括：
Bam文件生成模块，用于分别获取N个待测样本捕获测序的FASTQ文件，并分别与参考基因组进行比对生成Bam文件，所述待测样本为血浆样本；
位点甲基化水平提取模块，用于依次计算各待测样本Bam文件目标区域内
Figure 705843DEST_PATH_IMAGE001
个甲基化位 点上的甲基化水平和覆盖深度，并合并得到甲基化水平矩阵和位点深度矩阵；
甲基化连锁区域合并模块，用于基于所述甲基化水平矩阵和位点深度矩阵，针对待测 样本中的每一个甲基化位点，分别计算其与下一甲基化位点之间的距离及甲基化水平的线 性相关系数，并根据计算结果依次对相邻甲基化位点进行合并得到甲基化连锁区域，将包 含有预设数量甲基化位点的
Figure 646904DEST_PATH_IMAGE002
个甲基化连锁区域作为甲基化标志物候选区域输出；
区域甲基化平均水平提取模块，用于依次计算各待测样本中
Figure 317663DEST_PATH_IMAGE002
个甲基化连锁区域内的 甲基化水平均值和位点深度均值，生成连锁区域甲基化水平均值矩阵和位点深度均值矩 阵；
差异区域筛选模块，用于根据所述连锁区域甲基化水平均值矩阵和位点深度均值矩阵于合并得到的甲基化连锁区域中进一步筛选出与正常人群组存在设定差异的特异连锁区域，得到甲基化标志物；
基线构建与得分计算模块，用于根据筛选得到的特异连锁区域分别计算各待测样本的甲基化得分，并根据所述甲基化得分对甲基化标志物进行评价。
11.一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器运行所述计算机程序时实现如权利要求1-9中任一项所述基于靶向捕获测序的甲基化标志物筛选与评价方法的步骤。
12.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1-9中任一项所述基于靶向捕获测序的甲基化标志物筛选与评价方法的步骤。
CN202110078570.2A 2021-01-20 2021-01-21 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置 Active CN112397151B (zh)

Priority Applications (4)

Application Number Priority Date Filing Date Title
CN202110078570.2A CN112397151B (zh) 2021-01-21 2021-01-21 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置
EP21920475.7A EP4268231A1 (en) 2021-01-20 2021-04-30 Dna methylation sequencing analysis methods
PCT/CN2021/091761 WO2022156089A1 (en) 2021-01-20 2021-04-30 Dna methylation sequencing analysis methods
US17/490,549 US20220228209A1 (en) 2021-01-20 2021-09-30 Dna methylation sequencing analysis methods

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110078570.2A CN112397151B (zh) 2021-01-21 2021-01-21 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置

Publications (2)

Publication Number Publication Date
CN112397151A CN112397151A (zh) 2021-02-23
CN112397151B true CN112397151B (zh) 2021-04-20

Family

ID=74625106

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110078570.2A Active CN112397151B (zh) 2021-01-20 2021-01-21 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置

Country Status (1)

Country Link
CN (1) CN112397151B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2022156089A1 (en) * 2021-01-20 2022-07-28 Genecast Biotechnology Co., Ltd Dna methylation sequencing analysis methods
CN112951418B (zh) * 2021-05-17 2021-08-06 臻和(北京)生物科技有限公司 基于液体活检的连锁区域甲基化评估方法和装置、终端设备及存储介质
WO2023184330A1 (zh) * 2022-03-31 2023-10-05 京东方科技集团股份有限公司 基因组甲基化测序数据的处理方法、装置、设备和介质
CN115064211B (zh) * 2022-08-15 2023-01-24 臻和(北京)生物科技有限公司 一种基于全基因组甲基化测序的ctDNA预测方法及装置
CN115497561B (zh) * 2022-09-01 2023-08-29 北京吉因加医学检验实验室有限公司 一种甲基化标志物分层筛选的方法及装置
CN116287279B (zh) * 2023-05-25 2023-08-04 臻和(北京)生物科技有限公司 用于检测胰腺癌的生物标志物及其应用

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
PT3543356T (pt) * 2014-07-18 2021-10-04 Univ Hong Kong Chinese Análise dos padrões de metilação de tecidos em mistura de adn
CN108949970A (zh) * 2017-05-23 2018-12-07 中国科学院深圳先进技术研究院 基于多组学的宫颈癌特征获取方法和系统
CN107190076B (zh) * 2017-06-28 2019-12-27 中国科学院苏州生物医学工程技术研究所 一种人肿瘤相关的甲基化位点及其筛选方法和用途
JP7320067B2 (ja) * 2019-01-18 2023-08-02 ザ リージェンツ オブ ザ ユニバーシティ オブ カリフォルニア 保存された遺伝子座に基づく、哺乳動物に関するｄｎａメチル化測定
CN110982907B (zh) * 2020-02-27 2020-07-03 上海鹍远生物技术有限公司 甲状腺结节相关rDNA甲基化标志物及其应用

Also Published As

Publication number Publication date
CN112397151A (zh) 2021-02-23

Similar Documents

Publication Publication Date Title
CN112397151B (zh) 基于靶向捕获测序的甲基化标志物筛选与评价方法及装置
CN112029861B (zh) 基于捕获测序技术的肿瘤突变负荷检测装置及方法
CN112397150B (zh) 基于目标区域捕获测序的ctDNA甲基化水平预测装置及方法
CN109767810B (zh) 高通量测序数据分析方法及装置
CN110211633B (zh) Mgmt基因启动子甲基化的检测方法、测序数据的处理方法及处理装置
CN111304303B (zh) 微卫星不稳定的预测方法及其应用
CN106544407A (zh) 确定受体cfDNA样本中供体来源cfDNA比例的方法
CN104232777A (zh) 同时确定胎儿核酸含量和染色体非整倍性的方法及装置
CN108319813A (zh) 循环肿瘤dna拷贝数变异的检测方法和装置
CN112735531B (zh) 循环无细胞核小体活性区域的甲基化分析方法和装置、终端设备及存储介质
CN113096728B (zh) 一种微小残余病灶的检测方法、装置、存储介质及设备
WO2020224159A1 (zh) 基于二代测序用于脑胶质瘤的检测panel、检测试剂盒、检测方法及其应用
CN111968701A (zh) 检测指定基因组区域体细胞拷贝数变异的方法和装置
KR20200035427A (ko) 세포-무함유 바이러스 핵산을 사용하는 암 스크리닝의 증강
CN110106063B (zh) 基于二代测序的用于神经胶质瘤1p/19q联合缺失检测的系统
CN116064819A (zh) 一种检测ctDNA中肿瘤特异基因的突变和甲基化的方法
CN109461473B (zh) 胎儿游离dna浓度获取方法和装置
CN108570496A (zh) 一种遗传性骨病的分子诊断方法及试剂盒
CN109439741B (zh) 检测特发性癫痫病基因探针组合物、试剂盒及应用
CN103261442A (zh) Hpv 精确分型的生物信息学分析的方法及系统
CN116312779A (zh) 检测样本污染和识别样本错配的方法和装置
CN106874710A (zh) 一种用于利用肿瘤ffpe样本检测体细胞突变的装置
WO2019132010A1 (ja) 塩基配列における塩基種を推定する方法、装置及びプログラム
WO2020124625A1 (zh) 基于ctDNA的基因检测方法、装置、存储介质及计算机系统
CN110993024B (zh) 建立胎儿浓度校正模型的方法及装置与胎儿浓度定量的方法及装置

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
CP01 Change in the name or title of a patent holder

Address after: 100191 903, 9 / F, healthsmart Valley Building, 35 Huayuan North Road, Haidian District, Beijing

Patentee after: Zhenhe (Beijing) Biotechnology Co.,Ltd.

Patentee after: Wuxi Zhenhe Biotechnology Co.,Ltd.

Address before: 100191 903, 9 / F, healthsmart Valley Building, 35 Huayuan North Road, Haidian District, Beijing

Patentee before: Zhenhe (Beijing) Biotechnology Co.,Ltd.

Patentee before: Wuxi Zhenhe Biotechnology Co.,Ltd.

CP01 Change in the name or title of a patent holder