Identification and evaluation of candidate COVID-19 critical genes and medicinal drugs related to plasma cells

Abstract

The ongoing COVID-19 pandemic, caused by the SARS-CoV-2 virus, represents one of the most significant global health crises in recent history. Despite extensive research into the immune mechanisms and therapeutic options for COVID-19, there remains a paucity of studies focusing on plasma cells.

In this study, we utilized the DESeq2 package to identify differentially expressed genes (DEGs) between COVID-19 patients and controls using datasets GSE157103 and GSE152641. We employed the xCell algorithm to perform immune infiltration analyses, revealing notably elevated levels of plasma cells in COVID-19 patients compared to healthy individuals. Subsequently, we applied the Weighted Gene Co-expression Network Analysis (WGCNA) algorithm to identify COVID-19 related plasma cell module genes. Further, positive cluster biomarker genes for plasma cells were extracted from single-cell RNA sequencing data (GSE171524), leading to the identification of 122 shared genes implicated in critical biological processes such as cell cycle regulation and viral infection pathways.

We constructed a robust protein-protein interaction (PPI) network comprising 89 genes using Cytoscape, and identified 20 hub genes through cytoHubba. These genes were validated in external datasets (GSE152418 and GSE179627). Additionally, we identified three potential small molecules (GSK-1070916, BRD-K89997465, and idarubicin) that target key hub genes in the network, suggesting a novel therapeutic approach. These compounds were characterized by their ability to down-regulate AURKB, KIF11, and TOP2A effectively, as evidenced by their low free binding energies determined through computational analyses using cMAP and AutoDock.

This study marks the first comprehensive exploration of plasma cells’ role in COVID-19, offering new insights and potential therapeutic targets. It underscores the importance of a systematic approach to understanding and treating COVID-19, expanding the current body of knowledge and providing a foundation for future research.

Introduction

