A signaling pathway-driven bioinformatics pipeline for predicting therapeutics against emerging infectious diseases

Included datasets

A search of the Gene Expression Omnibus (GEO) database, hosted at the National Center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov/geo/)²⁰, was performed in mid-2020 to find RNA-sequencing datasets for various viruses including “MERS”, “SARS”, and “coronavirus”. The corresponding sequencing data for four GEO series were retrieved from the Sequence Read Archive (SRA; https://www.ncbi.nlm.nih.gov/sra) at NCBI²¹: GSE122876, GSE56192, GSE147507, GSE139516^22–24. These datasets were generated from cell cultures or patients infected with one of: respiratory syncytial virus (RSV), influenza A virus (IAV), MERS, SARS, and SARS-CoV-2. The SRA files were downloaded and converted to fastq format using version 2.10.1 of the NCBI sratools software package (https://github.com/ncbi/sra-tools).

Differential expression analysis

The Automated Reproducible MOdular Workflow for Preprocessing and Differential Analysis of RNA-seq Data (ARMOR) workflow was used to preprocess and analyze the fastq files against the Ensembl reference transcriptome for Homo sapiens build GRCh38, release 98 [GCA_000001405.15]²⁵. Briefly, this automated snakemake-based workflow performs quality control of the reads with fastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc/), trims the adapters and poor-quality regions with TrimGalore! (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), maps & quantifies reads to the human transcriptome using Salmon²⁶, performs differential expression using edgeR²⁷, identifies differential transcript usage with DRIMseq²⁸, and enriches on Gene Ontology terms²⁹ and Hallmark gene sets³⁰ using the Camera algorithm to adjust for inter-gene correlation³¹. Differential expression was performed by calculating log2 fold-change and the associated p-values from infected samples versus mock-infected samples. The significance threshold for differential gene expression and Gene Ontology enrichment was defined as an false discovery rate (FDR)-corrected p-value < 0.05. Gene Ontology enrichment results from all datasets were then combined to identify shared terms.

Pathway enrichment analysis

Once the differentially expressed gene lists were constructed, the Ensembl identifiers for these human gene lists were mapped to the corresponding NCBI Entrez Gene identifiers using BioMart and Bioconductor^32,33 prior to pathway enrichment using Signaling Pathway Impact Analysis (SPIA)³⁴. Briefly, this pathway analysis combines robust statistics and bootstrapping to identify enriched pathways from lists of genes by incorporating the directionality of expression. Public pathway databases in the Graphite package in R Bioconductor were used by SPIA for enrichment with those surpassing a Bonferroni-corrected p-value < 0.05 being retained³⁵. These databases include: KEGG³⁶, Reactome³⁷, Panther³⁸, NCI³⁹, and BioCarta. All pathways from each dataset were then compared to identify those that were shared among the various datasets, as well as those that were unique only to infection with SARS-CoV-2.

Prediction of relevant drug targets

The output from the SPIA pathway enrichment were then used as input for a custom bioinformatics pipeline to identify existing drugs and biologics that could be repurposed to reduce signs, symptoms, and/or replication of coronaviruses. Specifically, this pipeline iteratively 1) retrieves the genes that participate in each statistically significant signaling pathway, 2) maps the gene identifiers to the corresponding UniProt protein identifier⁴⁰, 3) searches the opentargets resource (www.opentargets.org⁴¹) to identify known drug targets and therapeutics, and 4) generates a table with various attributes of the target and small molecule/biologic treatment. These data were then integrated with the pathway comparisons across the various datasets and analyzed to determine what drugs would affect the pathways dysregulated in the greatest number of viruses and could be used as potential therapeutics for SARS-CoV-2 infection. The results were ranked by how many times a given protein target appeared in the results across the relevant virus taxa. The top ranked small molecules/biologic treatments were analyzed to determine which were predicted to reverse the effects of the viral infection on a given pathway. Manual review of high-ranking hits was then performed to determine the existing indication(s) for each treatment, followed by a literature search to determine which, if any, of the therapeutics identified had previously been used or considered for the treatment of COVID-19. Code for this workflow can be found on GitHub: https://github.com/bpickett/Pathway2Targets.

Differential expression analysis

Each of the algorithmic components of our computational workflow was dependent on first calculating differentially expressed (DE) genes for each of the datasets. We consequently began by using the same computational workflow to generate the DE genes from the raw transcriptomic data for human cells infected with one of: RSV, IAV, MERS, SARS, and SARS-CoV-2. Each set of results was specific to the comparison that was performed (e.g. mock-infected vs. infected) and included both log2 fold-change values and FDR-corrected p-values. We then integrated and compared the differential expression results across all virus comparisons to facilitate downstream comparisons.

