CN116153418B - 校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质 - Google Patents
校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质 Download PDFInfo
- Publication number
- CN116153418B CN116153418B CN202310413887.6A CN202310413887A CN116153418B CN 116153418 B CN116153418 B CN 116153418B CN 202310413887 A CN202310413887 A CN 202310413887A CN 116153418 B CN116153418 B CN 116153418B
- Authority
- CN
- China
- Prior art keywords
- methylation
- interval
- sample
- value
- window
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 230000000694 effects Effects 0.000 title claims abstract description 44
- 238000012164 methylation sequencing Methods 0.000 title claims abstract description 40
- 230000011987 methylation Effects 0.000 claims abstract description 224
- 238000007069 methylation reaction Methods 0.000 claims abstract description 224
- 206010028980 Neoplasm Diseases 0.000 claims abstract description 70
- 238000012937 correction Methods 0.000 claims abstract description 34
- 239000000523 sample Substances 0.000 claims description 180
- 238000012216 screening Methods 0.000 claims description 86
- 239000012634 fragment Substances 0.000 claims description 81
- 238000006243 chemical reaction Methods 0.000 claims description 25
- 238000004458 analytical method Methods 0.000 claims description 22
- 239000013610 patient sample Substances 0.000 claims description 22
- 238000012417 linear regression Methods 0.000 claims description 19
- 238000012360 testing method Methods 0.000 claims description 15
- 210000000349 chromosome Anatomy 0.000 claims description 14
- 238000003379 elimination reaction Methods 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 13
- 230000008030 elimination Effects 0.000 claims description 12
- 238000012163 sequencing technique Methods 0.000 claims description 12
- 108091029430 CpG site Proteins 0.000 claims description 8
- 230000001419 dependent effect Effects 0.000 claims description 8
- 238000004590 computer program Methods 0.000 claims description 6
- 238000001514 detection method Methods 0.000 abstract description 7
- 108090000623 proteins and genes Proteins 0.000 abstract description 3
- LFQSCWFLJHTTHZ-UHFFFAOYSA-N Ethanol Chemical compound CCO LFQSCWFLJHTTHZ-UHFFFAOYSA-N 0.000 description 32
- 210000002381 plasma Anatomy 0.000 description 32
- 239000008280 blood Substances 0.000 description 15
- 210000004369 blood Anatomy 0.000 description 14
- 108020004414 DNA Proteins 0.000 description 12
- 201000011510 cancer Diseases 0.000 description 12
- 239000007788 liquid Substances 0.000 description 12
- 239000000203 mixture Substances 0.000 description 7
- 238000000746 purification Methods 0.000 description 7
- 238000011282 treatment Methods 0.000 description 7
- HEMHJVSKTPXQMS-UHFFFAOYSA-M Sodium hydroxide Chemical compound [OH-].[Na+] HEMHJVSKTPXQMS-UHFFFAOYSA-M 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000010276 construction Methods 0.000 description 5
- OPTASPLRGRRNAP-UHFFFAOYSA-N cytosine Chemical compound NC=1C=CNC(=O)N=1 OPTASPLRGRRNAP-UHFFFAOYSA-N 0.000 description 5
- 239000000243 solution Substances 0.000 description 5
- 239000006228 supernatant Substances 0.000 description 5
- 239000011324 bead Substances 0.000 description 4
- 239000000872 buffer Substances 0.000 description 4
- 238000005352 clarification Methods 0.000 description 4
- 238000012546 transfer Methods 0.000 description 4
- 238000005406 washing Methods 0.000 description 4
- LRSASMSXMSNRBT-UHFFFAOYSA-N 5-methylcytosine Chemical compound CC1=CNC(=O)N=C1N LRSASMSXMSNRBT-UHFFFAOYSA-N 0.000 description 3
- 108091029523 CpG island Proteins 0.000 description 3
- KCXVZYZYPLLWCC-UHFFFAOYSA-N EDTA Chemical group OC(=O)CN(CC(O)=O)CCN(CC(O)=O)CC(O)=O KCXVZYZYPLLWCC-UHFFFAOYSA-N 0.000 description 3
- 101000653374 Homo sapiens Methylcytosine dioxygenase TET2 Proteins 0.000 description 3
- 102100030803 Methylcytosine dioxygenase TET2 Human genes 0.000 description 3
- 239000003153 chemical reaction reagent Substances 0.000 description 3
- 230000009615 deamination Effects 0.000 description 3
- 238000006481 deamination reaction Methods 0.000 description 3
- 238000002156 mixing Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 239000013598 vector Substances 0.000 description 3
- RYVNIFSIEDRLSJ-UHFFFAOYSA-N 5-(hydroxymethyl)cytosine Chemical compound NC=1NC(=O)N=CC=1CO RYVNIFSIEDRLSJ-UHFFFAOYSA-N 0.000 description 2
- FHSISDGOVSHJRW-UHFFFAOYSA-N 5-formylcytosine Chemical compound NC1=NC(=O)NC=C1C=O FHSISDGOVSHJRW-UHFFFAOYSA-N 0.000 description 2
- 108090000790 Enzymes Proteins 0.000 description 2
- 102000004190 Enzymes Human genes 0.000 description 2
- ISAKRJDGNUQOIC-UHFFFAOYSA-N Uracil Chemical compound O=C1C=CNC(=O)N1 ISAKRJDGNUQOIC-UHFFFAOYSA-N 0.000 description 2
- 229940104302 cytosine Drugs 0.000 description 2
- 230000002255 enzymatic effect Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- UYTPUPDQBNUYGX-UHFFFAOYSA-N guanine Chemical compound O=C1NC(N)=NC2=C1N=CN2 UYTPUPDQBNUYGX-UHFFFAOYSA-N 0.000 description 2
- 238000011534 incubation Methods 0.000 description 2
- 208000002154 non-small cell lung carcinoma Diseases 0.000 description 2
- 238000010606 normalization Methods 0.000 description 2
- 230000003647 oxidation Effects 0.000 description 2
- 238000007254 oxidation reaction Methods 0.000 description 2
- 238000002360 preparation method Methods 0.000 description 2
- 239000000047 product Substances 0.000 description 2
- 239000011535 reaction buffer Substances 0.000 description 2
- 239000011541 reaction mixture Substances 0.000 description 2
- 230000008439 repair process Effects 0.000 description 2
- 208000029729 tumor suppressor gene on chromosome 11 Diseases 0.000 description 2
- 238000002759 z-score normalization Methods 0.000 description 2
- HLHSUNWAPXINQU-GQCTYLIASA-N (E)-3-(3,4-dihydroxyphenyl)-N-prop-2-ynylprop-2-enamide Chemical compound OC=1C=C(C=CC=1O)/C=C/C(=O)NCC#C HLHSUNWAPXINQU-GQCTYLIASA-N 0.000 description 1
- BLQMCTXZEMGOJM-UHFFFAOYSA-N 5-carboxycytosine Chemical compound NC=1NC(=O)N=CC=1C(O)=O BLQMCTXZEMGOJM-UHFFFAOYSA-N 0.000 description 1
- 206010006187 Breast cancer Diseases 0.000 description 1
- 208000026310 Breast neoplasm Diseases 0.000 description 1
- 238000007399 DNA isolation Methods 0.000 description 1
- 108010067770 Endopeptidase K Proteins 0.000 description 1
- 208000000461 Esophageal Neoplasms Diseases 0.000 description 1
- 208000005016 Intestinal Neoplasms Diseases 0.000 description 1
- 208000035823 Non-specific autoimmune cerebellar ataxia without characteristic antibodies Diseases 0.000 description 1
- 108091028043 Nucleic acid sequence Proteins 0.000 description 1
- 206010030155 Oesophageal carcinoma Diseases 0.000 description 1
- 238000012408 PCR amplification Methods 0.000 description 1
- 229910019142 PO4 Inorganic materials 0.000 description 1
- 206010061902 Pancreatic neoplasm Diseases 0.000 description 1
- 208000005718 Stomach Neoplasms Diseases 0.000 description 1
- 238000000692 Student's t-test Methods 0.000 description 1
- 238000004833 X-ray photoelectron spectroscopy Methods 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 239000003146 anticoagulant agent Substances 0.000 description 1
- 229940127219 anticoagulant drug Drugs 0.000 description 1
- 238000013476 bayesian approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 210000000601 blood cell Anatomy 0.000 description 1
- 238000010241 blood sampling Methods 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 230000036425 denaturation Effects 0.000 description 1
- 238000004925 denaturation Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 201000004101 esophageal cancer Diseases 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 206010017758 gastric cancer Diseases 0.000 description 1
- 230000006607 hypermethylation Effects 0.000 description 1
- 201000002313 intestinal cancer Diseases 0.000 description 1
- 201000007270 liver cancer Diseases 0.000 description 1
- 208000014018 liver neoplasm Diseases 0.000 description 1
- 208000015486 malignant pancreatic neoplasm Diseases 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 201000002528 pancreatic cancer Diseases 0.000 description 1
- 208000008443 pancreatic carcinoma Diseases 0.000 description 1
- XEBWQGVWTUSTLN-UHFFFAOYSA-M phenylmercury acetate Chemical compound CC(=O)O[Hg]C1=CC=CC=C1 XEBWQGVWTUSTLN-UHFFFAOYSA-M 0.000 description 1
- NBIIXXVUZAFLBC-UHFFFAOYSA-K phosphate Chemical compound [O-]P([O-])([O-])=O NBIIXXVUZAFLBC-UHFFFAOYSA-K 0.000 description 1
- 239000010452 phosphate Substances 0.000 description 1
- 210000004180 plasmocyte Anatomy 0.000 description 1
- 229920000771 poly (alkylcyanoacrylate) Polymers 0.000 description 1
- 238000011176 pooling Methods 0.000 description 1
- 239000000843 powder Substances 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 239000012264 purified product Substances 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 239000002096 quantum dot Substances 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 201000011549 stomach cancer Diseases 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 238000012353 t test Methods 0.000 description 1
- 238000010257 thawing Methods 0.000 description 1
- 210000001519 tissue Anatomy 0.000 description 1
- 229940035893 uracil Drugs 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
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
- 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
-
- 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
- G16B35/00—ICT specially adapted for in silico combinatorial libraries of nucleic acids, proteins or peptides
- G16B35/20—Screening of libraries
-
- 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
- G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Medical Informatics (AREA)
- Biophysics (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Theoretical Computer Science (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Health & Medical Sciences (AREA)
- Biotechnology (AREA)
- Molecular Biology (AREA)
- Chemical & Material Sciences (AREA)
- Library & Information Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Software Systems (AREA)
- Public Health (AREA)
- Evolutionary Computation (AREA)
- Epidemiology (AREA)
- Databases & Information Systems (AREA)
- Data Mining & Analysis (AREA)
- Physiology (AREA)
- Bioethics (AREA)
- Analytical Chemistry (AREA)
- Genetics & Genomics (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Artificial Intelligence (AREA)
- Biochemistry (AREA)
- Investigating Or Analysing Biological Materials (AREA)
Abstract
本申请公开了一种校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质，属于基因检测技术领域。该方法通过预设数量的健康人血浆样本和肿瘤患者血浆样本的甲基化测序数据，筛选在健康人和肿瘤患者中甲基化水平相对稳定的区间作为甲基化稳定区间（house‑keeping window）；利用健康人血浆样本和待测样本的甲基化稳定区间的甲基化水平建立校正模型；利用校正模型校正待测样本全基因组内的甲基化水平。该方法在去除批次效应的同时保留了重要区域的甲基化特征，亦可在批次分组未知情况下去除批次效应。
Description
技术领域
本申请属于基因测序技术领域，具体涉及校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质。
背景技术
研究发现，癌症患者相比健康人群，甲基化水平在全基因组范围内会出现大范围的低甲基化，在部分CpG岛（CpG island）区域会出现高甲基化。当前的癌症研究，常基于全基因组甲基化测序数据，计算一定区域范围内的甲基化水平，并根据甲基化水平进行相关分析。由于建库批次和测序批次的存在，甲基化水平存在明显的批次效应，批次效应会影响最终的分析结果，因此需要对甲基化水平进行标准化处理以去除批次效应的影响。目前存在的去除批次效应的方法通常包括如下几种：
（1）Z-score标准化：样本中每个检测值减去检测值均值，再除以检测值的标准差，将检测值的均值变为0，标准差变为1；
（2）Quantile算法：根据分位数进行标准化，其所基于的统计思想是如果2个n维数据向量X、Y具有相同的数据分布，那么基于这2个向量的Quantile-Quantile图应该是平行于对角线的直线，可以将数据向量投影至对角线上来获得相同的数据分布；
（3）ComBat算法：使用经验贝叶斯方法来校正批次效应，其原理是基于估计参数的先验分布，独立估算并调整任一检测值每个批次的均值和方差。
上述去除批次效应的方法均是基于芯片数据开发，适用于对芯片数据做标准化处理。对于全基因组甲基化测序数据而言，Z-score标准化处理使得全基因组范围内甲基化水平的分布发生了改变，均值变为0，标准差变为1；Quantile算法同样也会改变甲基化水平的整体分布。
由于在肿瘤组织或患者血浆样本中，甲基化水平在全基因组范围内普遍会出现甲基化水平降低的情况，Z-score处理或Quantile算法会导致这些样本中低甲基化信号的丢失。ComBat算法对小样本数据更加有效，适用于批次分组已知情况下去除批次效应。
发明内容
1. 发明目的
本申请的目的在于提供一种校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质，该方法通过预设数量的健康人血浆样本和肿瘤患者血浆样本的甲基化测序数据，筛选在健康人和肿瘤患者中甲基化水平相对稳定的区间作为甲基化稳定区间（house-keeping window）；利用健康人血浆样本和待测样本的甲基化稳定区间的甲基化水平建立校正模型；利用校正模型校正待测样本全基因组内的甲基化水平。该方法在去除批次效应的同时保留了重要区域的甲基化特征，亦可在批次分组未知情况下去除批次效应。
2. 技术方案
为了解决上述问题，本申请所采用的技术方案如下：
本申请提供了一种校正全基因组甲基化测序数据批次效应的方法，该方法是通过预设数量的健康人血浆样本和肿瘤患者血浆样本的甲基化测序数据，筛选在健康人和肿瘤患者中甲基化水平相对稳定的区间作为甲基化稳定区间（house-keeping window）；利用健康人血浆样本和待测样本的甲基化稳定区间的甲基化水平建立校正模型；利用校正模型校正待测样本全基因组内的甲基化水平。
进一步地，上述预设数量的健康人血浆样本和肿瘤患者血浆样本尽可能覆盖健康人和各种肿瘤患者样本，健康人的年龄、性别尽可能多样，肿瘤患者与健康人之间的甲基化差异尽可能大，比如肿瘤晚期患者样本。
进一步地，上述甲基化水平包括全甲基化片段占比（methylated fragmentratio，MFR），MFR指落在某一个区间内的全甲基化片段数量/全部片段数量，全甲基化片段指该片段上所有的CpG位点全部甲基化。
进一步地，上述甲基化水平包括甲基化密度（methylation density，MD），MD指落在某一个区间内的甲基化CpG数量/全部CpG数量。
进一步地，上述甲基化稳定区间（house-keeping window）的筛选包括如下步骤：
S1，筛选唯一比对区间：从hg19参考基因组中提取每条染色体长度和碱基序列信息，根据碱基序列信息提取CpG的位置；使用GEM-library创建hg19参考基因组、CT转换后的基因组（C碱基转换成T碱基）和GA转换后的基因组（G碱基转换成A碱基）的mappability文件（对于PE100测序的数据，将kmer参数设置成100）；根据每个mappability文件，筛选不存在多重比对可能性的区间；取三组区间的交集，合并出唯一比对区间的bed文件；
S2，划分window区间并进行初步筛选：从hg19参考基因组中提取每条染色体长度和碱基序列信息，根据染色体长度信息、预设的window区间长度和步长划分window区间，同时统计window区间的唯一比对率，根据唯一比对率对划分的window区间进行初步筛选，唯一比对率指window区间上唯一比对区间占window区间总长度的比例，筛选规则如下：
，
其中，是第i个window区间内唯一比对率，/>是预设的唯一比对率阈值；
S3，甲基化水平计算：根据甲基化测序数据，筛选位于初步筛选的window区间的片段并统计片段信息，计算每个健康人血浆样本和肿瘤患者血浆样本在初步筛选的window区间的甲基化水平，甲基化水平包括全甲基化片段占比（MFR）和/或甲基化密度（MD），其计算方式如下：
区间全甲基化片段占比（MFR）计算为：
，
其中，是样本n第i个区间的MFR值，/>是样本n在第i个区间内全甲基化的片段数，/>是样本n在第i个区间片段数；
区间甲基化密度（MD）计算：
，
其中，是样本n第i个区间的MD值，/>是样本n在第i个区间内甲基化的CpG个数，/>是样本n在第i个区间的所有CpG个数；
S4，筛选甲基化稳定区间：将健康人血浆样本和肿瘤患者血浆样本在初步筛选的window区间的甲基化水平进行差异化分析，筛选满足预设条件的window区间为甲基化稳定区间，具体为：
当甲基化水平用MFR表示时，使用limma包对初步筛选的window区间的MFR进行差异分析（t-test），获得p-value后再计算FDR（使用R校正函数p.adjust计算FDR）；分别统计肿瘤患者和健康人的初步筛选的window区间的MFR的最大值、最小值和均值，按照如下条件筛选甲基化稳定区间（house-keeping window）：
，
，
，
，
其中，是第i个window区间的FDR值，/>是预设的FDR的阈值；是肿瘤患者样本在第i个window区间的最大MFR值，/>是健康人在第i个window区间的最大MFR值，/>是肿瘤患者样本在第i个window区间的最小MFR值，/>是健康人在第i个window区间的最小MFR值，/>是肿瘤患者样本在第i个window区间的MFR均值，/>是健康人在第i个window区间的MFR均值，/>是预设的house-keeping window在肿瘤患者和健康人之间甲基化水平差异的阈值；
当甲基化水平用MD表示时，用limma包对初步筛选的window区间的MD进行差异分析，获得p-value后再计算FDR；分别统计肿瘤患者和健康人的初步筛选的window区间的MD的最大值、最小值和均值，按照如下条件筛选甲基化稳定区间：
，
，
，
，
其中，是第i个window区间的FDR值，/>是预设的FDR的阈值；是肿瘤患者样本在第i个window区间的最大MD值，/>是健康人在第i个window区间的最大MD值，/>是肿瘤患者样本在第i个window区间的最小MD值，是健康人在第i个window区间的最小MD值，/>是肿瘤患者样本在第i个window区间的MD均值，/>是健康人在第i个window区间的MD均值，/>是预设的house-keeping window在肿瘤患者和健康人之间甲基化水平差异的阈值。
进一步地，上述S2中预设的window区间长度为5000 bp；步长（step）为5000 bp。
进一步地，上述S2中预设的唯一比对率阈值为0.95。
进一步地，上述S4中预设的FDR的阈值为0.05。
进一步地，上述S4中是预设的house-keeping window在肿瘤患者和健康人之间甲基化水平差异的阈值为0.2。
进一步地，上述S3中计算前还包括：筛选位于初步筛选的window区间的片段并统计片段信息，片段信息包括片段长度、CpG个数和依据none CpG位点计算出的转化率（转化率=甲基化的none CpG/所有none CpG），并根据片段长度、CpG个数和转化率这三个条件筛选片段，筛选条件如下：
，
，
，
其中，是样本在第i个区间上第j条片段的长度，是预设的片段长度的阈值下限，/>是预设的片段长度的阈值上限；/>是样本在第i个区间上第j条片段含有的CpG个数，/>是预设的片段上CpG个数的阈值；/>是样本在第i个区间上第j条片段的转化率，/>是预设片段转化率的阈值。
进一步地，上述预设的片段长度的阈值下限为80 bp。
进一步地，上述预设的片段长度的阈值上限为250 bp。
进一步地，上述预设的片段上CpG个数的阈值为3。
进一步地，上述预设的片段转化率的阈值为0.95。
进一步地，上述S4中，将健康人血浆样本和肿瘤患者血浆样本在初步筛选的window区间的甲基化水平进行差异化分析包括将健康人血浆样本和不同肿瘤患者样本、所有肿瘤患者样本分别在初步筛选的window区间的甲基化水平进行差异化分析。
进一步地，上述方法筛选的甲基化稳定区间，若两个甲基化稳定区间存在至少一个碱基的位置重叠，则直接合并这两个甲基化稳定区间，将其作为一个window区间重新计算筛选、甲基化水平计算步骤。
进一步地，上述利用健康人血浆样本和待测样本的甲基化稳定区间的甲基化水平建立校正模型，包括如下步骤：
以健康人血浆样本作为基线样本，提取基线样本的每个甲基化稳定区间的甲基化水平值，并计算平均值；
根据待测样本的测序数据，计算待测样本的每个甲基化稳定区间的甲基化水平值；
以基线样本的每个甲基化稳定区间的甲基化水平平均值作为因变量，以待测样本的每个甲基化稳定区间的甲基化水平值作为自变量，构建如下线性回归模型：
，
其中，是基线样本第i个甲基化稳定区间的甲基化水平平均值，/>是该待测样本第i个甲基化稳定区间的甲基化水平值，b是回归系数，a是常数项。
进一步地，上述构建线性回归模型包括使用R包stats的lm函数构建Linearregression (LM)模型。
进一步地，上述构建线性回归模型包括使用R包robustbase的lmrob函数构建Robust Linear regression(RLM)模型；该方法构建Robust Linear regression(RLM)模型时，使用MM-Estimate方法来估计权重，对不同的数据点赋予不同的权重，削减极端值的权重。
进一步地，当甲基化水平用MFR表示时，上述建立校正模型则包括如下步骤：
以健康人血浆样本作为基线样本，提取基线样本的每个甲基化稳定区间的MFR值，并计算平均值；
根据待测样本的测序数据，计算待测样本的每个甲基化稳定区间的MFR值；
以基线样本的每个甲基化稳定区间的MFR平均值作为因变量，以待测样本的每个甲基化稳定区间的MFR值作为自变量，构建如下线性回归模型：
，
其中，是基线样本第i个甲基化稳定区间的MFR平均值，/>是该待测样本第i个甲基化稳定区间的MFR值，b是回归系数，a是常数项。
进一步地，当甲基化水平用MFR表示时，为排除极端值的影响，针对每个基线样本甲基化稳定区间的MFR值，剔除MFR值最高与最低的5%的样本值，然后再求平均值。
进一步地，上述构建线性回归模型前还包括甲基化稳定区间筛选，具体为采用backward elimination（向后剔除法）按照预设比例剔除存在前后相关性差异大的house-keeping window，使用保留下来的house-keeping window MFR构建校正模型，以减小校正模型受到极端值的影响，所述backward elimination（向后剔除法）包括：
B1，用所有I个house-keeping window计算待测样本MFR和基线样本MFR均值的相关性：
，
其中，n表示第n个需要校正的样本，是样本n在所有house-keeping window的MFR值，/>是基线样本在所有house-keeping window的MFR均值；
B2，从I个house-keeping window剔除一个window，重新计算I-1个window上，待测样本MFR和基线样本MFR均值的相关性，每次剔除的window不同，可计算获得共I个相关性值：
，
是样本n和基线样本均剔除第i个house-keeping window后，剩余house-keeping window MFR的相关性；
B3，计算剔除前和剔除后相关性差异最大的house-keeping window。根据下列公式计算剔除每一个house-keeping window时的相关性差值，选取相关性差值最大的house-keeping window，作为本轮剔除的window，
，
其中，是样本n第i个house-keeping window剔除后相关性的差值；
B4，在剔除第i个house-keeping window的基础上，重复B2和B3，找到下一个剔除前后相关性差异最大的house-keeping window，进行剔除，重复步骤直到删除的window区间个数满足给定的比例后结束。
进一步地，上述预设比例可设定为0-100%之间的值。更进一步地，上述预设比例可设定为5%、10%或20%。
进一步地，上述采用backward elimination按照预设比例剔除存在前后相关性差异大的甲基化稳定区间前还包括甲基化稳定区间筛选，因为不同甲基化稳定区间的coverage（测序片段覆盖度）差异较大，coverage偏低的区间计算出的MFR值偏离真实值的可能性越大，参与模型构建会引入较大误差，所以在建模前，会依据待测样本甲基化稳定区间的片段数进行筛选，筛选规则如下：
其中，是待测样本n第i个甲基化稳定区间的片段数，/>是预设的甲基化稳定区间片段数的阈值；
当甲基化水平用MD表示时，根据待测样本甲基化稳定区间上CpG的总覆盖深度（测到该CpG的reads数）进行筛选，筛选公式为：
，
其中，是待测样本n第i个甲基化稳定区间的CpG总深度，/>是预设的甲基化稳定区间CpG总深度的阈值。
进一步地，上述预设的甲基化稳定区间片段数的阈值为10。
进一步地，上述利用校正模型校正全基因组内的甲基化水平包括：
M1，划分Bin区间并进行初步筛选：从hg19参考基因组中提取每条常染色体和碱基序列信息，根据染色体长度信息、预设的Bin区间长度L和步长S划分window区间，同时统计window区间的唯一比对率，根据唯一比对率对划分的Bin区间进行初步筛选，唯一比对率指Bin区间上唯一比对区间占Bin区间总长度的比例，筛选规则如下：
，
其中，是第i个Bin区间内唯一比对率，/>是预设的唯一比对率阈值；
M2，根据划分Bin区间，根据上述甲基化水平计算方法计算待测样本Bin区间的甲基化水平，并作为自变量，计算得到校正后的甲基化水平。
进一步地，上述预设的Bin区间长度L=1 Mb。
进一步地，上述预设的Bin区间步长S=0 bp。
进一步地，上述预设的唯一比对率阈值为0.95。
本申请还提供了一种校正全基因组甲基化测序数据批次效应的装置，包括：
输入模块，其被配置为输入健康人血浆样本、肿瘤患者血浆样本和待测样本的甲基化测序数据；
唯一比对区间筛选模块，其被配置为筛选唯一比对区间；
window区间划分模块，其被配置为根据染色体长度信息、预设的window区间长度和步长划分window区间，同时根据唯一比对率对划分的window区间进行初步筛选；
甲基化水平计算模块，其被配置为用于计算区间甲基化水平，如window区间、甲基化稳定区间、Bin区间；
甲基化稳定区间筛选模块，其被配置为筛选满足预设条件的window区间为甲基化稳定区间；
校正模型建立模块，其被配置为以基线样本的每个甲基化稳定区间的甲基化水平平均值作为因变量，以待测样本的每个甲基化稳定区间的甲基化水平值作为自变量，构建线性回归模型；
全基因组甲基化测序数据校正模块，其被配置为利用校正模型建立模块建立的校正模型校正待测样本的甲基化水平。
本申请还提供了一种电子设备，包括：一个或多个处理器；存储装置，其上存储有一个或多个程序，当一个或多个程序被一个或多个处理器执行，使得一个或多个处理器实现上述校正全基因组甲基化测序数据批次效应的方法。
本申请还提供了一种计算机存储介质，其上存储有计算机程序，其中，计算机程序被处理器执行时实现上述校正全基因组甲基化测序数据批次效应的方法。
3. 有益效果
本申请与现有技术相比，其有益效果在于：
（1）本申请提供的校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质，该方法通过建立基线样本甲基化稳定的区间（house-keeping window）的甲基化水平和单个样本甲基化稳定的区间（house-keeping window）的线性模型，校正全基因组范围内的甲基化水平，甲基化稳定的区间（house-keeping window）指在健康人血浆和患者血浆中甲基化水平相对稳定的区间。该方法在去除批次效应的同时保留了重要区域的甲基化特征，亦可在批次分组未知情况下去除批次效应。
（2）本申请提供的校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质，针对每个需要校正样本的house-keeping window甲基化水平建立线性校正模型来校正全基因组范围内的区域甲基化水平，不会大幅改变校正样本的甲基化水平分布，不会消除重要区域的甲基化特征。
（3）本申请提供的校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质，其校正模型只考虑筛选出的house-keeping window的甲基化水平，对于未知批次效应的样本亦可采用本方法进行校正，缩小因批次效应带来的系统误差，保证后续分析结果的准确性。
附图说明
图1是本申请筛选出2673个house-keeping window在参考基因组的分布。
图2是运输时间测试样本的原始bin-level MFR和校正后bin-level MFR的t-SNE分析结果，1-20表示受试者，A表示24 h处理，B表示72 h处理。
图3是采血管测试样本的原始bin-level MFR和校正后bin-level MFR的t-SNE分析结果，1-20表示受试者，A表示康为采血管，B表示EDTA采血管。
具体实施方式
下面结合具体实施例对本申请进一步进行描述。
除非另有定义，本文所使用的所有的技术和科学术语与属于本申请的技术领域的技术人员通常理解的含义相同；本文所使用的术语“和/或”包括一个或多个相关的所列项目的任意的和所有的组合。
实施例中未注明具体条件者，按照常规条件或制造商提供的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市售购买获得的常规产品。
如本文所使用，术语“待测样本”、“待校正样本”均是指需要进行矫正的样本，可以相互通用。
如本文所使用，术语“CpG”指由胞嘧啶(C)和鸟嘌呤(G)组成的一个二核苷酸链，当中的p是C和G之间的磷酸。
如本文所使用，术语“CpG island”指长度至少500 bp、GC含量超过55%、CpG比值大于0.65的 DNA 序列区域，主要位于基因的启动子（promoter）和第一外显子区域。
如本文所使用，术语“癌症患者”、“肿瘤患者”均非健康人群，可以相互通用。
浓度、量和其他数值数据可以在本文中以范围格式呈现。应当理解，这样的范围格式仅是为了方便和简洁而使用，并且应当灵活地解释为不仅包括明确叙述为范围极限的数值，而且还包括涵盖在所述范围内的所有单独的数值或子范围，就如同每个数值和子范围都被明确叙述一样。例如，约1至约4.5的数值范围应当被解释为不仅包括明确叙述的1至约4.5的极限值，而且还包括单独的数字（诸如2、3、4）和子范围（诸如1至3、2至4等）。相同的原理适用于仅叙述一个数值的范围，诸如“小于约4.5”，应当将其解释为包括所有上述的值和范围。此外，无论所描述的范围或特征的广度如何，都应当适用这种解释。
实施例1
本实施例提供校正全基因组甲基化测序数据批次效应的方法及其验证，具体包括如下步骤：
（1）用于分析的样本
筛选house-keeping window的样本涵盖了健康人（HEALTHY）和癌症患者，其中共有7类癌症患者，分别是ESCA（食管癌）、BRCA（乳腺癌）、PACA（胰腺癌）、LIHC（肝癌）、STAD（胃癌）、NSCLC（非小细胞肺癌）和COREAD（肠癌），健康人和癌症患者的样本数量如表1所示：
表1
使用上述352例健康人（HEALTHY）样本作为基线样本，参与校正模型的建立。
本实施例中，模拟了一批相同受试者不同运输时间和不同采血管处理的样本，这两批样本存在明显的批次效应，所以使用这些样本进行本申请校正全基因组甲基化测序数据批次效应的方法的验证。其中不同运输时间处理的有20位受试者，40例样本；不同采血管处理的有20位受试者，40例样本。具体为：
采血管测试：同一个志愿者分别使用康为游离DNA样本保存管（10 mL）和EDTA抗凝管（10 mL）采血，一共20位志愿者献血，得到20对样本测试采血管影响。
不同运输时间测试：同一个志愿者使用康为游离DNA样本保存管（10 mL）采2管血，分别运输24 h和48 h，一共20位志愿者献血，得到20对样本测试运输时间影响。
（2）实验过程
血浆cfDNA提取：每位受试者10 mL全血在4℃以1600 g转速离心10 min使血浆、血细胞分层。将上层血浆转移至新离心管，再次以12000 rpm转速4℃离心15 min取上清以去除细胞碎屑，得到约4 mL血浆，-80℃冻存备用。血浆样本融化后，每1 mL样本中加入15 μLProteinase K(20 mg/mL，thermoscientific cat#EO0492)和50 μL SDS(20%)（血浆量不足4 mL时，用PBS补足）。翻转混匀，60℃孵育20 min，然后冰浴5 min。使用MagMAX Cell-FreeDNA Isolation试剂盒（thermoscientific cat#A29319）提取cfDNA。使用Bioanalyzer2100 (Agilent Technologies)检测cfDNA的提取浓度和质量。
cfDNA文库构建：使用甲基化文库构建试剂盒NEBNext Enzymatic Methyl-seqKit（NEB，cat# E7120），以5-30 ng cfDNA起始量，通过TET2酶使5-甲基胞嘧啶（5-mC）转化为5-甲酰胞嘧啶（5-fC）和5-羧基胞嘧啶（5-caC），并且通过APOBEC酶，使非甲基化胞嘧啶（C）脱氨转化为尿嘧啶（U），然后进行扩增建库。具体文库构建过程如下：
内参准备：取50 μL CpG全甲基化的pUC19 DNA和50 μL CpG全非甲基化的LamdbaDNA混匀后加入100 μL打断管中，使用M220打断仪（Covaris）打断。建库时，向待测cfDNA加入0.001 ng的pUC19 DNA和0.02 ng 的Lambda DNA。
cfDNA样本的准备：cfDNA样本起始量为5-30 ng，不需要打断。
末端修复：在冰上混合表2的反应体系；反应体系置于PCR仪上，按表3进行末端修复反应。
表2
表3
连接接头（Adaptor）：在冰上操作，将表4组分加入上步的60 μL反应体系中，20℃孵育15 min。
表4
连接后纯化：上一步反应结束后，取出样本，加入110 μLNEBNext SamplePurification Beads，立即使用移液器吹打混匀。室温孵育5 min。离心管置于磁力架上5min待液体澄清，弃去上清。加入200 μL现配80％乙醇，孵育30 s后弃去。重复一次200 μL80％乙醇清洗步骤。用10 μL移液器吸尽离心管底部的残留乙醇，室温干燥3-5 min至乙醇完全挥发。从磁力架取下离心管，加入29 μL Elution Buffer （NEB），震荡混匀。室温孵育1min。短暂离心，离心管置于磁力架上3 min待液体澄清，取28 μL放进新的PCR管中。
5-甲基胞嘧啶和5-羟甲基胞嘧啶氧化反应：使用NEBNext Enzymatic Methyl-seqKit（NEB，cat# E7120）进行以下反应操作，TET2 Reaction Buffer Supplement干粉加入400 μL TET2 Reaction Buffer，充分混合。在冰上将表5组分加入上述28 μL已连接adapter的DNA。将500 mM Fe(II) 溶液按1：1250比例稀释。按表6往上步混匀的产物中加入已配好的Fe(II)。充分混合并在37℃孵育1 h。反应结束后移至冰上并按表7加入1 μL StopReagent。充分混合。按表8所示37℃孵育30 min。
表5
表6
表7
表8
氧化后纯化：上一步反应结束后，取出样本，加入90 μLNEBNext SamplePurification Beads，立即使用移液器吹打混匀。室温孵育5 min。离心管置于磁力架上5min待液体澄清，弃去上清。加入200 μL现配80％乙醇，孵育30 s后弃去。重复一次200 μL80％乙醇清洗步骤。用10 μL移液器吸尽离心管底部的残留乙醇，室温干燥3-5 min至乙醇完全挥发。从磁力架取下离心管，加入17 μL Elution Buffer，震荡混匀。室温孵育1 min。短暂离心，离心管置于磁力架上3 min待液体澄清，取16 μL放进新的PCR管中。
DNA变性：配制新鲜的0.1N NaOH。提前预热PCR仪到50℃。加入4 μL 0.1N NaOH到上步16 μL纯化产物中，充分混合。50℃孵育10 min。反应结束后立刻放入冰上。
胞嘧啶脱氨基：在冰上将表9组分加入上步20 μL变性DNA。充分混合。在PCR仪上37℃孵育3 h后转为4℃终止反应。
表9
/>
脱氨后纯化：上一步反应结束后，取出样本，加入100 μLNEBNext SamplePurification Beads，立即使用移液器吹打混匀。室温孵育5 min。离心管置于磁力架上5min待液体澄清，弃去上清。加入200 μL 现配80％乙醇，孵育30 s后弃去。重复一次200 μL80％乙醇清洗步骤。用10 μL移液器吸尽离心管底部的残留乙醇，室温干燥3-5 min至乙醇完全挥发。从磁力架取下离心管，加入21 μL Elution Buffer，震荡混匀。室温孵育1 min。短暂离心，离心管置于磁力架上3 min待液体澄清，取20 μL放进新的PCR管中。
文库PCR扩增：在冰上将表10组分加入上步脱氨后的20 μL DNA。充分混合后在PCR以上进行表11的程序进行PCR反应。
表10
表11
PCR后纯化：上一步反应结束后，取出样本，加入45 μLNEBNext SamplePurification Beads，立即使用移液器吹打混匀。室温孵育5 min。离心管置于磁力架上5min待液体澄清，弃去上清。加入200 μL 现配80％乙醇，孵育30 s后弃去。重复一次200 μL80％乙醇清洗步骤。用10 μL移液器吸尽离心管底部的残留乙醇，室温干燥3-5 min至乙醇完全挥发。从磁力架取下离心管，加入21 μL Elution Buffer，震荡混匀。室温孵育1 min。短暂离心，离心管置于磁力架上3 min待液体澄清，取20 μL放进新的PCR管中。
文库定量：使用Qubit高灵敏试剂（thermoscientific cat#Q32854）对所构建的文库进行定量，文库产量大于400 ng进行后续上机测序。
文库测序：
取100 ng上述文库加入10% PhiX DNA（Illumina cat#FC-110-3001）混合成上机样品，在Novaseq 6000（Illumina）平台进行PE100测序。
（3）测序数据预处理
处理下机FASTQ数据为各模块可使用的Bam文件，包括：
去接头：调用Trimmomatic-0.36将每一对FASTQ文件都作为配对的读段（pairedreads）比对到hg19人类参考基因组序列，除M参数与指定Reads Group的ID外，不使用其余参数选项，生成初始bam文件。
比对：调用Bismark-v0.19.0将去接头后的每一对FASTQ文件都作为配对读段比对到hg19人类参考基因组序列和Lambda DNA参考基因组序列，生成初始Bam文件。
去重：调用Bismark- v0.19.0的deduplicate模块，对初始Bam文件进行去重复处理，生成去重后的Bam文件。
排序标记：调用SAMtools-1.3的sort模块，对去重后的Bam文件进行排序，生成排序后的Bam文件。然后，调用Picard-2.1.0的AddOrReplaceReadGroups模块，对排序后的Bam文件进行标记分组。
筛选：调用BamUtil-1.0.14的clipOverlap模块对标记分组后的Bam文件进行筛选，去除重叠的配对读段，生成Bam文件。并调用SAMtools-1.3 view对去除重叠的Bam文件的比对质量进行过滤，采用“-q 20”作为参数，生成最终Bam文件。
建立索引：调用SAMtools-1.3的index模块对最终生成的Bam文件建立索引，生成与最终Bam文件配对的bai文件。
（4）全基因组甲基化测序数据批次效应校正
（a）筛选唯一比对区间，从ucsc网站（https://genome.ucsc.edu/index.html）下载hg19参考基因组，根据参考基因组碱基序列信息，获取CpG位点位置信息。使用GEM-library创建hg19参考基因组、CT转换后的基因组（C碱基转换成T碱基）和GA转换后的基因组（G碱基转换成A碱基）的mappability文件，针对PE100测序数据，kmer参数设置成100。使用gem-2-wig、wigToBigWig、bigWigToBedGraph和bedGraphTobed工具将mappability文件转换成bed文件，挑选不存在多重比对可能性的区间，使用bedtools coverage对3份bed文件取交集，合并出唯一比对区间的bed文件。
（b）划分window区间并进行初步筛选，从hg19参考基因组中提取每条染色体长度和碱基序列信息，设定参数window=5000 bp，step=5000 bp，使用bedtools makewindows对常染色体chr1, chr2, ..., chr21和chr22划分window区间，获得window区间个数如表12所示：
表12
使用bedtools coverage统计window区间的mappability，按照5000 bp大小的window mappability大于等于0.95这个条件筛选window区间，筛选条件如下：
，
其中，是第i个window区间内唯一比对率，唯一比对率为window区间上唯一比对区间占window总长度的比例，各染色体上window区间筛选结果如表13所示：
表13
（c）甲基化水平计算，根据甲基化测序数据，筛选位于初步筛选的window区间的片段并统计片段信息，计算每个健康人血浆样本和肿瘤患者血浆样本在初步筛选的window区间的甲基化水平，甲基化水平包括全甲基化片段占比（MFR）和/或甲基化密度（MD），在本实施例中，采用全甲基化片段占比（MFR）来表示甲基化水平。
首先，针对548例筛选house-keeping window样本的bam文件，统计每条片段（fragment）的长度、CpG个数、依据none CpG位点计算出的转化率，对每一个window区间上的片段，按照fragment长度大于等于80 bp，小于等于250 bp，每个fragment上至少3个CpG位点，转化率大于等于0.95来筛选，具体过滤条件如下：
，
，
，
其中，是样本在第i个window区间上第j条片段的长度，/>是样本在第i个window区间上第j条片段含有的CpG个数，/>是样本在第i个window区间上第j条片段的转化率。
然后，根据MFR公式计算每一个样本每一个window区间的MFR值，计算公式如下：
，
其中，是样本n第i个window区间的MFR值，/>是样本n在第i个window区间内全甲基化的片段数，/>是样本n在第i个window区间的片段数。
（d）筛选甲基化稳定区间（house-keeping window），将上述548例样本按照不同癌症 vs 健康人和所有癌症样本 vs 健康人，一共8组进行比较，对所有window区间MFR使用limma包进行差异分析，获得统计p-value，然后转化成FDR，再统计各组所有window区间MFR的最小值、最大值和均值。每组比较均按照如下条件来筛选：
，
，
，/>
，
其中，是第i个window区间的FDR值，/>是癌症样本在第i个window区间的最大MFR值，/>是健康人在第i个window区间的最大MFR值，是癌症样本在第i个window区间的最小MFR值，/>是健康人在第i个window区间的最小MFR值，/>是癌症样本在第i个window区间的MFR均值，是健康人在第i个window区间的MFR均值。
所有window区间在所有分组比较中均满足上述条件的被认定是甲基化稳定区间（house-keeping window），最后一共筛选出2673个house-keeping window，在参考基因组的分布如图1所示。
（e）RLM线性模型建立
提取352个基线样本在2673个house-keeping window的单样本MFR。对每个house-keeping window，剔除基线样本中MFR值最高与最低的5%的样本值，计算window MFR在基线样本内的均值。
计算80例（20对运输时间测试+20对采血管测试）待校正样本2673个house-keeping window的单样本MFR，同时统计各区间的单样本fragments数量。按照下列条件进行house-keeping window过滤：
，
其中，是待校正样本n第i个house-keeping window的片段数。
使用过滤后的house-keeping window进行建模，基线样本的house-keepingwindow MFR均值作为因变量，待校正样本的house-keeping window MFR值作为自变量，使用MM-Estimate方法来估计权重，建立Robust Linear regression(RLM)模型，得到回归方程。
，
其中，是基线样本第i个house-keeping window的MFR均值，/>是该样本第i个house-keeping window MFR，b是回归系数，a是常数项。
（f）全基因组甲基化测序数据批次效应校正（bin-level MFR校正）
划分Bin区间并进行初步筛选：bin区间的划分方法类比window区间，其中设定参数bin区间长度L=1 Mb，bin区间步长S=0 bp，使用bedtools makewindows对常染色体chr1,chr2, ..., chr21和chr22划分bin，按照bin区间内的唯一比对率大于等于0.95来筛选，最终获得1846个bin。筛选公式如下：
，
其中，是第i个bin的唯一比对率。
bin-level MFR校正，计算80例（20对运输时间测试+20对采血管测试）待校正样本1846个bin区间的单样本MFR，将1846个bin-level MFR作为自变量，代入确定了回归系数和常数项的回归方程，得到因变量的值，作为校正后的MFR值。
结果分析：
运输时间测试样本如图2所示，原始bin-level MFR的t-SNE分析结果如左图，校正后bin-level MFR的t-SNE分析结果如右图。图中数字1-20表示受试者，A表示24h处理，B表示72h处理。
采血管测试样本的结果如图3所示，原始bin-level MFR的t-SNE分析结果如左图，校正后bin-level MFR的t-SNE分析结果如右图。图中数字1-20表示受试者，A表示康为采血管，B表示EDTA采血管。
校正后，大部分的相同受试者的样本聚集在一起了，说明该校正方法可有效去除批次效应。
实施例2
本实施例提供一种校正全基因组甲基化测序数据批次效应的装置，包括：
输入模块，其被配置为输入健康人血浆样本、肿瘤患者血浆样本和待测样本的甲基化测序数据；
唯一比对区间筛选模块，其被配置为筛选唯一比对区间；
window区间划分模块，其被配置为根据染色体长度信息、预设的window区间长度和步长划分window区间，同时根据唯一比对率对划分的window区间进行初步筛选；
甲基化水平计算模块，其被配置为用于计算区间甲基化水平，如window区间、甲基化稳定区间、Bin区间；
甲基化稳定区间筛选模块，其被配置为筛选满足预设条件的window区间为甲基化稳定区间；
校正模型建立模块，其被配置为以基线样本的每个甲基化稳定区间的甲基化水平平均值作为因变量，以待测样本的每个甲基化稳定区间的甲基化水平平均值作为自变量，构建线性回归模型；
全基因组甲基化测序数据校正模块，其被配置为利用校正模型建立模块建立的校正模型校正待测样本的甲基化水平。
实施例3
本实施例提供一种电子设备，包括：一个或多个处理器；存储装置，其上存储有一个或多个程序，当一个或多个程序被一个或多个处理器执行，使得一个或多个处理器实现实施例1所述的校正全基因组甲基化测序数据批次效应的方法。
实施例4
本实施例提供一种计算机存储介质，其上存储有计算机程序，其中，计算机程序被处理器执行时实现实施例1所述的校正全基因组甲基化测序数据批次效应的方法。
Claims (10)
1.一种校正全基因组甲基化测序数据批次效应的方法，其特征在于，所述方法是通过预设数量的健康人血浆样本和肿瘤患者血浆样本的甲基化测序数据，筛选在健康人和肿瘤患者中甲基化水平相对稳定的区间作为甲基化稳定区间；利用健康人血浆样本和待测样本的甲基化稳定区间的甲基化水平建立校正模型；利用校正模型校正待测样本全基因组内的甲基化水平；
所述甲基化水平包括全甲基化片段占比MFR，MFR指落在某一个区间内的全甲基化片段数量/全部片段数量，全甲基化片段指该片段上所有的CpG位点全部甲基化；或甲基化水平包括甲基化密度MD，MD指落在某一个区间内的甲基化CpG数量/全部CpG数量；
所述甲基化稳定区间的筛选包括如下步骤：
S1，筛选唯一比对区间：从hg19参考基因组中提取每条染色体长度和碱基序列信息，根据碱基序列信息提取CpG的位置；使用GEM-library创建hg19参考基因组、CT转换后的基因组和GA转换后的基因组的mappability文件；根据每个mappability文件，筛选不存在多重比对可能性的区间；取三组区间的交集，合并出唯一比对区间的bed文件；
S2，划分window区间并进行初步筛选：从hg19参考基因组中提取每条染色体长度和碱基序列信息，根据染色体长度信息、预设的window区间长度和步长划分window区间，同时统计window区间的唯一比对率，根据唯一比对率对划分的window区间进行初步筛选，唯一比对率指window区间上唯一比对区间占window区间总长度的比例，筛选规则如下：
，
其中，是第i个window区间内唯一比对率，/>是预设的唯一比对率阈值；
S3，甲基化水平计算：根据甲基化测序数据，筛选位于初步筛选的window区间的片段并统计片段信息，计算每个健康人血浆样本和肿瘤患者血浆样本在初步筛选的window区间的甲基化水平，计算方式如下：
区间全甲基化片段占比MFR计算为：
，
其中，是样本n第i个区间的MFR值，/>是样本n在第i个区间内全甲基化的片段数，/>是样本n在第i个区间片段数；
区间甲基化密度MD计算为：
其中，是样本n第i个区间的MD值，/>是样本n在第i个区间内甲基化的CpG个数，/>是样本n在第i个区间的所有CpG个数；
S4，筛选甲基化稳定区间：将健康人血浆样本和肿瘤患者血浆样本在初步筛选的window区间的甲基化水平进行差异化分析，筛选满足预设条件的window区间为甲基化稳定区间，具体为：
当甲基化水平用MFR表示时，使用limma包对初步筛选的window区间的MFR进行差异分析，获得p-value后再计算FDR；分别统计肿瘤患者和健康人的初步筛选的window区间的MFR的最大值、最小值和均值，按照如下条件筛选甲基化稳定区间：
，
，
，
，
其中，是第i个window区间的FDR值，/>是预设的FDR的阈值；是肿瘤患者样本在第i个window区间的最大MFR值，/>是健康人在第i个window区间的最大MFR值，/>是肿瘤患者样本在第i个window区间的最小MFR值，/>是健康人在第i个window区间的最小MFR值，/>是肿瘤患者样本在第i个window区间的MFR均值，/>是健康人在第i个window区间的MFR均值，/>是预设的甲基化稳定区间在肿瘤患者和健康人之间甲基化水平差异的阈值；
当甲基化水平用MD表示时，使用limma包对初步筛选的window区间的MD进行差异分析，获得p-value后再计算FDR；分别统计肿瘤患者和健康人的初步筛选的window区间的MD的最大值、最小值和均值，按照如下条件筛选甲基化稳定区间：
，
，
，
，
其中，是第i个window区间的FDR值，/>是预设的FDR的阈值；是肿瘤患者样本在第i个window区间的最大MD值，/>是健康人在第i个window区间的最大MD值，/>是肿瘤患者样本在第i个window区间的最小MD值，是健康人在第i个window区间的最小MD值，/>是肿瘤患者样本在第i个window区间的MD均值，/>是健康人在第i个window区间的MD均值，/>是预设的甲基化稳定区间在肿瘤患者和健康人之间甲基化水平差异的阈值。
2.根据权利要求1所述的一种校正全基因组甲基化测序数据批次效应的方法，其特征在于，所述利用健康人血浆样本和待测样本的甲基化稳定区间的甲基化水平建立校正模型，包括如下步骤：
以健康人血浆样本作为基线样本，提取基线样本的每个甲基化稳定区间的甲基化水平值，并计算平均值；
根据待测样本的测序数据，计算待测样本的每个甲基化稳定区间的甲基化水平值；
以基线样本的每个甲基化稳定区间的甲基化水平平均值作为因变量，以待测样本的每个甲基化稳定区间的甲基化水平值作为自变量，构建线性回归模型：
，
其中，是基线样本第i个甲基化稳定区间的甲基化水平平均值，/>是该待测样本第i个甲基化稳定区间的甲基化水平值，b是回归系数，a是常数项。
3.根据权利要求2所述的一种校正全基因组甲基化测序数据批次效应的方法，其特征在于，所述利用校正模型校正待测样本全基因组内的甲基化水平包括：
M1，划分Bin区间并进行初步筛选：从hg19参考基因组中提取每条常染色体和碱基序列信息，根据染色体长度信息、预设的Bin区间长度L和步长S划分Bin区间，同时统计Bin区间的唯一比对率，根据唯一比对率对划分的Bin区间进行初步筛选，唯一比对率指Bin区间上唯一比对区间占Bin区间总长度的比例，筛选规则如下：
，
其中，是第i个Bin区间内唯一比对率，/>是预设的唯一比对率阈值；
M2，根据划分Bin区间，根据权利要求1中S3甲基化水平计算中所述计算方式计算待测样本Bin区间的甲基化水平，并作为自变量，计算得到校正后的甲基化水平。
4. 根据权利要求3所述的一种校正全基因组甲基化测序数据批次效应的方法，其特征在于，所述构建线性回归模型包括使用R包stats的lm函数构建Linear regression模型；或使用R包robustbase的lmrob函数构建Robust Linear regression模型。
5. 根据权利要求3或4所述的一种校正全基因组甲基化测序数据批次效应的方法，其特征在于，所述构建线性回归模型前还包括采用backward elimination按照预设比例剔除存在前后相关性差异大的甲基化稳定区间，使用保留下来的甲基化稳定区间构建线性回归模型，所述backward elimination包括：
B1，用所有I个甲基化稳定区间计算待测样本MFR和基线样本MFR均值的相关性：
，
其中，n表示第n个待测样本，是样本n在所有甲基化稳定区间的MFR值，是基线样本在所有甲基化稳定区间的MFR均值；
B2，从所有I个甲基化稳定区间剔除一个window区间，重新计算I-1个window区间上，待测样本MFR和基线样本MFR均值的相关性，每次剔除的window区间不同，可计算获得共I个相关性值：
，
其中，是样本n和基线样本均剔除第i个甲基化稳定区间后，剩余甲基化稳定区间MFR的相关性；
B3，计算剔除前和剔除后相关性差异最大的甲基化稳定区间，根据下列公式计算剔除每一个甲基化稳定区间时的相关性差值，选取相关性差值最大的甲基化稳定区间，作为本轮剔除的window区间，
，
其中，是待测样本n第i个甲基化稳定区间剔除后相关性的差值；
B4，在剔除第i个甲基化稳定区间的基础上，重复B2和B3，找到下一个剔除前后相关性差异最大的甲基化稳定区间进行剔除；重复步骤直到删除的window区间个数满足预设的比例。
6. 根据权利要求5所述的一种校正全基因组甲基化测序数据批次效应的方法，其特征在于，采用backward elimination按照预设比例剔除存在前后相关性差异大的甲基化稳定区间前还包括甲基化稳定区间筛选步骤：
当甲基化水平用MFR表示时，根据甲基化稳定区间的片段数进行筛选甲基化稳定区间，筛选规则如下：
，
其中，是待测样本n第i个甲基化稳定区间的片段数，/>是预设的甲基化稳定区间片段数的阈值；
当甲基化水平用MD表示时，根据甲基化稳定区间上CpG的总覆盖深度进行筛选甲基化稳定区间，筛选公式为：
，
其中，是待测样本n第i个甲基化稳定区间的CpG总深度，/>是预设的甲基化稳定区间CpG总深度的阈值。
7. 根据权利要求6所述的一种校正全基因组甲基化测序数据批次效应的方法，其特征在于，所述S3中甲基化水平计算前还包括片段筛选步骤，具体包括：筛选位于初步筛选的window区间的片段并统计片段信息，片段信息包括片段长度、CpG个数和依据none CpG位点计算出的转化率，并根据片段长度、CpG个数和转化率筛选片段，筛选条件如下：
，
，
，
其中，是样本在第i个window区间上第j条片段的长度，是预设的片段长度的阈值下限，/>是预设的片段长度的阈值上限；/>是样本在第i个window区间上第j条片段含有的CpG个数，/>是预设的片段上CpG个数的阈值；/>是样本在第i个window区间上第j条片段的转化率，是预设片段转化率的阈值。
8.一种校正全基因组甲基化测序数据批次效应的装置，其特征在于，包括：
输入模块，其被配置为输入健康人血浆样本、肿瘤患者血浆样本和待测样本的甲基化测序数据；
唯一比对区间筛选模块，其被配置为筛选唯一比对区间，所述筛选唯一比对区间包括：从hg19参考基因组中提取每条染色体长度和碱基序列信息，根据碱基序列信息提取CpG的位置；使用GEM-library创建hg19参考基因组、CT转换后的基因组和GA转换后的基因组的mappability文件；根据每个mappability文件，筛选不存在多重比对可能性的区间；取三组区间的交集，合并出唯一比对区间的bed文件；
window区间划分模块，其被配置为根据染色体长度信息、预设的window区间长度和步长划分window区间，同时根据唯一比对率对划分的window区间进行初步筛选，唯一比对率指window区间上唯一比对区间占window区间总长度的比例，筛选规则如下：
，
其中，是第i个window区间内唯一比对率，/>是预设的唯一比对率阈值；
甲基化水平计算模块，其被配置为用于计算区间甲基化水平，如window区间、甲基化稳定区间、Bin区间，
所述甲基化水平包括全甲基化片段占比MFR，MFR指落在某一个区间内的全甲基化片段数量/全部片段数量，全甲基化片段指该片段上所有的CpG位点全部甲基化；或甲基化水平包括甲基化密度MD，MD指落在某一个区间内的甲基化CpG数量/全部CpG数量；
window区间甲基化水平计算包括：根据甲基化测序数据，筛选位于初步筛选的window区间的片段并统计片段信息，计算每个健康人血浆样本和肿瘤患者血浆样本在初步筛选的window区间的甲基化水平，计算方式如下：
区间全甲基化片段占比MFR计算为：
，
其中，是样本n第i个区间的MFR值，/>是样本n在第i个区间内全甲基化的片段数，/>是样本n在第i个区间片段数；
区间甲基化密度MD计算为：
其中，是样本n第i个区间的MD值，/>是样本n在第i个区间内甲基化的CpG个数，/>是样本n在第i个区间的所有CpG个数；
甲基化稳定区间筛选模块，其被配置为筛选满足预设条件的window区间为甲基化稳定区间，包括将健康人血浆样本和肿瘤患者血浆样本在初步筛选的window区间的甲基化水平进行差异化分析，
当甲基化水平用MFR表示时，使用limma包对初步筛选的window区间的MFR进行差异分析，获得p-value后再计算FDR；分别统计肿瘤患者血浆样本和健康人血浆样本的初步筛选的window区间的MFR的最大值、最小值和均值，按照如下条件筛选甲基化稳定区间：
，
，
，
，
其中，是第i个window区间的FDR值，/>是预设的FDR的阈值；是肿瘤患者样本在第i个window区间的最大MFR值，/>是健康人在第i个window区间的最大MFR值，/>是肿瘤患者样本在第i个window区间的最小MFR值，/>是健康人在第i个window区间的最小MFR值，/>是肿瘤患者样本在第i个window区间的MFR均值，/>是健康人在第i个window区间的MFR均值，/>是预设的甲基化稳定区间在肿瘤患者和健康人之间甲基化水平差异的阈值；
当甲基化水平用MD表示时，使用limma包对初步筛选的window区间的MD进行差异分析，获得p-value后再计算FDR；分别统计肿瘤患者和健康人的初步筛选的window区间的MD的最大值、最小值和均值，按照如下条件筛选甲基化稳定区间：
，
，
，
，
其中，是第i个window区间的FDR值，/>是预设的FDR的阈值；是肿瘤患者样本在第i个window区间的最大MD值，/>是健康人在第i个window区间的最大MD值，/>是肿瘤患者样本在第i个window区间的最小MD值，是健康人在第i个window区间的最小MD值，/>是肿瘤患者样本在第i个window区间的MD均值，/>是健康人在第i个window区间的MD均值，/>是预设的甲基化稳定区间在肿瘤患者和健康人之间甲基化水平差异的阈值；
校正模型建立模块，其被配置为以健康人血浆样本的每个甲基化稳定区间的甲基化水平平均值作为因变量，以待测样本的每个甲基化稳定区间的甲基化水平值作为自变量，构建线性回归模型；
全基因组甲基化测序数据校正模块，其被配置为利用校正模型建立模块建立的校正模型校正待测样本的甲基化水平。
9.一种电子设备，包括：一个或多个处理器；存储装置，其上存储有一个或多个程序，当一个或多个程序被一个或多个处理器执行，使得一个或多个处理器实现权利要求1-7任一所述的校正全基因组甲基化测序数据批次效应的方法。
10.一种计算机存储介质，其上存储有计算机程序，其中，计算机程序被处理器执行时实现权利要求1-7任一所述的校正全基因组甲基化测序数据批次效应的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310413887.6A CN116153418B (zh) | 2023-04-18 | 2023-04-18 | 校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310413887.6A CN116153418B (zh) | 2023-04-18 | 2023-04-18 | 校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116153418A CN116153418A (zh) | 2023-05-23 |
| CN116153418B true CN116153418B (zh) | 2023-07-18 |
Family
ID=86360366
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310413887.6A Active CN116153418B (zh) | 2023-04-18 | 2023-04-18 | 校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116153418B (zh) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112735531A (zh) * | 2021-03-30 | 2021-04-30 | 臻和(北京)生物科技有限公司 | 循环无细胞核小体活性区域的甲基化分析方法和装置、终端设备及存储介质 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2020224159A1 (zh) * | 2019-05-06 | 2020-11-12 | 臻和精准医学检验实验室无锡有限公司 | 基于二代测序用于脑胶质瘤的检测panel、检测试剂盒、检测方法及其应用 |
| CN112951418B (zh) * | 2021-05-17 | 2021-08-06 | 臻和(北京)生物科技有限公司 | 基于液体活检的连锁区域甲基化评估方法和装置、终端设备及存储介质 |
| CN114974417A (zh) * | 2021-06-03 | 2022-08-30 | 广州燃石医学检验所有限公司 | 一种甲基化测序方法和装置 |
| KR20240073026A (ko) * | 2021-09-20 | 2024-05-24 | 그레일, 엘엘씨 | 노이즈 영역 필터링을 사용한 메틸화 단편 확률론적 노이즈 모델 |
| CN114171115B (zh) * | 2021-11-12 | 2022-07-29 | 深圳吉因加医学检验实验室 | 一种差异性甲基化区域筛选方法及其装置 |
| CN115064211B (zh) * | 2022-08-15 | 2023-01-24 | 臻和(北京)生物科技有限公司 | 一种基于全基因组甲基化测序的ctDNA预测方法及装置 |
-
2023
- 2023-04-18 CN CN202310413887.6A patent/CN116153418B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112735531A (zh) * | 2021-03-30 | 2021-04-30 | 臻和(北京)生物科技有限公司 | 循环无细胞核小体活性区域的甲基化分析方法和装置、终端设备及存储介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116153418A (zh) | 2023-05-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10741270B2 (en) | Size-based analysis of cell-free tumor DNA for classifying level of cancer | |
| US9965585B2 (en) | Detection of genetic or molecular aberrations associated with cancer | |
| CN106834275A (zh) | ctDNA超低频突变检测文库的构建方法、试剂盒及文库检测数据的分析方法 | |
| CN111647648A (zh) | 一种用于检测乳腺癌基因突变的基因panel及其检测方法与应用 | |
| CN108229103B (zh) | 循环肿瘤dna重复序列的处理方法及装置 | |
| CN104846089B (zh) | 一种孕妇外周血中胎儿游离dna比例的定量方法 | |
| CN115064211B (zh) | 一种基于全基因组甲基化测序的ctDNA预测方法及装置 | |
| He et al. | Assessing the impact of data preprocessing on analyzing next generation sequencing data | |
| CN106845155B (zh) | 一种用于检测内部串联重复的装置 | |
| CN108595918B (zh) | 循环肿瘤dna重复序列的处理方法及装置 | |
| CN108070658B (zh) | 检测msi的非诊断方法 | |
| CN107893116A (zh) | 用于检测基因突变的引物对组合、试剂盒以及构建文库的方法 | |
| EP3608452A1 (en) | Method for constructing amplicon library through one-step process | |
| Ma et al. | The analysis of ChIP-Seq data | |
| CN116356001B (zh) | 一种基于血液循环肿瘤dna的双重背景噪声突变去除方法 | |
| DiGuardo et al. | RNA-seq reveals differences in expressed tumor mutation burden in colorectal and endometrial cancers with and without defective DNA-mismatch repair | |
| CN112980961A (zh) | 联合检测snv、cnv和fusion变异的方法和装置 | |
| CN117253539B (zh) | 基于胚系突变检测高通量测序中样本污染的方法和系统 | |
| CN116153418B (zh) | 校正全基因组甲基化测序数据批次效应的方法、装置、设备和存储介质 | |
| CN110729025B (zh) | 基于二代测序的石蜡切片样本体细胞突变检测方法和装置 | |
| CN109680054A (zh) | 一种低频dna突变的检测方法 | |
| CN110373458B (zh) | 一种地中海贫血检测的试剂盒及分析系统 | |
| CN106815491B (zh) | 一种用于检测ffpe样本基因融合的装置 | |
| CN112410422A (zh) | 基于片段化模式预测肿瘤风险值的方法 | |
| CN115831355A (zh) | 多癌种wgs的肿瘤早期筛查方法 |
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 |