Integrated transcriptomics of multiple sclerosis peripheral blood mononuclear cells explored potential biomarkers for the disease

Abstract

Background: Despite their importance, blood RNAs have not been comprehensively studied as potential diagnostic markers for multiple sclerosis (MS). Herein, by the integration of GSE21942 and GSE203241 microarray profiles of peripheral blood mononuclear cells, this study explored potential biomarkers for the disease.

Methods: After identification of differentially expressed genes (DEGs), functional enrichment analyses were performed, and PPI and miRNA-mRNA regulatory networks were constructed. After implementing weighted gene co-expression network analysis (WGCNA) and discovering MS-specific modules, the converging results of differential expression analysis and WGCNA were subjected to machine learning methods. Lastly, the diagnostic performance of the prominent genes was evaluated by receiver operating characteristic (ROC) analysis.

Results: COPG1, RPN1, and KDM3B were initially highlighted as potential biomarkers based on their acceptable diagnostic efficacy in the integrated data, as well as in both GSE141804 and GSE146383 datasets as external validation sets. However, given that they were downregulated in the integrated data while they were upregulated in the validation sets, they could not be considered as potential biomarkers for the disease. In addition to this inconsistency, evaluating their diagnostic performance in other external datasets (GSE247181, GSE59085, and GSE17393) did not reveal their diagnostic efficacy.

Conclusions: This study could not unveil promising blood biomarkers for MS, possibly due to a small sample size and unaccounted confounding factors. Considering PBMCs and blood specimens as valuable sources for the identification of biomarkers, further transcriptomic analyses are needed to discover potential biomarkers for the disease.

Keywords: Biomarkers; Multiple sclerosis; PBMCs.

Fig. 1

Differential expression and functional enrichment analyses. (a) Datasets before batch-effect correction. (b) Datasets after batch-effect correction. (c) Samples are clustered according to the disease state after batch-effect correction. (d) A total of 106 DEGs were identified (|logFC| > 0.5 and adjusted p-value <0.05). (e) Gene ontology classifications and KEGG pathway enrichment analyses.

Fig. 2

Networks of the DEGs in MS. (a) PPI network. The intensity of the border color reflects the logFC values. (b) The miRNA-mRNA network illustrates the key regulators of the PPI network. Blue label reflects downregulation, while red label represents upregulation. The size of the nodes positively correlates with their number of PPIs in the networks.

Fig. 3

WGCNA detection of gene modules correlating with the disease state. (a) Clustering of samples after removing outliers. (b) A soft-thresholding power of 11 was selected as the optimal power. (c) Correlation of module eigengenes with disease state and their corresponding p-values (∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001). (d) The pie chart illustrates the proportion of each module. The most prominent KEGG pathways are presented below each module name based on the combined score in the Enrichr database. (e) The overlap of genes between DEGs and WGCNA analyses. Purple indicates positively correlated modules, whereas green indicates negatively correlated modules.

Fig. 4

In-depth analysis aiming at the identification of hub genes in module seven. (a) High-confidence (0.7) PPIs of hub genes for the co-expressed genes in module 7. (b) Gene significance against module membership for this module. The 23 central genes in the PPI network are colored. (c) Gene ontology analysis of the 23 central genes in the PPI network. (d) KEGG pathway analysis of the 23 central genes in the PPI network. Pathways with the subcategory of neurodegenerative disorders are represented by the green color.

Fig. 5

Employing machine learning methods for the identification of potential biomarkers. (a) Error against the number of trees used in the RF classifier. (b) Importance of genes based on the mean decrease in accuracy. (c) Importance of genes based on the mean decrease Gini. (d and e) Coefficient paths and binomial deviance plotted against the logarithm of lambda used in LASSO regression, respectively. (f) Genes with a non-zero coefficient were ordered according to the absolute value of the coefficient.

Fig. 6

ROC analysis evaluating the diagnostic efficacy of (a) COPG1, (b) RPN1, (c) KDM3B, and (d) ADAM28 across the datasets. Colors represent the diagnostic performance of the gene of interest across different datasets.

Fig. 7

Expression level of COPG1, RPN1, and KDM3B in (a) GSE21942, (b) GSE203241, (c) GSE141804, and (d) GSE146383. Blue and red boxes demonstrate statistically significant downregulation and upregulation, respectively, whereas the alterations of genes with black boxes were not statistically significant.

Fig. 8

Cell-type enrichment of DEGs in MS. Plots showing enriched cell types based on (a) the CellMarker 2024 and (b) the PanglaoDB Augmented 2021 databases.