key: cord-0897091-dytcppg2
authors: Hoque, M. Nazmul; Khan, Md. Arif; Hossain, Md. Arju; Hasan, Md Imran; Rahman, Md Habibur; Soliman, Mahmoud E.; Araf, Yusha; Zheng, Chunfu; Islam, Tofazzal
title: Differential gene expression profiling reveals potential biomarkers and pharmacological compounds against SARS-CoV-2: insights from machine learning and bioinformatics approaches
date: 2022-03-30
journal: bioRxiv
DOI: 10.1101/2022.03.30.486356
sha: 00d2c868091ce235365a4fd38a1eee2aa2555e80
doc_id: 897091
cord_uid: dytcppg2

SARS-CoV-2 continues to spread and evolve worldwide, despite intense efforts to develop multiple vaccines and therapeutic options against COVID-19. Moreover, the precise role of SARS-CoV-2 in the pathophysiology of the nasopharyngeal tract (NT) is still unfathomable. Therefore, we used the machine learning methods to analyze 22 RNA-seq datasets from COVID-19 patients (n=8), recovered individuals (n=7), and healthy individuals (n=7) to find disease-related differentially expressed genes (DEGs). In comparison to healthy controls, we found 1960 and 153 DEG signatures in COVID-19 patients and recovered individuals, respectively. We compared dysregulated DEGs to detect critical pathways and gene ontology (GO) connected to COVID-19 comorbidities. In COVID-19 patients, the DEG– miRNA and DEG–transcription factors (TFs) interactions network analysis revealed that E2F1, MAX, EGR1, YY1, and SRF were the most highly expressed TFs, whereas hsa-miR-19b, hsa-miR-495, hsa-miR-340, hsa-miR-101, and hsa-miR-19a were the overexpressed miRNAs. Three chemical agents (Valproic Acid, Alfatoxin B1, and Cyclosporine) were abundant in COVID-19 patients and recovered individuals. Mental retardation, mental deficit, intellectual disability, muscle hypotonia, micrognathism, and cleft palate were the significant diseases associated with COVID-19 by sharing DEGs. Finally, we detected DEGs impacted by SARS-CoV-2 infection and mediated by TFs and miRNA expression, indicating that SARS-CoV-2 infection may contribute to various comorbidities. These pathogenetic findings can provide some crucial insights into the complex interplay between COVID-19 and the recovery stage and support its importance in the therapeutic development strategy to combat against COVID-19 pandemic. IMPORTANCE Despite it has now been over two years since the beginning of the COVID-19 pandemic, many crucial questions about SARS-CoV-2 infection and the different COVID-19 symptoms it causes remain unresolved. An intriguing question about COVID-19 is how SARS-CoV-2 interplays with the host during infection and how SARS-CoV-2 infection can cause so many disease symptoms. Our analysis of three different datasets (COVID-19, recovered, and healthy) revealed significantly higher DEGs in COVID-19 patients than recovered humans and healthy controls. Some of these DEGs were found to be co-expressed in both COVID-19 patients. They recovered humans supporting the notion that DEGs level is directly correlated with the viral load, disease progression, and different comorbidities. The protein-protein interaction consisting of 24 nodes and 72 edges recognized eight hub-nodes as potential hub-proteins (i.e., RPL4, RPS4X, RPL19, RPS12, RPL19, EIF3E, MT-CYB, and MT-ATP6). Protein–chemical interaction analysis identified three chemical agents (e.g., Valproic Acid, Alfatoxin B1, and Cyclosporine) enriched in COVID-19 patients and recovered individuals. Mental retardation, mental deficiency, intellectual disability, muscle hypotonia, micrognathism, and cleft palate were the significant diseases associated with COVID-19 by sharing DEGs.

