Discovering methylation markers and development of a sense-antisense and dual-MGB probe PCR assay in plasma for colorectal cancer early detection

2.1. Data preparation

Thirteen methylation datasets, which were GSE77954, GSE101764, GSE131013, GSE199057, GSE164811, GSE193535, GSE129364, GSE139404, GSE107352, GSE75546, GSE77965, GSE68060 and EMTAB6450 from public databases were selected based on three criteria: 1) were generated by Illumina HumanMethylation 450k BeadChip and had the raw IDAT files, 2) the sample size was greater than 10, and 3) consisted of CRC or adenoma or normal adjacent tissue (NAT). All the IDAT files were then processed using the minfi tool [16] to obtain methylation β values. They were integrated as a single dataset (n= 1165), which we defined as Phase I discovery set for candidate marker identification.

Level 3 methylation data for 31 cancer types were retrieved from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Data for 8258 primary tumor tissues corresponding to 31 cancer types and their 710 NATs were retained. The 31 cancer types included ACC (n=79), BLCA (n=412), BRCA (n=778), CESC (n=306), CHOL (n=36), CRC (n=379), DLBC (n=48), ESCA (n=183), GBM (n=137), HNSC (n=523), KICH (n=65), KIRC (n=312), KIRP (n=271), LGG (n=513), LIHC (n=374), LUAD (n=456), LUSC (n=364), MESO (n=87), OV (n=10), PAAD (n=183), PCPG (n=178), PRAD (n=495), SARC (n=257), SKCM (n=104), STAD (n=393), TGCT (n=133), THCA (n=503), THYM (n=124), UCEC (n=418), UCS (n=57), UVM (n=80). The TCGA dataset was defined as Phase II discovery dataset. DMCs were deemed eligible if they exhibited methylation levels <0.2 on other cancer types (non-CRC), >=0.55 on CRC, and <0.15 on 710 NATs.

The other three datasets, GSE48684 [17], GSE40279 [18], and GSE122126 [19], generated by the same platform of Illumina HumanMethylation BeadChip in Gene Expression Omnibus (GEO) were used as validation sets to verify the methylation status of candidate differentially methylated CpGs (DMCs). The GSE48684 cohort consisted of 105 qualified samples, of which 41 were NAT and 64 were CRC. The GSE40279 cohort comprised whole blood cell (WBC) samples collected from 656 healthy individuals. The GSE122126 cohort included three CRC and 12 normal plasma samples.

2.2. Sample collection

This case-control study enrolled 772 cases, including 60 tissue samples (30 CRC and 30 NATs) and 712 plasma samples from the First Affiliated Hospital of Zhengzhou University between April 2022 and June 2022. Tissue samples were obtained from the preserved formalin-fixed paraffin-embedded (FFPE) sections collected from surgical patients. The plasma cohort was randomly divided into a training set and a validation set in a 2:1 ratio. The training set comprised 474 participants, including 115 healthy blood donors, 123 non-digestive disease patients (NDD), 65 intestinal disease patients (ID), 14 polyps patients, 10 non-advanced adenoma patients (Non-AA), 16 AA, 125 CRC patients, and 6 patients with other cancers. In the validation set, there were 235 participants, including 57 healthy blood donors, 60 NDD, 30 intestinal ID, 6 polyps patients, 5 Non-AA, 9 AA, 66 CRC patients, and 2 patients with other cancers (Supplemental Table 1). This study is a sub-project of the Clinical Study of Pan-cancer DNA Methylation Test in plasma (Clinical Trials ID: NCT05685524). Patients with polyps, adenomas, CRC and other cancers were confirmed by histopathological examinations. The included CRC patients were required to meet the following criteria: 1) did not receive any radiotherapy, chemotherapy or surgery, 2) ages older than 18. All CRC patients were classified as I, II, III, and IV stages according to the American Joint Cancer Committee (AJCC) staging system.

2.3. Differential methylation analysis

