To identify robust DNA methylation-based biomarkers, we searched the public datasets for the DNA methylation status in CRC cases and collected 1,104 CRC tumors, 268 para-tumor samples, and 54 adenomas for further analysis (See section “Materials and Methods,” Supplementary Table 1). We also obtained the complete list of genes belonging to the ZFP family (Supplementary Table 2). Based on the feature selection procedures described in the section “Materials and Methods,” we finally identified five candidate genes: ESR1, ZNF132, ZNF229, ZNF542, and ZNF677 (Figure 1).

These five genes were significantly hypermethylated in both CRC and adenoma tissues compared to the para-tumors (Supplementary Figure 1). Consistently, the expression levels of these genes were also significantly down-regulated in CRC tumors compared to para-tumors in the TCGA dataset (Figure 2). To characterize the abilities of these biomarkers quantitatively in the combined discovery dataset, we constructed a univariate logistic regression model for each gene and obtained robust discrimination between CRC tumors and para-tumor tissues (sensitivity = 0.82–0.90, specificity = 0.88–0.97, AUC = 0.93–0.97).

To further verify these DNA methylation-based biomarkers, we recruited 104 CRC patients of the Han Chinese population. The characteristics of these CRC patients (replication cohort 1) are shown in Table 1. The CRC tumors and their matched para-tumors were obtained for targeted bisulfite sequencing. The methylation profiles of the candidates and the known biomarker (SEPT9) were examined using targeted bisulfite sequencing (See section “Materials and Methods”). The bisulfite conversion rate (C to T) was high (>99%) in both tumors and para-tumors, and no significant difference in the read mapping rate was found between groups (Supplementary Figure 2). After quality control, 187 samples were retained for further analysis, including 98 CRC tumors and 89 para-tumors.

Based on the methylation profiles of these five ZFGs, we performed principal component analysis (PCA) and found a significant distinction between CRC tumors and para-tumors (Supplementary Figure 3). A differential methylation analysis was also conducted for these candidates, and we found that all candidate genes were significantly hypermethylated in CRC tumors (Figure 3 and Supplementary Figure 4). Furthermore, we performed a univariate logistic regression analysis without adjusting for covariates, and created the receiver operating characteristic (ROC) curve to reveal the diagnostic ability of each candidate gene. As shown in Supplementary Table 3, the area under the curve (AUCs) of these candidates ranged from 0.85 to 0.93. In particular, we found that the diagnostic capability of ZNF132 (AUC = 0.91) and ESR1 (AUC = 0.93) was equal to or better than that of SEPT9 (AUC = 0.91), indicating that these ZFGs might have great potential for CRC diagnosis.

In addition to the overall differential methylation analysis, we also evaluated the effects of age, gender, tumor stage, location of the tumor (colon or rectum), and the mutation status (KRAS + vs. KRAS-) of the CRC samples (Supplementary Tables 4–8). No significant differences were found in the diagnostic capability between CRC patients in the young/old, male/female, early/late stage, and colon/rectum subgroups (Supplementary Table 9, see section “Materials and Methods”). However, the diagnostic capability for CRC patients carrying the KRAS mutation (KRAS +) was significantly better than that for KRAS- samples (Supplementary Tables 8, 9). In the KRAS + group, the sensitivity of each gene ranged from 0.86 to 0.98, the specificity range was 0.89 to 1.00, and the AUC range was 0.97 to 1.00, which is significantly higher than the sensitivity (0.52 to 0.82), specificity (0.80 to 0.98) and AUC (0.74 to 0.86) in the KRAS- group (Supplementary Table 8). Moreover, the hierarchical clustering analysis revealed that all CRC tumors misclassified into the control group were KRAS- samples (Figure 4). In addition, we further examined if other factors may be associated with these misclassified CRC tumors (n = 17). However, none of the age, gender, tumor stage, and tumor locations showed significant differences in our study’s overall CRC tumors. Furthermore, we examined the diagnostic abilities of the ZFGs in the KRAS + and KRAS- subgroup from the TCGA dataset for verification. Consistently, the diagnostic capability of the ZFGs in the KRAS + subgroup (AUC = 0.96–1) was significantly higher than in the KRAS- subgroup (AUC = 0.93–0.97), suggesting the significant methylation changes between KRAS + and KRAS- subgroup (Supplementary Table 10).

Due to the limited sample size in replication cohort 1, we recruited another sample group (replication cohort 2), including 114 CRC tumors and matched para-tumors from patients of Han Chinese descent (Table 1) for further validation. The methylation profiles of these ZFGs were measured in replication cohort 2 (Supplementary Figures 5, 6). Consistently, we found that ESR1 (AUC = 0.90) and ZNF132 (AUC = 0.94) still showed the best diagnostic capability for CRC (Supplementary Table 11).