In late December 2019, a novel respiratory disease, now popularly termed as "COVID-19", caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), emerged in Wuhan, China [1] [2] [3] . Immediately after its first outbreak in China, this fearsome virus has emerged as one of the deadliest human pathogens [4] . Due to its worldwide spread and severity, COVID-19 has been declared a public health emergency of international concern by the World Health Organization (WHO) [5, 6] . As of January 10, 2022, COVID-19 disease affected 217 countries and territories, and more than 308 million cases have been confirmed around the globe, with about 5.6 million deaths [7] . In the early stage of the outbreak, the spectrum of clinical manifestations of COVID-19 ranges from the common cold to respiratory failure depending on the demography and environments [2, 5, 8] . However, recent data show that the clinical episodes of COVID-19 may range from asymptomatic infection to critical illness, with a dysregulated inflammatory response to infection a hallmark of severe cases and lifethreatening multi-system organ failure [8] [9] [10] [11] . In most cases (~ 80%), patients exhibit mild symptoms, while the remaining ∼20% may develop severe lung injury and death from respiratory failure [12, 13] . Some of the clinically infected patients may suffer from acute respiratory distress syndrome (ARDS) and multiple organ failure requiring intensive care unit (ICU) facilities for life support and medication [13] . Risk factors for severe SARS-CoV-2 include age, smoking status, ethnicity, and male sex [10, 14] . Notably, the persistence and prognosis of COVID-19 are greatly influenced by the patients' underlying health conditions and age [9, 15] . Given no effective antiviral treatment and slow vaccine rollout, COVID-19 continues to be a serious threat to public health worldwide [16] .

Despite increasing global threats of COVID-19, the host immune response against SARS-CoV-2 infection remains poorly understood, and the perturbations result in a severe outcome [12, 17] . The nasal epithelium serves as a portal for initial infection and transmission of the SARS-CoV-2 [5] . SARS-CoV-2 employs ACE2 (Angiotensin-converting enzyme 2) as a receptor for cellular entry [18] , and the binding affinity of the S protein and ACE2 was found to be a major determinant of SARS-CoV-2 replication rate and disease severity [17, 19] . After the entrance into the susceptible host, SARS-CoV-2 infects cells of the respiratory epithelium and mucous membranes, such as those of the nose or eyes [18, 20] . The host immune response to SARS-CoV-2 infection involves activating cellular and humoral arms. The innate immune system recognizes the SARS-CoV-2 RNAs through three major classes of cytoplasmic pattern recognition receptors: Toll-like receptors (TLRs), RIG-I-like receptors (RLRs), and NOD-like receptors (NLRs) [17, 21] . This response involves the release of interferons (IFNs) and inflammatory cytokines, including the IL-1 family, IL-6, and TNF, that activates a local and systemic response to infection [5, 17] . This inflammatory response cascade involves the recruitment, activation, and differentiation of innate and adaptive immune cells, including neutrophils, inflammatory myeloid cells, CD8+T cells, and natural killer (NK) cells [12] . The resolution of infection is largely dependent on the cytotoxic activity of CD8+T cells and NK cells, which enable the clearance of virus-infected cells [5, 17] . It is believed that dysregulated host immune response leads to the persistence of virus-infected cells and may facilitate a hyperinflammatory state termed Macrophage (MΦ) activation syndrome (MAS) or "cytokine storm", and ultimately damage to the infected lung [12, 17] . However, the underlying molecular mechanisms of the aberrant inflammatory responses in serology and histopathology under SARS-CoV-2 infection are still not clear.

The ongoing pandemic of SARS-CoV2 and lack of comprehensive knowledge regarding the progression of COVID-19 has constrained our ability to develop effective treatments for infected patients. To obtain a complete understanding of the host response to SARS-CoV2, one means to examine gene expression in relevant tissues. Until now, the scant amount of gene expression profiles are available from patients with COVID-19 and have yielded some insights into the pathogenic processes triggered by infection with SARS-CoV-2 [12, 17, 22] .

Transcriptomic analyses of cells upon viral infections are extremely useful to identify the host immune response dynamics and gene regulatory networks [12, 23] . However, because of the limited number of samples and preliminary analysis, a full picture of the biological state of SARS-CoV2-affected tissues has not emerged. To address this, we have employed RNA-Seq techniques to investigate the upper airway (nasopharyngeal tract) gene expression profile in 22 specimens of COVID-19 patients (n = 8), recovered (n = 7), and healthy (n = 7) individuals using several orthogonal bioinformatic tools to provide a complete view of the nature of the COVID-19 inflammatory response and the potential points of therapeutic intervention.

Through differentially expressed genes (DEGs) analyses in these datasets, we identified several genes coding for translational activities (e. g. RPL4, RPS4X, RPL19, RPS12, RPL19, EIF3E), ATP-synthesis (MT-CYB, MT-ATP6), transcription factors (e. g. E2F1, MAX, EGR1, YY1, SRF), hub-proteins (e. g. KIAA0355, DCUN1D3, FEM1C, ARHGEF12, THBS1), and mi-RNA (e. g. hsa-miR-19b, hsa-miR-495, hsa-miR-340, hsa-miR-101, and hsa-miR-19a) evidencing a sustained inflammation and cytokine storm in the COVID-19 patients.

