Deciphering SARS CoV-2-associated pathways from RNA sequencing data of COVID-19-infected A549 cells and potential therapeutics using in silico methods.

BACKGROUND: Coronavirus (CoV) disease (COVID-19) identified in Wuhan, China, in 2019, is mainly characterized by atypical pneumonia and severe acute respiratory syndrome (SARS) and is caused by SARS CoV-2, which belongs to the Coronaviridae family. Determining the underlying disease mechanisms is central to the identification and development of COVID-19-specific drugs for effective treatment and prevention of human-to-human transmission, disease complications, and deaths.
METHODS: Here, next-generation RNA sequencing (RNA Seq) data were obtained using Illumina Next Seq 500 from SARS CoV-infected A549 cells and mock-treated A549 cells from the Gene Expression Omnibus (GEO) (GSE147507), and quality control (QC) was assessed before RNA Seq analysis using CLC Genomics Workbench 20.0. Differentially expressed genes (DEGs) were imported into BioJupies to decipher COVID-19 induced signaling pathways and small molecules derived from chemical synthesis or natural sources to mimic or reverse COVID -19 specific gene signatures. In addition, iPathwayGuide was used to identify COVID-19-specific signaling pathways, as well as drugs and natural products with anti-COVID-19 potential.
RESULTS: Here, we identified the potential activation of upstream regulators such as signal transducer and activator of transcription 2 (STAT2), interferon regulatory factor 9 (IRF9), and interferon beta (IFNβ), interleukin-1 beta (IL-1β), and interferon regulatory factor 3 (IRF3). COVID-19 infection activated key infectious disease-specific immune-related signaling pathways such as influenza A, viral protein interaction with cytokine and cytokine receptors, measles, Epstein-Barr virus infection, and IL-17 signaling pathway. Besides, we identified drugs such as prednisolone, methylprednisolone, diclofenac, compound JQ1, and natural products such as Withaferin-A and JinFuKang as candidates for further experimental validation of COVID-19 therapy.
CONCLUSIONS: In conclusion, we have used the in silico next-generation knowledge discovery (NGKD) methods to discover COVID-19-associated pathways and specific therapeutics that have the potential to ameliorate the disease pathologies associated with COVID-19.

1. Introduction

