The whole-blood DNA methylation data for the 574 used samples were obtained from the GEO database, with the accession number: 179325 [26]. The samples were divided into three groups according to the COVID-19 severity and included 101 negative controls, 360 mildly infected individuals, and 113 severely infected individuals. The DNA methylation sites in the samples were considered as features, and each sample contained 768,067 methylation sites.

The number of DNA methylation sites involved in the dataset was large, although only a very small number of sites were significantly methylated during COVID-19 development. All sites were screened and ranked using the least absolute shrinkage and selection operator (LASSO) [27], light gradient-boosting machine (LightGBM) [28], and Monte Carlo feature selection (MCFS) [29]. A higher site was ranked corresponding to its higher degree of association with the target variable. These methods are extensively accepted in the life sciences [30, 31, 32, 33].

The use of IFS is well-established in machine learning research [33, 34]. When building a classification model, computational efficiency needs to be considered. For this purpose, we performed the IFS method to determine the optimum number of features required to build the classification model. It converted the list of methylation sites generated by the ranking method into a number of feature subsets. Depending on the set step size, the number of features contained in these subsets grew equally according to the ranking order. These subsets were used to train the downstream classification algorithm, where in reference to the performance of the obtained models, the finalized feature subset was the best feature subset, and the model at this point was the best classification model.

Looking at the sample dataset, a difference existed in the number of samples within the three classes, whereby there were 3.6 times more samples from the mildly infected individuals than from the negative controls. Moreover, directly training the model using an uneven dataset biased the results toward the majority classes. The SMOTE method generated new samples based on known samples. All samples were projected into the high-dimensional space, and one sample was randomly selected for the minority classes. Using Euclidean distance as a metric, SMOTE determined the k-nearest neighbors to that sample in the same class. Moreover, any point on the concatenation of this sample and any of its nearest neighbors can be selected as a new sample. The above process was repeated and the new samples were added into minority classes until the number of samples in each class was balanced.

Based on several feature subsets generated by IFS, this study used DT [35], kNN [36], RF [37], and SVM [38] algorithms to construct classification models. These algorithms have been approved by previous publications [34, 39, 40, 41, 42, 43, 44, 45].

Using the IFS method and above four classification algorithms, we obtained several classification models and evaluated the performance of these models using 10-fold cross-validation [46]. The F1 measure is often used as a key metric to evaluate the performance of classification models, although it does not perform fairly enough in multiclassification problems. The current study required weighting the F1 measure.

(2)$${\mathit{\text{Precision}}}_i=\frac{T{P}_i}{T{P}_i+F{P}_i},$$

(3)$${\mathit{\text{Precision}}}_{\mathit{\text{weighted}}}={\sum }_{i=1}^L{\mathit{\text{Precision}}}_i\times w_i,$$

(4)$${\text{𝑅𝑒𝑐𝑎𝑙𝑙}}_i=\frac{T{P}_i}{T{P}_i+F{N}_i},$$

(5)$${\text{𝑅𝑒𝑐𝑎𝑙𝑙}}_{\mathit{\text{weighted}}}={\sum }_{i=1}^L{\text{𝑅𝑒𝑐𝑎𝑙𝑙}}_i\times w_i,$$

(6)$$\mathit{\text{Weight}}F1=\frac{2\cdot {\mathit{\text{Precision}}}_{\mathit{\text{weighted}}}\cdot {\mathit{\text{Recall}}}_{\mathit{\text{weighted}}}}{{\mathit{\text{Precision}}}_{\mathit{\text{weighted}}}+{\mathit{\text{Recall}}}_{\mathit{\text{weighted}}}},$$

where $TP$ represents true positives, $FP$ represents false positives, $FN$ represents false negatives, $i$ represents the category, $L$ represents the number of classes, and $w_i$ represents the proportion of samples categorized to the overall sample. We also used accuracy (ACC) and Matthew’s correlation coefficient (MCC) [47] as references. A higher value corresponded to better model performance.

