Using Human Differentially Expressed Gene Lists to Perform Downstream Pathway Enrichment Analysis and Target Prioritization

Summary

The current work describes a protocol for running the Pathway2Targets algorithm, an R script that predicts and prioritizes therapeutic targets based on the profile of intracellular signaling pathways generated by comparing case versus control samples from a bulk RNA-sequencing experiment.

Abstract

This protocol outlines a multi-step computational pipeline for identifying potential therapeutic targets from RNA-sequencing data, including installation of relevant software, setup verification, and differential expression analysis using edgeR. We then show how to utilize the signaling pathway impact analysis (SPIA) algorithm to predict statistically significant pathways. To ensure confidence in the results, we focus on significant pathways (p < 0.05) to reduce false positive results. Unlike traditional gene sets, these pathways reflect protein-protein interaction networks, offering mechanistic insights into cellular processes like the cell cycle, immune response, and metabolism. These pathways are then analyzed using the Pathway2Targets algorithm, which interfaces with the OpenTargets.org database via an application programming interface (API). This algorithm incorporates a novel weighting approach that scores known drug targets within the identified pathways, while providing progress in real-time. Runtime depends on pathway complexity and target density. The output consists of two ranked files. The first file contains a list of predicted drug targets and their weighted scores, while the second consists of various details for the associated therapeutics. Together, this pipeline facilitates the prioritization of druggable targets and treatments in the context of disease-specific gene expression profiles.

Introduction

Bulk RNA-sequencing enables the comparison of expression levels for thousands of genes across a population of case cells versus a population of control cells. Experiments are typically designed to include at least triplicate samples, ideally biological replicates, though technical replicates can suffice. This design accounts for biological variability and reduces the impact of outlier samples. Analysis of these expression patterns provides detailed insight into the effect of the disease of interest on normal cellular processes and can potentially enable the prediction of relevant therapeutics.

Preprocessing of bulk RNA-sequencing data typically includes: quality control of sequencing reads (for repeats, sequencing adapters, GC%, etc.), read trimming and adapter removal, read mapping/quantitation^1,2,3, and differential expression analysis^4,5,6. Fortunately, a variety of analytical processes have been automated to reduce manual work associated with these steps^7,8,9. After preprocessing is complete, commonly performed downstream analyses include functional overrepresentation analysis with Gene Ontologies, signaling pathways enrichment, and variation in splicing. These downstream analyses summarize and facilitate interpretation of the differential expression results at a higher level of granularity than the gene lists alone.

Various tools have been developed with the goal of repurposing existing therapeutics for a well-defined type or subtype of disease. This is achieved by training the algorithm on multiple-omics data types for the intended disease. Unfortunately, such efforts to improve specificity and sensitivity in an intended disease often make using the tools in more general contexts suboptimal^10,11. Another set of tools is more broadly applicable to cases where gene expression profiles are matched to either existing signatures of gene expression^12,13, or to the quantified effects of current therapeutics^14,15. However, these more broadly applicable tools often achieve reduced specificity and sensitivity across a wide range of diseases and/or were trained with outdated data.

In contrast, the Pathway2Targets algorithm has previously been applied to predict potential therapeutic targets in B-cell lymphoma, periodontitis, estrogen-positive breast cancer, triple-negative breast cancer, and chikungunya virus^{16,17,18,19,20,21}. The results from these studies demonstrate that this tool is capable of predicting robust and biologically relevant targets. Impressively, Pathway2Targets predicted 392 potential drug targets for triple-negative breast cancer, among which 60 were tested in clinical trials; as well as 828 individual drugs for TBNC, with 37 being tested¹⁷. In the lymphoma study, this algorithm predicted 915 drugs, 461 of which are FDA-approved¹⁹.

The aim of the current work is to describe a computational protocol that will enable more researchers, who may benefit from access to more descriptive instructions on running programs on the command-line, to effectively use the recently developed Pathway2Targets algorithm (Figure 1). Pathway2Targets predicts targets for a given condition by combining differential expression data, gene-disease associations, clinical trial information, public target data²², pathway information, and other metrics. Importantly, this algorithm incorporates a unique and customizable weighting scheme, which allows users to determine the ~20 target-related metrics that they prefer to emphasize in their analysis, such as number of disease associations, number of signaling pathways, number of unique drugs, number of therapeutics in each phase of clinical trials, etc.²³. As an example use case for this protocol, we will re-analyze an existing colorectal cancer dataset²⁴.

Protocol

The bulk RNA-sequencing data analyzed in this study were acquired from publicly available databases (the NCBI Gene Expression Omnibus and Sequence Read Archive)^25,26. As such, the original data collectors ensured the ethical and appropriate collection of these samples from informed and consenting human subjects.

1. Download and install R software

