In total, 86 DLBCL patients were enrolled from the multi center studies including Peking University Third Hospital, Fifth Medical Center of PLA General Hospital, and Cancer Hospital Chinese Academy of Medical Sciences, from 2017 to 2019. All patients had signed the patient consent form. In all cases, the diagnosis of DLBCL was made using appropriate diagnostic criteria from the 2016 WHO classification of lymphoid tumors with combinations of histologic, immunohistochemical, and cell of origin (Coo) defined according to the Hans algorithm (Fang et al., 2010). This study was conducted in accordance with the Declaration of Helsinki.

Peripheral blood samples (8–10 mL) from DLBCL patients were obtained through routine intravenous blood sampling and collected into a Cell-Free DNA collection tube (Roche). Plasma separation was performed within 24 h, Whole blood samples were centrifugation twice at 4°C 1350×g for 12 min and once at 4°C 13,500×g for 12 min. The plasma layers were then transferred to a new tube. Then the plasma samples were immediately stored at −80°C for future use. The plasma cfDNA was extracted from 2–4 mL plasma using the Quick-cfDNA Serum & Plasma Kit (ZYMO) and then stored at −80 °C. The concentration and quality of cfDNA were quantified by Qubit fluorometer and nucleic acid electrophoresis before library preparation.

5hmC libraries for all cfDNA samples were constructed using the high-efficiency hmC-Seal technology described previously (Song et al., 2011). Given the high sensitivity of the chemical labeling method, we assigned low values of 1–10 ng for the input cfDNA. The plasma-derived cfDNA was subjected to end-repairing and then ligated with sequencing compatible adaptors. 5hmC bases were biotinylated via two-step chemistry and purified by the DNA Clean & Concentrator Kit (ZYMO), and subsequently enriched by binding to Streptavidin beads (Life Technologies). Then, the beads were re-suspended in RNase-free water and amplified with 14–16 cycles of PCR amplification. Finally, the PCR products were purified using AMPure XP beads (Beckman). All libraries were quantified with a Qubit 3.0 fluorometer (Life Technologies). 5hmC sequencing was performed on the Next-Seq 500 platform according to pairedend 38-bp high-throughput sequencing. Finally, we used the Agilent 2,100 bioanalyzer for quality control of the 5hmC library and based on the strip size to determine the presence of adaptor dimers (120-150bp).

FastQC (version 0.11.5) was used to check the sequence quality. Additionally, bowtie2 (version 2.2.9) was adopted for aligning raw reads to the human genome (version hg19) (Langmead and Salzberg, 2012), and then filtered with SAM tools (version 1.9) (parameters used: Sam tools view-f 2-F 1548 -q 30) (Li et al., 2009). Subsequently, the paired-end reads were converted into the Bed Graph format, and normalized them to the overall quantity of aligned reads by exploiting bed tools (version 2.19.1) (Langmead, 2010). Meanwhile, with the aid of UCSC BedGraphToBigWig, we also converted the paired-end reads into the BigWig format, so that the Integrated Genomics Viewer-assisted visualization could be achieved. Potential 5hmC enriched regions were identified using MACS2 (version 2.1.1) in each sample (Zhang et al., 2008). Peak regions that appeared in more than 10 samples and that were less than 1,000 bp were combined into one unified catalog for each patient. Genomic regions that tend to show artifact signals, according to ENCODE, were filtered out. The 5hmC enriched regions were generated by intersecting the individual peak call file with the merged peak file. We annotated the 5hmC-enriched region using the CHIP seeker package and used the genes closest to this region for subsequent analysis.

DLBCL patients were randomly divided into training and validation groups in a ratio of 2:1; using train_test_split in scikit-learn (version 0.22.1) package in Python (version 3.6.10), the logistic regression CV (LR) model was chosen to establish warning models (Abraham et al., 2014). In the training cohort, we identified deferentially 5hMc-enriched regions (DhMRs) using DESeq2 package (version 1.30.0) in R (version 3.5.0), with the filtering threshold (p-value <0.05 and |log2FoldChange| ≥ 0.5). To avoid overfitting, 10-fold cross-validation was carried out for 5 rounds in the following manner. Subsequently, this study carried out 100 repeats for further filtering with the use of Scikit-Learn module’s recursive feature elimination (parameters adopted: estimator = LogisticRegressionCV (class weight = “balanced”, cv = 2, max_iter = 1,000), scoring = “accuracy”). Meanwhile, tenfold cross-validation was repeated 100 times in each round, and the final markers observed in at least three rounds were used to build the final warning model in the training cohort. Next, trained LR model was used to predict the treatment outcome for patients in the validation cohort. Receiver operating characteristics (ROC) analysis was used to evaluate model performance. Area under the curve (AUC), best cutoff point, sensitivity, and specificity were computed with sklearn metrics module. A weighted diagnostic score (wd-score) was then calculated as the sum of the gene-wise product of logistic model coefficients and corresponding 5hmC marker value for each individual: wd-score = ∑ k = 1 n β k × g e n e k , where β _k is the coefficient from the final multivariable logistic model for the kth marker gene, and gene _k is the 5hmC level of the kth marker gene. The area under curve (AUC) and 95% CIs were generated to evaluate the model performance.