Whole-blood DNA methylation data from 574 samples, each containing 768,067 methylation-site features, were analyzed using a machine-learning workflow. LASSO, LightGBM, and MCFS were used to rank all methylation sites and three ranked lists were obtained (Supplementary Table 1). LightGBM discarded the redundant features, meaning only 33,537 features were output in LightGBM for ranking. The lists involved numerous site features, although only a very small number were associated with COVID-19 levels. Thus, for subsequent analyses, only the top 10,000 site features in the list were used.

The first 10,000 methylation sites in the three lists were input into the IFS framework, and the step parameter was set to 10. Each list was transformed into a subset of 1000 sites, and these subsets of methylated sites were used to construct 4 classification models, namely, DT, kNN, RF, and SVM. Supplementary Table 2 describes the classification performance of these models. Based on the performance results, IFS curves were plotted using the number of features as the independent variable to depict the weighted F1 trend (Fig. 2). The IFS curve observation revealed that the SVM model always had the highest performance regardless of the list being used. The best SVM model used the top 3390, 9830, and 8090 methylation sites from the Lasso, LightGBM, and MCFS lists, respectively, to provide weighted F1 values of 0.922, 0.918, and 0.947, respectively. Their ACC and MCC results were also excellent, with ACC of 0.922, 0.916, and 0.946, respectively, and MCC of 0.862, 0.857, and 0.904, respectively (Table 1).

Fig. 2.

IFS curves for evaluating the performance of the three classification algorithms based on the weighted F1. Four classifier models were constructed under each algorithm. (A) IFS curves based on LASSO results. (B) IFS curves based on LightGBM results. (C) IFS curves based on MCFS results. IFS, incremental feature selection.

To extract the most critical sites among the methylation sites used by the best classifier, we observed the IFS curves of the SVM models and set the inflection points of the curves. The inflection points were 150, 80, and 380 for the LASSO, LightGBM, and MCFS lists, respectively. The performance evaluation of the model at the inflection points is displayed in Table 1. Recurring methylation sites tended to play a more important role. We calculated the intersection from the pre-inflection point sites and plotted a Wayne diagram (Fig. 3). We did not obtain methylation sites that appeared in all three lists but obtained 19 features that appeared in two lists. The specific intersection results are shown in Supplementary Table 3. A detailed analysis of these sites is presented later.

Fig. 3.

Venn diagram of the most critical subset of features obtained using LASSO, LightGBM, and MCFS. The overlapping circles indicate methylation sites that are identified as the most critical features by the different ranking algorithms.

The DT model demonstrated a transparent decision process, which helped us obtain specific classification cues. Based on these cues, we summarized three sets of quantitative classification rules, which had implications for our understanding of the biological context of these key methylation sites. Supplementary Table 4 shows the specific classification rules. Each rule contained several parameters that represented the methylation levels of these sites. For the different lists, Fig. 4 shows the number of rules representing the different categories.

Fig. 4.

The number of rules for the three classes as summarized by the DT classification model, constructed using the three ranking methods.

As shown in Fig. 3, using 2 algorithms identified 19 methylation sites, which were generally more essential than other methylation sites. Here, we selected four for detailed analysis and these are listed in Table 2 (Ref. [48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67]).