Differentially expression and distribution of DEGs. To elucidate whether differentially expressed genes (DEGs) contribute to the SAR-CoV-2 inflammatory response and the potential points of therapeutic intervention, we analyzed 22 RNA-seq data of nasopharyngeal epithelial tissue of COVID-19 patients, recovered humans, and healthy controls. To perform RNA-seq analysis, we retrieved datasets from the National Center for Biotechnology Information (NCBI) that belonged to previously published BioProject under accession number PRJNA720904 (https://www.ncbi.nlm.nih.gov/bioproject). We identified 1960 and 153 gene signatures in COVID-19 patients and recovered human NT epithelial tissues, which were differentially expressed compared with healthy controls. We particularly focused on the dysregulation (up or down-regulation) of the identified DEGs during SARS-CoV-2 pathophysiology and its overlap with the recovered or healthy states of the humans. The volcano plots presented in Fig. 2 show the DEGs for COVID-19 with the red dots. The number of shared DEGs between COVID-19 and recovered datasets is presented in the Venn diagrams ( Fig. 2 C, D) . Thirty-seven shared DEGs were identified between COVID-19 patients and recovered subjects. Of the detected DEGs, 1,510 (77.04%) genes were upregulated (Up) during SARS-CoV-2 pathogenesis, of which 1,489 (98.61%) genes had a sole association with COVID-19 patients. Likewise, 90 (58.82%) genes were upregulated in recovered humans, and of them, 69 (76.67%) genes had a sole association with the recovery phage of SARS-CoV-2 infection (Fig. 2C) . By comparing the upregulated genes between COVID-19 patients and recovered individuals, we found that 21 genes (i.e., RPL4, MT-ND2, SCD5, MT-CYB, EZR) were shared between the conditions (Fig. 2C) . On the other hand, 450 (22.96%) and 63 (41.18%) DEGs were downregulated (Down) in COVID-19 patients and recovered subjects, respectively, and of them, only 12 genes (i.e., MAFF, ARHGEF12, DCUN1D3, DR1, MT-CO1.) were found to be shared between COVID-19 and recovered cases (Fig. 2D) . The DEGs shared between COVID-19 positive and recovered people and their relationships from the perspective of adjusted P-value, and log2 fold-change is presented in the heatmaps, respectively ( Fig. 3 A, B) . Finally, 33 common dysregulated (Up or Down) genes were presented in a bubble plot to show relationships with 10 log fold-change values (Fig. 3C) . patients. The 33 shared DEGs were used to identify key pathways and GOs that may be linked to COVID-19 comorbidities. We combined all the KEGG and Reactome pathway databases with Enrichr tools to create a single pathway database. We looked at the pathways whose significance was determined by the P-value and plotted the top 20 pathways for each condition (Fig. 4) . Consideration was given to the paths having a higher logarithmic P-value. The most significant pathways were the ribosome signaling pathway, coronavirus signaling pathway, and c-type-lectin receptor signaling pathway for KEGG analysis (Fig. 4B) and forming a pool of free 40S subunits 3-UTR-mediated translational regulation, and eukaryotic translational initiation signaling pathways for the Reactome database (Fig. 4B) .

We used the Enrichr tool to identify significantly enriched cellular signaling pathways and functional GO terms (molecular function, biological process, and cellular component) with DEGs in the nasopharyngeal epithelial cells on COVID-19 patients. The 33 shared DEGs were used to identify key pathways and GOs that may be linked to COVID-19 comorbidities. We looked at the pathways whose significance was determined by the P-value (having a higher logarithmic P-value) and plotted the top 20 pathways for each condition (Fig. 4) .

Furthermore, we have also conducted GO term enrichment using the same set of common DEGs. We employed the GO biological process, the GO molecular function, and the GO cellular component databases obtained from Enrichr libraries. The significantly enriched GO terms were identified if the enrichment yields the high logarithmic value of the adjusted Pvalue. The top 20 cellular signaling pathways in the COVID-19 patient's nasopharyngeal epithelial cells were selected in this study (Fig. 5 ) in relevance to the recovered phase. The most significant GO pathways were the ceramide 1-phosphate transfer activity, and ceramide 1-phosphate binding pathways for the molecular functions ( Fig. 5A) , database nucleartranscribed mRNA catabolic process and regulation of epithelial cell differentiation pathways for biological process (Fig. 5B) , and membrane raft, and cytosolic sizeable ribosomal subunit pathways for cellular component (Fig. 5C ).

network was built from the common DEGs interactions, consisting of 24 nodes and 72 edges.

The PPI network clustering highlighted RPL4, RPL18A, EIF3E, EIF3D, RPS4X, RPL19, EIF3K, RPS12, MT-ND2, MT-CO1, MT-ATP6, and MT-CYB with high interaction activity ( Fig. 6) . The proteins with several connecting edges can be identified as hub proteins. Fig. 7 shows the top-10 hub-nodes within the PPI network. As anticipated by five different methods (i.e., maximum neighborhood component; MNC, betweenness, degree, edge percolated component; EPC and maximal clique centrality; MCC), we recognized eight hub-nodes as potential hub-proteins (i.e., RPL4, RPS4X, RPL19, RPS12, RPL19, EIF3E, MT-CYB, and MT-ATP6) ( Fig. 7 A-E). Interestingly, these eight hub proteins were common in all methods.

Only RPS4X was found from 4 methods except for betweenness (Fig. 7B) . Conversely, the betweenness method predicted only three proteins (i.e., SCD5, EZR, and RIMS1) from the shared DEGs as hubs that were not found by other methods (Fig. 7B) . represented by different lines linking them. Nodes in a network that connect multiple edges are considered significant nodes since they are more critical. Out of 21 miRNAs detected, hsa-let-7e-5p, hsa-mir-7977, hsa-mir-155-5p, hsa-mir-186-5p, and hsa-mir-1827 were the most expressed miRNAs and had a stronger association with DEGs ( Fig. 8A) . Likewise, among the DEGs, DMD, AHDC1, BAG4, EMP2, TIMM50, RPL7L1, and THBS1 were more significant since these DEGs have a higher degree (number of connecting edges) among the others and miRNAs (Fig. 8A) . We further studied the interactions between TF and DEGs and identified 14 TFs, of which FOXC1, FOXL1, NFIC, YY1, and PPARG were significantly enriched and showed more interactions with DEGs (Fig. 8B) .

Apart from these, the present study included TFs and miRNAs that were highly relevant to SARS-CoV-2 interactions. This analysis identified 19 hub proteins, 10 TFs, and 5 miRNAs (Fig. 9A) . In COVID-19 interaction, the TF-miRNA network showed that E2F1, MAX, EGR1, YY1, and SRF were the highly expressed TFs, and hsa-miR-19b, hsa-miR-495, hsa-miR-340, hsa-miR-101, and hsa-miR-19a were among significant miRNAs (Fig. 9A) .

possible pathogenesis mechanisms and drug interactions that may not be evident using conventional approaches. To disrupt the SARS-CoV-2 pervasiveness, we sought to find pharmaceutical compounds that interact with viral proteins (Methods). We detected nine pharmacological compounds (for example, famoxadone, ubiquinone-2, 2-nonyl-4hydroxyquinoline, 5-n-undecyl-6-hydroxy-4,7-dioxobenzothiazole) acting against one protein, the human mitochondrial cytochrome b (MT-CYB) (Fig. 9B) .

Protein-chemical interaction (PCI) is an important study to understand the functionality of proteins underpinning the molecular mechanisms within the cell, which may also help in drug discovery. For example, it has been discovered that SARS-CoV-2 infection causes PCI networks in the COVID-19 patients and recovered humans. Fig. 10A depicts a network of PCI among significant proteins. The significant proteins identified from this network include FEM1C, NCALD, THBS1, PCDH9, DMD, and PDGFA. Similarly, we identified three chemical agents such as Valproic Acid, Alfatoxin B1, and Cyclosporine enriched in this interaction analysis (Fig. 10A) .

hypothesizes that many conditions can be associated or connected with COVID-19 by sharing some common genes among themselves. Disorder-specific therapeutic interface strategies attempt to discover the link between genes and diseases. In this study, we found 14 other diseases associated with COVID-19 by sharing four DEGs (i.e., DMD, C2CD3, WNT3, and AHDC1) most prevalent in COVID-19. Of the detected diseases, mental retardation, mental deficiency, intellectual disability, muscle hypotonia, micrognathism, and cleft palate were the significant diseases interconnected with COVID-19 (Fig. 10B) . DISCUSSION The SARS-CoV-2 infection causes many illnesses, from minor respiratory sickness that is generally asymptomatic to severe pneumonia that results in multisystem failure and death. [25] , modulating cell proliferation and differentiation (SCD5) [26] , mitochondrial deficiencies, and associated disorders (MT-CYB) [27] , epithelial marker ezrin (EZR) associated with cell surface structure adhesion, migration, and organization of the SARS-CoV-2 [28] were found to be co-expressed in the nasopharyngeal epithelial cells of COVID-19 patients and recovered humans (Fig. 2C) . Conversely, SARS-CoV-2 infection suppressed the expression genes associated with transcription factors (MAFF) [29] , erythropoiesis (ARHGEF12) [30] , membrane neddylation (DCUN1D3) [31] , global regulator of transcription (DR1) [32] , and cytochrome-c oxidase activity (MT-CO1) [33] in both COVID-19 patients and recovered humans (Fig. 2D) .

We next investigated whether host gene expression during SARS-CoV-2 pathophysiology is associated with functional enrichment, for example, cell signaling pathways and gene ontology. Our results showed that DEGs related to ribosome signaling pathway, coronavirus signaling pathway, c-type-lectin receptor signaling pathway, forming a pool of free can augment immunity and control COVID-19 infection by enhancing autophagy, adaptive immunity (Th1 programming) MHC-I dependent cytotoxic T lymphocytes (CTL) response [34] . The epithelium lining the airways plays a key role in the defense against infections.

Several lines of evidence showed that SARS-CoV-2 infection induces epithelial barrier function, as documented by decreased trans-epithelial resistance, increased permeability, and altered distribution of the tight junction protein [35, 36] . However, this functional impairment remained transient, with signs of epithelial regeneration during the recovery phage of SARS-CoV-2 infection. Basal cell mobilization and replication can also be observed to exert a moderate effect on epithelial barrier integrity [35] . With these dysregulated genes, we have conducted the PPI network analyses. PPIs network analysis is the most prominent section of the study as hub gene detection, analysis of modules, and drug identification thoroughly depends on the PPIs network. According to the PPIs network (Fig.s 6 and 7) , ribosomal proteins (RPL4, RPS4X, RPL19, RPS12, RPL19), translation initiation factor 3 subunit E (EIF3E), mitochondrial deficiencies and associated disorders (MT-CYB) [27] , and cytochrome-c oxidase activity (MT-CO1) [33] , and mitochondrial oxidative phosphorylation (MT-ATP6) proteins were declared as hub proteins because of their high interaction rate or degree value.

MT-ATP6) involved in ATP synthesis and respiratory activity, oxidative stress, proinflammatory state, and cytokine production [37] . The increased expression of ribosomal proteins can be attributed to the fact that the virus hijacks the host's translational machinery for its survival by the mechanisms such as ribosome shunting and phosphorylation of ribosomal proteins [24]. These hyper-interactive proteins' PPI analysis and gene enrichment analysis showed significant biological functions connected to COVID-19 related to cell signaling pathway and the host response to SARS-CoV-2 infections [38] . As discussed earlier, these proteins are involved in several other disorders [22, 37, 38] .