We performed differential methylation analysis using the rank-sum test on three groups in the discovery set: cancer vs. normal (NAT), adenoma vs. normal, and cancer vs. adenoma. DMCs were defined as significant if they meet two criteria: an FDR < 0.05 and a fold change ranking in the top 1% of all probes. The DMCs were then classified as hyper- or hypo- DMCs based on whether they exhibited high or low methylation levels in the tumor or adenoma samples compared to the normal samples. LASSO regression was implemented in the R package ‘glmnet’ to reduce the number of features. LASSO regression was repeated 100 times, and the frequencies of probes with non-zero coefficients in the regressions were counted.

2.4. Tissue and plasma DNA extraction

Genomic DNA from tissue samples and plasma cfDNA was isolated using the FFPE DNA rapid extraction kit (TianGen, Beijing) and the plasma/serum cfDNA extraction kit (Ammunition Life-tech, Wuhan), respectively, according to the protocols. Briefly, cfDNA extraction was divided into two steps. First, approximately 10 ml of blood was drawn from participants at room temperature and then centrifuged twice with 1300g and 14000g for 10 min at 4°C within 2 hours. Then, about 2 ml of plasma was retained for cfDNA purification. Purified cfDNA was eluted in 45 ul TE buffer and then treated with bisulfite using a DNA bisulfite modification kit (Ammunition Life-tech, Wuhan) in accordance with the instructions described in the aforementioned study [20], or stored at -80 °C until required.

2.5. Sanger sequencing

Sanger sequencing (ABI 3730XL, GENEWIZ, Suzhou) was performed for the target regions of NTMT1 and MAP3K14-AS1 after bisulfite treatment to confirm their methylation status on 30 CRC tissue samples and 30 NATs. The sequencing primers are shown in Supplemental Table 2. Following bisulfite treatment, cytosine in methylated CpG sites remained cytosine, while unmethylated cytosine was converted to thymine. Sanger sequencing results thus allowed for the direct determination of the methylation status of target regions.

2.6. Developing the SADMP and methylation-specific PCR (MSP) technique

The location and sequences of MSP primers/ probes are displayed in Supplemental Table 3. For NTMT1, two pairs of methylation primers targeted the sense and antisense strands, and their corresponding MGB probes were designed. For MAP3K14-AS1, the MGB-probe 1 and MGB-probe 2 were designed according to the antisense strand sequence and the reverse complementary sequence of the antisense strand.

Gradient dilutions of synthetic plasmids or plasmid mixtures are employed as templates for the determination of test parameters. Specifically, standard curve experiments were carried out to assess the amplification efficiency of each pair of primers. The experiments consisted of five ten-fold dilutions of fully methylated plasmid DNA templates with five replicates at each dilution (1, 10, 100, 1000, and 10000 copies per run). Primer tolerance experiments were designed to assess the specific amplification ability of methylated primers against methylated templates. Eight concentration dilutions of fully methylated DNA templates (0, 1, 5, 10, 50, 100, 200, and 400 copies/μL) mixed with a high concentration of unmethylated plasmid (10⁷ copies/μL) were prepared to evaluate the ability of methylated primers selectively amplifying methylated templates under the background of unmethylated templates. Each concentration experiment was repeated five times.

The amplification system of MSP was shown in Supplemental Table S4. The test employed a triplex PCR with a FAM channel for the internal reference gene ACTB, a ROX channel for NTMT1 (dual-strands), and VIC channel for MAP3K14-AS1 (dual-MGB probes). The PCR procedure was pre-denaturation at 95°C for 5 min (step 1), denaturation at 95°C for 15s (step 2), annealing at 60°C for 30s (step 3), repeating step 2 and 3 forty-five times in plasmid samples and fifity times in plasma samples on the 7500 Fast Real-Time PCR System (Applied Biosystems, USA), the cycle threshold (Ct) was assigned to 45 or 50 for the target without amplification.

