Similarity of immune-associated markers in COVID-19 and Kawasaki disease: analyses from bioinformatics and machine learning

Abstract

Background: Infection by the SARS-CoV-2 virus can cause coronavirus disease 2019 (COVID-19) and can also exacerbate the symptoms of Kawasaki disease (KD), an acute vasculitis that mostly affects children. This study used bioinformatics and machine learning to examine similarities in the molecular pathogenesis of COVID-19 and KD.

Methods: We first identified disease-associated modules in KD using weighted gene co-expression network analysis. Then, we determined shared differentially expressed genes (DEGs) in training datasets for KD (GSE100154) and COVID-19 (GSE225220), performed functional annotation of these shared DEGs, and used Cytoscape plug-ins (MCODE and Cytohubba) to characterize the protein-protein interaction (PPI) network and identify the hub genes. We performed Least Absolute Shrinkage and Selection Operator(LASSO) regression and receiver operating characteristic (ROC) curve analysis to identify the most robust markers, validated these results by analysis of two other datasets (GSE73461 and GSE18606), and then calculated the correlations of these key genes with immune cells.

Results: This analysis identified 26 shared DEGs in COVID-19 and KD. The results from functional annotation showed that the shared DEGs primarily functioned in immune responses, the formation of neutrophil extracellular traps, and NOD-like receptor signaling pathways. There were three key genes (PGLYRP1, DEFA4, RETN), and they had positive correlations with monocytes, M0 macrophages, and dendritic cells, which function as immune infiltrating cells in KD.

Conclusion: The potential immune-associated biomarkers (PGLYRP1, DEFA4, RETN) along with their shared pathways, hold promise for advancing investigations into the underlying pathogenesis of KD and COVID-19.

Keywords: Bioinformatics analysis; COVID-19; Immune cell infiltration; Kawasaki disease (KD); Weighted gene co-expression network analysis (WGCNA).

Fig. 1

Flow chart

Fig. 2

WGCNA analysis of the GSE100154 dataset (KD vs. control). (a) Network topology for different values of the soft-threshold power, and selection of a β value of 18. (b) Gene clustering dendrogram, with different colors indicating distinct clusters. Thirteen modules had similar co-expression characteristics. (c) Heatmap for the correlation of different modules with disease phenotype. The pink and purple modules had positive correlations with KD, while the green-yellow and turquoise were negatively correlated

Fig. 3

Identification of DEGs and in the GSE100154 dataset and the GSE22520 dataset. (a, b) Data normalization and volcano plot of DEGs in GSE100154(KD vs. control). (c, d) Data normalization and volcano plot of DEGs in GSE225220(COVID-19 vs. control). (adj. P < 0.05,|logFC|>1). (e) Shared genes in comparison of the WGCNA modules (n = 1047) and the DEGs in GSE100154 (n = 117). (f) Shared genes in comparison of the WGCNA modules and DEGs in GSE100154 (n = 56) and the DEGs in GSE225220 (n = 2711). logFC: log₂(fold-change)

Fig. 4

Functional enrichment analysis of COVID-19-related KD hub genes. (a) The top 5 GO terms related to three categories of hub genes (BP, biological process; CC, cellular component; MF, molecular function). (b) KEGG pathways related to the hub genes. Only highly enriched terms and pathways (P < 0.05) were retained

Fig. 5

PPI network analysis of COVID-19-related KD hub genes. (a) PPI network of 26 hub genes. (b) Fifteen genes in the most representative module screened by the ‘MCODE’ plug-in. (c) Common hub genes identified from four algorithms of cytoHubba(MCC, EPC, Degree and MNC). (d) Network of the nine common hub genes. MCC: Maximal Clique Centrality, EPC: Edge Percolated Component, MNC: Maximum Neighborhood Component

Fig. 6

Identification of the key COVID-19-related KD hub genes. (a, b) Screening of DEFA4, PGLYRP1, and RETN by LASSO regression analysis.The optimal lambda (lambda.1se) was 0.1563205 based on a 10-fold cross-validation. (c) ROC curves of DEFA4, PGLYRP1, RETN in the training dataset (GSE100154, KD vs. control). (d) ROC curves of DEFA4, PGLYRP1, and RETN in the validation dataset (GSE73461, KD vs. control). AUC: area under the curve

Fig. 7

Validation of the expression of the key COVID-19-related KD hub genes. (a) Expression of DEFA4, PGLYRP1, and RETN in the GSE100154 dataset (KD vs. control). (b) Expression of DEFA4, PGLYRP1, and RETN in KD patients before and after IVIG treatment in the GSE18606 dataset. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001. logFC: log₂(fold-change)

Fig. 8

Relationship of immune cell infiltration and key COVID-19-related KD hub genes. (a) Proportions of 22 types of immune cells in each sample of the GSE100154 dataset (KD vs. control). (b) Proportions of 22 types of immune cells in each sample of the GSE225220 dataset (COVID-19 vs. control). (c) Infiltration status of immune cells in the two groups of the GSE100154 dataset. (d) Infiltration status of immune cells in the two groups of the GSE225220 dataset. (e) Correlation of DEFA4, PGLYRP1, and RETN with infiltration of different immune cells in the GSE100154 dataset. (f) Correlation of DEFA4, PGLYRP1, and RETN with infiltration of different immune cells in the GSE225220 dataset. *p < 0.05; **p < 0.01; ***p < 0.001