IGSF6 (cg02481950) is located at chr16:21650858-21663981 and belongs to the immunoglobulin superfamily member 6, which is expressed in dendritic and myeloid cells [48]. Studies have shown that IGSF6 is associated with changes in the ratio of M0 macrophages and $\gamma$-delta T cells, alongside plasma cells and monocytes in atherosclerotic plaques. IGSF6 may be related to immune-related IFN-$\gamma$ and PD-1 signaling pathways [49] and is reportedly a susceptibility gene for inflammatory bowel disease; therefore, may influence disease manifestations [50]. IGSF6 is potentially involved in the atrial fibrillation development mechanism and is associated with the course or maintenance of autoimmune and chronic inflammatory diseases [51]. Moreover, it can act as a potential molecular marker for antigen presentation by DCs before host vaccinations have been identified. After stimulation with human papillomavirus E7 peptide (p11–20), the immune response in immature DCs (iDCs) is upregulated by IGSF6 and other molecules. With prolonged stimulation time, the genes related to the immune response become more significantly upregulated [52]. DCs play important roles in preventing viral infection. Thus, the proportion of plasmacytoid and myeloid DC levels in serum and the IFN level decrease in patients with severe COVID-19. The impairment of DC function and interferon secretion correlates to the COVID-19 severity [53]. Our analysis showed that the methylation level of IGSF6 increased in severe COVID-19 patients. These results suggest that the upregulation of these genes may be a marker of antigen presentation, and IGSF6 may be a potential molecular marker of the extent of the immune response in patients with COVID-19.

CD38 (cg04332373) encodes a type II transmembrane glycoprotein that synthesizes and hydrolyzes calcium-ion-mobilized messengers. The cg04332373 methylation site is located in the shore region of the CpG island. Studies have shown that CD38 levels in white adipose tissue (WAT) and the liver increase with age. Senescent cell signals promote the accumulation of CD38${}^{+}$ cells in WAT [54]. CD38 regulates nicotinamide dinucleotide (NAD${}^{+}$) metabolism and extracellular nucleotide homeostasis, while CD38 inhibition and “NAD${}^{+}$ enhancement” can help in metabolic disorders related to aging, inflammation, and tumor immune hyperplasia [55]. Increased CD38 expression is a consequence of aging [56], which is a major factor associated with the risk of SARS-CoV-2 infection. CD38 plays an important role in viral infections, including AIDS and COVID-19 [57]. Thus, in COVID-19 patients, CD38 mediates thrombosis and bacterial phagocytosis dependent on NAD${}^{+}$. CD38 promotes immune-cell migration into the site of infection through signaling [58, 59], meaning CD38 may aggravate SARS-CoV-2 infection and increase the risk of secondary bacterial infections [60, 61]. Activation of CD38 and decreased NAD${}^{+}$ can be considered features of aging and may be considered regulators of COVID-19 in old age. High inflammation in COVID-19 may lead to CD38 activation, especially causing severe reactions, such as tissue fibrosis and injury in the elderly [57]. Our analysis showed that CD38 methylation levels were lower in patients with severe COVID-19 than in patients with moderate COVID-19. This evidence suggests that SARS-CoV-2 infection may reduce CD38 methylation levels, promote CD38 activation, promote immune-cell migration, and induce NAD${}^{+}$-dependent bacterial phagocytosis, ultimately, leading to local tissue damage and disease risk.

cg03753191 is located at chr13:43566902 and in the Shore region of CpG islands. It corresponds to the EPSTI1 gene and encodes a protein that promotes tumor invasion and metastasis. Studies have reported that EPSTI1 expression levels were significantly upregulated in leukocytes and nasopharyngeal tissues of COVID-19 patients, compared to normal tissues. The downregulation of EPSTI1 also predicted poorer clinical outcomes in patients with COVID-19, such as intensive care unit hospitalization and increased viral load. Therefore, EPSTI1 may play an important role in antiviral immune regulation [62]. The expression profile of EPSTI1 can classify COVID-19 patients into different groups, whereby the younger patient population exhibited a stronger antiviral immune response, higher EPSTI1 expression level, and better clinical treatment effect. The expression level or methylation pattern of EPSTI1 may significantly affect the COVID-19 severity by modulating antiviral immune responses.