2.7. Statistical analysis

The data processing and analysis in this study were conducted using R software (version 4.1.0). The ‘glm’ function was employed to fit the logistic regression model with a parameter of ‘family=binomial’. Receiver operating characteristic (ROC) curve analysis was performed using the ‘pROC’ package, and the area under the ROC curve (AUC) was subsequently calculated to assess the test classification performance. The optimal sensitivity and specificity of the test were estimated when Youden’s index reached its maximal value. Sensitivity and specificity were calculated using the following formulas: And the Youden index = sensitivity+specificity-1.

The rank-sum test and Kruskal test were employed for the comparisons between two groups and the comparisons between multiple groups, respectively. The Chi-square test was utilized for the comparisons between categorical variables. The Other statistical methods used in this study were described in the corresponding results.

3.1. Study design and participant characteristics

The study is divided into four steps, as illustrated in Figure 1. In the initial step, candidate markers were identified in a set of 1165 samples from 13 datasets. Of these, 452 were normal, 168 were adenomas, and 545 were CRC. These markers were then narrowed down using data from 31 cancer types in the TCGA dataset, DMCs that were hypermethylated in cancer other than CRC were filtered out. The available DMCs were further validated in three independent datesets. In the second step, the methylation status of candidate CpGs was confirmed by Sanger sequencing using 30 CRC tissues and 30 NATs. In the third step, the SADMP technique and MSP system were established. This involved designing appropriate primers, optimizing the qPCR amplification system, and evaluating the technical parameters of the assay. In the last step, the developed assay was tested on training and validation sets of plasma. Its performance was assessed by estimating indicators such as sensitivity, specificity, and AUC.

• Download figure
• Open in new tab

Figure 1. The flowchart of this study. The four steps were illustrated by different colored blocks.

3.2. Landscape of the methylation patterns of the 13 discovery datasets

The methylation levels of adenoma and cancer samples were found to be lower than normal, as shown in Figure 2A. The methylation density curves for the three groups displayed bimodal distributions with hyper- and hypomethylation peaks, respectively (Figure 2B). The hyper-methylation peaks of adenoma and tumor samples were lower than those of normal samples, indicating higher methylation levels in normal samples.

• Download figure
• Open in new tab

Figure 2. Landscape of the methylation patterns of the phase I discovery set.

(A) boxplot showing the overall methylation levels of normal, adenoma and cancer samples. The average β value of all probes for each sample was calculated as the sample overall methylation level. P-values were estimated by Kruskal test. (B) Density curves of probe methylation levels in normal, adenoma and cancer samples. (C) t-SNE visualizing normal, adenoma and cancer samples in the discovery set. (D) Heatmap showing the most variable probes between normal, adenoma and cancer samples.

We employed t-SEN to analyze and visualize the structure of the discovery set and identified significant differences between normal and cancer samples. However, adenomas exhibited overlap with both normal and cancer samples (Figure 2C). We then selected the top 1% of most variable probes to cluster the discovery set using the K-means algorithm. The results showed that both normal and cancer samples clustered together, while adenomas were separated into two subgroups with a remarkable difference. They showed both high (Methy-H) and low (Methy-L) methylation status and were closer to cancer and normal samples, respectively (Figure 2D).

Further analysis revealed that tubular adenomas constituted the largest proportion of Methy-L adenomas (23/57 or 40.35%), while villous adenomas constituted the largest proportion of Methy-H adenomas (47/73 or 64.38%). To eliminate any potential bias between the datasets, we conducted Fisher’s test separately for the Methy-H and Methy-L adenomas, with the datasets serving as the stratification factor. However, we did not identify any significant differences (P > 0.05).

3.3. Identification of candidate DMCs