Furthermore, the significant differences of methylation profiles between the KRAS + and KRAS- subgroups were also validated (Supplementary Figure 7). In the KRAS + subgroup, the sensitivity (0.90 to 1.00), specificity (0.91 to 0.98), and the AUC (0.92 to 1.00) of each gene was significantly higher than that in the KRAS- subgroup (sensitivity: 0.58 to 0.85, specificity: 0.78 to 0.96, AUC: 0.71 to 0.88) (Supplementary Table 12). Similarly, the CRC tumors misclassified as the para-tumors were in the KRAS- subgroup, confirming that the KRAS + CRC samples were more epigenetically homogeneous than the KRAS- CRC samples (Figure 5).

To further identify the robust diagnostic biomarker for CRC diagnosis, we combined the datasets from replication cohort 1 and replication cohort 2. As shown in Figure 6, ZNF132 had the highest diagnostic capability (sensitivity = 0.83, specificity = 0.97, AUC = 0.93) than other genes, including SEPT9 (sensitivity = 0.83, specificity = 0.87, AUC = 0.91). In addition, ESR1 also achieved comparable diagnostic capability (sensitivity = 0.78, specificity = 0.97, AUC = 0.91) to SETP9 (Table 2), suggesting that these two genes have great potential in the liquid biopsy of CRC. Simultaneously, we assessed the diagnostic abilities of these ZFGs in the KRAS + subgroups and KRAS- group using all samples. Concordantly, all candidate biomarkers achieved better performance (AUC ≥ 0.95) in the KRAS + group than that in the KRAS- group, further indicating that KRAS- CRC patients are more epigenetically homogeneous.

To further improve the robustness of the diagnostic abilities of these ZFGs, we performed multiple machine-learning algorithms following the fivefold cross-validation method, which split the dataset into training and validation dataset randomly to obtain an unbiased estimation (see section “Materials and Methods,” Supplementary Table 13). As shown in Supplementary Table 13, the random forest (RF) model achieved the best accuracy (accuracy = 0.89) in the test dataset. In the KRAS + subgroup, the Naïve Bayes algorithm performed best on the test dataset (sensitivity = 0.99, specificity = 0.97, accuracy = 0.98), while the Neural Network model had the best performance in the KRAS- subgroup (sensitivity = 0.75, specificity = 0.89, accuracy = 0.82). Therefore, we suggested that the diagnostic abilities of these ZFGs are reliable for CRC diagnosis in test dataset.

Public high-throughput DNA methylation microarray datasets (Illumina HumanMethylation 450K) were searched exhaustively from TCGA and GEO databases. Two datasets from TCGA and nine datasets from GEO were included, yielding a total of 1,104 CRC tumor samples, 268 control samples, and 54 adenoma samples (Supplementary Table 1). The comprehensive list of genes in the zinc finger family (n = 1,594) was obtained from HGNC (Supplementary Table 2).

As shown in a previous study, the methylation linkage equilibrium decreased significantly when the block was longer than 1,000 bp (Guo et al., 2017a). Therefore, we defined the methylation region (MR) to have at least 6 CpG sites with less than 1,000 bp. We then arranged all the CpG sites in the high-throughput microarray according to their genomic coordinates and applied the sliding window method to identify all MRs. In total, we identified 6,166 MRs and performed differential methylation analyses (Supplementary Table 14). The 85 DMRs (McaM > 0.50, MadM > 0.50 and McoM < 0.30) were found for further analysis. To correct the noise from DNA originating from mixed-in peripheral blood, we suggested that the methylation rate of the DMRs need to be extremely low in the peripheral blood. Thus, we integrated the public high-throughput microarray datasets of the whole blood (WB, n = 1438), peripheral blood mononuclear cells (PBMC, n = 111), and peripheral blood leukocytes (PBL, n = 529) as references for DMR identification (Supplementary Table 15). The 32 DMRs were retained due to their low methylation rates in WB, PBMC, or PBL (mean methylation rate < 0.10).

Several DMRs were located at the same gene, and we selected the DMR with the most significant differences between CRC and control tissues for each gene (N = 10). To obtain the DMRs that may regulate the expression of neighboring genes, we further selected the DMRs with transcription factor binding sites (TFBSs) that correlate significantly with the expression of neighboring genes (Supplementary Tables 16, 17). In total, seven out of the ten candidate DMRs were selected for validation. However, due to the difficulties in the primer design due to CG percent, PolyT, and the number of SNPs, two candidate DMRs were removed (SALL1 and ZSCAN23). Finally, we obtained the top five candidate DMRs for further validation (ESR1, ZNF132, ZNF229, ZNF542, and ZNF677).