We began by identifying genes that were either up- or down-regulated across a large number of studies. Specifically, we calculated the number of DE genes (FDR-corrected p-value < 0.05) that significantly changed during infection with influenza A (4205 genes), RSV (3661 genes), MERS (range across studies: 1222-13006 genes), SARS-CoV (range across studies: 2130-5557 genes), or SARS-CoV-2 (range across studies: 427-6933 genes). Out of the 15 studies evaluated, we observed 11 genes that displayed increased expression across 10-13 study comparisons including: TRIM25, C3, NCOA7, PTAFR, TNFAIP3, EIF2AK2, HELZ2, HBEGF, NFKB2, REL and VEGFC. We found an additional 18 genes that were upregulated across 9 of the 15 comparisons including: IFIT2, EGFR, and FOSL1. In contrast, we identified 14 genes that showed decreased expression across 11-12 comparisons including: EIF4H, APH1A, NME2, SEC61B, TUBB4B, CHCHD2, TUBA1C, TP53I3, H2AFZ, PEBP1, HOXB6, TPM2, CBR1, and SIVA1.

We next analyzed the DE lists to identify significant genes that were shared across infections involving only the Betacoronaviruses including SARS-CoV, MERS-CoV, and SARS-CoV-2. This analysis revealed 12 genes that were upregulated in at least nine of thirteen comparisons including: TRIM25, REL, DNASE1, C3, NCOA7, PTAFR, EIF2AK2, HELZ2, NFKB2, ZNF385C, ZC3H12A, and OR7E122P. In contrast, the 10 genes that were found to be downregulated in ten of the thirteen comparisons involving any of the Betacoronaviruses included: CHCHD2, SMIM3, EIF4H, NME2, CBR1, SPARC, DSTN, CDC123, TIMM17A, and PDAP1.

We then performed a similar analysis to identify the five statistically significant DE genes (FDR-corrected p-value < 0.05) that displayed the highest and lowest average fold-change values during infection across all comparisons of SARS-CoV-2, Betacoronaviruses, or all viruses (Figure 1). This analysis showed that a subset of the genes such as CXCL11, which is induced by interferon and is involved in T-cell signaling, displayed similar fold-change values across multiple comparisons. This finding was somewhat expected since certain genes involved in the host innate immune response are expected to be modulated during virus infection. However, we also observed that genes such as MUC3A, PCSK5, MRC1, and CLEC3B displayed somewhat different average fold-change values across the included comparisons. This observation was also expected given the diverse virus replication requirements and the resulting host intracellular transcriptional response that occurs during infections with human viral pathogens.

Figure 1. Heatmap of human transcriptional biomarkers associated with infection by various virus taxa.

The five genes with the highest and lowest log2 fold-change values for each virus comparison were identified and compared against expression in the other datasets. Squares in red or blue represent genes that are upregulated or downregulated (respectively) in the relevant comparison.

Enrichment of annotated terms

Given the number of DE genes involved in the analysis, manual interpretation of the results would be intractable. We therefore performed an analysis to determine which Gene Ontology (GO) terms and Hallmark gene sets were specifically enriched to better understand which biological processes and molecular functions were represented by the genes in each of the DE gene lists. Overall, we observed a superset of 580 terms that were enriched in any of the included viral infections. The enrichment results for each taxon of virus were then reviewed to identify annotated terms that were shared between multiple types of viral infections, shared amongst coronaviruses, or specific to SARS-CoV-2.

We found several notable terms that were shared among the results obtained from infection of multiple unrelated taxa of viruses. These terms referred to functions such as the host interferon response, regulation of virus response, chemokine activity, and immune cell migration that were positively enriched in both the RSV infection of A549 cells and the SARS-CoV-2 infection of either A549 or NHBE cells. As expected, the statistical significance of these terms supports the important role that the human immune system plays in response to virus infection, which is an expected result and validates the upstream DE gene analysis.

Our method identified no significant GO terms that were shared across all Betacoronavirus comparisons. Interestingly, we found that each coronavirus species had its own set of uniquely enriched terms. In MERS-CoV, “negative regulation of inclusion body assembly” was positively enriched in four of the seven comparisons, while “aggrephagy” and “regulation of nucleotide binding oligomerization domain containing 2 signaling pathway” were positively enriched in three of the seven comparisons. The SARS-CoV samples had very few significantly enriched annotation terms including the “signal peptidase complex” and “cyclin dependent protein serine threonine kinase activator activity” terms that were each negatively enriched in only one of the four SARS-CoV comparisons.

We also identified 290 terms that were uniquely present among the SARS-CoV-2 comparisons, although no significant terms were identified in comparisons involving either the human post-mortem biopsies or the infected Calu-3 cells. The shared significant terms that we identified in the NHBE- and A549-infected cells included “response to chemokine”, “antimicrobial humoral response”, and “humoral immune response.” Enriched terms in the NHBE-infected cells alluded to the role of the inflammatory response, interleukin-1, interleukin-6, neuroinflammation, arachidonic acid binding, and T-cell activation.

Identification of affected signaling pathways

We then wanted to calculate which intracellular signaling pathways were significantly modulated during viral infection, based on the DE genes identified above. The pathway information used for this analysis was obtained from public databases that manually curate the proteins that participate in conveying a signal from receptors on the cell membrane to transcription factors within the nucleus in order to respond to a stimulus.