Differentially hydroxymethylated genes (DhMGs) were annotated using R Package’s ChIPseeker 1.20.0 (Schmitz et al., 2018), and further functional assessments were accomplished on the genes situated nearest to the marker zones. The enrichment analysis of the GO biological process (BP) was completed by the ClueGO (version 2.5.8) and CluePedia (version 1.5.8) plug-in from Cytoscape software (version 3.7.1). Additionally, medium network specificity, Bonferroni adjusted p < 0.01, and enriched gene number >5 were chosen as the criteria for significance. We used the Search Tool for the Retrieval of Interacting Genes (STRING) database (version 10.0, https://string-db.org) to find protein–protein interactions for 5hmC markers. Then, the Cytoscape software was used to construct the PPI network.

For our survival analysis, we utilized the gdc-client (version 1.5.0) to download mRNA HTseq FPKM data from 48 patients with Diffuse Large B-Cell Lymphoma (DLBCL) as part of the TCGA-DLBC dataset (Li T. et al., 2017) from the GDC Data Portal. Concurrently, we manually retrieved curated clinical data, encompassing overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI), from the UCSC Xena platform (Goldman et al., 2020). The survival analysis in this study utilized the Survminer package (version 0.4.6) and Surviva packages (version 2.44–1.1) in R. Forty-eight patients were categorized into either the high-expression group or low-expression group based on cutoff points determined by the maximally selected rank statistics algorithm. Survival analysis for each gene was conducted using Kaplan-Meier curves (Barakat et al., 2019) and the log-rank test (Torbicki et al., 2019). For the survival analysis, p-value <0.05 was considered statistically significant. For gene expression correlation analysis, we used a web tool called TIMER2.0 (Li et al., 2019), which incorporated all TCGA expression data, to explore the mRNA expression relationship between 5hmC markers and other genes of interests in the TCGA DLBC dataset. The correlation analysis was done using Spearman rank correlation.

With the use of GraphPad Prism 8, data were statistically processed as detailed in Table 1. For data showing normal distribution, two-tailed t-tests (either paired or unpaired) were employed. With the purpose of calculating 95% confidence intervals, the percentile method was used. Differences were thought to be of significance with p < 0.05.

Plasma samples from 86 DLBCL patients (Male 46, Female 40), and 90 healthy donors were collected. Clinical data were collected from all samples, and detailed information is listed in Table 1. The mean age of all patients was 54.6 years. Besides, there were 50 BCL2-positive patients, 18 BCL2-negative patients and 18 unknown patients (2.3%). Finally, the mean values of LDH and β2MG in all patients were 364.33 U.L⁻¹ and 2.84 mg.L⁻¹, respectively.

86 DLBCL patients and 90 healthy volunteers were randomly divided into a training cohort (DLBCL = 56, healthy = 60) and validation cohort (DLBCL = 30, healthy = 30) (Figure 1). First, 5hmC-Seal was performed with extracted cfDNA to map the genome-wide 5hmC profiles for all samples. In training cohort, sequencing data showed that 5hmC was mainly enriched at transcription start sites (TSS) and transcription end sites (TTS) (Figure 2A), which was consistent with previous reports (Quail and Joyce, 2013), suggesting that the accumulation of 5hmC is related to transcriptional activity. 5hmC was mostly distributed (75%) on the gene body (GB) in the four groups, and the relative enrichment of 5hmC on GBs was the highest in the DLBCL group (Figure 2B). With the aim of increasing the significance of the findings, our study employed highly stringent peak selection criteria and selected peak base pairs that overlapped in the biological replicates (Figures 2C,D; Supplementary Figure S1A). As a result, during disease development, 5hmC loss is attributed to the intergenic regions and tends to accumulate slowly on GBs. Differential analysis between healthy volunteers and DLBCL patients showed that there were 972 genes with high hydroxymethylation and 160 genes with low hydroxymethylation in DLBCL (Figure 2E, Supplementary Table S1). For instance, DDI1 (Figure 2F) was highly enriched in hydroxymethylation for DLBCL (p = 9.3e-07), and GPR12 (Figure 2G) was highly enriched in hydroxymethylation for healthy volunteers (p = 4.6e-08). Finally, using the default clustering methods, the heat map results showed that Top50 DhMRs in these 1,132 DhMRs could effectively separate DLBCL patients from healthy controls (Figure 2H).