First, we conducted a comparative analysis of the DMCs between cancer, adenoma, and normal samples using the 13 discovery datasets. The three comparisons (cancer vs. normal, adenoma vs. normal, and cancer vs. adenoma) yielded 3000, 3051, and 1545 DMCs, respectively. Notably, these DMCs were predominantly hyper-DMCs (Figure 3A). Further investigations revealed that the majority of DMCs were distributed in upstream regions of genes (5’UTR, TSS1500, TSS500, and 1st Exon), except for cancer vs. adenoma where DMCs were mainly located in gene bodies (Figure 3B). We observed an extremely high proportion of overlapping DMCs between cancer vs. normal and adenoma vs. normal, which was significantly higher than cancer vs. adenoma (Figure 3C). To identify the most appropriate DMCs, we focused on those that were overlapped between cancer vs. normal and adenoma vs. normal (2237 probes) and evaluated their methylation levels on 31 cancer types of TCGA. This yielded 75 accessible probes (Figure 3D). LASSO regression was employed to reduce the number of DMCs, with eight presented in 100 replicates with non-zero coefficients (Figure 3E).

• Download figure
• Open in new tab

Figure 3. Identification of candidate markers.

(A) Percentage of hyper-DMC and hypo-DMC between the three comparisons. (B) Percentage of DMC at different genomic regions between the three comparisons. (C) Venn diagram showing DMCs between the three comparisons. (D) Methylation values of 75 probes meeting the criteria on TCGA 31 cancer types. The “other” refers to 30 non-CRC cancer samples. “Normal” refers to 710 NATs. (E) Frequency of probes with non-zero coefficient in 100 LASSO regressions.

3.4. The methylation levels of NTMT1 and MAP3K14-AS1 in validation sets

The methylation levels of the eight probes were then analyzed in 656 WBC samples, and seven DMCs exhibited relatively low methylation levels (less than 0.1) with the exception of cg17892556 (Figure 4A). Therefore, the seven DMCs are recognized as the most promising markers.

• Download figure
• Open in new tab

Figure 4. Methylation profiles of the two candidate probes in different datasets.

(A) The methylation profiles of cg14015706 and cg08247376 in 656 healthy WBC in GSE40279. (B) The methylation β values of cg14015706 and cg08247376 between normal and cancer samples in GSE48684 dataset. (C) The cfDNA methylation profiles of cg14015706 and cg08247376 in normal and cancer plasma samples from GSE122126. Numbers in the heatmap indicated the methylation β values.

Two specific DNA methylation sites, cg14015706 and cg08247376, were selected for further analysis as they were suitable for primer/probe design. These sites are located in the first exon of NTMT1 and 200 bp upstream of the MAP3K14-AS1 transcription start site, respectively. The GSE48684 dataset revealed significantly higher methylation levels for both probes in cancer samples than in normal samples (Figure 4B). Furthermore, the two probes demonstrated hypermethylation in CRC plasma and hypomethylation in healthy plasma (Figure 4C), indicating a high consistency of the methylation status between tissue and plasma.

3.5. The SADMP technique

The Sanger sequencing of the target regions revealed that the CpGs sites among the amplicons were more frequently methylated in CRC than NATs, which provided the foundation for the development of SADMP technique. For NTMT1, two pairs of methylation primers targeted the sense and antisense strands, and their corresponding MGB probes were designed (Figure 5A). For MAP3K14-AS1, the MGB-probe1 and MGB-porbe2 were designed according to the antisense strand sequence and the reverse complementary sequence of antisense strand (Figure 5F). The estimated amplification efficiencies of NTMT1 and MAP3K14-AS1 were 98.14% and 110.08%, respectively (Figure 5B and 5G). The sense and anti-sense strands were transformed into completely different DNA sequences after bisulfite treatment and both could be used as PCR templates. It can be postulated that the copy numbers detected by any single pair of primers or single MGB-probe assay should be half of the diluted copy numbers calculated as the sum of sense and anti-sense strand templates, and that the SADMP assay should be equal to the diluted copy numbers. For NTMT1, the detected copy numbers by sense and antisense strand assays were slightly less than half of the theoretical copy numbers with slopes of 0.39 and 0.31, respectively. The detected copy numbers by dual-strand assay were approximately 0.73 times of theoretical copy numbers, which was approximately two-fold compared to any single-strand assay (Figure 5C), as expected. Meanwhile, Ct value of dual-strand assay was almost one cycle earlier than that of any single-strand assay (average △Ct=1.21, Figure 5D), which was consistent with the theoretical △Ct value (should be one). The MAP3K14-AS1 SADMP showed similar results, with the number of copies detected was doubled compared to the single probe assay (Figure 5H), and the Ct values were shifted forward by 0.99 cycles (Figure 5I).

