CN115132274B - 循环无细胞dna转录因子结合位点的甲基化水平分析方法及装置 - Google Patents

循环无细胞dna转录因子结合位点的甲基化水平分析方法及装置 Download PDF

Info

Publication number
CN115132274B
CN115132274B CN202211059714.0A CN202211059714A CN115132274B CN 115132274 B CN115132274 B CN 115132274B CN 202211059714 A CN202211059714 A CN 202211059714A CN 115132274 B CN115132274 B CN 115132274B
Authority
CN
China
Prior art keywords
transcription factor
methylation
binding site
factor binding
free dna
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
CN202211059714.0A
Other languages
English (en)
Other versions
CN115132274A (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 CN202211059714.0A priority Critical patent/CN115132274B/zh
Publication of CN115132274A publication Critical patent/CN115132274A/zh
Application granted granted Critical
Publication of CN115132274B publication Critical patent/CN115132274B/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
    • 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/20—Sequence assembly

Landscapes

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

Abstract

本发明提供了一种循环无细胞DNA转录因子结合位点的甲基化水平分析方法及装置，包括：接收待分析血浆样本的甲基化测序数据并从中提取循环无细胞DNA分子片段；提取转录因子结合位点上下游区域CpG位点的甲基化状态；以转录因子结合位点的基因组位置作为参照，将其上下游区域的CpG位点的坐标对齐；针对每个转录因子，分别统计各相对坐标上甲基化的CpG分子片段占比；计算转录因子结合位点中心点与侧翼区域的甲基化占比差值；将甲基化占比差值输入甲基化水平分析模型，根据甲基化水平分析模型的输出结果完成对待分析血浆样本循环无细胞DNA转录因子结合位点甲基化水平的分析，实现通过甲基化数据对转录因子的集合状态进行评估的目的。

Description

循环无细胞DNA转录因子结合位点的甲基化水平分析方法及 装置
技术领域
本发明涉及生物医学技术领域，尤其涉及一种循环无细胞DNA转录因子结合位点的甲基化水平分析方法及装置。
背景技术
癌症是身体的部分细胞不受控制地生长和繁殖的一种情况。癌细胞可以入侵并破坏周围的健康组织，包括器官。癌症的形成需要一个漫长的过程，在临床前期，从正常细胞经过异常增生、癌前病变再到形成可见的癌症通常需要10-20年的潜伏期。极早期的癌症通常没有明显的症状，但当出现明显的症状再进行就医，通常已经到了中晚期。当患者处于癌症潜伏期，只要有意识的进行干预，改善生存环境，改变生活习惯，增强免疫力，就可以避免癌症的发生。
目前，我国的常规体检通常缺乏早期癌症的筛查项目。日常生活中，普通群众也缺乏对癌症早期预防诊断的认知。但是，在早期对癌症进行治疗，平均治愈率达80%，癌症的早期治疗不仅可以降低癌症患者的费用支持，还可以减轻患者的身心痛苦。因此，随着我国人口老龄化愈发严重，恶性肿瘤在老年人中发病率更高，癌症的早期筛查急需进一步加强。
转录因子可以调节靶基因的表达，这些基因通常在发育和分化中起关键作用。转录因子结合通常与核小体占位相关。近来，有研究表明，核小体定位可以从循环无细胞DNA推断，这表明从肿瘤循环无细胞DNA中不仅可以推断基因表达，还可以推断肿瘤样本中转录因子的结合。从循环无细胞DNA推断转录因子结合在癌症和其他潜在疾病中的诊断潜力很大，但迄今为止，转录因子的集合状态还不能通过甲基化数据进行无创评估。
发明内容
针对上述问题，本发明提供了一种循环无细胞DNA转录因子结合位点的甲基化水平分析方法及装置，通过甲基化数据对转录因子的集合状态进行评估。
本发明提供的技术方案如下：
一方面，本发明提供了一种循环无细胞DNA转录因子结合位点的甲基化水平分析方法，包括：
S100 接收待分析血浆样本的全基因组甲基化测序数据并从中提取循环无细胞DNA分子片段；
S200 针对提取的循环无细胞DNA分子片段，进一步提取转录因子结合位点上下游区域胞嘧啶-磷酸-鸟嘌呤位点的甲基化状态；
S300 以转录因子结合位点的基因组位置作为参照，将其上下游区域的胞嘧啶-磷酸-鸟嘌呤位点的坐标对齐，得到各胞嘧啶-磷酸-鸟嘌呤位点与其对应的转录因子结合位点的相对坐标；
S400 针对每个转录因子，分别统计各相对坐标上甲基化的胞嘧啶-磷酸-鸟嘌呤分子片段占比；
S500 针对每个转录因子，计算转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta，所述侧翼区域对应转录因子结合位点的上下游区域；
S600 将所有转录因子的甲基化占比差值Delta输入预先训练的甲基化水平分析模型，并根据所述甲基化水平分析模型的输出结果完成对待分析血浆样本循环无细胞DNA转录因子结合位点甲基化水平的分析。
另一方面，本发明提供了一种循环无细胞DNA转录因子结合位点的甲基化水平分析装置，包括：
循环无细胞DNA分子片段提取模块，用于接收待分析血浆样本的全基因组甲基化测序数据并从中提取循环无细胞DNA分子片段；
甲基化状态提取模块，用于针对提取的循环无细胞DNA分子片段，进一步提取转录因子结合位点上下游区域胞嘧啶-磷酸-鸟嘌呤位点的甲基化状态；
坐标对齐模块，用于以转录因子结合位点的基因组位置作为参照，将其上下游区域的胞嘧啶-磷酸-鸟嘌呤位点的坐标对齐，得到各胞嘧啶-磷酸-鸟嘌呤位点与其对应的转录因子结合位点的相对坐标；
计算模块，用于针对每个转录因子，分别统计各相对坐标上甲基化的胞嘧啶-磷酸-鸟嘌呤分子片段占比；及用于针对每个转录因子，计算转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta，所述侧翼区域对应转录因子结合位点的上下游区域；
甲基化水平分析模块，用于将所有转录因子的甲基化占比差值Delta输入预先训练的甲基化水平分析模型，并根据所述甲基化水平分析模型的输出结果完成对待分析血浆样本循环无细胞DNA转录因子结合位点甲基化水平的分析。
本发明提供的循环无细胞DNA转录因子结合位点的甲基化水平分析方法及装置，能够通过甲基化数据对转录因子的集合状态进行分析评估，该方法简单易实现，在甲基化水平分析评价中具备高灵敏度和特异性，为后续区分待测血浆样本是否来源于癌症组织提供依据，尤其能够提高某些良性结节、早期癌症患者的检测灵敏度，从而有效辅助癌症的早期诊断以及癌症的早期筛查，提高筛查效率和精度。
附图说明
下面将以明确易懂的方式，结合附图说明优选实施方式，对上述特性、技术特征、优点及其实现方式予以进一步说明。
图1为本发明的实施例中循环无细胞DNA转录因子结合位点的甲基化水平分析方法流程图；
图2为本发明的实施例中提取循环无细胞DNA中分子片段长度与不同长度对应甲基化水平对应图；
图3为本发明实施例中转录因子CTCF结合位点上甲基化水平分布图；
图4为本发明实施例中不同长度循环无细胞DNA分子片段在CTCF位点的甲基化水平分布图；
图5为本发明实施例中转录因子Lyl1的甲基化占比差值Delta在健康人和不同分期癌症患者中的分布箱示图；
图6为本发明实施例中终端设备示意图。
具体实施方式
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对照附图说明本发明的具体实施方式。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，并获得其他的实施方式。
在对本发明实施例进行说明之前，对其中涉及的部分英文缩写进行说明：
FASTQ是一种保存生物序列，通常是核酸序列以及其测序质量得分信息的文本格式。序列与质量得分皆由单个ASCII字符表示。是目前保存高通量测序结果的标准。
Bam是一种用于储存序列数据二进制格式，其中包括序列的比对信息、质量信息以及其他由比对软件保留的序列相关信息。
BisMark是一种比对方法软件，用于查找测序序列在基因参考序列中的位置，可输出Bam格式结果文件。
caret包是一个用于机器学习的集成R包。
本发明的第一种实施例，一种循环无细胞DNA转录因子结合位点的甲基化水平分析方法，如图1所示，该甲基化分析方法包括：
S100 接收待分析血浆样本的全基因组甲基化测序数据并从中提取循环无细胞DNA分子片段。
可以基于一代、二代、三代等测序技术对待分析血浆样本进行测序，得到测序数据。完成测序之后，从中提取循环无细胞DNA分子片段的步骤中包括：S101 对待分析血浆样本的全基因组甲基化测序数据进行质控，去除低质量的测序序列；S102 将质控后余下的测序序列比对至人类参考基因组，并根据比对结果对重复测序序列进行去除；S103 针对去重后的序列文件，提取成对测序读段（reads）的起始位置和终止位置，及测序读段上甲基化位点的甲基化状态，完成对循环无细胞DNA分子片段的提取，即该循环无细胞DNA分子片段中包括甲基化信号及分子片段信息。
具体来说，进行甲基化测序生成FASTQ文件之后，对其中的低质量测序序列进行去除：切去剩余部分开头和结尾处碱基质量低于20的碱基，从reads的5’端开始，以大小为5的窗口进行划窗计算平均质量，如果窗口内平均碱基质量低于20，则切除该窗口，并要求切除后剩余碱基数量超过75。之后，调用BisMark将每一对FASTQ文件作为成对reads与hg19人类参考基因组序列进行比对去重，生成初始Bam文件和比对报告。最后，根据去重后的序列文件，对循环无细胞DNA分子片段进行提取，包括提取成对reads的起始/终止位置及序列上甲基化位点的甲基化状态。
S200 针对提取的循环无细胞DNA分子片段，进一步提取转录因子结合位点上下游区域胞嘧啶-磷酸-鸟嘌呤位点（CpG位点）的甲基化状态。
在CpG提取位点的甲基化状态之前，在现有的已知数据库（如GTRD等数据库）中下载已知的人类转录因子结合位点的基因组坐标，并基于下载的信息对筛选规则进行配置，以对步骤S100提取的循环无细胞DNA分子片段的甲基化测序数据进一步进行筛选得到本实施例中需要的甲基化测序数据。这里，配置于转录因子结合位点基因位置上下分别延伸预设长度基因组距离（对于每一个转录因子结合位点来说，以其中心位点为标志，上下游分别延伸相应数量的碱基）对相应转录因子结合位点上下游区域的甲基化测序数据进行筛选。基于此，对每条循环无细胞DNA进行遍历，得到筛选序列对应基因组位置上CpG位点的甲基化状态。
应当清楚，在实际应用中，可以根据需求对转录因子及延伸的基因组距离的长度进行配置，如，一实例中，为了进一步提升甲基化水平分析的精度，可选定数据库中有1000以上结合位点的转录因子进行分析，并于选定的转录因子结合位点基因位置上下分别延伸3000bp基因组距离（对于每一个转录因子结合位点来说，以其中心位点为标志，上下游分别延伸3000个碱基）。在其他实例中，还可以选定数据库中有500、800等以上结合位点的转录因子进行分析，于选定的转录因子结合位点基因位置上下分别延伸1000bp、2000bp等基因组距离进行筛选，这里不进行具体限定。
S300 以转录因子结合位点的基因组位置作为参照，将其上下游区域的CpG位点的坐标对齐，得到各CpG位点与其对应的转录因子结合位点的相对坐标。
对于同一转录因子来说，其转录因子结合位点的基因组坐标各不相同，是以为了便于数据统计，转录因子结合位点的基因组位置作为参照，将其上下游区域的CpG位点的坐标对齐。该坐标对齐过程为：S301 将所有转录因子结合位点按照其所属转录因子进行分类；S302 针对覆盖同一个转录因子的每个转录因子结合位点，以转录因子结合位点的基因组坐标l bind 为参照，将其上下游区域的各CpG位点的物理坐标l分别减去该转录因子结合位点的基因组坐标l bind 得到CpG位点相对坐标l’，完成各CpG位点相对其所属转录因子结合位点的坐标对齐。
可具体选定转录因子结合位点中心位点的基因组坐标作为参照，将其上下游区域CpG位点中心位点的基因组坐标减去该转录因子结合位点中心位点的基因组坐标完成坐标对齐。在如步骤S200中将转录因子结合位点基因位置上下分别延伸3000bp基因组距离作为筛选规则的实例中，每个转录因子上，坐标对齐后CpG位点的相对位置均在-3000到3000之间。
S400 针对每个转录因子，分别统计各相对坐标上甲基化的胞嘧啶-磷酸-鸟嘌呤分子片段占比。
本实施例将所有转录因子结合位点按照其所属转录因子分类后，后续针对每个转录因子的相关数据进行统计，是以，在坐标对齐后，针对同一转录因子的转录因子结合位点来说，统计过程包括：S401 对每个相对坐标上的CpG位点，计算其甲基化的循环无细胞DNA分子片段数量N m ；S402 对每个相对坐标上的CpG位点，计算其覆盖的循环无细胞DNA分子数量N；S403 根据统计的循环无细胞DNA分子数量N和甲基化的循环无细胞DNA分子片段数量N m ，对各相对坐标上的甲基化占比F m 进行计算。
这一过程中，将同一转录因子所有转录因子结合位点上下游区域的CpG位点坐标对齐后，将其叠加在一起针对各新的坐标点进行叠加统计，即步骤S401中统计的是相对于各转录因子结合位点位置相同的所有CpG位点（坐标对齐后同一相对坐标上的所有CpG位点，相对坐标为0的位点为核心位点）上甲基化的循环无细胞DNA分子片段数量；同样的，步骤S402中统计的也是坐标对齐后同一相对坐标上所有CpG位点覆盖的循环无细胞DNA分子数量。基于此，计算得到每个相对坐标上的甲基化占比，其中，第i个相对坐标上的甲基化占比F im 如式(1)：
Figure DEST_PATH_IMAGE001
(1)
其中，
Figure 14638DEST_PATH_IMAGE002
为第i个相对坐标上胞嘧啶-磷酸-鸟嘌呤位点甲基化的循环无细胞 DNA分子片段数量，
Figure DEST_PATH_IMAGE003
为第i个相对坐标上胞嘧啶-磷酸-鸟嘌呤位点覆盖的循环无细胞 DNA分子片段数量。
S500 针对每个转录因子，计算转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta，侧翼区域对应转录因子结合位点的上下游区域。
针对一转录因子的甲基化占比差值Delta计算过程包括：S501 对该转录因子所有CpG位点的甲基化占比F m 进行平滑处理；S502 对转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta进行计算，如式(2)：
Figure 435254DEST_PATH_IMAGE004
(2)
其中，
Figure DEST_PATH_IMAGE005
表示转录因子结合位点中心点的甲基化占比，
Figure 729270DEST_PATH_IMAGE006
表示转 录因子结合区域侧翼区域的甲基化占比。
在步骤S501中，针对每个转录因子上每个相对位置上CpG位点的甲基化值得到曲 线S使用如GAM等方法进行平滑处理得到曲线S’。在步骤S502中，针对每个转录因子，分别计 算坐标对齐后核心位点和侧翼区域相对位置上的甲基化占比差值Delta，这里的核心位点 为相对位置为0的位点，侧翼区域即坐标对齐前基于转录因子结合点位上下延伸的基因组 距离对应的区域，即转录因子结合区域侧翼区域的甲基化占比
Figure DEST_PATH_IMAGE007
为侧翼区域上各 相对位点上CpG位点的甲基化占比F m ，针对一转录因子得到的甲基化占比差值Delta为转录 因子结合位点中心点的甲基化占比
Figure 354942DEST_PATH_IMAGE008
分别减去侧翼区域各相对坐标的甲基化占比
Figure DEST_PATH_IMAGE009
生成的一组数据。
S600 将所有转录因子的甲基化占比差值Delta输入预先训练的甲基化水平分析模型，并根据甲基化水平分析模型的输出结果完成对待分析血浆样本循环无细胞DNA转录因子结合位点甲基化水平的分析。
基于待分析血浆样本中转录因子的数量，步骤S500中得到相应数量的甲基化占比差值Delta，以此将其输入预先训练的甲基化水平分析模型中进行甲基化水平分析。在对其甲基化水平进行分析中，将甲基化水平分析模型输出的概率与预先设定的概率阈值进行比较，完成对所述待测血浆样本的转录因子结合状态进行评价。当模型输出的概率大于等于预先设定的概率阈值，表明待分析待测血浆样本可能来源于肿瘤患者血浆样本；当计算得到的概率小于预先设定的概率阈值，表明该待测血浆样本可能来源于健康人血浆样本。在此基础上，在后续诊断过程中可以辅助医生进行综合判断，为诊断结果提供部分依据，辅助癌症筛查工作，尤其是早期癌症的诊断和筛查。概率阈值的确定，由阈值确定过程中区分不同类型患者过程中的特异性确定。
转录因子区域甲基化水平模型可以采用支持向量机（SVM）、决策树（Randomforest）、随机森林（DT-tree）等方式进行建模，为了提高模型的分类精度，在将所有转录因子甲基化占比差值Delta输入该甲基化水平分析模型之前，还可采用非参数检验等方法进行特征筛选。
在使用本实施例中的循环无细胞DNA转录因子结合位点的甲基化水平分析方法对待分析血浆样本进行分析之前，还应包括对甲基化水平分析模型进行建模及训练的步骤。训练过程中，从数据库中采集大量的健康人血浆样本和癌症患者血浆样本，按照7:3的比例分成训练集和验证集；之后，将这些样本进行步骤S100~S500的处理步骤之后，将训练集中样本的甲基化占比差值Delta输入创建的甲基化水平分析模型中对其进行分类，并使用验证集中样本的甲基化占比差值Delta对模型进行验证。这一过程中，根据训练结果中对区分不同类型患者过程中的特异性对概率阈值进行确定，如一实例中，根据训练集样本的甲基化水平分析模型预测得到每一个样本属于癌症样本的概率值，基于训练集样本已有的样本类型（癌症或健康人）与预测得到的概率值进行计算，选定98%特异性下的概率阈值作为判断标准。
本发明的另一实施例，一种循环无细胞DNA转录因子结合位点的甲基化水平分析装置，包括：循环无细胞DNA分子片段提取模块，用于接收待分析血浆样本的全基因组甲基化测序数据并从中提取循环无细胞DNA分子片段；甲基化状态提取模块，用于针对提取的循环无细胞DNA分子片段，进一步提取转录因子结合位点上下游区域CpG位点的甲基化状态；坐标对齐模块，用于以转录因子结合位点的基因组位置作为参照，将其上下游区域的CpG位点的坐标对齐，得到各CpG位点与其对应的转录因子结合位点的相对坐标；计算模块，用于针对每个转录因子，分别统计各相对坐标上甲基化的胞嘧啶-磷酸-鸟嘌呤分子片段占比；及用于针对每个转录因子，计算转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta，侧翼区域对应转录因子结合位点的上下游区域；甲基化水平分析模块，用于将所有转录因子的甲基化占比差值Delta输入预先训练的甲基化水平分析模型，并根据甲基化水平分析模型的输出结果完成对待分析血浆样本循环无细胞DNA转录因子结合位点甲基化水平的分析。
循环无细胞DNA分子片段提取模块中可以基于一代、二代、三代等测序技术对待分析血浆样本进行测序，得到测序数据。完成测序之后，进一步从中提取循环无细胞DNA分子片段，提取的过程包括：对待分析血浆样本的全基因组甲基化测序数据进行质控，去除低质量的测序序列；将质控后余下的测序序列比对至人类参考基因组，并根据比对结果对重复测序序列进行去除；针对去重后的序列文件，提取成对测序读段（reads）的起始位置和终止位置，及测序读段上甲基化位点的甲基化状态，完成对循环无细胞DNA分子片段的提取，即该循环无细胞DNA分子片段中包括甲基化信号及分子片段信息。
具体来说，进行甲基化测序生成FASTQ文件之后，对其中的低质量测序序列进行去除：切去剩余部分开头和结尾处碱基质量低于20的碱基，从reads的5’端开始，以大小为5的窗口进行划窗计算平均质量，如果窗口内平均碱基质量低于20，则切除该窗口，并要求切除后剩余碱基数量超过75。之后，调用BisMark将每一对FASTQ文件作为成对reads与hg19人类参考基因组序列进行比对去重，生成初始Bam文件和比对报告。最后，根据去重后的序列文件，对循环无细胞DNA分子片段进行提取，包括提取成对reads的起始/终止位置及序列上甲基化位点的甲基化状态。
甲基化状态提取模块中包括：筛选单元，用于根据选定的转录因子及预先设定的筛选规则遍历每条循环无细胞DNA分子片段的甲基化测序数据进行筛选得到筛选序列，其中，选定的转录因子为结合位点数量在1000以上的转录因子，预先设定的筛选规则为：基于选定的转录因子结合位点的基因位置筛选其上下分别延伸预设长度基因组距离的测序序列；甲基化状态提取单元，用于根据筛选得到的筛选序列，进一步提取对应基因组位置上CpG位点的甲基化状态。
在CpG提取位点的甲基化状态之前，在现有的已知数据库（如GTRD等数据库）中下载已知的人类转录因子结合位点的基因组坐标，并基于下载的信息对筛选规则进行配置，以提取的循环无细胞DNA分子片段的甲基化测序数据进一步进行筛选得到本实施例中需要的甲基化测序数据。这里，配置于转录因子结合位点基因位置上下分别延伸预设长度基因组距离对相应转录因子结合位点上下游区域的甲基化测序数据进行筛选。基于此，对每条循环无细胞DNA进行遍历，得到筛选序列对应基因组位置上CpG位点的甲基化状态。
应当清楚，在实际应用中，可以根据需求对转录因子及延伸的基因组距离的长度进行配置，如，一实例中，为了进一步提升甲基化水平分析的精度，可选定数据库中有1000以上结合位点的转录因子进行分析，并于选定的转录因子结合位点基因位置上下分别延伸3000bp基因组距离（对于每一个转录因子结合位点来说，以其中心位点为标志，上下游分别延伸3000个碱基）。在其他实例中，还可以选定数据库中有500、800等以上结合位点的转录因子进行分析，于选定的转录因子结合位点基因位置上下分别延伸1000bp、2000bp等基因组距离进行筛选，这里不进行具体限定。
坐标对齐模块中包括转录因子分类单元，用于将所有转录因子结合位点按照其所属转录因子进行分类；坐标对齐单元，用于针对覆盖同一个转录因子的每个转录因子结合位点，以转录因子结合位点的基因组坐标l bind 为参照，将其上下游区域的各CpG位点的物理坐标l分别减去该转录因子结合位点的基因组坐标l bind 得到CpG位点相对坐标l’，完成各CpG位点相对其所属转录因子结合位点的坐标对齐。
对于同一转录因子来说，其转录因子结合位点的基因组坐标各不相同，是以为了便于数据统计，转录因子结合位点的基因组位置作为参照，将其上下游区域的CpG位点的坐标对齐。在坐标对齐中，首先，转录因子分类单元将所有转录因子结合位点按照其所属转录因子进行分类；之后，坐标对齐单元针对覆盖同一个转录因子的每个转录因子结合位点，以转录因子结合位点的基因组坐标l bind 为参照，将其上下游区域的各CpG位点的物理坐标l分别减去该转录因子结合位点的基因组坐标l bind 得到CpG位点相对坐标l’，完成各CpG位点相对其所属转录因子结合位点的坐标对齐。
可具体选定转录因子结合位点中心位点的基因组坐标作为参照，将其上下游区域CpG位点中心位点的基因组坐标减去该转录因子结合位点中心位点的基因组坐标完成坐标对齐。坐标对齐单元将转录因子结合位点基因位置上下分别延伸3000bp基因组距离作为筛选规则的实施例中，每个转录因子上，坐标对齐后CpG位点的相对位置均在-3000到3000之间。
计算模块包括统计单元，用于针对同一转录因子，对每个相对坐标上的CpG位点，计算其甲基化的循环无细胞DNA分子片段数量N m ；及对每个相对坐标上的CpG位点，计算其覆盖的循环无细胞DNA分子数量N；第一计算单元，用于针对同一转录因子，根据统计的循环无细胞DNA分子数量N和甲基化的循环无细胞DNA分子片段数量N m ，对各相对坐标上的甲基化占比F m 进行计算。
这一过程中，将同一转录因子所有转录因子结合位点上下游区域的CpG位点坐标对齐后，将其叠加在一起针对各新的坐标点进行叠加统计，即统计单元统计的是相对于各转录因子结合位点位置相同的所有CpG位点（坐标对齐后同一相对坐标上的所有CpG位点，相对坐标为0的位点为核心位点）上甲基化的循环无细胞DNA分子片段数量；同样的，之后统计的也是坐标对齐后同一相对坐标上所有CpG位点覆盖的循环无细胞DNA分子数量。基于此，计算得到每个相对坐标上的甲基化占比，其中，第i个相对坐标上的甲基化占比F im 如式(1)。
计算模块还包括平滑处理单元，用于针对同一转录因子的所有CpG位点的甲基化占比F m 进行平滑处理；第二计算单元，用于针对同一转录因子的转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta进行计算，如式(2)。
在平滑处理单元中，针对每个转录因子上每个相对位置上CpG位点的甲基化值得 到曲线S使用如GAM等方法进行平滑处理得到曲线S’。第二计算单元针对每个转录因子，分 别计算坐标对齐后核心位点和侧翼区域相对位置上的甲基化占比差值Delta，这里的核心 位点为相对位置为0的位点，侧翼区域即坐标对齐前基于转录因子结合点位上下延伸的基 因组距离对应的区域，即转录因子结合区域侧翼区域的甲基化占比
Figure 720937DEST_PATH_IMAGE010
为侧翼区域 上各相对位点上CpG位点的甲基化占比F m ，针对一转录因子得到的甲基化占比差值Delta为 转录因子结合位点中心点的甲基化占比
Figure DEST_PATH_IMAGE011
分别减去侧翼区域各相对坐标的甲基化 占比
Figure 355312DEST_PATH_IMAGE012
生成的一组数据。
基于待分析血浆样本中转录因子的数量，计算模块得到相应数量的甲基化占比差值Delta，以此将其输入预先训练的甲基化水平分析模型中进行甲基化水平分析。在对其甲基化水平进行分析中，将甲基化水平分析模型输出的概率与预先设定的概率阈值进行比较，完成对所述待测血浆样本的转录因子结合状态进行评价。当模型输出的概率大于等于预先设定的概率阈值，表明待分析血浆样本的甲基化程度高，该待测血浆样本可能来源于肿瘤患者血浆样本；当计算得到的概率小于预先设定的概率阈值，表明待分析血浆样本的甲基化程度低，该待测血浆样本可能来源于健康人血浆样本。在此基础上，在后续诊断过程中可以辅助医生进行综合判断，为诊断结果提供部分依据，辅助癌症筛查工作，尤其是早期癌症的诊断和筛查。概率阈值的确定，由阈值确定过程中区分不同类型患者过程中的特异性确定。
甲基化水平分析模型可以采用支持向量机（SVM）、决策树（Randomforest）、随机森林（DT-tree）等方式进行建模，为了提高模型的分类精度，在将甲基化占比差值Delta输入该甲基化水平分析模型之前，还可采用非参数检验等方法进行特征筛选。
在使用本实施例中的循环无细胞DNA转录因子结合位点的甲基化水平分析方法对待分析血浆样本进行分析之前，还应包括对甲基化水平分析模型进行建模及训练的步骤。训练过程中，从数据库中采集大量的健康人血浆样本和癌症患者血浆样本，按照7:3的比例分成训练集和验证集；之后，将这些样本进行步骤S100~S500的处理步骤之后，将训练集中样本的甲基化占比差值Delta输入创建的甲基化水平分析模型中对其进行分类，并使用验证集中样本的甲基化占比差值Delta对模型进行验证。这一过程中，根据训练结果中对区分不同类型患者过程中的特异性对概率阈值进行确定，如一实例中，根据训练集样本的甲基化水平分析模型预测得到每一个样本属于癌症样本的概率值，基于训练集样本已有的样本类型（癌症或健康人）与预测得到的概率值进行计算，选定98%特异性下的概率阈值作为判断标准。
以下通过一实例对本发明的循环无细胞DNA转录因子结合位点的甲基化水平分析方法/装置及其有益效果进行说明。
1.1 样本收集
参照表1为收集到的样本，其中训练集的样本总数为913，包括健康人血浆样本352份，癌症血浆样本561份，癌症血浆样本中包含的癌种为BRCA、COREAD、ESCA、GAC、HCC、NSCLC、OV及PAAD；测试集的样本总数为385，包括健康人血浆样本145份，癌症血浆样本240份。
表1 样本信息
Figure DEST_PATH_IMAGE013
1.2 测序处理
1.2.1 cfDNA提取
血浆cfDNA的分离血浆样品从采集的10 ml外周血中分离，储存在cfDNA储存管(Cwbiotech)中。血液在4℃下1600 g离心10 min，血浆转移到新管中。第二个离心机在4℃下进行15 min，12000 rpm，以清除剩余的细胞碎片，获得4 ml血浆，并在-80℃保存直到使用。使用MagMAX Cell-Free DNA分离试剂盒(Thermo Fisher Scientific)根据制造商说明书提取cfDNA。采用生物分析仪2100 (Agilent)检测cfDNA的数量和质量。
1.2.2 循环无细胞DNA甲基化测序
对5~30 ng的cfDNA使用NEBNext建立DNA甲基化建库，按照(New EnglandBiolabs)根据制造商说明。9个周期PCR扩增文库，使用Qubit dsDNA HS分析试剂盒(ThermoFisher Scientific)进行定量分析。文库测序使用Illumina公司的NovaSeq 6000(Illumina)进行，测序为双端测序，reads长度为100 bp。
1.3 数据预处理
对于测序得到FASTQ文件，首先使用Trimmomatic-0.36对测序数据中低质量序列以及残余的接头进行处理。其次再使用Bismark(v0.23.1)将序列比对到人类参考基因组上，以及对比对后的文件进行去重。再次对序列进行质控处理，去掉比对质量较低的序列。
1.4 甲基化状态计算
对于上一步获得的序列比对文件，首先，对其中每对reads根据其序列名称、起始终止位置以及是否比对到基因组正链上，将reads组装成分子片段，即循环无细胞DNA分子片段。之后，根据比对文件中的XM:Z:标签对位置上的CpG位点甲基化状态进行判断，其中，Z代表此位点甲基化，z代表此位点未甲基化。对于全基因组范围内的分子片段，得到的分子片段长度与其甲基化水平的对应关系如图2所示，其中，横轴为分子片段长度；下方曲线a为分子片段长度分布密度图，对应右侧坐标轴（分子片段长度分布密度）；上方曲线b为不同长度的分子片段序列中甲基化CpG位点占序列上所有CpG位点数量比例，对应左侧坐标轴（不同长度分子片段长度甲基化CpG占比）。正常包裹在组蛋白八聚体上的分子片段长度为167bp，分子片段在120bp-350bp之间的分子片段上甲基化CpG位点占比相对较高，说明包裹在核小体上的cfDNA相对于较短的cfDNA有更高的甲基化水平。循环无细胞DNA测序数据中甲基化水平和核小体占位相关，而核小体占位会影响转录因子的结合。说明甲基化测序可以评估转录因子的可及性。
1.5 单个转录因子结合位点甲基化状态
对于覆盖同一个转录因子的转录因子结合位点，对其以每个转录因子结合位点的基因组坐标l bind 为参照，用每个CpG位点的物理坐标l减去对应转录因子结合位点的基因组坐标l bind 得到相对坐标l’。一实例中，转录因子CTCF结合位点及其上下游区域甲基化水平分布图如图3所示。
对每个相对坐标上的CpG位点，计算其甲基化的循环无细胞DNA分子片段数量Nm， 及计算其覆盖的循环无细胞DNA分子数量N；进而计算各相对坐标上的甲基化占比
Figure 441735DEST_PATH_IMAGE014
。
分别对各转录因子所有位点上的甲基化占比F m进行平滑处理，得到平滑后每个位点上的F m ’。之后，对转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta进行计算。
1.6 使用甲基化占比差值Delta区分样本来源
对于步骤1.5中计算所得的504个转录因子的甲基化占比差值Delta。使用不同长度的分子片段长度计算转录因子结合位点上下游3000bp区间相对位置CpG位点的甲基化水平。图4为不同长度循环无细胞DNA分子片段在CTCF位点的甲基化水平分布图，图(a)~(d)中循环无细胞DNA分子片段的长度分别为120bp、170bp、250bp和400bp，从图中可以看出，当用于计算的分子片段长度变短时，侧翼序列甲基化水平有较为明显的下降，计算的甲基化占比差值Delta也会相应变小，可见，甲基化占比差值Delta与用于计算的循环无细胞DNA分子片段长度之间具备相关性。
对于队列中801例癌症患者样本以及497例健康人样本，使用秩和检验以及AUC(受试者曲线下面积)对每一个转录因子的甲基化占比差值Delta的分类性能进行评估。其中，有434个转录因子秩和检验p值小于0.05，414个转录因子秩和检验p值小于0.01，387个转录因子秩和检验p值小于0.001。504个转录因子中甲基化占比差值Delta在癌症和健康人血浆样本之间的分类性能较好的10个转录因子的AUC以及秩和检验P值如表2所示。转录因子Lyl1的甲基化占比差值Delta在健康人和晚期癌症患者中的分布箱示图如图5所示，其中，横坐标分别表示健康人样本、癌症I期患者样本、癌症II期患者样本、癌症III期患者样本、癌症IV期患者样本及分期不明样本，从图中可以看出，健康人样本的甲基化占比差值Delta较癌症患者样本的高，随着癌症从早期到晚期（对应图示中的I期到IV期），甲基化占比差值Delta逐步下降，晚期癌症患者相比早期癌症患者血液中有更多来自肿瘤组织的ctDNA，信号也越明显。说明通过计算甲基化占比差值Delta对循环无细胞DNA转录因子结合位点的甲基化水平进行分析能够提高某些良性结节、早期癌症患者的检测灵敏度，从而有效辅助癌症的早期诊断以及癌症的早期筛查，提高筛查效率和精度。
表2 10个转录因子的AUC以及秩和检验P值
Figure DEST_PATH_IMAGE015
1.7 训练模型
对于1.5中获得504个转录因子的甲基化占比差值Delta。对于训练集合的数据，首先使用caret包中的preprocess模块（预处理模块，用于对数据进行标准化处理）对数据进行处理，将每个转录因子的甲基化占比差值Delta缩放到相同范围内。
分别选择cforest（条件随机森林）、svmLinear（线性核支持向量机）、LogitBoost（一种使用逻辑回归方法作为损失函数的梯度下降算法）、gbm（Gradient BoostingMachine，梯度提升机）及Ensemble（集成学习）模型进行建模分析，使用召回率作为参数选择标准，在训练集中使用重复的10次交叉验证方法进重采样，获得每个建模方法的最优参数以及训练集合中每个样本属于癌症或者健康人的概率。对于测试集合，使用和训练集合相同的参数进行数据的缩放，得到与训练集合相同范围的数据，使用上述训练集合中得到模型，对测试集合预处理得到的数据进行预测，得到每个样本属于癌症或者健康人来源的概率值。对于测试集合和训练集合，按照队列的临床信息以及模型预测出的概率值，计算如表3所示的模型AUC（受试者曲线下面积），评估模型分类性能。以训练集中特异性等于0.98时的概率值作为判断样本来源于癌症还是健康人的阈值，在cforest模型中概率值大于0.583836084795837为来源于癌症的样本；在svmLinear模型中概率值大于0.420882507176859为来源于癌症的样本；在LogitBoost模型中概率值大于0.952574126822433为来源于癌症的样本；在gbm模型中概率值大于0.760563706382831为来源于癌症的样本；在Ensemble模型中概率值大于0.648726690600909为来源于癌症的样本。根据分析结果可知，该方法具备高灵敏度和特异性。
表3 不同模型训练得到AUC
Figure 372346DEST_PATH_IMAGE016
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各程序模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的程序模块完成，即将装置的内部结构划分成不同的程序单元或模块，以完成以上描述的全部或者部分功能。实施例中的各程序模块可以集成在一个处理单元中，也可是各个单元单独物理存在，也可以两个或两个以上单元集成在一个处理单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件程序单元的形式实现。另外，各程序模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。
图6是本发明一个实施例中提供的终端设备的结构示意图，如所示，该终端设备200包括：处理器220、存储器210以及存储在存储器210中并可在处理器220上运行的计算机程序211，例如：循环无细胞DNA转录因子结合位点的甲基化水平分析关联程序。处理器220执行计算机程序211时实现上述各循环无细胞DNA转录因子结合位点的甲基化水平分析方法实施例中的步骤，或者，处理器220执行计算机程序211时实现上述循环无细胞DNA转录因子结合位点的甲基化水平分析装置实施例中各模块的功能。
终端设备200可以为笔记本、掌上电脑、平板型计算机、手机等设备。终端设备200可包括，但不仅限于处理器220、存储器210。本领域技术人员可以理解，图6仅仅是终端设备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 (8)

1.一种循环无细胞DNA转录因子结合位点的甲基化水平分析方法，其特征在于，包括：
S100 接收待分析血浆样本的全基因组甲基化测序数据并从中提取循环无细胞DNA分子片段；
S200 针对提取的循环无细胞DNA分子片段，进一步提取转录因子结合位点上下游区域胞嘧啶-磷酸-鸟嘌呤位点的甲基化状态；
S300 以转录因子结合位点的基因组位置作为参照，将其上下游区域的胞嘧啶-磷酸-鸟嘌呤位点的坐标对齐，得到各胞嘧啶-磷酸-鸟嘌呤位点与其对应的转录因子结合位点的相对坐标；
S400 针对每个转录因子，分别统计各相对坐标上甲基化的胞嘧啶-磷酸-鸟嘌呤分子片段占比；
S500 针对每个转录因子，计算转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta，所述侧翼区域对应转录因子结合位点的上下游区域；
S600 将所有转录因子的甲基化占比差值Delta输入预先训练的甲基化水平分析模型，并根据所述甲基化水平分析模型的输出结果完成对待分析血浆样本循环无细胞DNA转录因子结合位点甲基化水平的分析；
步骤S400针对每个转录因子，分别统计各相对坐标上甲基化的胞嘧啶-磷酸-鸟嘌呤分子片段占比中，针对一转录因子的统计过程包括：
S401 对每个相对坐标上的胞嘧啶-磷酸-鸟嘌呤位点，计算其甲基化的循环无细胞DNA分子片段数量N m ；
S402 对每个相对坐标上的胞嘧啶-磷酸-鸟嘌呤位点，计算其覆盖的循环无细胞DNA分子数量N；
S403 根据统计的循环无细胞DNA分子数量N和甲基化的循环无细胞DNA分子片段数量N m ，对各相对坐标上的甲基化占比F m 进行计算；其中，第i个相对坐标上的甲基化占比F im 为：
Figure 535748DEST_PATH_IMAGE001
其中，
Figure 7881DEST_PATH_IMAGE002
为第i个相对坐标上胞嘧啶-磷酸-鸟嘌呤位点甲基化的循环无细胞DNA分 子片段数量，
Figure 771306DEST_PATH_IMAGE003
为第i个相对坐标上胞嘧啶-磷酸-鸟嘌呤位点覆盖的循环无细胞DNA分子 片段数量；
步骤S500针对每个转录因子，计算转录因子结合位点中心点与周边侧翼区域各相对坐标的甲基化占比差值Delta中，针对一转录因子的甲基化占比差值Delta计算过程包括：
S501 对该转录因子所有胞嘧啶-磷酸-鸟嘌呤位点的甲基化占比F m 进行平滑处理；
S502 对转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta进行计算：
Figure 930630DEST_PATH_IMAGE004
其中，
Figure 78758DEST_PATH_IMAGE005
表示转录因子结合位点中心点的甲基化占比，
Figure 196493DEST_PATH_IMAGE006
表示转录 因子结合区域侧翼区域的甲基化占比。
2.如权利要求1所述的甲基化水平分析方法，其特征在于，步骤S200针对提取的循环无细胞DNA分子片段，进一步提取转录因子结合位点上下游区域胞嘧啶-磷酸-鸟嘌呤位点的甲基化状态中包括：
S201 根据选定的转录因子及预先设定的筛选规则遍历每条循环无细胞DNA分子片段的甲基化测序数据进行筛选得到筛选序列，所述预先设定的筛选规则为：基于选定的转录因子结合位点的基因位置筛选其上下分别延伸预设长度基因组距离的测序序列；
S202 根据筛选得到的筛选序列，进一步提取对应基因组位置上胞嘧啶-磷酸-鸟嘌呤位点的甲基化状态。
3.如权利要求1或2所述的甲基化水平分析方法，其特征在于，步骤S300以转录因子结合位点的基因组位置作为参照，将其上下游区域的胞嘧啶-磷酸-鸟嘌呤位点的坐标对齐，得到各胞嘧啶-磷酸-鸟嘌呤位点与其对应的转录因子结合位点的相对坐标中包括：
S301 将所有转录因子结合位点按照其所属转录因子进行分类；
S302 针对覆盖同一个转录因子的每个转录因子结合位点，以转录因子结合位点的基因组坐标l bind 为参照，将其上下游区域的各胞嘧啶-磷酸-鸟嘌呤位点的物理坐标l分别减去该转录因子结合位点的基因组坐标l bind 得到胞嘧啶-磷酸-鸟嘌呤位点相对坐标l’，完成各胞嘧啶-磷酸-鸟嘌呤位点相对其所属转录因子结合位点的坐标对齐。
4.一种循环无细胞DNA转录因子结合位点的甲基化水平分析装置，其特征在于，包括：
循环无细胞DNA分子片段提取模块，用于接收待分析血浆样本的全基因组甲基化测序数据并从中提取循环无细胞DNA分子片段；
甲基化状态提取模块，用于针对提取的循环无细胞DNA分子片段，进一步提取转录因子结合位点上下游区域胞嘧啶-磷酸-鸟嘌呤位点的甲基化状态；
坐标对齐模块，用于以转录因子结合位点的基因组位置作为参照，将其上下游区域的胞嘧啶-磷酸-鸟嘌呤位点的坐标对齐，得到各胞嘧啶-磷酸-鸟嘌呤位点与其对应的转录因子结合位点的相对坐标；
计算模块，用于针对每个转录因子，分别统计各相对坐标上甲基化的胞嘧啶-磷酸-鸟嘌呤分子片段占比；及用于针对每个转录因子，计算转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta，所述侧翼区域对应转录因子结合位点的上下游区域；
甲基化水平分析模块，用于将所有转录因子的甲基化占比差值Delta输入预先训练的甲基化水平分析模型，并根据所述甲基化水平分析模型的输出结果完成对待分析血浆样本循环无细胞DNA转录因子结合位点甲基化水平的分析；
所述计算模块中包括：
统计单元，用于针对同一转录因子，对每个相对坐标上的胞嘧啶-磷酸-鸟嘌呤位点，计算其甲基化的循环无细胞DNA分子片段数量N m ；及对每个相对坐标上的胞嘧啶-磷酸-鸟嘌呤位点，计算其覆盖的循环无细胞DNA分子数量N；
第一计算单元，用于针对同一转录因子，根据统计的循环无细胞DNA分子数量N和甲基化的循环无细胞DNA分子片段数量N m ，对各相对坐标上的甲基化占比F m 进行计算；其中，第i个相对坐标上的甲基化占比F im 为：
Figure 635432DEST_PATH_IMAGE007
其中，
Figure 365490DEST_PATH_IMAGE008
为第i个相对坐标上胞嘧啶-磷酸-鸟嘌呤位点甲基化的循环无细胞DNA分 子片段数量，
Figure 532904DEST_PATH_IMAGE009
为第i个相对坐标上胞嘧啶-磷酸-鸟嘌呤位点覆盖的循环无细胞DNA分子 片段数量；
所述计算模块中还包括：
平滑处理单元，用于针对同一转录因子的所有胞嘧啶-磷酸-鸟嘌呤位点的甲基化占比F m 进行平滑处理；
第二计算单元，用于针对同一转录因子的转录因子结合位点中心点与侧翼区域各相对坐标的甲基化占比差值Delta进行计算：
Figure 952034DEST_PATH_IMAGE010
其中，
Figure 673740DEST_PATH_IMAGE011
表示转录因子结合位点中心点的甲基化占比，
Figure 795106DEST_PATH_IMAGE012
表示转录因 子结合区域侧翼区域的甲基化占比。
5.如权利要求4所述的甲基化水平分析装置，其特征在于，所述甲基化状态提取模块中包括：
筛选单元，用于根据选定的转录因子及预先设定的筛选规则遍历每条循环无细胞DNA分子片段的甲基化测序数据进行筛选得到筛选序列，其中，所述选定的转录因子为结合位点数量在1000以上的转录因子，所述预先设定的筛选规则为：基于选定的转录因子结合位点的基因位置筛选其上下分别延伸预设长度基因组距离的测序序列；
甲基化状态提取单元，用于根据筛选得到的筛选序列，进一步提取对应基因组位置上胞嘧啶-磷酸-鸟嘌呤位点的甲基化状态。
6.如权利要求4或5所述的甲基化水平分析装置，其特征在于，所述坐标对齐模块中包括：
转录因子分类单元，用于将所有转录因子结合位点按照其所属转录因子进行分类；
坐标对齐单元，用于针对覆盖同一个转录因子的每个转录因子结合位点，以转录因子结合位点的基因组坐标l bind 为参照，将其上下游区域的各胞嘧啶-磷酸-鸟嘌呤位点的物理坐标l分别减去该转录因子结合位点的基因组坐标l bind 得到胞嘧啶-磷酸-鸟嘌呤位点相对坐标l’，完成各胞嘧啶-磷酸-鸟嘌呤位点相对其所属转录因子结合位点的坐标对齐。
7.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质中存储有计算机指令，所述计算机指令能够被处理器执行以实现如权利要求1-3任意一项所述的循环无细胞DNA转录因子结合位点的甲基化水平分析方法的步骤。
8.一种终端设备，其特征在于，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器运行所述计算机程序时实现如权利要求1-3任意一项所述的循环无细胞DNA转录因子结合位点的甲基化水平分析方法的步骤。
CN202211059714.0A 2022-09-01 2022-09-01 循环无细胞dna转录因子结合位点的甲基化水平分析方法及装置 Active CN115132274B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211059714.0A CN115132274B (zh) 2022-09-01 2022-09-01 循环无细胞dna转录因子结合位点的甲基化水平分析方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211059714.0A CN115132274B (zh) 2022-09-01 2022-09-01 循环无细胞dna转录因子结合位点的甲基化水平分析方法及装置

Publications (2)

Publication Number Publication Date
CN115132274A CN115132274A (zh) 2022-09-30
CN115132274B true CN115132274B (zh) 2022-11-25

Family

ID=83387434

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211059714.0A Active CN115132274B (zh) 2022-09-01 2022-09-01 循环无细胞dna转录因子结合位点的甲基化水平分析方法及装置

Country Status (1)

Country Link
CN (1) CN115132274B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115323058A (zh) * 2022-10-17 2022-11-11 深圳市睿法生物科技有限公司 一种基于ctDNA甲基化模式的癌种定位方法
CN117831623A (zh) * 2024-03-04 2024-04-05 阿里巴巴（中国）有限公司 对象检测方法、对象检测模型训练方法、转录因子结合位点检测方法、目标对象处理方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112823213A (zh) * 2019-05-31 2021-05-18 福瑞诺姆控股公司 用于甲基化核酸的高深度测序的方法和系统
CN113215264A (zh) * 2021-07-07 2021-08-06 上海伯豪生物技术有限公司 用于子宫内膜癌早期筛查人外周血循环肿瘤dna中tmem101基因甲基化检测试剂盒

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AU2018266964A1 (en) * 2017-05-12 2019-12-05 Dana-Farber Cancer Institute, Inc. Universal early cancer diagnostics
CN108949970A (zh) * 2017-05-23 2018-12-07 中国科学院深圳先进技术研究院 基于多组学的宫颈癌特征获取方法和系统
US20220033880A1 (en) * 2018-10-17 2022-02-03 The University Of Queensland Epigenetic biomarker and uses therefor
CN110157804A (zh) * 2019-04-04 2019-08-23 广州优泽生物技术有限公司 用于肺癌诊断、疗效预测或预后的甲基化位点、检测引物及试剂盒
CN111575306B (zh) * 2020-06-11 2022-08-05 南方医科大学 一种激活红细胞中γ珠蛋白基因表达的方法
CN112735531B (zh) * 2021-03-30 2021-07-02 臻和(北京)生物科技有限公司 循环无细胞核小体活性区域的甲基化分析方法和装置、终端设备及存储介质

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112823213A (zh) * 2019-05-31 2021-05-18 福瑞诺姆控股公司 用于甲基化核酸的高深度测序的方法和系统
CN113215264A (zh) * 2021-07-07 2021-08-06 上海伯豪生物技术有限公司 用于子宫内膜癌早期筛查人外周血循环肿瘤dna中tmem101基因甲基化检测试剂盒

Also Published As

Publication number Publication date
CN115132274A (zh) 2022-09-30

Similar Documents

Publication Publication Date Title
CN115132274B (zh) 循环无细胞dna转录因子结合位点的甲基化水平分析方法及装置
CN113257350B (zh) 基于液体活检的ctDNA突变程度分析方法和装置、ctDNA性能分析装置
CN112951418B (zh) 基于液体活检的连锁区域甲基化评估方法和装置、终端设备及存储介质
CN112086129B (zh) 预测肿瘤组织cfDNA的方法及系统
CN109767810B (zh) 高通量测序数据分析方法及装置
CN112951327B (zh) 药物敏感预测方法、电子设备及计算机可读存储介质
CN113838533B (zh) 一种癌症检测模型及其构建方法和试剂盒
CN109830264B (zh) 肿瘤患者基于甲基化位点进行分类的方法
US20220336043A1 (en) cfDNA CLASSIFICATION METHOD, APPARATUS AND APPLICATION
Lee et al. A deep learning and similarity-based hierarchical clustering approach for pathological stage prediction of papillary renal cell carcinoma
CN117275585A (zh) 基于lp-wgs和dna甲基化的肺癌早筛模型构建方法及电子设备
CN112037863B (zh) 一种早期nsclc预后预测系统
CN111951893B (zh) 构建肿瘤突变负荷tmb面板的方法
CN112837748A (zh) 一种区分不同解剖学起源肿瘤的系统及其方法
CN110570951A (zh) 构建乳腺癌新辅助化疗疗效分类模型的方法
CN113793638B (zh) 一种同源重组修复基因变异的解读方法
WO2018209704A1 (zh) 基于dna测序数据的样本来源检测方法、装置和存储介质
CN116762132A (zh) 基于游离dna的疾病预测模型及其构建方法和应用
US11535896B2 (en) Method for analysing cell-free nucleic acids
CN113160895A (zh) 一种结直肠癌风险评估模型及系统
CN113159529A (zh) 一种肠道息肉的风险评估模型及相关系统
CN113643759B (zh) 基于液体活检的染色体稳定性评价方法和装置、终端设备及存储介质
Durmaz et al. A multimodal analysis of genomic and RNA splicing features in myeloid malignancies
CN117035113A (zh) 结肠癌预后风险模型构建方法、装置、设备及存储介质
CN117577182B (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