Pathway analysis of 1,132 differentially hydroxymethylated genes (DhMGs) (Supplementary Table S1) in DLBCL patients suggested functional enrichment in certain canonical pathways. The GO biological pathways is mainly concentrated in immune and inflammation related signaling pathways such as myeloid leukocyte activation, CD4-positive alpha-beta T cell activation, peptidyl-serine modification, T cell differentiation involved in immune response and cell activation involved in immune response (Figure 3A). Among these pathways, signaling by alpha-beta T cell differentiation was known to be relevant to tumor growth and apoptosis, which suggested that the DhMRs might be involved in the immunity system (Aucamp et al., 2018; Cioroianu et al., 2019; Roma-Rodrigues et al., 2019). Meanwhile, the hubs of the GO functional interaction networks (Figure 3B) showed that these genes, including Interleukin-6 (IL-6), BCL2 apoptosis regulator (BCL2), Interleukin-18 (IL-18), CD33, CD44, phosphodiesterase 4D (PDE4D) and Mitogen activated protein kinase 8 (MAPK8), participated in the immune and inflammation related signaling pathways.

First, using RFECV based on logistic regression CV estimator, we reduced the number of DhMRs (Supplementary Table S2) in the training cohort, which achieved the best cross-validation score. Then, Using a logistic regression method, we constructed a diagnostic model with these ten markers. Applying the model yielded a sensitivity of 93% and specificity of 98% for DLBCL in the training data set of 56 DLBCL and 60 normal samples (Figure 4A) and a sensitivity of 83% and specificity of 87% in the validation data set of 30 DLBCL and 30 normal samples (Figure 4B). Subsequently, We also demonstrated this model could differentiate DLBCL patients from normal controls both in the training data set (AUC = 0.96) and the validation data set (AUC = 0.94) (Figures 4C,D). Unsupervised hierarchical clustering of these 10 markers was able to distinguish DLBCL from normal controls with high specificity and sensitivity. (Figures 4E,F). Finally, we also calculated the individual AUC for each of the 10 5hmC markers in the training and validation cohorts (Supplementary Figure S2B). Among these, THRAP3 showed the best diagnostic performance, yielding an AUC of 0.711 in the validation cohort (Supplementary Figure S2). These results indicate that plasma cfDNA-based 5hmC markers are a promising diagnostic tool for DLBCL. Meantime, we also applied this model to calculate the wd-score for every single sample, and showed that the wd-score in DLBCL patients was significantly higher than those in healthy volunteers (Supplementary Figure S3A). In addition, we also used the wd-score value to draw the box plot of different stages of the disease and clinical diagnostic indicators. The results showed that each stage and clinical diagnostic indicators were significantly different from the wd-score value of the healthy volunteers, but there was no difference between the groups (Supplementary Figure S3B).

Next, we sought to investigate the potential association of those 10 markers with DLBCL. Since our previous data showed that significantly expressed hydroxymethylation genes were mainly concentrated in immune-related signaling pathways, it was believed that these genes were related to DLBCL immune response. Therefore, we first intersected these 10 genes with immune-related genes to obtain two genes, DENND1A (DENN domain containing 1A) and TSC22D1 (TSC22 domain family member 1) (Supplementary Figure S2A), and then among the two 5hmC-modified marker genes, DENND1A had the higher AUC in the validation cohort (Supplementary Figure S2B). In our study, DENND1A was highly enriched in hydroxymethylation in the DLBCL patients (p = 4.7e-08) (Figure 5A) and its mRNA expression level in the TCGA-DLBC dataset was consistent with the hydroxymethylation level in our data-set (Figure 5B).

In addition, from the PPI network constructed from the STRING database, we identifed several genes linked to DENND1A, including AKT serine/threonine kinase 1 (AKT1), AKT serine/threonine kinase 2 (AKT2), AKT serine/threonine kinase 3 (AKT3), Ras-related protein Rab-35 (RAB35), MAP kinase activating death domain (MADD) and CHM like Rab escort protein (CHML) (Figure 5C). Interestingly, we found that all of these gene expressions (AKT1 (rho = 0.485), AKT2 (rho = 0.523), AKT3 (rho = 0.342), RAB35 (rho = 0.585), CHML (rho = 0.453), MADD (rho = 0.612)) were positively associated with that of DENND1A (Figures 5D–I) (Supplementary Figure S4A). Moreover, from survival analysis results in the TCGA DLBC dataset, we found that the overall survival time (OS, days) of patients with high expression of DENND1A and AKT1 was significantly lower than that of patients with low expression in these two genes (Figures 5J,K).