We then investigated the protein-protein, gene-miRNA, TF-gene, protein-drug, and protein-chemical interactions of the common DEGs in COVID-19 patients and recovered humans. Our results showed that hsa-miR-19b, hsa-miR-495, hsa-miR-340, hsa-miR-101, and hsa-miR-19a were the mostly expressed miRNAs (Fig. 8A) , and E2F1, MAX, EGR1, YY1, and SRF were the highly expressed transcription factors (TFs) (Fig. 8B ). While host responses to infection are critical in differential outcomes of SARS-CoV-2 infection, the role of miRNAs in COVID-19 pathogenesis is poorly understood. We observed that most of these miRNAs were strongly upregulated in COVID-19 patients, which could be used as the circulating biomarkers for the diagnosis or prognosis of COVID-19 [39] . Circulating miRNAs are extracellular serum/plasma miRNAs that could be involved in cell-cell communication and thus, might contribute to disease progression. Besides their diagnostic value, miRNAs are well known for their therapeutic potential, especially in viral diseases. A recent report has compared the miRNA signature in the peripheral blood of COVID-19 patients versus healthy donors, and several miRNAs have been identified to be deregulated and interfered with the shaping of the immune responses [40] . Therefore, the upregulated levels of miRNAs could be involved in the inflammatory storm seen in COVID-19 patients via inhibiting the immunosuppressive and antiinflammatory role ensured by the transcription signaling pathway. Recent studies reported that transcription of mRNAs in epithelial cells is induced by TNF-α and triggers a negative feedback loop involving E-selectin to control inflammatory signaling [24, 41] . Although we identified the 14 TF-genes showing more interactions with DEGs, we further assess whether these genes have the potential causal effects on the COVID-19 development. Our networkbased approach identified TF hubs that likely regulate many cellular functions (e.g., cytokine storm) overexpressed i COVID-19 patients. Previous research identified 95 TFs in cytokines upregulated in the COVID-19 patients, and of them, 19 TFs are targets of FDA-approved drugs [42] . Targeting TFs associated with the cytokine releasing syndrome provides candidate drugs and targets to treat COVID-19 [42] . However, additional research is needed to determine whether these combinations elicit the same immunomodulatory response in the context of SARS-CoV-2 infection.