• Download figure
• Open in new tab

Figure 5. The development and validation of SADMP.

(A) SADMP strategy for NTMT1. Two pairs of primers were designed according to the sense and antisense BS-strand sequences. CBS indicated the complementary sequence of sense and antisense BS-strands. (B) The standard amplification curve of NTMT1. (C) The agreement between detected copies by SADMP and single-strand assays and diluted template DNA copies (theoretical copies) of NTMT1. The solid lines were fitted by a simple linear model. (D) and (E) The detected Ct values of SADMP NTMT1 and single-strand NTMT1 assays for different diluted concentrations. (F) SADMP strategy for NTMT1. Two MGB probes were designed. Probe 1 is located downstream of the forward primer, targeting the BS-strand template. Probe 2 is located downstream of the reverse primer, targeting the complementary sequence of BS-strand. (G) The standard amplification curve of MAP3K14-AS1. (H) The agreement between detected copies by SADMP and single-probe assays and diluted template DNA copies (theoretical copies) of MAP3K14-AS1. The solid lines were fitted by a simple linear model. (I) and (J) The detected Ct values of SADMP MAP3K14-AS1 and single-probe MAP3K14-AS1 assays for different diluted concentrations.

Additionally, in the serially diluted experiment, no amplification cures for both targets were observed when a high proportion of unmethylated DNA was used as templates (the methylated DNA copies were 0) (Figure 5E and 5G). These results suggested that the developed assays were explicitly targeted for methylated DNAs even at the high background of unmethylated DNAs. We also found that the NTMT1 SADMP assay was able to robustly detect methylated DNAs at a concentration of 10 copies/μL (Figure 5E), while it was 5 copies/μL for the MAP3K14-AS1 SADMP (Figure 5G).

3.6. Validation and evaluation of biomarkers using MSP in plasma samples

In the plasma sample set, the Ct values of ACTB exhibited a slight decrease in trend from the healthy control groups to the cancer group (p < 0.0001) (Figure 6A). Ct values of NTMT1 and MAP3K14-AS1 were much lower in CRC samples and AA samples than other control groups (Figure 6B and 6C).

• Download figure
• Open in new tab

Figure 6. dual-target assessment and validation in plasma samples.

(A-C) Ct values of ACTB (A), NTMT1 (B) and MAP3K14-AS1 (C) in various clinical groups. NDD: Non-digestive disease, ID: Intestinal disease, AA: advanced adenoma. The eight groups were compared by Kruskal-Wallis test. P<0.05 is considered significant. **** is p<0.0001. (D-E) ROC curves of two biomarkers in CRC vs Control groups (D) and CRC+AA vs Control groups (E) in training set. (F-G) Sensitivity and specificity of individual biomarker and dual-target in training set (F), validation set (G) and total set (H) using 1/2 algorithm. Error bars represent 95% CI.

The Ct values could discriminate the CRC and AA samples from healthy samples very well. Furthermore, the Ct values-based result determination exhibited favorable outcomes in a previous study [21,22]. Consequently, the Ct values were utilized as the result analysis indicators. In order to determine the algorithm and cutoff, ROC curve analysis was conducted in the training set. In the CRC vs control groups, the AUC values for NTMT1 and MAP3K14-AS1 were 0.882 (95% CI: 0.838-0.925) and 0.785 (95% CI: 0.729-0.840), respectively (Figure 6D). In the CRC+AA vs control groups, the AUC values for NTMT1 and MAP3K14-AS1 were 0.852 (95% CI: 0.807-0.898) and 0.761 (95% CI: 0.707-0.815), respectively (Figure 6E).

