The experimentally determined three-dimensional structures of SARS-CoV-2 [envelope protein (PDB ID: 7K3G), main protease (PDB ID: 5R7Y), spike glycoproteins (PDB ID: 6VSB, 6VXX), HR2 domain (PDB ID: 6LVN, 6M1V), NSP2 (PDB ID: 7MSW, 7MSX), NSP3 (PDB ID: 6W02, 6W6Y, 6W9C), NSP7 (PDB ID: 6WTC), NSP8 (6M5I), NSP9 (PDB ID: 6W9Q, 6W4B, 6WCI), NSP12 (PDB ID: 6XEZ, 6M71), NSP13 (PDB ID: 5RL6), NSP14 (PDB ID: 7DIY), NSP15 (PDB ID: 6VWW, 7K0R), NSP16 (PDB ID: 6W4H), ORF3a (PDB ID: 6XDC), ORF7a (PDB ID: 6W37), and ORF8 (PDB ID: 7JTL)] were obtained from the RCSB Database (https://www.rcsb.org/) (Rose et al., 2017). The remaining coronavirus structures [membrane protein (M), NSP4, and ORF6] that are not available in the Protein Data Bank (PDB) were modeled using Contact-guided Iterative Threading ASSEmbly Refinement (C-I-TASSER) (https://zhanggroup.org/C-I-TASSER/), an extended version of I-TASSER for high-accuracy protein structure and function predictions (Zheng et al., 2021). The modeled structures of M protein, NSP4, and ORF6 were verified using the Structural Analysis Verification Server (SAVES) program (https://saves.mbi.ucla.edu/) and protein structure analysis (ProSA)-web program (https://prosa.services.came.sbg.ac.at/prosa.php?pdbCode=1zyb) (Wiederstein and Sippl, 2007). The structural quality (PROCHECK) (Laskowski et al., 1993), non-bonded interactions (ERRAT) (Colovos and Yeates, 1993), and energy profile or Z-score (ProSA using molecular mechanics force field) were thoroughly checked and calibrated. The Z-scores obtained for the M protein, NSP4, and ORF6 were close to the related native conformations regarding their residue length, as determined by X-ray and NMR structures available in the RCSB database. Additionally, nearly 83%–88% of residues in these respective modeled structures were found to occupy favorable allowed regions in the Ramachandran plot (Hooft et al., 1997). It indicates that the modeled structures are high quality and likely biologically relevant.

The structurally similar human proteins of SARS-CoV-2 were obtained by submitting the PDB IDs of viral protein structures (either known or C-I-TASSER-generated) to the Dali webserver (http://ekhidna2.biocenter.helsinki.fi/dali/) (Holm, 2022). In this study, only the Homo sapiens proteins with structural similarity score (Z-score) ≥ 2 were considered structurally similar proteins of the respective SARS-CoV-2 protein and are referred to as “hCoV-2” proteins. The PDB codes of hCoV-2 obtained from Dali were mapped to their corresponding UniProt ID and gene names by the UniProt Retrieve/ID mapping tool (https://www.uniprot.org/id-mapping) (Pundir et al., 2016).

The Dali database may contain multiple structures of the same proteins in the PDB, which could result in repetitive interaction predictions. This issue can arise with certain SARS-CoV-2 proteins with multiple PDB structures, leading to repeating similar predictions. Hence, duplicate PDB IDs were removed before converting them to UniProt accession ID and gene names for target prediction. In addition, multiple hCoV-2 proteins (human proteins structurally similar to SARS-CoV-2) can have common cellular partners. Among these, unique pairs of interactions between human UniProt accession and SARS-CoV-2 proteins were considered.

Human endogenous protein interactors of SARS-CoV-2 were identified by obtaining the target protein interactors of hCoV-2 during various cellular processes. These cellular partners of hCoV-2 were obtained from HPRD (https://www.hprd.org/), BioGRID (https://thebiogrid.org/), HIPPIE (http://cbdm-01.zdv.uni-mainz.de/∼mschaefer/hippie/), and MINT (https://mint.bio.uniroma2.it/) databases and are referred to as T_hCoV-2 (target protein interactors of hCOV-2) (Baolin et al., 2007; Licata et al., 2012; Alanis-Lobato et al., 2017; Oughtred et al., 2019). These datasets are from literature-curated interactions established through in vitro and/or in vivo methods among human proteins. These cellular proteins, known to interact with hCoV-2 proteins, are presumed to interact with SARS-CoV-2 proteins due to their structural similarities with hCoV-2 proteins.

The target protein interactors of hCoV-2, referred to as T_hCoV-2, of each of the SARS-CoV-2 proteins, were merged into one dataset corresponding to their HUGO Gene Nomenclature Committee (HGNC) symbols (https://www.genenames.org/) by eliminating the duplicate gene IDs to avoid redundancy. Pre-processed data corresponding to COVID-19 mRNA expression profiling were obtained from the NCBI-GEO (https://www.ncbi.nlm.nih.gov/geo/) (Jha et al., 2022). The unique target genes of T_hCoV-2 are then compared with the COVID-19 patient gene sample dataset from which only the overlapping genes between targets and infected datasets were considered for further DEG analyses.

To identify DEGs between T-hCoV-2 and COVID-19 samples, a two-sample statistical t-test was used, followed by obtaining their log₂ (fold change) and Benjamin–Hochberg (BH) p − value through the limma R package (Ritchie et al., 2015). Genes with BH p < 0.05 and |log₂ (FC) | > 2 were considered DEGs. DEGS with log₂ (FC) > 2 and log₂ (FC) < −2 were categorized as up- and downregulated, respectively. The DEGs that overlapped between targets and infected groups were referred to as TIG-DEGs.

Pathway and GO term enrichment data for TIG-DEGs were compiled using various libraries, i.e., Kyoto Encyclopedia of genes and genomes (KEGG) and GO-biological process (BP), GO-molecular function (MF), and GO-cellular compartment (CC) terms available within the Enrichr database (https://maayanlab.cloud/Enrichr/) (Kuleshov et al., 2016). All the pathways and GO terms corresponding to p < 0.01 were determined as statistically significant.

The genes corresponding to the enrichment analysis were then merged into one dataset (eliminating duplicate gene names) and subjected to PPI network construction using the STRING v11.5 web-based tool (https://string-db.org/) (Szklarczyk et al., 2021). The construction of this PPI network was based on the highest confidence score >0.9, and it was visualized using Cytoscape v3.9.1 (Shannon et al., 2003). The topological properties of the PPI network were analyzed using the NetworkAnalyzer plugin in Cytoscape. Furthermore, the Molecular Complex Detection (MCODE) plugin was used to identify highly correlated gene clusters/modules within the PPI network. The parameters set in MCODE for cluster detection were as follows: “degree cutoff = 2,” “node score cutoff = 0.2,” “k-score = 2,” “max. depth = 100,” and “cut style = haircut.” The top-scoring genes in the PPI cluster were considered the hub genes.

The RPS3 gene has the most significant interaction partners of the top-scoring PPI cluster genes. Ribosomal protein subunit 3 (RPS3) is a protein-coding gene involved in the viral mRNA translation pathway. Since the viruses exploited the ribosomal proteins and the ribosomal biogenic processes to facilitate their replication, the RPs have been considered effective targets for developing antiviral agents (Dong et al., 2021). The three-dimensional electron microscopy structure of the human ribosomal subunit, PDB ID: 6ZLW (resolution: 2.60 Å), bounded to the SARS-CoV-2 NSP1 protein, was downloaded from the Protein Data Bank (https://www.rcsb.org/) (Rose et al., 2017). The two α-helices of NSP1 (residues: 154–179) of SARS-CoV-2 directly interact with RPS3, RPS2, and the phosphate backbone of rRNA helix 18 (h18) of 40S ribosomal protein. This interaction with RPS3, RPS2, and h18 rigidly anchors NSP1, obstructing the mRNA entry channel. Therefore, for this study, we took a subcomplex comprising RPS3, RPS2, and h18 as the targets for the molecular docking studies.

An in silico high-throughput virtual screening of known therapeutic COVID-19 drugs from the PubChem library was obtained (https://pubchem.ncbi.nlm.nih.gov/#query=covid-19) (Kim et al., 2019). The total search results were filtered thrice for final docking studies. Following are the filters for the selection of final screened compounds: filter 1, database filters (molecular weights should be from 100 to 500 g/mol, rotatable bond count from 0 to 7, h-bond acceptor count from 0 to 10, the polar area from 0 to 150 Å², and XlogP from −6.3 to 5); filter 2, calculation of physicochemical properties (no Lipinski’s rule violations, no lead likeness violations, and high GI absorption) using the SwissADME web server (http://www.swissadme.ch/) (Daina et al., 2017); and filter 3, calculation of pharmacokinetics and toxicological properties (cLogP, solubility, TPSA, drug-likeness, drug score, mutagenic, tumorigenic, irritant, and reproductive effective) using the Osiris Property Explorer (Sander, 2001).

The 3D structures of the final screened compounds were generated using the OpenBabel program (O Boyle et al., 2011). The 2D SDF files were used as input files for OpenBabel for converting them to 3D structures, followed by energy minimization in MMFF94 force field (Halgren, 1996). Finally, all the structures of ligands were converted to the PDBQT format with appropriate rotatable bonds.

The PDBQT file of the subcomplex receptor was prepared using Kollman united atom charges for molecular docking study using AutoDockTools v 1.5.6 (Singh and Kollman, 1984; Morris et al., 2009). As the subcomplex receptor contained two proteins and one nucleic acid, the tentative binding site for the inhibitor is unknown. Hence, a blind docking study using AutoDock Vina (v1.1.2) was performed (Trott and Olson, 2009). Grid point spacing was set at 1 Å, grid point 56 was taken in each direction of the grid box, and the concerned box was centered at the subcomplex receptor. The ligands bind to two probable sites of the subcomplex receptor—one at the rRNA h18 (ligand-binding site-I) and another at the RPS3 protein (II). The best three conformations of each ligand at two binding locations were identified after cluster analysis according to their lowest binding energy. The receptor–ligand hydrogen bond interactions were analyzed using the Swiss PDB viewer (Guex and Peitsch, 1997) and visualized using PyMOL v 2.5 (the PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger, LLC).

We used the Dali Lite webserver to identify human host proteins that exhibit structural similarity to the proteins of SARS-CoV-2. Twenty SARS-CoV-2 proteins were selected for this study (17 proteins from known sources, i.e., available via the RCSB PDB, and 3 proteins predicted or modeled using the C-I-TASSER program) (Table 1). The viral protein structures were then submitted to the Dali webserver to compare 3D structural coordinates of virus and human host proteins. To make the study more comprehensive and avoid redundancy, all the structurally similar host proteins for each viral protein were merged into one dataset, and duplicates were removed. For example, for NSP3 of SARS-CoV-2, the human proteins structurally similar to these three PDB IDs: 6W02, 6W6Y, and 6W9C were merged into one column, and duplicate human PDB IDs were then removed. A similar protocol was followed for the other viral proteins, with more than one PDB ID obtained for this study. It resulted in the identification of 6,116 human proteins (hCoV-2) similar to 20 SARS-CoV-2 proteins (Table 2 and Supplementary Table S1).

Afterward, we identified all possible interacting partners of hCoV-2 proteins of the 20 proteins of SARS-CoV-2 from four different databases (HPRD, MINT, HIPPIE, and BioGRID). The target protein interactors of hCoV-2 from each database corresponding to 165,051 genes were merged into one unique dataset without duplicates. A total of 19,051 unique interactors were identified for these 20 SARS-CoV-2 proteins. These unique interactors are termed T_hCoV-2 proteins (Table 2 and Supplementary Table S2).

The pre-processed SARS-CoV-2 mRNA expression profile (GSE164805) datasets comprised five controls, and 10 COVID-19 (five mild and five severe) samples were compared with the T_hCoV-2 genes. The overlapping genes between targets and mRNA datasets amounted to 8,336 genes (Figure 1A). These genes were then differentially expressed corresponding to p − value < 0.05 and log 2 fold change > 2 . The DEGs overlapping between targets and infection severity groups were considered TIG-DEGs, which amount to 1,120 unique genes. Amongst these genes, 663 and 457 were up- and downregulated, respectively (Figure 1B). The most highly up- and downregulated DEGs were HBD ⁡ log 2 F C = 6.00 and TEX101 ⁡ log 2 F C = − 10.34 , respectively (Figure 1C).

All the 1,120 TIG-DEGs participated in the Gene Ontology and pathway enrichment analyses. All the significant pathways and GO-BP, GO-MF, and GO-CC terms ( p − v a l u e < 0.01 ) are shown in Supplementary Table S3. The most significant pathways and GO-BP, GO-MF, and GO-CC terms were coronavirus disease ( p − value = 8.49 × 10 − 08 ), cellular macromolecule biosynthetic process ( p − value = 1.88 × 10 − 08 ), RNA binding ( p − value = 5.64 × 10 − 05 ), and ribosome ( p − value = 1.18 × 10 − 06 ) (Supplementary Table S3). For a more comprehensive representation, only the top 10 most significant pathways and GO terms are shown in Figure 2.

A total of 606 genes comprising all the pathways and GO terms p − v a l u e < 0.01 ) were inputted into the STRING database. The PPI network constructed using Cytoscape software included 302 nodes, and 688 edges are shown in Figure 3A. A total of 302 genes out of 1,120 TIG-DEGs participated in the PPI network at the selected highest confidence level, i.e., > 0.9. Moreover, Molecular Complex Detection (MCODE) revealed that the most highly significant cluster with a top score = 18 was selected amongst all the 17 identified clusters, consisting 19 nodes and 171 edges (Figure 3B). Essential centrality measures such as node degree, betweenness, closeness, average shortest path length, clustering coefficient, radiality, and topological coefficient of the PPI network are provided in Supplementary Table S4. Based on the centrality measures, the significant hub gene RPS3 (ribosomal protein subunit 3, degree = 25), with the highest number of interacting partners, was identified.

The PubChem library was searched with the keyword “COVID-19,” resulting in 1,709 compounds. The compounds were filtered out through three filtering processes (Martín-Sánchez et al., 2023): database filters (Ali et al., 2020), pharmacokinetic properties, and (Morse et al., 2020) toxicological properties. Through database filtering, only 886 candidates were selected that have a molecular weight between 100 and 500 g/mol, rotatable bond count between 0 and 7, h-bond acceptor count between 0 and 10, polar area between 0 and 150 Å², and XlogP between −6.3 and 5. Government organizations, research and development, journal publishers, and NIH initiatives were selected as the data source category. The SwissADME webserver was accessed to calculate pharmacokinetic parameters. After the pharmacokinetic filtration process, 42 compounds that show zero violations of Lipinski’s rule, zero lead likeness violations, and high gastrointestinal (GI) were selected. To calculate toxicological parameters, the OSIRIS Property Explorer was used to predict the toxicity risks and drug scores of the 42 compounds. The program evaluated the risks of different side effects, such as mutagenic, irritant, tumorigenic, reproductive effects, and drug-related properties. Furthermore, the overall drug score was calculated by summing up various parameters such as cLogP, solubility (logS), molecular weight, drug-likeness, and toxicity risk. The final filtration process resulted in 28 compounds showing no toxicity risks and drug score values between 0.40 and 1.00, making it a good library for further docking studies (Table 3).

In this study, blind docking was performed since the subcomplex receptor of ribosomal protein has no ligand-binding site (Figures 4A, B). The rotatable bonds of the 28 ligands were kept flexible, while the subcomplex macromolecule was adopted as a rigid structure. The analysis of the molecular docking results shows that most of the ligands screened bind into two active sites: (I) the center of RPS3 protein and (II) the head of rRNA helix 18 (Figure 4C). The top three binding poses with the best binding energies at each ligand-binding site are given in Table 4. The best binding energies of compound baicalein, vadadustat, and estetrol are observed at the head of the rRNA helix 18 site with binding energies −8.24, −8.01, and −7.98 kcal/mol, respectively. Similarly, baicalein, estetrol, and 3-{[(3-methyl-1,2,4-oxadiazol-5-yl)methyl]carbamoyl}benzoic acid show the best binding energies at the center of RPS3 protein with binding energies −7.54, −7.33, and −7.24 kcal/mol, respectively. Baicalein has been shown to have the highest binding pose for both ligand-binding sites. All the compounds form hydrogen bonds with the rRNA helix18 site that contains the following nucleotides: U607, C608, U609, G610, G611, G613, U630, A629, and U631, whereas the RPS3-binding site contains the following residues: Ala30, Gly33, Arg54, Arg94, Ala100, Gln101, Glu103, and Ser104 (Table 4). The interactions of baicalein, estetrol, and vadadustat (rRNA h18) are shown in Figure 5A, whereas the interactions of baicalein, estetrol, and 3-{[(3-methyl-1,2,4-oxadiazol-5-yl)methyl]carbamoyl}benzoic acid (RPS3 protein) are shown in Figure 5B. In Figure 5A, (i) baicalein binds in the region of rRNA h18 through h-bond interactions with U607, G609, G610, U630, and U631; (ii) estetrol shows h-bond interactions with U607, U609, G611, G613, and A629; and (iii) vadadustat forms h-bond interactions with U608, G610, G611, and U630. Meanwhile, in Figure 5B, (iv) baicalein binds in the region of RPS3 protein through h-bond interactions with residues Ala30, Gly33, Arg54, Arg94, and Ala100; (v) estetrol forms h-bond interactions with residues Arg94, Ala100, and Gln101; and (vi) 3-{[(3-methyl-1,2,4-oxadiazol-5-yl)methyl]carbamoyl}benzoic acid forms h-bond with residues Arg54, Arg94, Gln101, Glu103, and Ser104.