Nine pharmacological compounds were found to be effective against SARS-CoV-2, and of them, fungicide (famoxadone), blood pressure controlling coenzyme (ubiquinone-2), VPA can reduce the SARS-CoV-2 receptor ACE-2 expression level and be used as a potential drug candidate for the prevention strategy against COVID-19 [44] . Aflatoxin B1 (AFB1), which alters immune responses to mammals, is one of the most common mycotoxins in feeds and food and a potential aggravating risk factor in COVID-19 patients [45] . The effect of cyclosporine on coronaviruses, including the new SARS-CoV2, has been extensively studied [46] . Several earlier studies showed that cyclosporine could prevent uncontrolled inflammatory response, SARS-CoV-2 replication, and acute lung injury [46, 47] . Therefore, there is an urgent anxiety, and depression [48] . These complications are not necessarily short-lived and can cause long-term effects of multi-organ injury following SARS-CoV-2 infections. COVID-19 appears to present a greater risk to people with intellectual and developmental disabilities, especially at younger ages, and recent evidence suggests that mental health problems significantly increased worldwide during this pandemic [49, 50] . Since muscle possesses the ACE2 receptor to which SARS-Cov-2 binds, it follows that the involvement of the muscle could be due not only to the secondary effects of the infection (e.g., reduced oxygen supply from persistent lung disease, perfusion defects from cardiovascular defects, and vascular damage) but also to the direct action of virus (SARS-Cov-2 myositis) [51] .