Coronavirus disease 2019 (COVID-19) is caused by a type of coronavirus (CoV), severe acute respiratory syndrome (SARS) virus 2 (SARS CoV-2). COVID-19 is characterized by symptoms ranging from a mild cold to more severe illnesses, such as SARS, sudden stroke, gastrointestinal complications, and multiple organ failure, even leading to death in some patients.[ Coronaviruses belong to the Coronaviridae family, and the presence of viral spike proteins in the virus gives it a halo or corona-like appearance under the electron microscope (Fig. 1A). A novel coronavirus (nCoV) was discovered in Wuhan, China in 2019 as the cause of a human respiratory outbreak that resulted in severe atypical pneumonia.[ and is the source of the current global pandemic affecting all levels of society.[

Figure 1.

Structure of SARS-CoV 2. (A) The structure depicted based on electron microscopic observations of coronavirus showing the surface protein particles S, N, M, and E and shows a corona-like shape. (B) The host tissues expressing ACE2 receptors. (C) The mechanism of entry of SARS CoV2 into the host cells (this figure was created using the graphic tools offered by BioRender.com with an academic license).

The World Health Organization (WHO) has renamed this nCoV as SARS-CoV 2, which is the causative agent of COVID-19.[ COVID-19 is highly transmissible and pathogenic compared with other viral infections, and the exact mortality rate has yet to be determined because the pandemic is not yet under control in several countries, resulting in unprecedented protective measures, partial or complete lockdowns, travel restrictions, etc.[ As of March 7, 2022, COVID-19 had already infected more than 446 million people in 195 countries and territories around the world and killed approximately 6 million people, according to data from the Johns Hopkins Coronavirus Dashboard.[ However, the exact mortality rate will not be calculated or determined until the COVID-19 epidemic reaches a plateau. The United States of America and WHO have declared the SARS-CoV-2 outbreak a public health emergency because it is more contagious than the severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV).[ SARS-CoV-2 possesses a nucleocapsid with a positive-sense RNA genome. Host cells express SARS-CoV-2 nucleoproteins and the nucleocapsid protein (N protein), which is the most abundant, highly immunogenic protein, and is required for CoV RNA synthesis. The N protein is a structural protein that binds to the CoV RNA genome and forms a capsid around viral RNA. However, the spike protein (S protein) is critical for binding between SARS CoV-2 and angiotensin-converting enzyme 2 (ACE2) surface receptors on host cells (Fig. 1B), thus facilitating coronavirus entry into host cells.,[ respectively (Fig. 1C).

Although COVID-19 vaccines are currently available as preventive measures and many are in the research and development phase,[ deciphering the underlying pathological mechanisms is central to identifying and developing COVID-19 specific drugs to effectively treat and prevent human-to-human transmission, COVID-19 complications, and deaths. In silico methodologies can be successfully used to identify potential drugs and natural products based on high-dimensional RNA-seq datasets derived from various disease pathologies.[ We have recently shown that the next-generation knowledge discovery (NGKD) platforms can effectively be used to uncover the gene signatures regulated by COVID-19 and the potential therapeutics using RNAseq datasets derived from normal human primary bronchial epithelial (NHBE) cells.[ However, in the present study, the raw RNA Seq reads (single-end) (FASTQ files) in quadruplicate obtained from SARS CoV-infected A549 cells and mock-treated A549 cells using Illumina Next Seq 500 were obtained from the Gene Expression Omnibus (GEO) (accession number: GSE147507) and quality control (QC) was evaluated before RNA Seq analysis using CLC Genomics Workbench 20.0 (Qiagen, USA). After the initial QC, the RNA Seq reads were imported into the CLC Genomics Workbench 20.0 (Qiagen, USA) before RNA Seq analysis and evaluated using NGKD platforms such as BioJupies[ and iPathwayGuide (Advaita Bioinformatics, USA) to decipher the disease-specific molecular signatures and a series of small molecules derived from either synthetic or natural sources to mimic or reverse the COVID-19 gene signatures.

2. Materials and Methods

2.1. Ethical statement

Animal models and human subjects were not used in this study. This study was performed using RNA-seq datasets from next-generation sequencing experiments with A549 cells. The raw data were obtained from the Gene Expression Omnibus (GEO), as indicated in the Data Source section below. Therefore, they were exempt from institutional review board (IRB) approval.[

2.2. Next-Generation Sequencing (NGS) data source

Raw RNA Seq reads (single-end) (FASTQ format) in quadruplicate obtained with Illumina Next Seq 500 from A549 cells infected with SARS CoV-2 and mock-treated A549 cells were obtained from the Gene Expression Omnibus (GEO) (accession number: GSE147507)[ and were used for subsequent downstream analysis with high-throughput NGKD platforms.

2.3. COVID-19 RNA Seq data from A549 cells – quality control

Raw RNA Seq reads (single-end) in quadruplicates (FASTQ files) derived from SARS CoV-infected A549 cells and mock-treated A549 cells using Illumina Next Seq 500 were derived from the GEO, and quality control (QC) was evaluated using CLC Genomics Workbench 20.0 (Qiagen)[ to obtain the differentially expressed genes (DEGs) before RNA Seq analysis.

2.4. COVID-19 RNA Seq data from A549 cells – differential gene and transcript expression analysis

RNA Seq reads were imported into CLC Genomics Workbench 20.0 (Qiagen) after the QC step. The RNA Seq Analysis Tool in the Biomedical Genomics Analysis plugin of the CLC Genomics Workbench was used to extract all annotated transcripts using both Homo sapiens (hg38) _genes (Gene track) and Homo sapiens (hg38) _mRNA (mRNA track) and mapped to the human reference genome (GRCh38). A gene expression track (GE) was generated for A549 cells infected with SARS CoV-2 and corresponding mock reads (test vs. control). In addition, the differential expression tool was used in the two groups in the CLC Genomics Workbench to perform a statistical test for differential expression for a set of expression tracks (test vs. control). Here, a multifactorial statistic based on a negative binomial generalized linear model (GLM) is used, and the differential expression in the two groups tool deals with one factor and two groups. In this analysis, “Total Exon Read” values were used for GE. Differentially expressed genes (DEGs) were generated for the test compared to the corresponding control and used for further downstream analysis using the NGKD platform.

2.5. BioJupies analysis of RNASeq data

BioJupies was used to analyze the DEGs generated using the CLC Genomics Workbench to identify novel signaling pathways, disease-specific gene networks, and a range of drugs and small molecules derived from natural sources to mimic or reverse disease-specific gene signatures.[ In Biojupies, RNASeq datasets were compressed into an HDF5 data package and uploaded to Google Cloud. Raw data were normalized to log10 counts per million (log CPM) and differentially expressed genes between the control and experimental groups were determined using the R package limma.[ The principal component analysis function in the Python module of sklearn was used to transform log CPM based on the Z-score method to generate the PCA plot, and Clustergrammer[ was used to generate interactive heat maps, and the DEGs were then sent to Enrichr.[ In the volcano plot, DEGs were plotted on the x-axis and P values were corrected using the Benjamini-Hochberg method, transformed (−log10), and plotted on the y-axis.[ However, average gene expression is shown on the x-axis in the MA plot, and P values were corrected, transformed (−log10), and plotted on the y-axis using the Benjamini-Hochberg method.[ Gene Ontology (GO) and pathway enrichment analyses were performed with both upregulated and downregulated genes in Enrichr. Significant GO terms and pathways (KEGG, WiKiPathways, and Reactome) were calculated using a cut-off value of < 0.1 after applying the Benjamini–Hochberg correction.[

2.6. L1000CDS2 and L1000FWD analyses

The L1000CDS2 analysis was performed by submitting the best 2000 DEGs to the L1000CDS2 signature search API.[ Similarly, the L1000FWD analysis was performed by submitting the top 2000 DEGs to the L1000FWD signature search API.[

2.7. In silico analysis of RNASeq expression data using iPathwayGuide

The impact analysis method (IAM)[ in the iPathwayGuide was used to determine the differentially regulated signaling pathways, gene ontologies, and upstream drugs or natural products, as previously described.[ Briefly, the pathway score was calculated based on the p-value obtained using Fisher’s method. The P value was corrected using multiple testing corrections for the false discovery rate (FDR) and Bonferroni correction.[ The FDR has significant power, but it controls only the family-based false-positive rate.[ Pathways and gene interactions with DEGs were generated using the KEGG database.[ For each Gene Ontology term (GO), the number of DEGs annotated with the term was compared with the randomly expected DEGs.[ iPathwayGuide used an overrepresentation approach to calculate the statistical significance of observing at least the specified number of DEGs.[ The hypergeometric distribution was used to calculate P values in the iPathwayGuide analysis and corrected for multiple comparisons using FDR and Bonferroni.[

2.8. Prediction of upstream drugs or natural products with iPathwayGuide

The prediction of upstream chemicals, drugs, and toxicants (CDTs) was based on two types of information: (i) the enrichment of DEGs from experiments and (ii) a network of interactions from the Advaita Knowledge Base (AKB v2012).[

2.9. Upstream CDTs predicted to be present (or overabundant)

The research hypothesis refers to the presence of CDT. This hypothesis is useful for investigating whether a given phenotype is influenced by the presence of a particular chemical, drug, or toxicant.[ For each CDT u, the number of consistent DE genes after u, DTA(u), is compared to the number of measured target genes expected to be both consistent and DE. iPathwayGuide uses an over-representation approach to calculate the statistical significance of observing at least a given number of consistent DE genes. The P value Ppres was calculated using a hypergeometric distribution.[ The analysis uses Fisher’s standard method to combine the P values into a test statistic.[

2.10. Upstream CDTs predicted to be absent (or insufficient)

In parallel with the upstream CDTs predicted to be present, Pabs and Pz were used to predict upstream CDTs that were absent. This hypothesis is important when investigating whether a given phenotype is affected by the absence of a particular chemical necessary for the proper functioning of the organism or cell. The research hypothesis states that the upstream CDT is insufficient under the conditions under study. For each upstream CDT u, the number of consistent DE genes downstream of u, DTI(u), was compared to the number of measured target genes expected to be both consistent and DE by chance. Using Fisher’s method, the analysis combines Pabs and Pz, with Pz considered only for significantly negative z-scores (z ≤ −2).[

2.11. Swiss target prediction of potential anti- COVID-19 compounds

The isomeric simplified molecular-input line-entry system (SMILES) codes of prednisolone and withaferin-A were used in the SwissTarget Prediction tool to identify protein targets.[ Ligand-based target prediction for both prednisolone and withaferin-A was performed as previously described.[

2.12. The Open Targets Platform analysis of anti-COVID-19 compounds

The Open Targets Platform web tool was used to uncover the molecular targets of prednisolone and withaferin-A associated with COVID-19 disease pathology (date accessed: January 12, 2020)[ The Open Targets Platform uses scientific evidence to assess and rank associations between targets and disease and to help prioritize targets.[ The query list of approximately 100 candidate molecular targets of prednisolone and withaferin-A was used to discover protein targets significantly (P < .05) associated with COVID-19.

3. Results

The present study was done using RNA-seq datasets obtained from next-generation sequencing experiments with mock-treated and SARS CoV-2 infected A549 cells.[ The raw RNA Seq reads (Single-End) (FASTQ files) in quadruplicates derived using Illumina Next Seq 500 from SARS CoV-infectedA549 cells, and mock-treated A549 cells were obtained from the Gene Expression Omnibus (GEO) (GSE147507), and quality control (QC) was evaluated before RNA Seq analysis using the CLC Genomics Workbench 20.0 (Qiagen). The DEGs were further analyzed using BioJupies and iPathwayGuide (Advaita Bioinformatics, USA) to decipher disease-specific signatures and an array of drugs and small molecules derived from natural sources to mimic or reverse disease-specific gene signatures.

The global patterns in the high-dimensional RNA-seq datasets were uncovered using PCA analysis (Fig. 2A). The Clustergrammer web tool was used to generate interactive heatmaps for visualization and in-depth analysis of DEGs derived from high-dimensional RNASeq data of SARS CoV-infected A549 cells and mock-treated A549 cells (Fig. 2A–C). A volcano plot was generated using transformed gene fold changes using log2 and is shown on the x-axis (Fig. 2D). The MA plot was based on the average gene expression, which was calculated using the mean of the normalized gene expression values and is shown on the x-axis (Fig. 2E).

Figure 2.

(A) Principal component analysis (PCA) was applied to identify global patterns in high-dimensional RNASeq datasets. (B) The heatmaps were generated top 250 DEGs and (C) the top 50 DEGs using Clustergrammer web tool for visualizing and analyzing high-dimensional RNASeq data. (D) Volcano plot was generated using transformed gene fold changes using log2 and displayed on the x-axis. (E) MA plot was based on average gene expression which was calculated using mean of the normalized gene expression values and displayed on the x-axis.

The bar chart (Fig. 3A) shows the top small molecules identified by the L1000CDS2 query. The left panel displays small molecules such as calyculin A, emetine hydrochloride, narliclasine, NVP-TAE684, wiskostatin, NCGC00185684-02, and amsacrine, which mimic the observed gene expression signature, while the right panel displays small molecules such as trichostatin A, vorinostat, afatinib, DL-PDMP, withaferin-A, IMD 0354, and 2-[(chloroacetyl) (4-fluorophenyl] amino-N-cyclohexyl-2 pyridine 3, which reverse it. In addition, natural products and drugs with opposite (Table 1) and similar molecular signatures (Table 2) based on the L1000FWD tool, which contains gene signatures from an array of human cell lines administered with more than 20,000 drugs and natural products. Withaferin-A, an active ingredient of the medicinal plant (Fig. 3A), Withania somnifera was found to reverse the COVID-19 induced molecular signatures in both L1000CDS2 and L1000FWD analyses along with other small molecule drugs.

Figure 3.

(A) The bar chart displaying the top small molecules identified by the L1000CDS2 query. The left panel displays the small molecules which mimic the observed gene expression signature, while the right panel displays the small molecules which reverse it. (B) Bar graphs show the synthetic drugs and natural compounds with similar (mimic) and opposite (reverse) molecular signatures based on iPathwayGuide analysis.

Table 1

Natural products and drugs with opposite and similar molecular signatures based on L1000FWD web-based tool.

Opposite molecular signatures
Signature ID	Drugs or natural products	Similarity score	P value	q value	Z-score	Combined score
CPC019_VCAP_24H:BRD-K50234570-001-06-6:10	EMF-bca1-16	−0.0569	2.03E-10	8.67E-07	1.67	−16.15
ERG005_VCAP_6H:BRD-K88378636-001-02-8:20	Withaferin-a	−0.0544	1.54E-09	2.74E-06	1.65	−14.57
CPC006_HCC515_24H:BRD-A28105619-001-01-3:10	Cucurbitacin-i	−0.0531	2.51E-09	3.98E-06	1.81	−15.59
CPC006_HCC515_6H:BRD-K16406336-311-01-2:10	Methylene-blue	−0.0544	6.11E-09	8.44E-06	1.77	−14.52
CPC016_MCF7_24H:BRD-K08547377-003-03-2:10	Irinotecan	−0.0506	1.88E-08	2.12E-05	1.7	−13.16
CPC001_VCAP_24H:BRD-K12516989-001-01-9:10	Zaprinast	−0.0442	6.06E-08	5.64E-05	1.93	−13.91
CPC016_NPC_24H:BRD-A22783572-065-01-3:10	Vinblastine	−0.0493	1.28E-07	1.08E-04	1.69	−11.63
CPC004_PC3_6H:BRD-A69815203-001-05-0:10	Cyclosporin-a	−0.0455	1.90E-07	1.51E-04	1.84	−12.37
CPC008_PC3_6H:BRD-K66037923-001-04-4:10	BRD-K66037923	−0.048	1.96E-07	1.53E-04	1.76	−11.79
MUC.CP003_MCF7_24H:BRD-K02407574-001-04-8:0.3704	Parbendazole	−0.0468	2.45E-07	1.79E-04	1.63	−10.75
Similar Molecular Signatures
CPC013_SKB_24H:BRD-K61175124-001-01-0:10	BRD-K61175124	0.0556	2.63E-13	1.12E-08	−1.83	23.09
CPC016_SKB_24H:BRD-A06352508-001-02-9:10	SB-218078	0.0544	1.17E-12	1.67E-08	−1.87	22.28
CPC006_HT29_24H:BRD-A67788537-001-01-7:120	Salermide	0.0493	1.78E-12	1.91E-08	−1.85	21.75
CPC002_PC3_6H:BRD-A22684332-003-03-1:10	Procaterol	0.0582	3.07E-12	2.63E-08	−1.64	18.87
CPC007_HT29_24H:BRD-A09719808-001-02-3:10	BRD-A09719808	0.0506	8.82E-12	6.29E-08	−1.81	20.02
CPC019_VCAP_6H:BRD-K23282736-001-01-1:10	BRD-K23282736	0.0594	1.12E-11	6.83E-08	−1.78	19.53
CPC007_HT29_6H:BRD-A69470004-019-04-0:10	BRD-A69470004	0.0556	5.64E-11	3.02E-07	−1.69	17.3
CPC013_SKB_24H:BRD-K74623475-001-02-7:10	BRD-K74623475	0.048	1.67E-10	7.94E-07	−1.86	18.18
CPC006_A673_6H:BRD-K84924563-001-01-2:40	BRD-K84924563	0.0531	3.62E-10	1.29E-06	−1.68	15.85
CPC013_SKB_24H:BRD-K16541732-001-01-3:10	BRD-K16541732	0.0493	7.61E-10	1.92E-06	−1.81	16.48

Table 2

Top pathways and their associated P values are stated in the table.

Pathway name	Pathway Id	P value	P value (FDR)	P value (Bonferroni)
Influenza A	05164	7.635e-7	6.389e-5	9.849e-5
Viral protein interaction with cytokine and cytokine receptor	04061	1.533e-6	6.389e-5	1.978e-4
Measles	05162	1.665e-6	6.389e-5	2.148e-4
Epstein-Barr virus infection	05169	2.019e-6	6.389e-5	2.605e-4
IL-17 signaling pathway	04657	2.476e-6	6.389e-5	3.194e-4

The P value corresponding to the pathway was computed using only over-representation analysis.

The GO enrichment analysis for the biological processes, molecular function, and cellular components was generated using Enrichr (Fig. 4). The x-axis indicates the −log10(P value) for each term, and significant terms enriched in each GO category are highlighted in bold. Similarly, Figure 5 shows the results of pathway enrichment analysis using Enrichr. The x-axis indicates the −log10(P value) for each term, and the significantly enriched pathways (KEGG, Wiki pathways, and Reactome) are highlighted in bold.

Figure 4.

Gene Ontology Enrichment Analysis. The bar charts display the results of the Gene Ontology enrichment analysis generated using Enrichr. The x-axis indicates the −log10(P value) for each term. Significant terms are highlighted in bold.

Figure 5.

Pathway Enrichment Analysis. The bar charts displaying the results of the pathway enrichment analysis generated using Enrichr. The x-axis indicates the −log10(P-value) for each term. Significant terms are highlighted in bold.

In this experiment, 141 DEGs were identified from a total of 9665 DEGs obtained from BioJupies analysis of the RNASeq reads of SARS CoV-infected A549 cells and mock-treated A549 cells based on a P value cut-off (.05) and a fold change cut-off of 1.5. The DEGs were analyzed in the context of pathways obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, gene ontologies from the Gene Ontology Consortium database, and iPathwayGuide analysis, which further showed that 34 pathways were significantly affected in the SARS CoV2 infected A549 cells compared to the mock-treated A549 cells. In addition, 557 Gene Ontology (GO) terms, 224 gene upstream regulators, 451 chemical upstream regulators, and 31 diseases were found to be significantly (P < .05) enriched before the correction for multiple comparisons.

The top five upstream regulators identified after the Bonferroni correction for signal transducer and activator of transcription 2 (STAT2), interferon regulatory factor 9 (IRF9), interferon-beta (IFNβ), interleukin-1-beta (IL-1β), and interferon regulatory factor 3 (IRF3) were predicted to be activated (Table 7).

Table 7

Top upstream regulators after Bonferroni Correction are given in the table.

Upstream Regulator (u)	DTA(u)	DT(u)	P value	P value (FDR)	P value (Bonferroni)
STAT2	11	11	1.655e-14	6.848e-12	8.092e-12
IRF9	10	10	2.801e-14	6.848e-12	1.370e-11
IFNB1	6	7	1.526e-6	2.488e-4	7.464e-4
IL1B	7	8	1.188e-4	.014	.058
IRF3	3	3	1.464e-4	.014	.072

COVID-19 infection activates key infectious disease-specific immune-related signaling pathways such as influenza A, viral protein interaction with cytokine and cytokine receptors, measles, Epstein-Barr virus infection, and IL-17 signaling pathway (Table 3). Likewise, significantly enriched Gene Ontology (GO) terms such as biological, molecular, and cellular processes based on the false discovery rate (q value) were identified using iPathwayGuide. The top identified biological processes were innate immune response, response to external biotic stimulus, response to other organisms, response to biotic stimulus, and defense response to other organisms, including chemokine receptor binding, chemokine activity, CXCR chemokine receptor binding, receptor-ligand activity, signaling receptor activator activity. The top cellular components identified included blood microparticles, fibrinogen complexes, nuclear outer membranes, extracellular spaces, and extracellular regions for each pruning type (Tables 4–6).

Table 3

Top identified biological processes. The top-scoring biological process, molecular function, and cellular component for each pruning type are described below in the table.

Pruning type: None				Pruning type: High-specificity		Pruning type: Smallest common denominator
GO Term	P value	P value (FDR)	P value (Bonferroni)	GO Term	P value	GO Term	P value
Biological processes
Innate immune response	1.000e-24	4.427e-22	4.427e-22	Type I interferon signaling pathway	2.380e-12	Type I interferon signaling pathway	9.961e-14
Response to external biotic stimulus	1.000e-24	6.549e-22	1.965e-21	Defense response to virus	8.301e-12	Defense response to virus	4.012e-13
Response to other organism	1.000e-24	6.549e-22	1.965e-21	Negative regulation of viral genome replication	2.490e-8	Interferon-gamma-mediated signaling pathway	1.383e-10
Response to biotic stimulus	3.000e-24	2.075e-21	8.301e-21	Interferon-gamma-mediated signaling pathway	8.301e-6	Negative regulation of viral genome replication	1.868e-8
Defense response to other organism	4.900e-23	2.712e-20	1.356e-19	Innate immune response
Molecular functions
Chemokine receptor binding	1.300e-9	3.897e-7	5.629e-7	Chemokine activity	7.794e-7	Chemokine receptor binding	5.629e-7
Chemokine activity	1.800e-9	3.897e-7	7.794e-7	CXCR chemokine receptor binding	4.546e-5	2′-5′-Oligoadenylate synthetase activity	.040
CXCR chemokine receptor binding	1.300e-8	1.876e-6	5.629e-6	2′-5′-Oligoadenylate synthetase activity	.030	Double-stranded RNA binding	.040
Receptor ligand activity	1.700e-6	1.840e-4	7.361e-4	Double-stranded RNA binding	.030	Protein ADP-ribosylase activity	.086
Signaling receptor activator activity	2.400e-6	2.078e-4	.001	Protein ADP-ribosylase activity	.068	ISG15 Transferase activity	.109
Cellular components
Blood microparticle	9.800e-8	3.763e-5	3.763e-5	Blood microparticle	3.763e-5	Blood microparticle	3.763e-5
Fibrinogen complex	.002	.358	.768	Fibrinogen complex	.358	Fibrinogen complex	.358
Nuclear outer membrane	.003	.358	1.000	Nuclear outer membrane	.358	Nuclear outer membrane	.358
Extracellular space	.006	.425	1.000	Specific granule lumen	.553	Extracellular region	.553
Extracellular region	.008	.425	1.000	Costamere	.553	Specific granule lumen	.553

Table 4

Top upstream regulators after Bonferroni Correction is given in the table.

Upstream Regulator (u)	DTA(u)	DT(u)	P value	P value (FDR)	P value (Bonferroni)
STAT2	11	11	1.655e-14	6.848e-12	8.092e-12
IRF9	10	10	2.801e-14	6.848e-12	1.370e-11
IFNB1	6	7	1.526e-6	2.488e-4	7.464e-4
IL1B	7	8	1.188e-4	.014	.058
IRF3	3	3	1.464e-4	.014	.072

Table 6

Top identified cellular components. Only the top scoring cellular component for each pruning type is described below the table.

Pruning type: None				Pruning type: High-specificity		Pruning type: Smallest common denominator
GO Term	P value	P value (FDR)	P value (Bonferroni)	GO Term	P value	GO Term	P value
Blood microparticle	9.800e-8	3.763e-5	3.763e-5	Blood microparticle	3.763e-5	Blood microparticle	3.763e-5
Fibrinogen complex	.002	.358	.768	Fibrinogen complex	.358	Fibrinogen complex	.358
Nuclear outer membrane	.003	.358	1.000	Nuclear outer membrane	.358	Nuclear outer membrane	.358
Extracellular space	.006	.425	1.000	Specific granule lumen	.553	Extracellular region	.553
Extracellular region	.008	.425	1.000	Costamere	.553	Specific granule lumen	.553

The upstream regulator drugs obtained either based on chemical synthesis or natural sources with opposite molecular signatures were also identified based on iPathwayGuide Analysis The drugs that can significantly reverse the molecular impact of COVID-19 infection are Methyl Prednisolone, Prednisolone, Gold Sodium Thiomalate, Tofacitinib, Diclofenac, JQ1 Compound, Azathioprine, etc. (Fig. 3B). The upstream regulator drugs and natural products with opposite molecular signatures identified using iPathwayGuide sorted based on the Z score are listed in Supplementary Table 1, http://links.lww.com/MD/G901.

In the present study, SwissTargetPrediction was performed for prednisolone and withaferin A, using the canonical SMILES code. The Open Targets Platform was applied to uncover Withaferin-A molecular targets associated with COVID-19 disease pathology. Scientific evidence was used in the Open Targets Platform to assign a score and rank target-disease associations and help target prioritization. Among the molecular targets of prednisolone and withaferin-a, 40 and 36 targets, respectively, were significantly associated with COVID-19 pathology (Table 8).

Table 8

COVID-19 associated targets regulated by prednisolone and withaferin-A.

Disease	Drug or natural product	Number of associated targets	Therapeutic area	All targets
COVID-19	Prednisolone	40	Infectious disease	DPP4 JAK1 NR3C1 JAK2 AR PTGS1 CHRNA4 IL6 PDE10A SLC5A2 FLT3 MAPK14 TYK2 OPRM1 KIT PPARG KDR ABL1 NR3C2 ESR2 CNR1 ADORA3 MPEG1 PGR ADAM17 CD38 MTOR MPO EGFR SLC6A3 MAPK1 ALK NOS2 SLC5A1 BRD4 MAPK3 ADK LCK RORA SHBG
COVID-19	Withaferin-A	36	Infectious disease	NR3C1 PTGS2 AR HMGCR PTGS1 GSK3B F10 PDE4D GSK3A IMPDH1 PDE3A PDE3B PDE10A MAPK14 JAK3 NR3C2 IKBKB ADORA2A PGR REN PARP1 ERBB2 CCR1 MAPK1 ALK HDAC3 PRKCB BRAF IL6ST CXCR3 MAPK8 IARS1 BRD4 BCL2L1 MAPK3 MDM2

4. Discussion

COVID-19 is highly infectious and pathogenic compared to other viral infections, and the exact mortality rate has yet to be determined because the pandemic is not yet under control in several countries.[ Therefore, deciphering the underlying pathologic mechanisms is central to identifying and developing COVID-19-specific drugs to effectively treat and prevent person-to-person transmission, COVID-19 complications, and reduce mortality. COVID-19 is usually characterized by cough, breathing problems, high body temperature, diarrhea, and abdominal discomfort, and in severe cases, it causes atypical pneumonia, SARS, stroke, thrombosis, multiple organ failure, and in some cases, death.[ It was found that approximately 80% of COVID-19 cases had mild or asymptomatic symptoms, with the elderly and those with other comorbid conditions more likely to develop severe symptoms and succumb to the disease.[

Distinguishing COVID-19 from other influenza viruses, SARS, and MERS coronaviruses is essential in the clinical setting to develop effective or efficient treatment strategies for patients.[ Noninfectious diseases such as idiopathic interstitial pneumonia, cryptogenic organizing pneumonia, dermatomyositis, and vasculitis also need to be differentially diagnosed from COVID-19[

The COVID -19 infection of A549 cells activated upstream genes, such as STAT2, IRF9, IFNB1, IL1B, and IRF3. Biological processes such as the type I interferon signaling pathway, defense response to viruses, negative regulation of viral genome replication, and interferon-gamma-mediated signaling pathways were differentially regulated. Molecular functions such as chemokine activity, CXCR chemokine receptor binding, 2′-5′-oligoadenylate synthetase activity, double-stranded RNA binding, and protein ADP-ribosylase activity were enriched in the COVID-infected cells. Cytokines are hormones of the immune system that are important for innate and adaptive host responses, cell growth and differentiation, repair and development, cellular homeostasis, and cell death.[ Cytokines are glycoproteins that are released upon any external stimulus and bind to specific cell surface receptors on the plasma membrane of target cells to elicit their responses.[

The cytokine/chemokine storm seen in moderate to severe cases of COVID -19 is caused by a significant increase in the levels of several circulating cytokines and chemokines, such as interleukin-6 (IL-6), interleukin-8 (IL-8), tumor necrosis factor-alpha (TNF-α), C-X-C motif chemokine ligand 10 (CXCL-10), and interferon-gamma induced protein-10 (IP-10), and contributes to poor prognosis.[ In general, viruses evolve mechanisms to avoid detection and subsequent destruction in the host by remodeling and copying cytokine and cytokine receptor genes.[ Similarly, COVID-19 induced cytokines, cytokine receptors, chemokines, and other specific cytokine receptors and binding proteins to destabilize and alter host cytokine responses and immune networks.[ Here, COVID-19-induced chemokines and cytokines can either enhance or prevent cytokine signaling and significantly alter or attenuate various arms of the host immunity. In addition, cellular processes such as the blood microparticle-fibrinogen complex were activated in COVID-infected A549 cells. The increase in cellular processes, such as blood microparticles, observed in the present study was confirmed by a recent study showing an increase in circulating blood microparticles and activated platelets in COVID-19 patients.[

The COVID-19 pandemic is currently being addressed with vaccines, convalescent plasma, monoclonal antibodies, antiviral drugs such as remdesivir, and preventive measures such as wearing masks, hand hygiene, and social distancing.[ In the present study, withaferin-A was predicted to counteract the molecular signatures triggered by COVID-19. Using NGKD platforms, we recently found that withaferin-A reverses the gene signatures induced by COVID-19 in NHBE cells.[

Analysis of the open-target platform revealed that 36 targets played a role in COVID-19 pathology. Withaferin-A is a constituent of the medicinal plant W. somnifera (Indian ginseng or ashwagandha). Its active constituents include withanolides, saponins, alkaloids, and steroidal lactones. W. somnifera is used in herbal preparations in traditional medicine and has antioxidant, anti-anxiety, anti-inflammatory, antibacterial, and aphrodisiac properties, among others[ Ashwagandha has neuroprotective, cardioprotective, immunomodulatory, and anticancer properties.[ In a recent in silico screening study, ashwagandha was also found to contain natural compounds against COVID-19.[

Traditional Chinese Medicine (TCM) has also been used in the treatment of COVID-19.[ The traditional Chinese herbal formula, JinFuKang, consists of 12 medicinal plants, with each dose containing 10 mL.[ JinFuKang has anticancer properties and numerous medicinal benefits.[ Antiviral remdesivir reduces mortality only very slightly,[ the use of corticosteroids increases the possibility of secondary infections,[ and monoclonal antibody therapies are either expensive or difficult to obtain for COVID-19 therapy. However, oral antiviral drugs such as paxlovid and molnupiravir introduced by Pfizer and Merck, respectively, are authorized by the Food and Drug Administration (FDA), USA, for COVID-19 treatment.[ Nevertheless, it may also be valuable to explore the gene signatures triggered by COVID-19 and its variants in different experimental model systems to identify potential drugs or natural products for COVID-19 therapy.

5. Conclusions

The present study demonstrated the application of RNA sequencing technologies in conjunction with NGKD platforms to decipher specific compounds, either synthetic or derived from natural products, for the potential amelioration of COVID-19. However, further in-depth studies are needed to validate drugs such as prednisolone, methylprednisolone, diclofenac, and JQ1, and natural products such as Withaferin-A and JinFuKang in COVID-19 infection model systems, such as primary human alveolar epithelial cells and human small intestinal organoids (hSIOs)[ to determine mechanisms of action before preclinical and clinical trials for the potential treatment of COVID-19 and related pathologies. In conclusion, this study outlines a valuable method for applying NGKD platforms to discover precise drugs and natural products for the potential treatment of COVID-19-related disease pathologies.

Acknowledgments

The authors gratefully acknowledge the technical and financial support from King Abdulaziz University, DSR, Jeddah. Saudi Arabia

Author contributions

PNP, LAD, LD, SB, and MR designed experiments. PNP and MR conducted experiments. PNP, LAD, LD, SB, and MR analyzed the data. PNP and MR prepared the manuscript. PNP and MR revised the manuscript. All authors contributed to the editing of the manuscript and scientific discussions.

Table 5

Top identified molecular functions. Only the top scoring molecular function for each pruning type is described below the table.

Pruning type: None				Pruning type: High-specificity		Pruning type: Smallest common denominator
GO term	P value	P value (FDR)	P value (Bonferroni)	GO Term	P value	GO Term	P value
chemokine receptor binding	1.300e-9	3.897e-7	5.629e-7	Chemokine activity	7.794e-7	Chemokine receptor binding	5.629e-7
chemokine activity	1.800e-9	3.897e-7	7.794e-7	CXCR chemokine receptor binding	4.546e-5	2′-5′-Oligoadenylate synthetase activity	.040
CXCR chemokine receptor binding	1.300e-8	1.876e-6	5.629e-6	2′-5′-Oligoadenylate synthetase activity	.030	Double-stranded RNA binding	.040
receptor ligand activity	1.700e-6	1.840e-4	7.361e-4	Double-stranded RNA binding	.030	Protein ADP-ribosylase activity	.086
signaling receptor activator activity	2.400e-6	2.078e-4	.001	Protein ADP-ribosylase activity	.068	ISG15 transferase activity	.109