DTX3L (cg22930808) encodes Deltex E3 ubiquitin ligase 3L, which plays a role in DNA damage repair and interferon-mediated antiviral response [63, 64, 65, 66]. DTX3L plays a role in the antiviral response by mediating the “LS-48”-linked ubiquitination in the C3 protease of the cerebral myocarditis virus and human rhinovirus, thereby promoting their proteome-mediated degradation [65]. Dtx3l-related pathways include DNA damage and the innate immune system. The activation of the IFN response has been demonstrated to induce the adp nucleation of host proteins, which depend on PARP9 and its binding partner DTX3L. Expression of the large domain of SARS-CoV-2 nonstructural protein 3 (Nsp3) or deletion of PARP9 or DTX3L had no effect on IFN signaling. PARP9/DTX3L-dependent adp nucleation is a downstream effector of the host IFN response, and the SARS-CoV-2 Nsp3 macrodomain can hydrolyze the IFN signaling end product [67]. This study has added evidence to the relationship between DTX3L methylation patterns and disease severity in patients with SARS-CoV-2 infection, thereby providing more reference for the study of COVID-19 progression and pathogenesis.

In addition to identifying essential methylation sites by using various algorithms, we also obtained several classification rules, as listed in Supplementary Table 4. These rules can predict disease severity based on the methylation level of key DNA methylation sites. However, it is impossible to analyze them in detail because a huge number of rules were constructed. Herein, we focused on valuable features because they can identify important DNA methylation sites and demonstrate their important role as epigenetic-susceptibility sites in patients infected with SARS-CoV-2. We collected the scientific findings of other researchers and preliminarily summarized experimental evidence to conduct the discussion. The key methylation sites are listed in Table 3 (Ref. [68, 69, 70, 71]).

TBC1D4 (cg00197681) belongs to TBC1 domain family member 4, and cg00197681 is located at chr13:76056419 and CpG island region. TBC1D4 encodes a rabbit-GTPase activating protein, which is involved in regulating glucose homeostasis through transporting glucose transporter 4 (GLUT4). Human genetic variants may influence SARS-CoV-2 infection and COVID-19 pathology. Researchers have applied machine-learning algorithms to predict phosphorylation-associated missense single nucleotide variants (pSNVs) and found that these variants are associated with the evolution of host defense systems [68]. Proteins with pSNVs include TBC1D4 and TRIM28, which are involved in the regulation of the viral life cycle and host antiviral response. TBC1D4 may also be associated with the dysregulation of glucose homeostasis in SARS-CoV-2 infection. In the present study, we found that TBC1D4 methylation levels were higher in patients with severe COVID-19 than in the normal population. Combined with previous studies, we hypothesized that SARS-CoV-2 infection induced an increase in the methylation levels of TBC1D4, thereby affecting the host antiviral immunity.

The methylation site cg09623286 is located at chr4:154605468 and the island region of CpG islands. Its related gene is Toll-like receptor (TLR) 2, which encodes products belonging to the TLR family. As a cell surface protein, TLR2 can form dimers with other TLRs to recognize pathogen-associated molecular patterns and activate host immune responses. Activation of different TLR pathways leads to the secretion of proinflammatory factors, including IL, tumor necrosis factor-alpha, and interferon, and is involved in SARS-CoV-2 invasion and infection. However, some TLRs, such as TLR2, may play a dual role in COVID-19 infection [69]. SARS-CoV-2 infection can cause severe disease features (e.g., cytokine storm) and organ failure (e.g., testicular and germ-cell damage). Studies have shown that exposure to SARS-CoV-2 envelope proteins causes TLR2 receptor-dependent testicular cytopaplasia [70]. TLR2 mRNA expression levels also significantly increase in patients with moderate to severe COVID-19. Moreover, the mRNA expression levels of CK-MB, ACE2, and neuroproteinase 1 receptors were positively correlated with TLR2 expression in all patients. The mRNA expression of TLR2 was positively correlated with renal biomarkers and cardiac enzymes in severe and moderate patients. These results suggested that it may be related to COVID-19 severity [71]. The present study showed that patients with moderate COVID-19 had higher methylation levels than normal controls; however, this was not observed in severe patients. Therefore, the TLR2 methylation level may serve as an indicator of COVID-19 severity, although the specific mechanism needs further verification.