CONCLUSIONS Gene expression analysis may potentially reveal disease-pathogenesis pathways and point to novel targets for potential therapeutic approaches. This study examines the RNA-seq data of COVID-19 patients, recovered persons, and healthy individuals to find DEGs and biomarkers between the SARS-CoV-2 pathogenesis and recovery stage from a molecular and cellular standpoint. We found that COVID-19 patients had a much larger number of DEGs than recovered humans and healthy controls and that some of these DEGs were co-expressed in both COVID-19 patients and recovered humans. We used gene expression analysis with the biomarker to identify cellular signaling pathways and GO terms.

In the COVID-19 patients, we found several genes coding for translational activities, transcription factors, hub-proteins, and miRNA expressions, all of which indicated a persistent inflammation and cytokine storm. The signaling pathways, GO terms, and chemical compounds discovered in this study could help researchers Fig. out how genes are linked together to find possible therapeutic approaches. However, the DEGs' direct molecular biological functions and significant pathways discovered in this study should be investigated further to understand better the mechanisms underlying the host response to SARS-CoV-2 and identify potential therapeutic targets and drug candidates for COVID-19.

Overview of the proposed bioinformatics pipelines. Network-based approaches are common to identify and analyze the pathogenesis of SARS-CoV-2. Datasets required in this work were constructed and collected at the initial phase and detailed in the following subsections. Gene expression analysis was performed to identify the DEGs from each dataset ( Fig. 1) . Next, the common DEGs between two groups of COVID-19 datasets were identified.

These common DEGs were further used to discover their protein-protein interactions (PPIs) and to perform gene set enrichment analysis (GSEA) to identify enriched cell signaling pathways and functional gene ontology (GO) terms. Next, the same set of common DEGs was used to discover three types of GRNs: DEGs-micro RNAs (miRNA) network, DEGstranscription factors (TFs) network, and TF-miRNA network. Finally, protein-chemical compound and protein-drug interactions were also investigated for the common DEGs (Fig. 1) . The demographic information on the study subjects, COVID-19 diagnosis, sample collection, and details about library preparation and RNA sequencing techniques is available in the previously published articles of Hoque et al. [52] .