As AA was precancerous lesions, it was also benefit patients from detecting it. So we use the CRC+AA vs control groups of training set for further analysis. The performance of the combination of the two biomarker were evaluated by two strategies (Table 1). Strategy 1 involved the construction of a logistic regression model, with an estimated AUC value was 0.881 (95% CI: 0.840-0.922), the optimal sensitivity for CRC and specificity of 86.4% and 91.0%, respectively (Table 1). Strategy 2 was the 1/2 algorithm, which determined a positive measurement when the Ct of any single marker was less than its corresponding threshold. The cutoff was determined as Ct value corresponding to the maximum Youden index in the training dataset using CRC+ AA vs control data. The cutoff values were 48.2 and 48.4 for NTMT1 and MAP3K14-AS1, respectively. At these cutoff, the two target sensitivities were 80.0% and 59.2%, with specificities of 92.8% and 97.0%, respectively (Figure 6F). The sensitivity and specificity of the combination of the two biomarkers were 86.4% and 91.0% respectively. Notably, strategy 2 yielded an identical sensitivity and specificity to strategy 1. Given that strategy 2 is considerably more straightforward for physicians to interpret test results in clinical practice, we have adopted the 1/2 algorithm as the combination algorithm for the dual-target test.

After the algorithm and cutoff was set, the validation set and all sample set were analyzed. The validation set had a 81.8% sensitivity and 92.5% specificity (Figure 6G) and the all sample set had a 84.8% sensitivity and 91.5% specificity (Figure 6H). The dual-target test had a higher sensitivity than that of any single biomarker and a slightly lower specificity in the three data sets (Figure 6F-H). It indicated that the combination of the two biomarkers were better than single biomarker.

3.7. Performance in subgroups of plasma samples

The performance of the dual-target test was evaluated in subgroups of all plasma samples. The test demonstrated a detection rate of 26.7% in non-AA samples and a detection rate of 32.0% in AA samples (Figure 7A). It exhibited sensitivities of 75.0%, 81.2% and 90.3% in stages I, II and III-IV of CRC samples, respectively (Figure 7B). The test exhibited specificities of 95.9%, 92.3%, 84.2% and 90.9% in healthy, NDD, ID and polyps samples, respectively (Figure 7C). The positive predictive value (PPV), negative predictive value (NPV) and accuracy for the dual-target test were 79.4%, 94.0% and 89.6%, respectively (Figure 7D).

• Download figure
• Open in new tab

Figure 7. Sensitivity, specificity in various clinical groups and other performance indicators of single biomarker and dual-target in total plasma sample set.

(A-B) Sensitivity in Non-AA and AA groups (A), stages subgroups of CRC (B). (C) Specificity in healthy, NDD, ID and polyps groups. (D) PPV, NPV and accuracy of individual biomarkers and dual-target in CRC vs control groups. Control groups include healthy, NDD, ID, polyps, Non-AA and other cancers. Error bars represent 95% CI. NDD: Non-digestive disease, ID: Intestinal disease, AA: advanced adenoma.

With regard to the detection rates in adenoma and different stages of CRC, the dual-target test demonstrated superior performance compared to that of a single biomarker (Figure 7A-B). With respect to specificities in subgroups, the dual-target test exhibited a marginal decline in value compared to that of a single biomarker (Figure 7C). However, with regard to NPV and accuracy, the dual-target test demonstrated superior performance compared to that of a single biomarker (Figure 7D). With regard to PPV, the dual-target test exhibited a slight decline in performance relative to that of a single biomarker (Figure 7D).