Infectious disease COVID-19, primarily attributed to the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), continues its relentless spread globally [1, 2]. The infection typically manifests as severe flu-like symptoms, progressing to acute respiratory distress syndrome (ARDS), pneumonia, renal failure, and, in severe cases, death [3,4,5,6]. Epidemiological evidence indicates fever, cough, and dyspnea as the predominant symptoms among affected individuals [7]. As of May 12^th, 2024, the World Health Organization (WHO) reported 775,481,326 confirmed cases and 7,049,376 fatalities due to COVID-19 (https://data.who.int/dashboards/covid19/). This underscores the necessity of comprehending the disease’s characteristics and the role of pharmacological interventions in its management.

Viral invasion, immune system failure, microcirculation dysfunction, endothelial cell damage, and dysregulation of the renin-angiotensin-aldosterone system are recognized as pathophysiological hallmarks of SARS-CoV-2 infection [8, 9]. Since the emergence of COVID-19, extensive studies have focused on elucidating the molecular mechanisms and host-pathogen interactions pertinent to COVID-19 pathogenesis by analyzing biomolecular profiles in patients [10,11,12,13]. The spike protein of SARS-CoV-2 or SARS-CoV, for example, targets the ACE2 membrane protein, which possesses an enzymatic domain on human cell surfaces, facilitating viral entry [14]. Chua et al. observed that epithelial cells in COVID-19 cases exhibit a threefold increase in ACE2 expression, linked to immune cell interferon signaling [15]. Similarly, Winkler et al. reported significant inflammation, immune cell infiltration, and impaired respiratory function in infected human ACE2 transgenic mice [16]. In COVID-19 patients, a notable decrease in peripheral blood lymphocytes correlates with disease severity [17]. Improvement in lymphocyte count during treatment indicates a favorable prognosis, while a continual decrease suggests a negative outcome [18,19,20,21,22,23]. Severely affected patients exhibit a low T cell count but high inflammatory cytokine levels, primarily secreted by monocytes, macrophages, or other innate immune cells, leading to adverse patient prognoses [24, 25]. Despite the pivotal role of the immune system in COVID-19 progression, the transcriptomic alterations and key genes affecting plasma cells in the disease remain underexplored.

Weighted gene co-expression network analysis (WGCNA) is a system biology method used to analyze gene expression profiles across multiple samples [24, 25]. It clusters genes with similar expression patterns and assesses correlations between modules and specific traits [24, 25], and is widely used in studying conditions such as acute lung damage, chronic obstructive pulmonary disease, and heart failure [26,27,28]. Highly connected gene clusters, potential biomarkers, and therapeutic targets for a variety of illnesses, including viruses, can all be found using it [25, 29, 30]. Using the WGCNA algorithm on high-throughput data, we identified genes with similar expression changes. WGCNA’s criteria, rooted in biological significance, make its results more reliable than traditional methods [24, 25]. These findings enable further research, including trait associations, metabolic pathway modeling, and gene interaction network creation [26, 28, 31]. Single-cell RNA sequencing (scRNA-seq) has greatly advanced our understanding of COVID-19’s etiology and immunological characteristics [32]. It profiles transcriptomes from individual cells within a sample, enabling detailed gene expression analysis [32]. scRNA-seq has been instrumental in studying SARS-CoV-2-infected cells, compositional changes in cell subpopulations, and profiling immune signaling pathways and transcriptomes during disease progression [33, 34]. WGCNA and high-throughput technologies together offer a fresh chance to comprehend the molecular causes of COVID-19 and other infectious disorders [35, 36]. This evidence suggests a significant potential to identify crucial COVID-19 genes through integrated bioinformatics analysis.

In the face of such a highly infectious respiratory condition, developing anti-COVID-19 medications is critical due to the absence of specific treatments for COVID-19. A range of existing medications can mitigate COVID-19 symptoms. The WHO recommends nirmatrelvir and ritonavir for high-risk patients with mild to moderate COVID-19, at risk of hospitalization [37]. Nirmatrelvir inhibits SARS-CoV-2’s 3Cl protease, essential for viral replication [38,39,40], and is combined with ritonavir, a CYP3A inhibitor, to sustain effective drug concentrations [40,41,42]. However, side effects, drug interactions, and cost considerations are important. Common side effects include disturbances to the nervous and cardiovascular systems, gastrointestinal issues, and drug interactions, especially in patients with chronic conditions [43,44,45,46,47]. Hence, research into refining therapeutic strategies to alleviate symptoms and enhance effectiveness is highly promising.

Plasma cells, derived from B lymphocytes, can secrete immunoglobulins or antibodies [48, 49]. Antibodies produced against SARS-CoV-2 can neutralize the virus, reducing its pathogenicity [50]. Because plasma cells produce antibodies that specifically target the SARS-CoV-2 virus; they are essential to the immune response to COVID-19. Convalescent plasma has historically been used in treating viral diseases like Ebola and various fevers, including SARS and influenza strains [51]. The successful application of plasma in pulmonary diseases offers feasibility for COVID-19 treatment [51, 52]. Studies pertaining to COVID-19 and plasma cells have been published. Following COVID-19 or immunization against SARS-CoV-2, SARS-CoV-2 specific antibody-secreting plasma cells promoting specific normal immunity have been found in the human bone marrow. Following basic mRNA vaccination, Schulz et al. demonstrated the production of phenotypically varied SARS-CoV-2-specific plasma cells in the human bone marrow, including the development of memory characteristics [18, 53]. They also came to the conclusion that the basic mRNA vaccine against SARS-CoV-2 induces long-lasting humoral protection [53]. Long-lived bone marrow plasma cells are a reliable and necessary source of antibodies that provide protection. Turner reported that aspirates from 11 healthy people without a history of SARS-CoV-2 infection did not contain S-specific bone marrow plasma cells, indicating their possible membership in a stable compartment [54]. To sum up, in order to effectively mount an antibody response against COVID-19, plasma cells are necessary. Their function changes according to the disease’s severity, and their continued existence is essential for maintaining immunity against both vaccination and spontaneous infection.

Large-scale cytokine and chemokine release causes excessive inflammation and immune cell mobilization, which is what causes the pathogenic phenomena known as “cytokine storm” [55]. Patients with COVID-19 may have a poor prognosis and a persistent infection as a result of plasma cell malfunction [56,57,58]. As a result, the discovery of potential biomarkers linked to plasma cells during COVID-19 infection may help clarify the infection’s progression, enhance our knowledge of the molecular underpinnings of SARS-Cov-2 infection, and enhance COVID-19 diagnosis and therapy.

Herein, the purpose of this study is further to explore the role of plasma cells in COVID-19, providing new insights and novel therapeutic targets, and to highlight the significance of plasma cells in COVID-19 with regard to prognosis and treatment. It emphasizes how crucial it is to approach comprehending and treating COVID-19 methodically, adding to the body of information already in existence and laying the groundwork for further research. Therefore, we developed a workflow for identifying and evaluating candidate COVID-19 critical genes and medicinal drugs related to plasma cells in Fig. 1. In this study, we identified 122 shared genes among differentially expressed genes (DEGs), co-expressed genes, and positive cluster biomarker genes. Functional enrichment analysis annotated these genes to processes like nuclear division, organelle fission, DNA helicase and ATP hydrolysis activities, cell cycle, p53 signaling pathway, human T-cell leukemia virus 1 infection, and protein processing in the endoplasmic reticulum. We constructed a protein-protein interaction (PPI) network of 89 genes and identified 20 hub genes using the cytoHubba plugin. Finally, three potential small compounds (GSK-1070916, BRD-K89997465, and idarubicin) were identified that down-regulate key genes (AURKB, KIF11, and TOP2A) with the lowest free binding energy, using the connectivity map (cMAP) database and AutoDock modeling. Our study offers a novel perspective for identifying and evaluating plasma cell-related candidate targets and their drugs for COVID-19 patients.

Materials and methods

Data source and pre-operation

Our study compiled five independent datasets (Table 1), categorized into discovery and validation sets. The discovery datasets comprised GSE157103, GSE152641, and GSE171524, while GSE152418 and GSE179627 served as validation datasets.

GSE157103 encompasses 126 RNA-seq samples, consisting of 100 from COVID-19 patients and 26 from non-COVID-19 hospitalized individuals [59]. GSE152641, another publicly available RNA-seq dataset, includes 86 samples: 62 from patients and 24 from healthy controls [60]. Additionally, we incorporated a single-nucleus RNA-seq (snRNA-seq) dataset, GSE171524, comprising 20 samples from 19 COVID-19 decedents and 7 control donors [61].

For external validation, GSE152418 and GSE179627 were utilized. GSE152418 contains RNA-seq data of 17 COVID-19 patients and 17 healthy controls, further classified by severity status (17 healthy, 1 convalescent, 4 moderate, 8 severe, and 4 ICU) [62]. GSE179627 includes 70 RNA-seq samples from COVID-19 patients, differentiated by five immune response stages (22 uninfected, 9 asymptomatic, 13 symptomatic, 14 recovering, and 12 re-detectable positive patients) [35].

For bulk RNA-seq datasets, raw read counts and sample labels were retrieved and prepared for analysis. Gene annotation was performed using the bitr() function from the clusterProfiler (version 4.4.4) R package, converting EntrezIDs to gene symbols [63, 64]. Similarly, gene names, cell names, sparse matrices, and metadata were extracted from snRNA-seq data for subsequent analyses.

Fig. 1

The flowchart figure for our overall methodology. The workflow of this study is mainly divided into the following two steps: Step (1) Database-based screening for hub genes associated with COVID-19 symptoms. Step (2) Verify the reliability of hub genes based on external datasets and further analyze docking results

Table 1 Summary of the datasets used in this study

Identification of DEGs

We combined the bulk RNA-seq datasets GSE157103 and GSE152641 as our discovery dataset. This amalgamation resulted in a raw read count matrix comprising 19,184 genes across 212 samples. FactoMineR (version:2.11) and factoextra (version: 1.0.7) R packages were used for performing PCA analysis and visualization. The PCA() function in FactoMineR R package was used to perform PCA analysis. All features (genes) were included in the PCA analysis. The fviz_pca_ind() function in the factoextra R package was applied for visualization. The top two principal components were selected for visualization. Given the disparate origins of the datasets, which may affect gene expression levels, we incorporated batch effect correction in our analysis (Supplementary Fig. 1). The detailed correction method is that the sample label (COVID-19 or non-COVID-19) was designated as the primary variable factor, and the batch effect was considered an additional variable factor by utilizing the “design” parameter in DESeq2 (version 1.36.0) R package for the differential gene expression (DEG) analysis [65].

The identification of DEGs is crucial in analyzing RNA-seq count data. DESeq2 (version 1.36.0), a prevalent R package, facilitates the quantification and statistical inference of RNA-seq data by constructing a negative binomial generalized linear model [65]. For count data transformation analysis, given the sample size exceeding 100, we opted for the vst() function over the rlog() function. The resulting normalized expression matrix was used as the input for heatmap generation.

Significant up-regulation was attributed to genes with an adjusted P-value less than 0.05 and a log2 fold change (log2FC) greater than 1. Conversely, significant down-regulation was assigned to genes with an adjusted P-value less than 0.05 and log2FC less than − 1. We utilized the ggplot2 (version 3.6.6) R package in RStudio for volcano plot visualization and the pheatmap (version 1.0.12) R package for depicting the top 20 up- and down-regulated genes, sorted by fold change magnitude.

Cell type enrichment analysis

Human tissue comprises various cell types. Accurately estimating the proportion of these cell types, especially immune cell types, is critical for elucidating disease mechanisms and suggesting potential treatments in the realm of precision medicine [66, 67]. The xCell algorithm is frequently used for cell type enrichment analysis from gene expression profiles [68]. It utilizes 1,822 human cell-type transcriptome datasets to generate 64 gene signatures, which are categorized into five groups: lymphoid, stem cells, myeloid, stromal cells, and others [68].

For this analysis, we converted raw read counts to fragments per kilobase million (FPKM) by normalizing for sequence depth and gene length, adhering to the xCell requirements. We employed the xCellAnalysis() function from the xCell package (version 1.1.0) on our RNA-seq dataset, yielding an xCell score matrix for 64 cell types across 212 samples.

Subsequently, we conducted comparative analyses to discern differences in xCell scores between COVID-19 patients and healthy controls. To visualize and statistically analyze these differences, we employed the ggpubr (version 0.4.0) and rstatix (version 0.7.0) R packages for generating boxplots and conducting statistical tests. We deemed cell type enrichments as significantly different if the P-value between COVID-19 and healthy controls was less than 0.05.

Weighted gene co-expression network analysis (WGCNA)

The module associated with COVID-19 and plasma cells was identified using WGCNA [24, 25]. WGCNA is a systematic biological approach designed to identify clusters of genes with similar expression patterns, often indicative of shared roles in biological processes. The WGCNA computational framework includes building a gene co-expression network, identifying modules, associating modules with clinical traits, exploring module interactions, and identifying key genes within the modules of interest.

Using WGCNA (version 1.71) in R, we sought to discover modules related to COVID-19 and plasma cells. Our analysis incorporated 162 COVID-19 and 50 healthy control samples. We categorized all samples into low and high plasma cell cohorts based on the median value of the xCell score. We input an expression matrix consisting of 19,060 genes across 212 samples into WGCNA. Only genes within the top 75% of median absolute deviation (MAD) and with a minimum MAD of 0.01 were considered. After this gene-level filtering, we processed 14,295 genes for 212 samples. Outlier samples were identified using the hclust() method. We only contained 209 samples from 212 samples by removing three outliers for the further analysis.

Many data in the relevant references prove that the network already conforms to the distribution of a scale-free network when the signed R² is greater than 0.85 [24, 25]. Expecting to consider the recommended soft power by the script “sft$powerEstimate”, we added two criteria: (1) Draw the relationship between the soft threshold and the scale-free topology model fit, and take the soft threshold that makes scale-free topology model fit, signed R² > 0.85 as the soft threshold. (2) Draw the relationship between the soft threshold and the mean connectivity, and take the soft threshold that makes mean connectivity < = 500 as the soft threshold. Based on the above discussion, 6 was set to the soft power. The minModuleSize parameter represents the minimum number of genes in a module. The smaller the value, the smaller the module will be retained. The mergeCutHeight indicates the distance of merging similar modules. The smaller the value, the less likely they will be merged, and the more modules are retained. The minModelSize was set to 100, and the mergeCutHeight was set to 0.25 based on the defaulted setting consistent with the WGCNA tutorial (https://rpubs.com/natmurad/WGCNA) and related published studies [69,70,71]. A module was deemed related to COVID-19 and plasma cells if it had a correlation coefficient greater than 0.3 in both the COVID-19 and high-plasma cell-enriched groups and a P-value less than 0.05. Furthermore, genes within this module were defined as module genes related to both COVID-19 and plasma cells. The relationship between module membership and gene significance was visualized using the verboseScatterplot() function.

Single-cell RNA sequencing (scRNA-seq) analysis

We analyzed positive cluster biomarker genes of plasma cells in COVID-19 samples as they are indicative of key roles in the immune microenvironment. Using Seurat (version 4.2.0), we initiated objects by introducing sparse matrices, gene names, and cell names [72,73,74,75]. Notably, the sparse matrix was first converted to a dgCMatrix format as per Seurat’s default requirement. A total of 34,546 genes across 116,313 cells were integrated into the Seurat object. In this study, we utilized the cleaned data download from the link: https://singlecell.broadinstitute.org/single_cell/study/SCP1219/columbia-university-nyp-covid-19-lung-atlas#study-download. Based on the description from Melms et al., the data was cleaned [61].

For uniformity in cell type annotation, we accessed and downloaded processed data, including metadata, from their portal website (https://singlecell.broadinstitute.org/single_cell/study/SCP1219/columbia-university-nyp-covid-19-lung-atlas). The “cell_type_intermediate” column was used to assign cell type identities. Subsequently, UMAP was employed for non-linear dimension reduction and two-dimensional visualization.

A stacked barplot visualized the distribution of 19 cell types across 26 COVID-19 samples. To compare cell ratios between COVID-19 and healthy control groups across each cell type, an inter-group t-test was conducted, setting a significance cutoff at P-value < 0.05. Positive cluster biomarker genes in plasma cells were identified using the FindMarkers() method with the parameters: only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, and test.use = “wilcox”. The only.pos parameter was set to TRUE because we only focus on positive markers. At the same time, 0.25 is a common parameter that can be used at the min.pct and logfc.threshold based on the associated published research [76,77,78]. The “wilcox” was set to the test method for identifying differentially expressed genes between two groups of cells using the Wilcoxon Rank Sum test by default.

Functional annotation of significant genes

Functional annotation analysis is a critical tool for unraveling the biological significance of genes [79,80,81]. It aids in understanding the functionality of gene groups, discerning gene-disease relationships, and inferring gene products and their roles in biological processes [82,83,84]. For this purpose, various gene annotation databases have been developed to annotate and store gene functions across different biological paradigms.

In this study, our gene lists were derived from differentially expressed genes (DEGs) and overlapping shared genes. We utilized the clusterProfiler package (version 4.4.4) in R to conduct functional enrichment analysis of both downregulated and upregulated genes between COVID-19 and non-COVID-19 patients. Initially, the bitr() function was employed for gene name transformation from the official gene symbol to Entrez ID. Subsequently, the top 30 Gene Ontology (GO) biological process (BP) category terms, sorted by significance level in ascending order, were visualized.

We integrated GO cellular component (CC), biological process (BP), molecular function (MF), and KEGG pathway analysis to deduce functional pathways of shared genes implicated in COVID-19, with a particular focus on plasma cells. Shared genes were defined as those appearing at least twice among DEGs, module genes, and positive cluster biomarker genes. The clusterProfiler package (version 4.4.4) facilitated the functional enrichment analysis of these shared genes [63, 64]. The Database for Annotation, Visualization, and Integrated Discovery (DAVID) was employed to provide enriched KEGG pathways of shared genes [85, 86]. The top 10 GO terms in the BP, CC, and MF categories were individually visualized, and KEGG pathways with a P-value < 0.05 were highlighted. This comprehensive analysis enables a deeper understanding of the functional implications of genes in the context of COVID-19 and plasma cells.

Construction of the PPI network

Protein-protein interactions (PPIs) are fundamental elements of cellular biochemical networks and are crucial for regulating cellular functions and signaling pathways [87,88,89]. Understanding the intricate relationships among multiple proteins via systematic analysis facilitates insights into biological signal response mechanisms and energy substance metabolism in specific physiological conditions, including disease states, as well as elucidates the functional interplay among proteins [87, 90].

To analyze PPIs, shared genes identified in our study were input into the STRING database (version 11.5, https://string-db.org/) [91]. In total, there are 19,566 proteins (nodes) in H. sapiens (https://version-11-5.string-db.org/organism/9606). We focused on retaining PPIs with a false discovery rate (FDR) of < 0.05 and a combined score of > 0.8 to ensure the relevance and reliability of the interactions based on the recommended results [92, 93]. The resulting PPI network was then visualized using Cytoscape software (version 3.9.0) [94], which is adept at depicting complex interaction networks. Within the network, node border thickness corresponds to the degree of connectivity, with thicker borders indicating higher degrees of interaction. Similarly, the edge thickness represents the interaction score, with thicker lines suggesting stronger associations. This graphical representation allows for a more intuitive understanding of the extensive and nuanced interplay between proteins in the context of COVID-19 and plasma cells, thereby providing a foundational framework for further biological and therapeutic explorations.

Identification of sub-networks and hub genes

In the construction of a protein-protein interaction (PPI) network using the STRING database, often hundreds of PPI pairs are identified. To manage this complexity and distill the network down to its most influential components, it is necessary to employ algorithms that can identify key sub-networks and central genes, known as hub genes, based on network topology.

CytoHubba, a widely acclaimed plugin for Cytoscape, is particularly adept at identifying pivotal nodes within biological networks and is thus a valuable tool for experimental biologists seeking to uncover critical regulatory networks and potential drug targets [95]. Among the various methods available in CytoHubba, the Matthews correlation coefficient (MCC) method is recommended for its superior predictive accuracy compared to other topological analysis algorithms [95].

Utilizing the MCC method, a scoring system is applied to each gene within the PPI network, facilitating the identification of the most centrally involved genes, or hub genes. The resultant sub-network, comprising the top 20 hub genes, is visually distinguished in the network diagram, with the genes colored according to their rank. This prioritization not only highlights the genes most central to the network’s structure and function but also provides a focused list of candidates for subsequent biological or therapeutic investigation. This analytical approach is particularly crucial in understanding the complex biological interactions at play in diseases such as COVID-19 and can significantly aid in identifying potential targets for further research and drug development.

Validation of hub genes

Validation of hub genes is crucial for substantiating their relevance and reliability in differentiating COVID-19 patients from non-COVID-19 individuals and across various severity stages of the disease. Initially, the expression patterns of the 20 hub genes identified in the previous step were extracted from the discovery dataset to establish a baseline comparison between COVID-19 and non-COVID-19 groups.

To ensure the robustness and applicability of these hub genes, we utilized an independent external dataset, GSE152418, which includes detailed severity clinical information. Samples within this dataset were categorized into five distinct groups: convalescent, healthy, ICU (Intensive Care Unit), moderate, and severe. We computed and compared the expression levels of the 20 hub genes across these varied groups to observe any consistent trends or significant differences that could correlate with disease severity.

Further validation was performed using another independent external dataset, GSE179627, which provided a range of immune response classifications, including asymptomatic, recovering, re-detectable positive patients, uninfected, and symptomatic groups. This diversity in immune response categories allowed for a comprehensive evaluation of the hub genes as potential biomarkers across different stages and responses to COVID-19 infection.

To quantify the predictive power of the hub genes, the R package pROC (version 1.18.4) was employed to calculate the area under the receiver operating characteristic curve (AUC) and the 95% confidence intervals (CIs). Genes with AUC values exceeding 0.65 were considered to have satisfactory discriminatory ability and were retained for further analysis based on the previous studies [96, 97]. This step is vital to ensuring that the hub genes not only show differential expression between COVID-19 and non-COVID-19 cases but also possess the ability to distinguish between various disease severities and immune responses, thereby affirming their potential utility in clinical and research settings.

Screening of candidate small molecular compounds

We incorporated the 20 hub genes into the Connectivity Map (cMAP) database [98] to explore potential small molecular compounds that may influence these genes’ expression. The cMAP database, hosting over 1.5 million expression profiles from more than 5,000 small molecule compounds and approximately 3,000 genetic perturbations across diverse cell types, provided a normalized connectivity score reflecting each compound’s effect on specific cell types and perturbations [98]. These proteins, encoded by the hub genes, were considered potential targets for novel therapeutic drugs. cMAP facilitated the prediction of small molecule medicines that might down-regulate these proteins [98]. We particularly emphasized the top ten promising small molecule substances or medications that demonstrated potential in reducing the expression of these critical genes.

In order to make the research more in-depth, we sourced structural and descriptive information on these therapeutic agents from the PubChem database [99], a leading repository for chemical data. Proteins encoded by genes with Area Under Curve (AUC) values greater than 0.75 were categorized as significant proteins associated with COVID-19, aiding the prediction of immune responses [99].

Molecular modeling simulation

Molecular docking employs a theoretical modeling approach to design drugs by analyzing the properties of receptors and how they interact with drug molecules, aiming to predict their binding affinities. Autodock Vina, version 1.1.2, is a widely recognized application in this field [100]. We collected the crystal structures of three proteins, AURKB (4AF3) [101], KIF11 (3K5E) [102], and TOP2A (4FM9), from the Protein Data Bank (PDB) in PDB format for our study [103, 104]. Correspondingly, the chemical structures for the potential inhibitory small molecules—GSK-1070916, BRD-K89997465 (chlorpromazine), and idarubicin—were sourced from the PubChem database in SDF format [99]. PyMOL, version 1.7.X, was utilized for the removal of extraneous molecules and data format transformation, while Autodock tools, version 1.5.7, were used for hydrogen addition and charge editing. Finally, LigPlus, version 2.2.8, facilitated the detailed analysis of interactions between the proteins and ligands [105, 106]. This comprehensive simulation provided a deeper insight into the potential efficacy of these compounds as COVID-19 therapeutics.

Molecular dynamics simulation

We used GROMACS (version: 2022.4) software to perform an all-atom molecular dynamics simulation analysis on the protein-ligand complex obtained by molecular docking. The protein part was parameterized using the Amber14SB force field, while the topology file of the small molecule part was generated by ACPYPE and Antechamber programs. We selected a cubic solvation box and set the shortest distance from the edge of the system to the complex to 1 nm. We selected the TIP3P water model and added appropriate amounts of sodium and chloride ions to balance the overall charge of the system.

Furthermore, the system used the steepest descent method for energy minimization, and the temperature was regulated by the NVT ensemble, and the pressure was controlled by the NPT ensemble to ensure that the temperature of the entire system was stable at 300 K and the pressure was maintained at 101.325 kPa. For each balanced system, we performed a 100ns molecular dynamics simulation at 300 K and obtained a total of 10,000 frames of simulation trajectories. We used the trajectory files obtained from the simulation to conduct an in-depth analysis of key parameters such as the root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), and the total number of hydrogen bonds between the protein and the ligand. The free energy landscape (FEL) reflects the relationship between the Gibbs free energy of the protein and the RMSD and Rg. (The presence of two low energy points (i.e., local minima) in the free energy landscape (FEL) usually means that the system has two stable conformational states.)

In addition, we take a 2 ns trajectory (i.e., 200 frames) during the equilibrium period and calculate the binding free energy using the MM-PBSA method of gmx_mmpbsa. By definition, in solution, we can write the binding free energy as

$$\:\varDelta\:{G}_{bind}={G}_{complex}-\left({G}_{free-protein}+{G}_{free-ligand}\right)$$

Then, the solvation free energy can be further decomposed into polar and nonpolar parts.

$$\:{G}_{bind}={E}_{gas}-{TS}_{gas}+{G}_{solvation}$$

In which, the solvation free energy can be further decomposed into polar and nonpolar parts.

$$\:{G}_{solvation}={G}_{polar}+{G}_{nonpolar}$$

In MM-PBSA method, the energy and entropy contributions of the gas phase are calculated according to the MM method as follows

$$\:{E}_{gas}={E}_{MM}={E}_{bond}+{E}_{angle}+{E}_{dihedral}+{E}_{vdw}+{E}_{coulomb}$$

$$\:{S}_{gas}={S}_{MM}$$

Among them, $\:{E}_{bond}$, $\:{E}_{angle}$, and $\:{E}_{dihedral}$ correspond to bond, bond angle, and dihedral interactions, respectively, and $\:{E}_{vdw}$ and$\:\:{E}_{coulomb}$ represent van der Waals forces and Coulomb electrostatic interactions, respectively.

In the MM-PBSA method, the solvation energy contains two parts: acute solvation energy and nonpolar solvation energy. The acute solvation energy comes from the electrostatic interaction between solute and solvent molecules and is calculated using an implicit solvent model, in which the solvent is considered as a continuous medium and the corresponding Poisson-Boltzmann equation is linearized and numerically solved:

$$\:{G}_{polar}={G}_{GB}$$

The nonpolar solvation free energy can be calculated based on the empirical surface area method, so it is also called the surface solvation energy. The calculation requires knowledge of the solvent accessible surface area (A) of the molecule and uses two empirical parameters $\:\gamma\:$ and b:

$$\:{G}_{nonpolar}={G}_{surface}=\gamma\:A+b$$

Combining the above items, the free energy equation of MM-PBSA is obtained

$$\:{G}_{bind}={E}_{MM}-{TS}_{MM}+{G}_{GB}+{G}_{surface}$$

This analysis process is similar to published analysis [107, 108].

Other statistical analysis

In this study, bioinformatics and statistical analyses were conducted using R and RStudio, version 4.2.1. We employed the VennDiagram package, version 1.7.3, to calculate and illustrate the overlap among DEGs, module genes, and positive plasma cell biomarker genes. The Student’s t-test was utilized to ascertain differences in cell ratios between COVID-19 patients and healthy controls across various cell types. Additionally, the Kruskal-Wallis test was applied to assess variations in the expression of hub genes among convalescent, healthy, ICU, moderate, and severe COVID-19 patient groups in our validation cohorts. A P-value threshold of 0.05 was adopted to determine statistical significance in all tests. This comprehensive statistical approach ensures rigorous and robust analysis of the data, aiding in reliable and insightful conclusions.

Results

DEGs between COVID-19 and non-COVID-19 patients

Upon integrating the datasets GSE157103 and GSE152641, we assembled a matrix comprising 19,184 genes across 212 samples, including 162 from COVID-19 patients and 50 from non-COVID-19 individuals. The expression profiles were adjusted for batch effects and designated as the discovery dataset.

Analysis revealed 936 genes significantly upregulated and 125 significantly downregulated, adhering to the criteria: log2FC > 1, adjusted P-value < 0.05 for upregulation, and log2FC < -1, adjusted P-value < 0.05 for downregulation (Supplementary Table 1). We highlighted DEGs with log2FC > 2 and adjusted P-value < 0.01 in the volcano plot (Fig. 2A). The top 20 upregulated and downregulated genes were extracted from the normalized expression matrix for detailed visualization in a heatmap (Fig. 2B), indicating a predominance of upregulated genes.

Separate functional enrichment analyses were conducted for up- and down-regulated genes. Up-regulated genes are predominantly associated with organ or tissue-specific immune response, mucosal immune response, and viral life cycle regulation, among others (Fig. 2C). Conversely, down-regulated genes corresponded to pathways like leukocyte chemotaxis and various metabolic processes (Fig. 2D). Noteworthy is the overlap in immune response and signaling pathways across both gene sets, despite distinct pathway enrichments unique to each. Specifically, virus invasion pathways were exclusive to up-regulated genes, while metabolic pathways were more characteristic of down-regulated genes in GO-BP terms. This distinct pathway involvement underscores the complex interplay of biological processes in COVID-19 pathogenesis.

Fig. 2

Visualization of DEGs and functional enrichment GO-BP terms of significantly up- and down-regulated genes between COVID-19 and non-COVID-19 patients. A Volcano plot of DEGs. The red color dots indicate up-regulated genes, the blue dots indicate down-regulated genes, and the grey dots indicate the non-significantly expressed genes in COVID-19 compared to non-COVID-19. B The heatmap shows the expression pattern of the top 20 up- and down-regulated genes. The functional enrichment in GO-BP of up- (C) and down-regulated (D) genes. The x-axis represents the gene ratio value, while the y-axis represents the GO term

The difference in plasma cell infiltration between COVID-19 and non-COVID-19 patients

Gene expression profiles normalized to gene length (FPKM) were analyzed using the xCell algorithm [68]. Cell infiltration scores for 64 cell types across 212 samples were calculated and compared, as shown in Fig. 3 and Supplementary Table 2. Notably, a substantial number of cell types exhibited significant differences between the COVID-19 and non-COVID-19 cohorts. Within the lymphoid category, B cells, various T cell subsets (including CD4 + and CD8 + populations), memory B cells, naïve B cells, NKT cells, plasma cells, Pro B cells, Th1 cells, Th2 cells, and Tregs all demonstrated significant differential enrichment between COVID-19 and non-COVID-19 patients (Fig. 3A). Among myeloid cells, aDC, basophils, cDCs, DCs, iDCs, macrophages, and neutrophils were notably differentiated between the two patient groups (Fig. 3B). Other cell types such as astrocytes and neurons also showed significant differences (Fig. 3C), as did certain stem cells including CMP, GMP, megakaryocytes, and MEP (Fig. 3D), and stromal cells like endothelial cells, MSCs, myocytes, pericytes, and smooth muscle (Fig. 3E).

Such cellular disparities align with the known pathological progression of COVID-19, where viral infection triggers an exaggerated immune response, leading to hyperinflammation or a cytokine storm, further escalating multi-organ damage and increasing mortality risk [109].

Specifically, plasma cells were markedly more prevalent in COVID-19 patients compared to non-COVID-19 individuals (P-value ~ 2e-10), highlighting them as crucial for subsequent focused analyses. This significant increase underscores the potential role of plasma cells in the immune response against SARS-CoV-2, warranting further investigation into their involvement in disease progression and potential as targets for therapeutic intervention.

Fig. 3

Immune enumerate levels of 64 cell types between COVID-19 and non-COVID-19 patients. The 64 cell types were divided into five main categories. A Lymphoid. B Myeloids. C Others. D Stem cells. E Stromal cells. The highlighted red cell type is plasma cells. Student’s t-test was used for statistical analysis, and P < 0.05 was set as the cutoff. The boxplot shows the results from the analysis of xCell score values of 64 cell types in COVID-19 and non-COVID-19 cohorts based on the t-test statistical method. The ns, *, **, ***, and **** correspond separately to 1, 0.05, 0.01, 0.001, and 0.0001, respectively

Identification of modules and genes related to both COVID-19 and plasma cells

We constructed a co-expression network pertinent to plasma cells and COVID-19 using Weighted Gene Co-expression Network Analysis (WGCNA). This method groups genes into modules based on similar expression patterns and examines the correlation between these modules and specific traits [24, 25].

For this analysis, we established a network with a chosen soft power of 6, ensuring an R² > 0.9 for optimal network connectivity (Fig. 4A). The gene dendrogram and module color plot are illustrated in Fig. 4B. The upper panel shows the gene hierarchical cluster dendrogram, while the lower panel denotes the corresponding gene modules. This combination indicates that genes closely located in the dendrogram are grouped into the same module. A total of ten modules were identified: black (197 genes), blue (1,583 genes), brown (1,173 genes), green (472 genes), gray (657 genes), magenta (131 genes), pink (146 genes), red (214 genes), turquoise (8,928 genes), and yellow (794 genes). Notably, the black module exhibited a positive correlation with both plasma cell presence and COVID-19 status (Fig. 4C), with statistical significance (P-value < 0.05). Figure 4D depicts the correlation between module membership and gene significance within the black module. The findings indicate a strong association of the black module with plasma cells and COVID-19 (correlation = 0.32 and P-value = 4.6e-6). Consequently, genes within the black module were selected for subsequent detailed analysis due to their potential relevance to both COVID-19 and plasma cell activity (Supplementary Table 3).

Fig. 4

WGCNA for identifying the modules and genes. A The soft thresholding power’s scale-free fit index. The correlation between soft power and R² is demonstrated in the left panel, while the correlation between soft power and mean connectivity is shown in the right panel. B Gene and module color dendrogram. The y-axis indicates the distance between clusters. C Heatmap of the relationship between modules and clinical features of interest (plasma cells and COVID-19). D The scatter plot of the relationship between module membership and gene significance in the black module

Positive cluster biomarker genes of plasma cells in COVID-19 samples

Advancements in single-cell technology have allowed for an extensive exploration of cellular heterogeneity [110]. The cell type distribution in control and COVID-19 samples is depicted in Fig. 5A. We identified 19 distinct cell types: Airway epithelial cells, AT1, AT2, B cells, CD4 + T cells, CD8 + T cells, Cycling NK/T cells, Dendritic cells, Endothelial cells, Fibroblasts, Macrophages, Mast cells, Monocytes, Neuronal cells, NK cells, Other epithelial cells, Plasma cells, Smooth muscle, and Tregs. The cell counts for each cell type in both control and COVID-19 groups are specified in the supplementary materials.

The proportion of each cell type in the respective samples is visualized in Fig. 5B. Notably, there was a marked increase in the percentage of immune cells, particularly plasma cells (Fig. 5C), macrophages (Fig. 5D), and monocytes (Fig. 5E), in COVID-19 patients compared to controls (Supplementary Fig. 2). These alterations corroborate existing literature on the pathogenic mechanisms of COVID-19, highlighting the virus’s impact on the immune system and the subsequent hyperactive immune response, often referred to as a cytokine storm, leading to aggravated disease progression and potentially fatal outcomes [111, 112]. Consequently, 296 positive cluster biomarker genes specific to plasma cells were identified for subsequent analysis (Supplementary Table 4).

Fig. 5

The cellular landscape of COVID-19 patients. A Major clusters and respective cell-type assignments in UMAP deduce dimension visualization of control and COVID-19 samples. B The fraction of 19 cell types in 26 individual lungs. The difference in cell compositions from two groups of patients for (C) plasma cells, (D) macrophages, and (E) monocytes. The suffix “cov” represents COVID-19, and “ctl” represents control samples

Shared genes among DEGs, module genes, and positive cluster biomarker genes

We conducted an intersection analysis to determine commonalities among DEGs, module genes, and positive cluster biomarker genes, illustrated in Fig. 6A. We prioritized genes implicated in multiple dysregulated pathways due to their potential significance in pathogenesis, designating genes appearing in at least two of the three categories as “shared genes.” Consequently, 122 shared genes were identified, including ACOXL, AURKB, BHLHA15, and others listed in Supplementary Table 5.

The shared genes were postulated to contribute to at least two of the following three dysregulated mechanisms: (1) differential expression between COVID-19 and non-COVID-19 samples; (2) association with plasma cells and COVID-19 phenotype through co-expression patterns; and (3) regulation as positive cluster biomarker genes in plasma cells.

To investigate the biological relevance of these shared genes, we conducted GO and KEGG pathway enrichment analyses, the results of which are depicted in Fig. 6B-C. The analyses revealed significant enrichment in nuclear division, organelle fission, mitotic cell cycle phase transition, DNA helicase activity, ATP hydrolysis activity, cell cycle, p53 signaling pathway, human T-cell leukemia virus 1 infection, and protein processing in the endoplasmic reticulum. These findings align with the heightened cellular and immune responses observed in combating SARS-CoV-2 infection, indicating a vigorous immune response to viral invasion [113,114,115].

Fig. 6

Shared genes and enriched pathways. A The Venn diagram and UpSet plots show the overlapping genes among DEGs, module genes, and positive cluster biomarker genes. B The enriched GO terms of shared genes and the top 10 GO terms were selected to visualize each GO category. CC stands for cellular component, BP stands for biological process, and MF stands for molecular function. C The enriched KEGG pathways of shared genes

Analysis of PPI networks, sub-networks, and hub genes

Protein-protein interactions (PPIs) are essential for various biological functions including signaling, gene regulation, metabolism, and cell cycle control. In biomedicine, PPI network analysis is instrumental for identifying new therapeutic targets and unraveling disease mechanisms systematically. Generally, PPI network analysis employs databases like STRING for initial analysis, followed by tools like Cytoscape for network visualization and refinement [31, 116, 117].

Our analysis yielded an initial undirected network comprising 111 nodes and 4,886 edges. By applying stringent criteria, we refined this to a robust network consisting of 89 nodes and 1,283 edges, as illustrated in Fig. 7A. Using the cytoHubba plugin, we identified the top 20 hub genes: CCNA2, TOP2A, CCNB2, TPX2, BUB1, CDK1, KIF2C, RRM2, AURKB, CDCA8, KIF11, CCNB1, CDC20, UBE2C, KIF20A, PTTG1, TTK, DLGAP5, BUB1B, and PBK. A sub-network featuring these genes was subsequently delineated (Fig. 7B). An examination of the gene expression profiles of these hub genes revealed a consistent and significant upregulation in COVID-19 patients compared to non-COVID-19 individuals, as depicted in Fig. 7C. This observation underscores their potential relevance in the pathophysiology of COVID-19.

Fig. 7

PPI network and sub-network. A The PPI network. The node indicates the shared genes, while the edge indicates the interaction obtained from the STRING database. The depth of the node and edge represents the score of interaction. The larger the interaction combined score is, the darker the color. B The sub-network is mined from network analysis. The depth of the node and edge represents the degree of the node. The larger the degree is, the darker the color. C The expression pattern of 20 hub genes of sub-network. Up-regulated illustrates that the gene is significantly up-regulated in COVID-19 compared to non-COVID-19 samples

Validation and identification of candidate targets

Identifying clinical biomarkers and developing diagnostic and therapeutic tools are critical for advancing clinical applications [118, 119]. One of the key challenges in this process is extracting potential biomarkers with high sensitivity, stability, and accuracy from the vast amounts of multi-omic data available [120, 121].We employed two external RNA-seq datasets, GSE152418 and GSE179627, to validate the robustness and potential of our identified hub genes in terms of expression levels and classification capabilities. As depicted in Fig. 8A-E and Supplementary Fig. 3, a discernible difference in expression levels among various COVID-19 patient groups was noted for the 20 hub genes (P-value < 0.05), as per external dataset 1. Moreover, these genes exhibited higher expression in ICU, moderate, and severe COVID-19 patient cohorts compared to convalescent and healthy groups, aligning with inferences drawn from the discovery dataset.

Additionally, these 20 hub genes were evaluated for their ability to classify immune responses in external dataset 2, which encompassed diverse immune statuses including asymptomatic, recovering, re-detectable positive, uninfected, and symptomatic cases (Table 2). Notably, the genes TOP2A, BUB1, KIF2C, AURKB, CDCA8, KIF11, PTTG1, and BUB1B demonstrated AUC values exceeding 0.65, thereby qualifying them as promising candidate targets (Fig. 8F). These findings suggest that these eight genes, due to their substantial classification performance, are worthy of further exploration in the context of COVID-19 treatment.

Fig. 8

The validation and selection of target genes. The normalized expression levels of the CDCA8 (A), KIF11 (B), AURKB (C), BUB1 (D), and PTTG1 (E) genes in external dataset 1. F The ROC curves of the CDCA8, KIF11, and PTTG1 genes in by an external dataset 2

Table 2 Summary of the datasets used in this study

Docking model of candidate compounds and target proteins

Molecular docking has become an increasingly vital tool in computer-aided drug discovery, offering insights into the potential interactions between drugs and targets [122, 123]. While docking scores may not precisely predict drug inhibitory effects on proteins, the conformations derived from docking often closely align with crystal structures, making it a cost-effective approach for medicinal chemists [124, 125].

In our study, we utilized the cMAP database to identify candidate small compounds capable of down-regulating the expression of hub genes. We focused on compounds with normalized connectivity scores below zero and further refined our selection to those targeting any of the 20 hub genes. Ultimately, this process yielded candidate compounds for AURKB, KIF11, and TOP2A genes (Supplementary Tables 6–8). The leading candidates for down-regulating AURKB include GSK-1070916, KW-2449, and danusertib. For KIF11, the top candidates are BRD-K89997465, SB-743921, and ispinesib, and for TOP2A, the leading compounds are idarubicin, etoposide, and mitoxantrone.

AURKB, also known as aurora kinase B, is a serine/threonine kinase implicated in cell division. GSK-1070916 is noted as an Aurora kinase inhibitor [126,127,128]. BRD-K89997465 is known to down-regulate KIF11 expression through dopamine receptor antagonism, as indicated in the PubChem database [99]. TOP2A, or DNA topoisomerase II alpha, plays a role in DNA topology during transcription and is a known target of the anticancer drug idarubicin [129,130,131].

To assess the binding affinity between these compounds and their respective protein targets, we conducted molecular docking studies (Fig. 9A-C). The binding affinity of GSK-1070916 to AURKB is -6 kcal/mol, of BRD-K89997465 to KIF11 is -6.4 kcal/mol, and of idarubicin to TOP2A is -10 kcal/mol. Compounds with a free binding energy above − 5.5 kcal/mol are generally considered inactive [132,133,134]. Consequently, the most effective drug-protein interaction complexes are illustrated in Fig. 9.

A molecular dynamic simulation is a tool used to assess the protein-ligand system at the atomistic level and articulate the stability of the protein-ligand complex in a dynamic environment. In this study, a 100 ns simulation was performed to evaluate the complex’s stability. The RMSD of the protein and small compound drug were used to calculate the stability of the complex (Fig. 10A). The structural dynamics of the protein residues were investigated by calculating the RMSF values, which show the flexibility of protein residues in response to small compound drug binding during simulation (Fig. 10B).

Fig. 9

The interactions between small compound drugs and proteins. A GSK-1070916 and AURKB. B BRD-K89997465 and KIF11. C Idarubicin and TOP2A

Fig. 10

Molecular dynamic simulation of protein (KIF11) and small compound drug (BRD-K89997465). A RMSD plot for protein-small compound drug complex during 100 ns of molecular dynamics simulation. B RMSF trajectories of residues for protein-small compound drug complex

Discussion

In a matter of weeks, the COVID-19 pandemic rapidly spread to over 100 nations and territories, attracting substantial attention from the academic community and resulting in numerous studies and reports [1, 2, 135]. While some studies have proposed potential pathological mechanisms, therapeutic drugs, and target genes related to COVID-19, there remain several unresolved issues, such as drug side effects and the disease’s pathogenesis [8, 9, 37]. Particularly, there is a paucity of research on plasma cell biomarkers and therapeutic agents.

In this study, we focused on identifying and assessing candidate key genes and their associated drugs to enhance COVID-19 treatment, with a specific emphasis on plasma cells. Using cell infiltration analysis conducted via the xCell method, we observed significantly higher enrichment scores for immune cells in COVID-19 samples compared to normal controls. This finding aligns with the established conclusion that as COVID-19 advances, the virus-infected immune system triggers hyperinflammation or a cytokine storm, which can lead to severe organ damage and, in some cases, death [109]. Notably, we found a remarkably elevated enrichment score for plasma cells in COVID-19, underscoring their pivotal role in the disease’s development. Consequently, we employed Weighted Gene Co-expression Network Analysis (WGCNA) to identify co-expressed genes associated with COVID-19 and plasma cells. WGCNA, with its soft power approach, enhances the biological significance of strong correlations while diminishing weak or negative correlations, thus approximating a scale-free network [24, 25]. Simultaneously, we identified positive cluster biomarker genes highly expressed in plasma cells compared to the other 18 cell-type positive cluster biomarker genes within the plasma cell sub-population in COVID-19 samples, using single-cell RNA sequencing (scRNA-seq) analysis. These positive cluster biomarker genes are known to play pivotal roles in cell type-specific biological processes, including antiviral activity and immune response [136,137,138]. In-depth investigation of these genes aids in comprehending the pathogenic mechanisms and specific therapeutic targets, fostering the development of clinically relevant strategies. For instance, we identified IGKC as a positive cluster biomarker gene of plasma cells with the most significant fold change. IGKC is associated with the humoral immune response and has links to non-small cell lung cancer (NSCLC), breast cancer, and colon cancer [139,140,141]. Another positive cluster biomarker gene, IFNG-AS1, is implicated in COVID-19 and its severity, with its expression corresponding to increased inflammation levels [142, 143]. Furthermore, IFNG-AS1 is known to promote inflammation in ulcerative colitis [143].

Our systematic bioinformatics analysis unveiled shared genes among differentially expressed genes (DEGs), module genes, and positive cluster biomarker genes. In addition, we collected and organized the clinical characteristics of these five datasets, including three discovery datasets and two validation datasets, to test whether there is a significant effect between gender and COVID-19. The chi-square test of independence evaluates whether there is an association between the categories of the two variables. In this study, chisq.test() function when the number of samples is larger than 5 and fisher.test() function when the number of samples is smaller than 5, which was performed in R language to test the effect. The results validated there is no significant impact between gender and COVID-19 (p > 0.01, Supplementary Table 9). Given that proteins seldom function in isolation but rather form specific functional structures through interactions, these functional groups and complexes can be elucidated through network analysis of these interaction relationships [144, 145]. Therefore, we constructed a protein-protein interaction (PPI) network using Cytoscape. Within this network, we identified a sub-network comprising 20 hub genes (CCNA2, TOP2A, CCNB2, TPX2, BUB1, CDK1, KIF2C, RRM2, AURKB, CDCA8, KIF11, CCNB1, CDC20, UBE2C, KIF20A, PTTG1, TTK, DLGAP5, BUB1B, and PBK) based on Maximal Clique Centrality (MCC) values from the cytoHubba plugin. To validate the robustness of these 20 hub genes, we included two external independent bulk RNA-seq datasets with detailed clinical information and immune response data. Our analysis revealed significantly distinct expression patterns for these genes in COVID-19 patients, with higher expression levels in moderate, ICU, and severe cases compared to convalescent and healthy individuals. Notably, genes such as TOP2A, BUB1, KIF2C, AURKB, CDCA8, KIF11, PTTG1, and BUB1B demonstrated satisfactory classification ability in predicting the immune response to COVID-19. For instance, BUB1, also known as budding uninhibited by benzimidazole 1, encodes a serine/threonine-protein kinase that plays a central role in mitosis and has been implicated in the pathophysiological processes of COVID-19 [117, 146].

Furthermore, we identified three essential genes (AURKB, KIF11, and TOP2A) as potential therapeutic targets for COVID-19. Aurora kinase B (AURKB) belongs to the aurora kinase subfamily of serine/threonine kinases [147, 148]. Studies have demonstrated AURKB’s excellent classification performance in COVID-19 using common algorithms [148]. Kinesin family member 11 (KIF11) encodes a motor protein from the kinesin-like protein family, known to participate in various aspects of spindle dynamics and lung-related diseases [149]. Published results indicate that KIF11 is significantly up-regulated in lung adenocarcinoma and influences the cell cycle, tumor microenvironment alteration, and tumor-infiltrating immune cell proportions [150, 151]. DNA topoisomerase II alpha (TOP2A) produces an enzyme known as DNA topoisomerase, which regulates and modifies DNA’s topological states during transcription [152]. A study by Stephenson et al. found that TOP2A expression increased with the severity of COVID-19 according to single-cell multi-omics analysis [153]. Additionally, Samy et al. proposed TOP2A as a potential drug target, with digitoxin and hydroxyquinone being effective inhibitors against COVID-19 [154]. The literature review thus validates the reliability of our results, suggesting that controlling SARS-CoV-2 replication and excessive inflammation can be achieved by inhibiting these three genes.

We also identified three potential small compounds (GSK-1070916, BRD-K89997465, and idarubicin) with the lowest free binding energy, which were evaluated and visualized using cMAP and AutoDock. GSK-1070916, an Aurora inhibitor, has been reported to affect cell cycle progression and viability [127, 155]. Given AURKB’s involvement in chromosomal segregation and cytokinesis during mitosis, its inhibition by GSK-1070916 is particularly relevant [156]. A proteomics study by Bock et al. indicated that SARS-CoV-2 infection may impact numerous proteins related to inflammatory responses and chromosomal segregation [157]. Furthermore, DEGs between mild and severe COVID-19 cases were found to be enriched in cytokinesis failure, interferon response, and host translation inhibition [158]. Moreover, cytokinesis is a significant pathological mechanism underlying liver injury in COVID-19 [159]. BRD-K89997465, as recorded in the PubChem database, down-regulates KIF11 expression through a dopamine receptor antagonist mechanism [99]. KIF11, identified as a biomarker of COVID-19 through integrated blood RNA-seq and PPI network analysis, further underscores the potential significance of BRD-K89997465 [160, 161]. Idarubicin, an anthracycline anti-cancer drug used in the treatment of adult patients with acute myeloid leukemia, has also been found to inhibit SARS-CoV-2 endoribonuclease, particularly NSP15, a hexameric endoribonuclease unique to single-stranded RNA viruses [130, 162,163,164]. Our findings thus suggest that these three candidate small compounds positively impact COVID-19 treatment, with molecular docking models supporting their efficacy.

However, it is important to acknowledge the limitations of our study. It primarily relied on retrospective analysis, necessitating further research to ascertain the precise roles of the identified genes and potential treatment drugs. Additionally, in vitro and in vivo experiments are warranted to explore the immunotherapeutic efficacy of these interventions.

In conclusion, our study identified 122 shared genes among DEGs in COVID-19 and non-COVID-19 samples, co-expressed genes associated with both COVID-19 and plasma cell infiltration, and positive cluster biomarker genes for plasma cells in COVID-19 samples. We constructed a highly reliable PPI network comprising 89 genes and identified 20 hub genes with the highest MCC values, which were validated using external datasets. Furthermore, we proposed three potential small compounds (GSK-1070916, BRD-K89997465, and idarubicin) with promising therapeutic potential by targeting AURKB, KIF11, and TOP2A. These compounds offer a potential avenue for treating COVID-19 by modulating cell cycling and the hyperimmune response.