Dataset preparation and analysis of differentially expressed genes. This study prepared two datasets as COVID-19 patients versus recovered humans with the same healthy control group for analytical purposes. We performed several statistical operations on the datasets to determine the DEGs. Moreover, the Benjamini-Hochberg false discovery rate method was used to provide a good balance between the discovery of statistically significant genes and the limitation of false positives. The BioJupies generator online server (https://maayanlab.cloud/biojupies/) was used for RNA-seq raw data analysis [1] . In this study, genes with adjusted P-value < 0.05 and absolute value of log2 fold-change ≥ 1 were considered DEGs. Next, to determine the shared DEGs, we compared two COVID-19 datasets between each other using the Venny v2.1 web tool [2] . In this article, we use the term combined DEGs' to refer to the collection of these two sets of DEGs, which have been used in the downstream bioinformatics analyses.

Functional enrichment analysis. We utilized Enrichr [3] with Fisher's exact test to conduct the functional enrichment analysis with the combined DEGs. After performing an overrepresentation analysis, a collection of enriched cellular signaling pathways and functional GO keywords were discovered, revealing the biological importance of the previously detected DEGs. In Enrichr analysis, we combined the signaling pathways from two libraries, including KEGG and Reactome, to create a single route. Only the important paths for which the P-value was less than 0.05 were evaluated and considered after deleting duplicate pathways. For functional GO annotations, we looked at the GO biological process, GO molecular function, and GO cellular component datasets in Enrichr. We selected the most important GO terms based on a set of criteria and with an adjusted P-value < 0.05.

The shared DEGs' protein-protein interaction (PPI) was analyzed using the STRING database [4] . We applied different local-and global-based methods using the cytoHubba plugin [5] in Cytoscape v3.8.2 [6] to determine potential hubs proteins within the PPI network. While the local method ranked hubs based on the relationship between the node and its direct neighbor, the global method ranked hubs based on the interaction between the node and the whole network. In total, five different methods were considered, including three local rank methods, i.e., degree, maximum neighborhood component (MNC), maximal clique centrality (MCC), and two global rank methods, i.e., edge percolated component (EPC) and betweenness. Next, we compared the results and identified the common nodes as the most potential hubs proteins. Finally, the protein networks were analyzed through Cytoscape v3.8.2.

Differential gene regulatory network (GRN) analysis. The finding of DEG-miRNA, TF-DEG, and TF-miRNA interaction networks is a part of the GRN analysis. Using the Network Analyst platform [7] , the commonly dysregulated genes were utilized to identify GRN networks. Discovering DEG-miRNA interaction networks was accomplished by using the miRTarBase [8] database. The JASPAR [9] database was used to identify the TF-DEG interaction network. Employing TF-miRNA coregulatory network database, the TF-miRNA interaction was analyzed. The networks were filtered with a betweenness value of 100 and degree centrality of 0 to 10 to remove unnecessary information.

Protein-chemical compound analysis. Analyses of protein-chemical compounds can identify the chemical molecules responsible for the interaction of proteins in comorbidities [53] . For example, this study found protein-chemical interactions using the enriched gene (common DEGs) that COVID-19 patients developed several digestive problems. Furthermore, using the Comparative Toxicogenomics Database [10] , we have identified the proteinchemical interactions through Network Analyst [7] .

Protein-drug interaction network. One of the key goals of this study is to identify potential therapeutic compounds which could be effective in mitigating SARS-CoV-2 pervasiveness. Using the shared DEGs, we constructed the protein-drug interaction (PDI) network through the NetworkAnalyst v3.0 web server [7] in conjunction with the DrugBank v5.0 database (https://go.drugbank.com/docs/drugbank_v5.0.xsd). We downloaded the network data and conFig.d the data with Cytoscape v3.8.2 to aid the analysis [6] .

Gene-disease association network. DisGeNET (https://www.uniprot.org/database/DB-0218) is a standardized gene-disease association database that incorporates correlations from various sources involving various biological features of disorders. It emphasizes the increasing understanding of human genetic illnesses. We examined the gene-disease connection using a network analyzer [11] to find diseases and chronic problems associated with common DEGs. 

Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China

Genomic diversity and evolution, diagnosis, prevention, and therapeutics of the pandemic COVID-19 disease

Probable pangolin origin of SARS-CoV-2 associated with the COVID-19 outbreak

Epitope-based chimeric peptide vaccine design against S, M and E proteins of SARS-CoV-2, the etiologic agent of COVID-19 pandemic: an in silico approach

Microbial co-infections in COVID-19: Associated microbiota and underlying mechanisms of pathogenesis

World Health Organization declares global emergency: A review of the 2019 novel coronavirus (COVID-19)

Johns Hopkins coronavirus resource center

Prevalence of comorbidities among individuals with COVID-19: A rapid review of current literature

Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study

Chronic lung diseases are associated with gene expression programs favoring SARS-CoV-2 entry and severity

Upper airway gene expression reveals suppressed immune responses to SARS-CoV-2 compared with other respiratory viruses

Comprehensive transcriptomic analysis of COVID-19 blood, lung, and airway

Differential expression of viral transcripts from single-cell RNA sequencing of moderate and severe COVID-19 patients and its implications for case severity

Risk factors of critical & mortal COVID-19 cases: A systematic literature review and meta-analysis

Transcriptomic studies revealed pathophysiological impact of COVID-19 to predominant health conditions

An interactive web-based dashboard to track COVID-19 in real time

Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in COVID-19 patients

Differentially expressed immune response genes in COVID-19 patients based on disease severity

SARS-CoV-2 entry factors are highly expressed in nasal epithelial cells together with innate immune genes

Respiratory virus infections: understanding COVID-19

Viral innate immune evasion and the pathogenesis of emerging RNA virus infections

Imbalanced host response to SARS-CoV-2 drives development of COVID-19

Gene expression profiling reveals the shared and distinct transcriptional signatures in human lung epithelial cells infected with SARS-CoV-2, MERS-CoV, or SARS-CoV: potential implications in cardiovascular complications of COVID-19

The expression characteristics of mt-ND2 gene in chicken

Stearoyl-CoA desaturase 5 (SCD5), a Δ-9 fatty acyl desaturase in search of a function

Circulating mitochondrial DNA is an early indicator of severe illness and mortality from COVID-19

SARS-CoV-2 infection and replication in human gastric organoids

The small MAF transcription factors MAFF, MAFG and MAFK: current knowledge and perspectives

ARHGEF12 regulates erythropoiesis and is involved in erythroid regeneration after chemotherapy in acute lymphoblastic leukemia patients

The human Dcn1-like protein DCNL3 promotes Cul3 neddylation at membranes

The Dr1/DRAP1 heterodimer is a global repressor of transcription in vivo

Cytochrome c oxidase deficiency

Host sphingolipids: Perspective immune adjuvant for controlling SARS-CoV-2 infection for managing COVID-19 disease

SARS-CoV-2 infection induces the dedifferentiation of multiciliated cells and impairs mucociliary clearance

A novel coronavirus from patients with pneumonia in China

Implications of oxidative stress and potential role of mitochondrial dysfunction in COVID-19: Therapeutic effects of vitamin D

Network-based identification genetic effect of SARS-CoV-2 infections to Idiopathic pulmonary fibrosis (IPF) patients

Circulating miR-29c, miR-30c, miR-193a-5p and miR-885-5p: Novel potential biomarkers for HTLV-1 infection diagnosis, Infection

Differential microRNA expression in the peripheral blood from human patients with COVID-19

miR-19a/b and miR-20a promote wound healing by regulating the inflammatory response of keratinocytes

Therapeutic Targeting of Transcription Factors to Control the Cytokine Release Syndrome in COVID-19

Medicinal significance of benzothiazole scaffold: an insight view

Immunopharmacological management of COVID-19: Potential therapeutic role of valproic acid

Herbal Immunity Booster-Associated Liver Injury During COVID-19 Pandemic and Aflatoxins

Cyclosporine and COVID-19: Risk or favorable?

Cyclosporine A: a valid candidate to treat COVID-19 patients with acute respiratory failure

COVID-19 and comorbidities: Deleterious impact on infected patients

Impact of COVID-19 pandemic on mental health among general Bangladeshi population: a cross-sectional study

Intellectual and developmental disability and COVID-19 case-fatality trends: TriNetX analysis

ACE2 receptor polymorphism: Susceptibility to SARS-CoV-2, hypertension, multi-organ failure, and COVID-19 disease outcome

SARS-CoV-2 Infection Reduces Human Nasopharyngeal Commensal Microbiome With Inclusion of Pathobionts

Network-based analysis of fatal comorbidities of COVID-19 and potential therapeutics

The authors declare no competing interests.

All data needed to evaluate the conclusions in the paper are present in the paper. The sequence data reported in this article has been deposited in the National Center for Biotechnology Information (NCBI) under BioProject accession number PRJNA720904.

This manuscript utilized RNA-seq data from a publicly available repository, i.e., NCBI, and has no ethical issue to be declared.