We used the results of the differential gene expression analysis as input to the SPIA algorithm. The robust bootstrap-based approach of this method identified 249 pathways that were significantly perturbed across at least one viral comparison. We subsequently analyzed the lists of significantly perturbed pathways for each comparison to determine those that were shared across virus taxa (e.g. influenza A, respiratory syncytial, MERS-CoV, SARS-CoV, and SARS-CoV-2), shared among coronaviruses, and unique to SARS-CoV-2. This analysis revealed a subset of the total number of pathways to be dysregulated across many of the viral taxa studied. These pathways included translation (affected in 12/17 comparisons and 4/5 viral taxa), cytokine-cytokine receptor interaction (modulated in 11/17 comparisons and 4/5 viral taxa), as well as rRNA processing (modulated in 10/17 comparisons and 3/5 viral taxa). We also observed that some of these pathways were predicted to be activated during infection with certain viral taxa in specific cell types, while inhibited in others.

We did not observe any signaling pathways that were significantly and consistently affected during infection by the individual coronavirus species examined. However, we did detect relevant pathways across multiple virus taxa. Interestingly, we found that the direction (i.e. activated or inhibited) was occasionally dependent on virus taxa or timepoint of infection. For example, the cytokine-cytokine receptor interaction pathway was predicted to be activated during infection with RSV, a subset of comparisons involving MERS-CoV, and many SARS-CoV-2 infections; while the same pathway was predicted to be inhibited during infection with SARS-CoV.

We predicted a total of 38 pathways that were uniquely affected during SARS-CoV-2 infection (Table 1). We noted that 35 of these pathways were only found to be significant in one SARS-CoV-2 comparison and included “NF-kB signaling pathway”, “Interleukin-1 signaling”, “IL6-mediated signaling events”, “PI3K-Akt signaling pathway”, “Jak-STAT signaling pathway”, “Apoptosis”, “Complement and coagulation cascades”, and other processes associated with either the immune system or infectious diseases. The remaining three pathways that were predicted as affected during at least two SARS-CoV-2 comparisons were “Cytokine signaling in Immune system”, “Tuberculosis”, and the “Innate immune system”. In short, these findings indicate a set of signaling pathways that are strongly associated with virus infection and/or immune activation in the host, some of which are uniquely detected during SARS-CoV-2 infection. The “Cytokine Signaling in Immune System” pathway is stored in the Reactome database and consists of interferon alpha/beta and gamma signaling, interleukin 1, 2, 6, and 10 signaling, and other components³⁷. To better understand the impact of SARS-CoV-2 on interleukin-6 signaling, we overlaid the differential expression data on a representation of the signaling pathway (Figure 2). This analysis revealed six of the eleven total proteins in the pathway were upregulated, while another protein was downregulated during SARS-CoV-2 infection.

Figure 2. Infection with SARS-CoV-2 shows differential gene expression in Interleukin-6 signaling.

This signaling pathway is a component of the larger Cytokine Signaling in Immune System pathway in Reactome. Each node represents a protein in the network, while each edge represents a characterized interaction between the proteins. Higher log2 fold-change values (up-regulation) are represented by increasing saturation of red, while more negative log2-fold-change values (down-regulation) are colored with increased saturation of blue. White nodes indicate no measured log2-fold change.

Prediction of relevant drug targets

We next wanted to determine whether any of the significantly affected signaling pathways contained known drug targets that could be modulated to reduce infection, virus replication, and/or clinical signs and symptoms associated with infection by a panel of Betacoronaviruses or by SARS-CoV-2 alone. To do so we cross-referenced the results from our significant pathways analysis with the drug-target information accessible through an application programming interface (API) to the opentargets.org database, which yielded 179 potential human drug targets (Table 2). We ranked the predicted drugs and their associated targets according to how many coronaviruses shared the same drug target across the various datasets included in our analysis. We performed a separate ranking based on the data obtained solely from the SARS-CoV-2 studies.

After reviewing the results of this analysis, we identified 38 potential human protein targets to counteract MERS-CoV, SARS-CoV, and/or SARS-CoV-2. Specifically, 38 targets were predicted to be relevant in at least 10 of the 15 coronavirus comparisons (Table 3). Seventeen of these targets were identified across 11 comparisons and included interferon (IFN)-A, IFN A receptors, IL1A and IL12B. Twenty-one additional targets were identified across 10 comparisons and included members of the tumor necrosis factor (TNF) superfamily, CCR4, CCR5, GHR, CXCR4, IL2R, IL4R, IL5, IL5R, IL6R, MPL, and IL20. We also predicted seven targets that were relevant specifically to SARS-CoV-2 including JAK1, FDPS, TLR7, TLR9, PRKCA, TYK2, and VDR.

Our analysis predicted twelve existing drugs that are predicted to be useful as repurposed therapeutics against 73.3% of all coronavirus taxa evaluated in this work and 100% of the SARS-CoV-2 comparisons included in our analysis (Table 4). Eight of these twelve drugs are used to treat common autoimmune disorders such as systemic lupus erythematosus (SLE), Crohn’s disease, and multiple sclerosis (MS). Five of the twelve drugs have been associated or used as antiviral measures largely against hepatitis C virus. Another 27 drugs were predicted to have potential therapeutic activity against 66.7% of all coronavirus infections and 30 were predicted as potential therapeutics against 75% of SARS-CoV-2 infections.