1. Install R (version 4.0 or later) by clicking an appropriate link from the Comprehensive R Archive Network (CRAN) https://cran.r-project.org/mirrors.html and use the option 0-Cloud, then follow the appropriate instructions for the computer's operating system. This process generally takes 5-10 min.
2. Install R Studio (version 2024 or later) from https://posit.co/download/rstudio-desktop/, then follow the instructions on the download page. Installing RStudio generally takes ~10 min.
   NOTE: Installing RStudio is optional (but highly recommended), as it provides an integrated development environment that simplifies code execution, offers syntax highlighting, facilitates package management, and helps users visualize outputs, which is especially beneficial for users who are less familiar with R.

2. Download and install R scripts for relevant tools

1. Download the following necessary R scripts from the GitHub repository https://github.com/bpickett/Pathway2Targets. This URL is for reference only.
2. Download the scripts using the below links: SPIA version 1.0 (download by clicking Download Raw File): https://github.com/bpickett/Pathway2Targets/blob/main/SPIA_Code.Rmd (commit ID: 60fcd46); Pathway2Targetsversion 3.1 (download by clicking Download Raw File): https://github.com/bpickett/Pathway2Targets/blob/main/Pathway2Targets.R (commit ID: 8e4c7c8)

3. Download R libraries for relevant tools

1. While running R (either in RStudio or in a terminal window), enter the following commands to download and install the required additional R libraries needed to run the software.
   1. Start RStudio program. By default, the console panel is located at the bottom left corner of RStudio. Click anywhere within the console panel window, and a typing cursor should appear at the bottom after the arrow ">" symbol.
   2. Copy and paste the following command in the console area and press the Enter key: 
      install.packages(c("RCurl", "stringr", "jsonlite", "httr")). 
      Once successfully installed, a status message will appear that says "The downloaded binary packages are in ....".
   3. Copy and paste the following command into the console panel area and press the Enter key.
      BiocManager::install(c("SummarizedExperiment", "EnrichmentBrowser", "biomaRt", "org.Hs.eg.db")).
      Once successfully installed, a similar message will appear that says "The downloaded binary packages are in ....".
      NOTE: The first set of libraries (step 3.1.1) consists of typical R libraries, while the second (2.1.2) consists of BioConductor libraries. As such, the commands need to be entered separately. The appropriate versions of these libraries should automatically download based on the R version that is installed on the computer. Each of these commands should take ~5 min to complete.

4. File processing

1. Download the previously computed edgeR-generated differential expression output file (in RDS format), generated by ARMOR software (or similar), to the local computer. The name of this file is generally edgeR_dge.rds.
2. Manually review the results of edgeR (or a similar differential expression) to begin biologically relevant interpretation of the results. To do this, filter by, at a minimum, a corrected p-value < 0.05, and potentially the absolute value of the log₂ fold-change value > 1.5. If testing out the software, a sample edgeR_dge.rds file can be found on Zenodo here: https://doi.org/10.5281/zenodo.15186609
   NOTE: Reviewing the genes that remain in these filtered results can begin to explain the underlying molecular mechanism(s) of the phenotype associated with the case samples (compared to the control samples). It is important to recognize that the ability to interpret lists of genes in an unbiased way is extremely difficult due to the relatively low number of gene symbols that can be quickly recalled. As such, signaling pathway analysis is a useful way to summarize the filtered genes based on how they interact and/or communicate with each other in the cell.
3. Pre-processing bulk RNA-seq data can take anywhere from several hours to days of compute time, depending on the size of the dataset being analyzed. Save this .rds file in the Downloads folder on the computer. Note that this .rds file type is not human-readable.

5. Run SPIA pathway enrichment algorithm

1. If running using R, use the R script either from GitHub or Supplementary Coding File 1. Type in the following command, assuming the edgeR_dge.rds file is in the Downloads folder
   Rscript --vanilla SPIA_Code.Rmd ~/Downloads/edgeR_dge.rds
   1. If the edgeR_dge.rds file is in a different folder (or directory), then replace this command with the following
      Rscript --vanilla SPIA_Code.Rmd <path_to_'edgeR_dge.rds'_file>
2. If running using RStudio, use the R script from GitHub or Supplementary Coding File 2.
   NOTE: In the language R, adding hash tag # symbol before a line of code temporarily disables it. The scripts were originally designed to run in a command-line environment instead of RStudio. Enabling or disabling certain lines of code is the easiest way to reconfigure the setup for the input files.
   1. Open the SPIA_Code.Rmd script in R Studio by clicking the Open File option in the File menu, then selecting the script name. In the R Studio code windowby default, this is located at the top left panel.
   2. Select all lines of code in the file and click the Run button (or Run Selected Line(s)) button, which is located above and to the right of the code window. A successful run will produce a file named similar to
      edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23_