Colorectal cancer patients in both replication cohort 1 and cohort 2 were recruited from the Affiliated Hospital of Nantong University between 2016 and 2018. The earlier recruited samples (N = 104) were included in replication cohort 1, while the others belonged to replication cohort 2. The recruited patients with available CRC tumors and matched para-tumors (the corresponding adjacent normal tissue at least 5 cm distant from the tumor tissue) were selected in this study. The patients recruited had not been treated with any neoadjuvant therapy before. Tumors were classified according to TNM (Tumor Node Metastasis)/UICC (Union for International Cancer Control) criteria following the histopathological examination. At least two professional pathologists evaluated all tumor samples carefully. All procedures performed in this study were in accordance with the ethical standards of the institutional research committee, as well as the 1964 Declaration of Helsinki and its later amendments. The study was approved by the institutional review boards of the Affiliated Hospital of Nantong University. Written informed consent was obtained from each participant of the study. All tumors and para-tumors were immediately frozen at −80°C after surgical resection.

DNA extraction was performed using the AIIperp DNA/RNA Mini Kit (Qiagen, Duesseldorf, Germany) according to the manufacturer’s protocols. The EpiTect Fast DNA Bisulfite Kit (Qiagen, Duesseldorf, Germany) was further used for bisulfite conversion. After carefully evaluating CG percent, PolyT, and the occurrence of SNPs in the targeted regions of the candidate DMRs, we designed primers to detect them in a panel for NGS sequencing (Supplementary Table 18). The PCR amplicons were then diluted and amplified using these primers, and the products (170–270 bp) were further separated and purified by QIAquick Gel Extraction kit (Qiagen, Duesseldorf, Germany). Finally, the sequencing libraries from different samples were pooled together and sequenced using the Illumina Hiseq 2000 platform according to the manufacturer’s protocols.

The BSseeker2 was applied for read mapping and methylation calling. After that, we removed the samples with a bisulfite conversion rate < 98%. The methylation level of each CpG site for each sample has been provided in Supplementary Tables 19, 20. The average coverage and missing rate for each CpG site were calculated and utilized for quality control (average coverage > 20X, missing rate < 20%). In addition, the sample whose missing rate > 30% were also filtered out. FastTarget next-generation sequencing was used to detect tumor DNA for the mutations in codons 12, 13, 59, 61, 117, and 146 of the KRAS and NRAS genes, as well as the mutation in codon 600 of the BRAF gene (Liao et al., 2019).

The Wilcoxon rank-sum test was performed in the discovery stage to identify the differential methylation sites and regions between CRC tumors, adenomas, and para-tumors. Moreover, the differential methylation status (odds ratios) between tumors and para-tumors of the DMRs were calculated with logistic regression. Benjamini-Hochburg correction was utilized for multiple test correction.

To identify the factors that may affect the diagnostic abilities of these ZFGs in CRC, we first split the dataset based on these covariates, including young (age ≤ median age) vs. old (age > median age), male vs. female, early (Stage I/II) vs. late (Stage III/IV), colon vs. rectum and KRAS + vs. KRAS- subgroups. Second, we performed the univariate logistic regression model using each gene in these subgroups separately. The cut-off value was determined as the value which maximized the sum of both sensitivity and specificity of the model. After that, the predicted sample type for each sample was obtained, and was used to calculate the number of accurate predictions (true positive + true negative) and inaccurate predictions (false positive + false negative) of the logistic model. Finally, we utilized the Fisher exact test to explore if there are significant differences between the prediction outcomes between these subgroups.

To further give an unbiased estimation of the diagnostic ability of these ZFGs, we performed fivefold cross-validation of our dataset using the logistic regression (Package stats), support vector machine (SVM, Package e1071), random forest (Package randomForest), Naïve Bayes (Package e1071), neural network (Package nnet), linear discriminant analysis (LDA, Package mda), mixture discriminant analysis (MDA, Package mda), flexible discriminant analysis (FDA, Package mda), gradient boosting machine (Package gbm), catboost (Package catboost), and XGBoost (Package xgboost) methods. We randomly selected 80% of our samples as the training dataset, and the remaining 20% of samples were used as the test dataset. After that, we optimized these machine learning methods using the training dataset and validate these optimized models in the test dataset. We repeated these procedures 1,000 times, and obtained the sensitivities, specificities and accuracies for both training and test datasets. The averaged sensitivities, specificities and accuracies for both training and test datasets were calculated and shown in Supplementary Table 13. All statistical analyses were conducted using R (v3.4.3).