10-56-45.12767_SPIA_Results.csv
      in the download directory. This file will contain the statistically significant signaling pathways
   3. Review the file with the statistically significant results manually by opening it as a spreadsheet. The contents of this file should help to summarize the underlying intracellular signaling cascades that are significantly represented by the differentially expressed genes.
      ​NOTE: Completing the calculation of significant pathways may take anywhere from ~30 min to several hours, depending on the strength of the signal in the dataset being analyzed. When the program is running, real-time progress messages will continuously update in the console window. The frequently updated messages show that the program is working successfully. More detailed descriptions of what occurs in this step can be found in the GitHub repository: https://github.com/bpickett/Pathway2Targets/tree/main
   4. If using the sample input file, the output file from this step will be found in the Downloads folder and will be named
      edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_"timestamp"_SPIA_Results.csv.
      This naming style reflects the input file name, the process being performed, and the output. This avoids confusion in identifying files if more than one is processed.
   5. Adjust other parameters as described below.
      1. The default parameters for the SPIA algorithm in this library is 1,000 permutations. Increase this to 2,000 permutations to improve confidence in the results. Adjust the number of permutations in lines 84 and 85 of this script by changing the perm = 2000 to the desired number of permutations. Adjust other parameters as described below.
      2. Adjust the p-value correction approach in lines 84 and 85 by removing the padj.method = 'BH'. This will not correct p-values, which may make false-positive results more likely.

6. Running Pathway2Targets target prioritization algorithm on SPIA output

1. If running R, use the R script from GitHub or Supplementary Coding File 3. Use the following command to invoke this algorithm
   Rscript --vanilla Pathway2Targets.R
2. If running using R Studio, use the R script from GitHub or Supplementary Coding File 4. Open the Pathway2Targets.R script in R Studio by clicking the Open File option in the File menu, then selecting the script name.
   1. In the RStudio code window (top-left panel), replace the file name on line 22 with the SPIA results file name, such as (from the example data)
      infile <- "edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23_
10-56-45.12767_SPIA_Results.csv"
   2. Select all lines of code in the file and click the Run button, which is located above and towards the right side of the code window. Real-time progress status messages will be continuously displayed in the bottom-right panel. A successful run will produce a file called (similar to)
      "edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23_
10-56-45.12767_SPIA_Results.csv-RankedTargets.tsv" in the download directory. The naming style of the files reflects the input, the process, and the output 
      edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_2025-04-23
_10-56-45.12767_SPIA_Results.csv
      is the input for such file.
      NOTE: This step can take anywhere from one to several hours, depending on the number of signaling pathways with a significant p-value, the number of gene products in those significant pathways, and the number of gene products that are known drug targets.
3. Some parameters for the Pathway2Targets algorithm can be adjusted. Specifically, adjust the multiplier values (lines 31-38 of the script) to customize the weighting scheme for each metric. More detailed descriptions of what occurs in this step can be found in the corresponding GitHub repository: https://github.com/bpickett/Pathway2Targets/tree/main
   NOTE: For reference, in the example, identifying potential drugs from 132 pathways, consisting of hundreds of individual targets, takes approximately 2 h to complete. Based on this metric, it is reasonable to estimate total computing time.

7. Open results files for prioritized targets and therapeutics

1. The file with prioritized targets and their metrics will be generated. For the sample input files, use the output file found in the Downloads folder, named:
   edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_"timestamp"_SPIA_Results.csv-RankedTargets.tsv
   1. This file will, by default, be sorted with the targets sorted in descending order based on the custom-weighed metric. Manually review the output file to ensure that the results are biologically relevant and that the targets are logical for the phenotype being evaluated.
2. The file with prioritized therapeutics and their metrics will also be generated. For the sample input files, use the output file in the Downloads folder, named:
   edgeR_dge.rds-TreatmentTumor-TreatmentNativeTissue_"timestamp"_SPIA_Results.csv-Treatments.tsv
   1. Similarly, this output file will also, by default, be sorted with the therapeutics for the various targets (in step 7.1.1) sorted in descending order based on the weighted metric. Review the file manually with sufficient background knowledge of the underlying biological system to determine whether additional experiments are justified. It is expected that multiple therapeutics may have the same weighted metric since many of the targets will be affected by more than one therapeutic on the market.

Results

The setup described in steps 1-3 of the protocol is necessary to enable subsequent execution of SPIA and the Pathway2Targets algorithm. At the end of each step, a message will be generated to confirm the successful installation of the software. Step 4 consists of either downloading an existing set of differential expression results, which could include the provided example file, a different existing file, or preprocessing a custom RNA-sequencing dataset. The principal requirement for step 4 is that the workflow uses edge...

Discussion

Steps 1-3 of the protocol specifically pertain to installing the underlying R software, scripts, and dependencies to successfully run the downstream software. A detailed list of required R libraries is available (Supplementary Table 4). Step 4 of the protocol involves retrieving an R data (.rds format) file, which contains the output from a differential expression analysis. Common software for this step involves edgeR⁶, DESeq2⁴, and limma^5...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Pathway2Targets R script	Brigham Young University (Pickett Lab)	Version 3.1	https://github.com/bpickett/Pathway2Targets/blob/main/Pathway2Targets.R
R Software	Comprehensive R Archive Network (CRAN)	Version: 4.4.3	https://cran.r-project.org
R Studio Desktop Software	posit	Version: 2024.12.1+563	https://posit.co/download/rstudio-desktop/
SPIA R script	Brigham Young University (Pickett Lab)	Version: 3.1	https://github.com/bpickett/Pathway2Targets/blob/main/SPIA_Code.Rmd