Pre-requisites

• Basic knowledge of R syntax
• Familiarity with RStudio
• Familiarity with pathway analysis
• Knowledge of Principal Components Analysis (PCA)⁸ is helpful but not required

Workshop participation

This will be a hands-on workshop, students will live-code with us. It would be helpful for participants to bring a laptop with RStudio and pathwayPCA package installed.⁹

R/Bioconductor packages used

To install the workshop dependencies, assuming you have already installed Bioconductor:¹⁰

install.packages(c("survival", "survminer", "tidyverse"))
BiocManager::install("rWikiPathways")
BiocManager::install("pathwayPCA")

We also strongly recommend the RStudio¹¹ integrated development environment as a user-friendly graphical wrapper for R. We require R version 3.6 or later.

library(survival)
library(survminer)
library(Bioc2019pathwayPCA)

devtools::install_github("gabrielodom/pathwayPCA")

Time outline

50m total

Learning goals

• Describe PCA-based pathway analysis
• Apply the pathwayPCA workflow to typical gene expression and copy number variations data
• Perform integrative pathway-based analysis for multiple types of -omics data jointly
• Apply pathwayPCA to experiments with complex designs, such as those with covariate and/or interaction effects

Learning objectives

• Identify pathways significantly associated with survival, binary, or continuous outcome using the pathwayPCA approach
• Estimate and visualize individual gene effects within a significant pathway
• Estimate and visualize sample-specific pathway activities
• Prioritize genes most likely to be driver (instead of passenger) genes using integrative analysis
• Identify pathways with significant sex-specific effects

Features of pathwayPCA

pathwayPCA allows users to:

1. Test pathway association with binary, continuous, or survival phenotypes.
2. Compute principal components (PCs) based on selected genes. These estimated latent variables represent pathway activities for individual subjects, which can then be used to perform integrative pathway analysis, such as multi-omics analysis.
3. Extract relevant genes that drive pathway significance (as well as data corresponding to these relevant genes) using the SuperPCA and AESPCA approaches for additional in-depth analysis.
4. Analyze studies with complex experimental designs, with multiple covariates, and with interaction effects, e.g., testing whether pathway association with clinical phenotype is different between male and female subjects.

How does pathwayPCA fit into the plethora of pathway analysis tools available?

• pathwayPCA tests the self-contained null hypothesis (Q1) (PMID: 21565265²²); that is, the features (e.g. genes) in a pathway are not associated with the disease phenotype.
• In contrast, many popular pathway analysis tools such as DAVID,²³ IPA,²⁴ Enrichr,²⁵ and goseq,²⁶ perform over-representation analysis and test for competitive null hypothesis (Q2); that is, the genes in a pathway show the same magnitude of associations with the disease phenotype compared with genes in the rest of the genome.
• When the “real” causal genes are fully contained in one particular pathway, testing Q1 and Q2 are approximately the same. However, when genes in multiple pathways are associated with the disease (as in many cancer studies) or when causal genes are shared by multiple gene sets, using competitive tests that compare pathway association signals with the rest of the genome may result in loss of power.
• pathwayPCA models correlations within pathways when constructing PCA scores.
• Many tools, including globaltest,²⁷ ignore expert-based biological information and fail to consider gene-gene correlations contained in biological pathways.
• Other tools, such as MOFA,²⁸ mixOmics,²⁹ mogsa³⁰, also perform integrative analysis using PCA. However, there is one important difference: these tools typically extract PCs on all the genes in dataset, and then map selected genes (e.g. those with nonzero loadings) to each pathway. In contrast, pathwayPCA subsets genes for each pathway first, then extracts PCs for each pathway: this helps to reduce noisy signals from irrelevant genes and improves sensitivity and specificity (on association with phenotype) of the PCs.
• pathwayPCA provides seamless multi-omics integration via estimated subject-specific pathway activity
• Some tools require additional software other than R packages to run, such as CoGAPS³¹ (requires OpenMP) or MOFA (requires Python).

Caveats for using pathwayPCA:

• Missing values in -omics data need to be imputed prior to analysis.
• In general, PCA is a poor choice for binary data. Therefore, pathwayPCA is a poor choice for GISTIC calls (copy number data) or mutation data.

Main functions

pathwayPCA uses four main groups of functions: data import and wrangling, object creation, object analysis, and analysis inspection. The main function groups are

1. Data:
   • read_gmt() imports a .gmt file as a pathway collection.
   • SE2Tidy() extracts an assay from a SummarizedExperiment object,³² and turns it into a “tidy” data frame.
   • TransposeAssay() is a variant of the base t() function designed specifcially for data frames and tibbles. It preserves row and column names after transposition.
2. Omics* Objects:
   • CreateOmics() takes in a collection of pathways, a single -omics assay, and a clinical response data frame and creates a data object of class Omics*.
   • SubsetPathwayData() can extract the pathway-specific assay values and responses for a given pathway from an Omics* object.
3. Methods:
   • AESPCA_pVals() takes in an Omics* object and calculates pathway p-values (parametrically or non-parametrically), principal components, and loadings via AESPCA. This returns an object of class aespcOut.
   • SuperPCA_pVals() takes in an Omics* object with valid response information and calculates pathway parametric p-values, principal components, and loadings via SuperPCA. This returns an object of class superpcOut.
4. Results:
   • getPathPCLs() takes in an object of class aespcOut or superpcOut and the TERMSname of a pathway. This function extracts 1) the data frame of principal components and subject IDs for the given pathway, and 2) a data frame of sparse loadings and feature names for the given pathway.
   • getPathpVals() takes in an object of class aespcOut or superpcOut and returns a table of the p-values and false discovery rates for each pathway.

pathwayPCA inputs

# DONT RUN
library(SummarizedExperiment)
data(airway, package = "airway")

airway_df <- SE2Tidy(airway, whichAssay = 1)

# DONT RUN
library(rWikiPathways)
# library(XML) # necessary if you encounter an error with readHTMLTable
downloadPathwayArchive(
  organism = "Homo sapiens", format = "gmt"
) 

trying URL 'http://data.wikipathways.org/current/gmt/wikipathways-20190110-gmt-Homo_sapiens.gmt'
Content type '' length 174868 bytes (170 KB)
downloaded 170 KB

#> [1] "wikipathways-20190110-gmt-Homo_sapiens.gmt"

dataDir_path <- system.file(
  "extdata", package = "pathwayPCA", mustWork = TRUE
)
wikipathways_PC <- read_gmt(
  paste0(dataDir_path, "/wikipathways_human_symbol.gmt"),
  description = TRUE
)

wikipathways_PC
#> Object with Class(es) 'pathwayCollection', 'list' [package 'pathwayPCA'] with 3 elements: 
#>  $ pathways   :List of 457
#>  $ TERMS      : chr [1:457] "WP23" ...
#>  $ description: chr [1:457] "B Cell Receptor Signaling Pathway" ...

Methods

Pathway-based Adaptive, Elastic-net, Sparse PCA (AESPCA)³⁹ combines the following methods into a single objective function: Elastic-Net⁴⁰, Adaptive Lasso⁴¹, and Sparse PCA⁴². AESPCA extracts principal components from pathway k which minimize this composite objective function. These extracted PCs represent activities within each pathway. The estimated latent variables are then tested against phenotypes using either a parametric or permutation test (permuting the sample responses). Note that the AESPCA approach does not use the response information to estimate pathway PCs, so it is an unsupervised approach.

Pathway-based Supervised PCA (SuperPCA):^{43 44} ranks each feature in pathway k by its univariate relationship with the outcome of interest (survival time, tumor size, cancer subtype, etc.), then extracts principal components from the features most related to the outcome. Because of this gene selection step, this method is “supervised”. Therefore, test statistics from the SuperPCA model can no longer be approximated using the Student’s t-distribution. To account for the gene selection step, SuperPCA as implemented in pathwayPCA estimates p-values from a two-component mixture of Gumbel extreme value distributions instead.

Ovarian cancer dataset

For this example, we will use the mass spectrometry based global proteomics data for ovarian cancer recently generated by the Clinical Proteomic Tumor Analysis Consortium (CPTAC).⁴⁵ The gene-level, log-ratio normalized protein abundance expression dataset for tumor samples can be obtained from the LinkedOmics database.⁴⁶ We used the dataset “Proteome (PNNL, Gene level)” which was generated by the Pacific Northwest National Laboratory (PNNL).⁴⁷ One subject was removed due to missing survival outcome. Proteins with greater than 20% missingness were removed from the assay. The remaining missing values were imputed using the Bioconductor package impute under default settings. The final dataset consisted of 5162 protein expression values for 83 samples.

Creating an Omics data object for pathway analysis

First, we need to create an Omics-class data object that stores

1. the expression dataset
2. phenotype information for the samples
3. a collection of pathways

data("ovProteome_df")
ovProteome_df[, 1:5]
#> # A tibble: 83 x 5
#>    Sample          A1BG    A2M    AAAS    AAK1
#>    <chr>          <dbl>  <dbl>   <dbl>   <dbl>
#>  1 TCGA-13-1489  0.0154 -0.419  0.0359  0.014 
#>  2 TCGA-42-2590 -0.508  -0.769 -0.0443  0.109 
#>  3 TCGA-36-2529 -0.148   0.209 -0.212   0.599 
#>  4 TCGA-24-1105  0.302  -0.343  0.116  -0.403 
#>  5 TCGA-29-1785 -0.584  -0.353 -0.192  -0.220 
#>  6 TCGA-24-2290 -0.954  -0.670 -0.0364  0.590 
#>  7 TCGA-24-1428 -0.352  -1.34   0.0283 -0.200 
#>  8 TCGA-24-1923  0.245  -0.471  0.0789 -0.0568
#>  9 TCGA-24-1563  0.0212  0.242  0.0636  0.785 
#> 10 TCGA-24-1430 -0.205  -0.154  0.0028  0.105 
#> # … with 73 more rows

data("ovPheno_df")
ovPheno_df
#> # A tibble: 565 x 3
#>    Sample       OS_time OS_status
#>    <chr>          <dbl>     <int>
#>  1 TCGA-3P-A9WA     420         0
#>  2 TCGA-59-A5PD     624         1
#>  3 TCGA-5X-AA5U     361         0
#>  4 TCGA-04-1331    1336         1
#>  5 TCGA-04-1332    1247         1
#>  6 TCGA-04-1335      55         1
#>  7 TCGA-04-1336    1495         0
#>  8 TCGA-04-1337      61         1
#>  9 TCGA-04-1338    1418         0
#> 10 TCGA-04-1341      33         0
#> # … with 555 more rows

ov_OmicsSurv <- CreateOmics(
  # protein expression data
  assayData_df = ovProteome_df, 
  # pathway collection
  pathwayCollection_ls = wikipathways_PC,
  # survival phenotypes
  response = ovPheno_df,
  # phenotype is survival data
  respType = "survival",
  # retain pathways with > 5 proteins
  minPathSize = 5  
)
There are 83 samples shared by the assay and phenotype data.

  ======  Creating object of class OmicsSurv  =======
The input pathway database included 5831 unique features.
The input assay dataset included 4763 features.
Only pathways with at least 5 or more features included in the assay dataset are
  tested (specified by minPathSize parameter). There are 312 pathways which meet
  this criterion.
Because pathwayPCA is a self-contained test (PMID: 17303618), only features in
  both assay data and pathway database are considered for analysis. There are 1936
  such features shared by the input assay and pathway database.

ov_OmicsSurv
#> Formal class 'OmicsSurv' [package "pathwayPCA"] with 6 slots
#>   ..@ eventTime            : num [1:83] 2279 3785 2154 2553 2856 ...
#>   ..@ eventObserved        : logi [1:83] TRUE TRUE TRUE TRUE TRUE TRUE ...
#>   ..@ assayData_df         :Classes 'tbl_df', 'tbl' and 'data.frame':    83 obs. of  4763 variables:
#>   ..@ sampleIDs_char       : chr [1:83] "TCGA-09-1664" "TCGA-13-1484" "TCGA-13-1488" "TCGA-13-1489" ...
#>   ..@ pathwayCollection    :List of 4
#>   .. ..- attr(*, "class")= chr [1:2] "pathwayCollection" "list"
#>   ..@ trimPathwayCollection:List of 5
#>   .. ..- attr(*, "class")= chr [1:2] "pathwayCollection" "list"

Testing pathway association with phenotypes

Once we have a valid Omics object, we can perform pathway analysis using the AESPCA or SuperPCA methodology as described above. Because the syntax for performing SuperPCA is nearly identical to the AESPCA syntax, we will illustrate only the AESPCA workflow below.

ovarian_aespcOut <- AESPCA_pVals(
  # The Omics data container
  object = ov_OmicsSurv,
  # One principal component per pathway
  numPCs = 1,
  # Use parallel computing with 2 cores
  parallel = TRUE,
  numCores = 2,
  # # Use serial computing
  # parallel = FALSE,
  # Estimate the p-values parametrically (AESPCA only)
  numReps = 0,
  # Control FDR via Benjamini-Hochberg
  adjustment = "BH"
)
Part 1: Calculate Pathway AES-PCs
Initializing Computing Cluster: DONE
Extracting Pathway PCs in Parallel: DONE

Part 2: Calculate Pathway p-Values
Initializing Computing Cluster: DONE
Extracting Pathway p-Values in Parallel: DONE

Part 3: Adjusting p-Values and Sorting Pathway p-Value Data Frame
DONE

names(ovarian_aespcOut)
#> [1] "pVals_df"    "PCs_ls"      "loadings_ls"

getPathpVals(ovarian_aespcOut)

ovOutGather_df <- getPathpVals(ovarian_aespcOut, score = TRUE, numPaths = 10)

ggplot(ovOutGather_df) +
  # set overall appearance of the plot
  theme_bw() +
  # Define the dependent and independent variables
  aes(x = reorder(terms, score), y = score) +
  # From the defined variables, create a vertical bar chart
  geom_col(position = "dodge", fill = "#005030", width = 0.7) +
  # Add a line showing the alpha = 0.0001 level
  geom_hline(yintercept = -log10(0.0001), size = 1, color = "#f47321") +
  # Add pathway labels
  geom_text(
    aes(x = reorder(terms, score), label = reorder(description, score), y = 0.1),
    color = "white",
    size = 4, 
    hjust = 0
  ) +
  # Set main and axis titles
  ggtitle("AESPCA Significant Pathways: Ovarian Cancer") +
  xlab("Pathways") +
  ylab("Negative Log10 (p-Value)") +
  # Flip the x and y axes
  coord_flip()

wp195PCLs_ls <- getPathPCLs(ovarian_aespcOut, "WP195")
wp195PCLs_ls
#> $PCs
#> # A tibble: 83 x 2
#>    sampleID          V1
#>    <chr>          <dbl>
#>  1 TCGA-09-1664 -2.90  
#>  2 TCGA-13-1484  0.0985
#>  3 TCGA-13-1488 -0.0724
#>  4 TCGA-13-1489  0.201 
#>  5 TCGA-13-1494 -0.506 
#>  6 TCGA-13-1495 -0.360 
#>  7 TCGA-13-1499 -0.529 
#>  8 TCGA-13-2071  2.11  
#>  9 TCGA-23-1123 -0.874 
#> 10 TCGA-23-1124  1.52  
#> # … with 73 more rows
#> 
#> $Loadings
#> # A tibble: 25 x 2
#>    featureID    PC1
#>    <chr>      <dbl>
#>  1 MAPK14    0     
#>  2 AKT1      0     
#>  3 IKBKB     0.531 
#>  4 NFKB1     0.528 
#>  5 PIK3R1    0.191 
#>  6 PLCG1     0     
#>  7 MAPK1     0.0882
#>  8 MAP2K1    0.300 
#>  9 MAP2K2    0     
#> 10 MAP2K6    0     
#> # … with 15 more rows
#> 
#> $pathway
#> [1] "path22"
#> 
#> $term
#> [1] "WP195"
#> 
#> $description
#> [1] "IL-1 signaling pathway"

wp195Loadings_df <- 
  wp195PCLs_ls$Loadings %>% 
  filter(PC1 != 0)
wp195Loadings_df
#> # A tibble: 8 x 2
#>   featureID    PC1
#>   <chr>      <dbl>
#> 1 IKBKB     0.531 
#> 2 NFKB1     0.528 
#> 3 PIK3R1    0.191 
#> 4 MAPK1     0.0882
#> 5 MAP2K1    0.300 
#> 6 TAB1      0.0992
#> 7 MAPK3     0.324 
#> 8 SQSTM1    0.436

wp195Loadings_df <- 
  wp195Loadings_df %>% 
  # Sort Loading from Strongest to Weakest
  arrange(desc(abs(PC1))) %>% 
  # Order the Genes by Loading Strength
  mutate(featureID = factor(featureID, levels = featureID)) %>% 
  # Add Directional Indicator (for Colour)
  mutate(Direction = factor(ifelse(PC1 > 0, "Up", "Down"))) 

ggplot(data = wp195Loadings_df) +
  # Set overall appearance
  theme_bw() +
  # Define the dependent and independent variables
  aes(x = featureID, y = PC1, fill = Direction) +
  # From the defined variables, create a vertical bar chart
  geom_col(width = 0.5, fill = "#005030", color = "#f47321") +
  # Set main and axis titles
  labs(
    title = "Gene Loadings on IL-1 Signaling Pathway",
    x = "Protein",
    y = "Loadings of PC1"
  ) +
  # Remove the legend
  guides(fill = FALSE)

ggplot(data = wp195PCLs_ls$PCs) +
  # Set overall appearance
  theme_classic() + 
  # Define the independent variable
  aes(x = V1) +
  # Add the histogram layer
  geom_histogram(bins = 10, color = "#005030", fill = "#f47321") +
  # Set main and axis titles
  labs(
    title = "Distribution of Sample-specific Estimate of Pathway Activities",
    subtitle = paste0(wp195PCLs_ls$pathway, ": ", wp195PCLs_ls$description),
    x = "PC1 Score for Each Sample",
    y = "Count"
  )

wp195Data_df <- SubsetPathwayData(ov_OmicsSurv, "WP195")
wp195Data_df
#> # A tibble: 83 x 28
#>    sampleID EventTime EventObs  MAPK14    AKT1   IKBKB   NFKB1  PIK3R1
#>    <chr>        <dbl> <lgl>      <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
#>  1 TCGA-09…      2279 TRUE     -0.995   0.958  -1.70   -0.730  -1.35  
#>  2 TCGA-13…      3785 TRUE      1.71    0.575  -0.744   0.783   1.18  
#>  3 TCGA-13…      2154 TRUE     -0.127  -0.542  -0.908   0.418   0.324 
#>  4 TCGA-13…      2553 TRUE     -0.0591 -0.291  -0.416  -0.488  -2.37  
#>  5 TCGA-13…      2856 TRUE      0.610   1.25   -0.943   0.394   0.309 
#>  6 TCGA-13…      2749 TRUE      0.265   1.79   -0.515  -0.0499 -0.0504
#>  7 TCGA-13…      3500 FALSE     0.516  -0.519  -0.487   0.104   1.54  
#>  8 TCGA-13…       773 TRUE     -0.107   0.323   1.05    1.27    1.35  
#>  9 TCGA-23…      1018 TRUE     -0.419  -3.63    0.0398 -0.114  -0.0345
#> 10 TCGA-23…      1768 TRUE     -1.81   -0.0215  0.817   1.17    0.567 
#> # … with 73 more rows, and 20 more variables: PLCG1 <dbl>, MAPK1 <dbl>,
#> #   MAP2K1 <dbl>, MAP2K2 <dbl>, MAP2K6 <dbl>, PTPN11 <dbl>, RELA <dbl>,
#> #   MAP3K7 <dbl>, IKBKG <dbl>, TAB1 <dbl>, IRAK4 <dbl>, ECSIT <dbl>,
#> #   TOLLIP <dbl>, MAPK3 <dbl>, MAP2K3 <dbl>, MAP2K4 <dbl>, TRAF6 <dbl>,
#> #   UBE2N <dbl>, SQSTM1 <dbl>, MAPKAPK2 <dbl>

library(survival)
NFKB1_df <- 
  wp195Data_df %>% 
  select(EventTime, EventObs, NFKB1)

wp195_mod <- coxph(
  Surv(EventTime, EventObs) ~ NFKB1,
  data = NFKB1_df
)

summary(wp195_mod)

NFKB1_df <- 
  NFKB1_df %>% 
  # Group subjects by gene expression
  mutate(NFKB1_Expr = ifelse(NFKB1 > median(NFKB1), "High", "Low")) %>% 
  # Re-code time to years
  mutate(EventTime = EventTime / 365.25) %>% 
  # Ignore any events past 10 years
  filter(EventTime <= 10)

# Fit the survival model
NFKB1_fit <- survfit(
  Surv(EventTime, EventObs) ~ NFKB1_Expr,
  data = NFKB1_df
)

ggsurvplot(
  NFKB1_fit,
  # No confidence intervals; add the p-value
  conf.int = FALSE, pval = TRUE,
  # Show times to median survival
  surv.median.line = "hv",
  xlab = "Time in Years",
  palette = c("#f47321", "#005030")
)

Creating copy number Omics data object

We can identify copy number (CNV) pathways significantly associated with survival in the same way as we did for protein expressions. This gene-level CNV data was downloaded from the same LinkedOmics database.

data("ovCNV_df")
ovCNV_df[, 1:5]
#> # A tibble: 579 x 5
#>    Sample        ACAP3 ACTRT2   AGRN ANKRD65
#>    <chr>         <dbl>  <dbl>  <dbl>   <dbl>
#>  1 TCGA-04-1331 -0.703 -0.703 -0.703  -0.703
#>  2 TCGA-04-1332  0.08   0.08   0.08    0.08 
#>  3 TCGA-04-1335 -0.807 -0.807 -0.807  -0.807
#>  4 TCGA-04-1336  0.101  0.101  0.101   0.101
#>  5 TCGA-04-1337  0.021  0.021  0.021   0.021
#>  6 TCGA-04-1338 -0.999 -0.999 -0.999  -0.999
#>  7 TCGA-04-1341 -0.421 -0.421 -0.421  -0.421
#>  8 TCGA-04-1342  0.089  0.089  0.089   0.089
#>  9 TCGA-04-1343  0.279  0.279  0.279   0.279
#> 10 TCGA-04-1346 -0.396 -0.396 -0.396  -0.396
#> # … with 569 more rows

And now we create an Omics data container. (Note: because these analysis steps take a little longer than the steps shown previously–3 minutes over 20 cores, we include the AESPCA output directly. You don’t need to execute the code in the next two chunks.)

# DONT RUN
ovCNV_Surv <- CreateOmics(
  assayData_df = ovCNV_df,
  pathwayCollection_ls = wikipathways_PC,
  response = ovPheno_df,
  respType = "survival",
  minPathSize = 5
)

1230 gene name(s) are invalid. Invalid name(s) are: ...

There are 549 samples shared by the assay and phenotype data.

  ======  Creating object of class OmicsSurv  =======
The input pathway database included 5831 unique features.
The input assay dataset included 24776 features.
Only pathways with at least 5 or more features included in the assay dataset are
  tested (specified by minPathSize parameter). There are 424 pathways which meet
  this criterion.
Because pathwayPCA is a self-contained test (PMID: 17303618), only features in
  both assay data and pathway database are considered for analysis. There are 5637
  such features shared by the input assay and pathway database.

# DONT RUN
ovCNV_aespcOut <- AESPCA_pVals(
  object = ovCNV_Surv,
  numPCs = 1,
  parallel = TRUE,
  numCores = 20,
  numReps = 0,
  adjustment = "BH"
)

Part 1: Calculate Pathway AES-PCs
Initializing Computing Cluster: DONE
Extracting Pathway PCs in Parallel: DONE

Part 2: Calculate Pathway p-Values
Initializing Computing Cluster: DONE
Extracting Pathway p-Values in Parallel: DONE

Part 3: Adjusting p-Values and Sorting Pathway p-Value Data Frame
DONE

data(ovCNV_aespcOut)

Combine significant pathways from CNV and protein analyses

Next, we identify the intersection of significant pathways based on both CNV and protein data. First, we will create a data frame of the pathway p-values from both CNV and proteomics.

# Copy Number
CNVpvals_df <- 
  getPathpVals(ovCNV_aespcOut, alpha = 0.05) %>% 
  mutate(rawp_CNV = rawp) %>% 
  select(description, rawp_CNV)

# Proteomics
PROTpvals_df <- 
  getPathpVals(ovarian_aespcOut, alpha = 0.05) %>% 
  mutate(rawp_PROT = rawp) %>% 
  select(description, rawp_PROT)


# Intersection
SigBoth_df <- inner_join(PROTpvals_df, CNVpvals_df, by = "description")
# WnT Signaling Pathway is listed as WP363 and WP428

The results showed there are 32 pathways significantly associated with survival in both CNV and protein data, which is significantly more than expected by chance (p-value = 0.00065; Fisher’s Exact Test; shown in multi_pathway_overlap_fishers.R). Here are the top-10 most significant pathways (sorted by protein data significance):

SigBoth_df[1:10, ]

# Copy Number Loadings
CNVwp195_ls <- getPathPCLs(ovCNV_aespcOut, "WP195")
CNV195load_df <- 
  CNVwp195_ls$Loadings %>% 
  filter(PC1 != 0) %>% 
  rename(PC1_CNV = PC1)

# Protein Loadings
PROTwp195_ls <- getPathPCLs(ovarian_aespcOut, "WP195")
PROT195load_df <- 
  PROTwp195_ls$Loadings %>% 
  filter(PC1 != 0) %>% 
  rename(PC1_PROT = PC1)

# Intersection
inner_join(CNV195load_df, PROT195load_df)
#> Joining, by = "featureID"
#> # A tibble: 5 x 3
#>   featureID PC1_CNV PC1_PROT
#>   <chr>       <dbl>    <dbl>
#> 1 IKBKB      0.0615   0.531 
#> 2 NFKB1      0.133    0.528 
#> 3 PIK3R1     0.0704   0.191 
#> 4 TAB1       0.0384   0.0992
#> 5 MAPK3      0.0560   0.324

(OPTIONAL) Integrating sample-specific pathway activities

Also in Section 1.5.3, we have seen that there can be considerable heterogeneity in pathway activities between patients. One possible reason could be that copy number changes might not directly result in changes in protein expression for some of the patients. pathwayPCA can be used to estimate pathway activities for each patient, for copy number, gene, and protein expressions separately. These estimates can then be viewed jointly using a Circos plot.

The accompanying Circos plot shown normalized copy number (inner circle) and protein expression (outer circle) pathway activities for the IL-1 signaling pathway in the ovarian cancer dataset samples. Each bar corresponds to a patient sample. Red bars indicate lower expression values and negative pathway activity for the sample, while blue color bars indicate higher expression values and stronger pathway activity for the sample. Note that only some patients have concordant changes in copy number and protein expression. The code to make this plot is included in inst/scripts/circle_plot.R.

Data setup and AESPCA analysis

For this example, we used IlluminaHiSeq gene expression RNAseq data from the TCGA KIRP cohort, which we downloaded from the Xena browser datahub.⁵² (This process is described in inst/scripts/explore_KIRP.R.) First, we load the KIRP data, create an Omics data container, and extract first AESPC from each pathway.

data("kirpPheno_df")
data("kirpRNAseq_df")

Once again, this will take a few moments (on my machine, this takes 1.9 minutes over 20 cores). (You can skip the next two chunks and load the results directly.)

# DONT RUN
kidney_Omics <- CreateOmics(
  assayData_df = kirpRNAseq_df,
  pathwayCollection_ls = wikipathways_PC,
  # Exclude the gender indicator for now
  response = kirpPheno_df[, 1:3],
  respType = "surv",
  minPathSize = 5
)

317 genes have variance < epsilon and will be removed. These gene(s) are: ...

30 gene name(s) are invalid. Invalid name(s) are: ...

There are 320 samples shared by the assay and phenotype data.

  ======  Creating object of class OmicsSurv  =======
The input pathway database included 5831 unique features.
The input assay dataset included 20213 features.
Only pathways with at least 5 or more features included in the assay dataset are
  tested (specified by minPathSize parameter). There are 423 pathways which meet
  this criterion.
Because pathwayPCA is a self-contained test (PMID: 17303618), only features in
  both assay data and pathway database are considered for analysis. There are 5566
  such features shared by the input assay and pathway database.

# DONT RUN
kirpRNAseq_aespcOut <- AESPCA_pVals(
  object = kidney_Omics,
  numPCs = 1,
  parallel = TRUE,
  numCores = 20,
  numReps = 0,
  adjustment = "BH"
)

Part 1: Calculate Pathway AES-PCs
Initializing Computing Cluster: DONE
Extracting Pathway PCs in Parallel: DONE

Part 2: Calculate Pathway p-Values
Initializing Computing Cluster: DONE
Extracting Pathway p-Values in Parallel: DONE

Part 3: Adjusting p-Values and Sorting Pathway p-Value Data Frame
DONE

Rather than execute this code yourself, we have included the output object with this package:

data("kirpRNAseq_aespcOut")

Test for sex interaction with first PC

Next we can test whether there is differential pathway association with survival outcome for males and females by the following model,
h(t) = h₀(t)exp [β₁PC₁+β₂male+β₃(PC₁×male)].

In this model, h(t) is expected hazard at time t, h₀(t) is baseline hazard when all predictors are zero, variable male is an indicator variable for male samples, and PC₁ is a pathway’s estimated first principal component based on AESPCA.

In order to test the sex interaction effect for all pathways, we will write a function which tests the interaction effect for one pathway.

TestIntxn <- function(pathway, pcaOut, resp_df){
  
  # For the given pathway, extract the PCs and loadings from the pcaOut list
  PCL_ls <- getPathPCLs(pcaOut, pathway)
  # Select and rename the PC
  PC_df <- PCL_ls$PCs %>% rename(PC1 = V1)
  # Join the phenotype data to this PC
  data_df <- inner_join(resp_df, PC_df, by = c("Sample" = "sampleID"))
  
  # Construct a survival model with sex interaction
  sex_mod <- coxph(
    Surv(time, status) ~ PC1 + male + PC1 * male, data = data_df
  )
  # Extract the model fit statistics for the interaction
  modStats_mat <- t(
    coef(summary(sex_mod))["PC1:maleTRUE", ]
  )
  colnames(modStats_mat) <- c("coef", "exp_coef", "se_coef", "z", "pVal")
 
  # Return a data frame of the pathway and model statistics
  list(
    statistics = data.frame(
      terms = pathway,
      description = PCL_ls$description,
      modStats_mat,
      stringsAsFactors = FALSE
    ),
    model = sex_mod
  )

}

As an example, we can test it on pathway WP195,

TestIntxn("WP195", kirpRNAseq_aespcOut, kirpPheno_df)$model
#> Call:
#> coxph(formula = Surv(time, status) ~ PC1 + male + PC1 * male, 
#>     data = data_df)
#> 
#>                  coef exp(coef) se(coef)      z     p
#> PC1           0.03278   1.03332  0.08206  0.399 0.690
#> maleTRUE     -0.38822   0.67826  0.32827 -1.183 0.237
#> PC1:maleTRUE  0.05649   1.05812  0.10438  0.541 0.588
#> 
#> Likelihood ratio test=3.97  on 3 df, p=0.2652
#> n= 320, number of events= 51

We can also apply this function to our list of pathways.

# List of pathways
paths_char <- kirpRNAseq_aespcOut$pVals_df$terms

# Apply over this list
interactions_ls <- sapply(
    paths_char,
    FUN = TestIntxn,
    pcaOut = kirpRNAseq_aespcOut,
    resp_df = kirpPheno_df,
    simplify = FALSE
)

We will tabulate this list and sort the pathways by significance:

interactions_df <- 
  # Take list of interactions
  interactions_ls %>% 
  # select the first element (the data frame of model stats)
  lapply(`[[`, 1) %>% 
  # stack these data frames
  bind_rows() %>% 
  as_tibble() %>% 
  # sort the rows by significance
  arrange(pVal)

Now, we can inspect the pathways with significant gender-PC1 interaction (note that we have not adjusted these pathway results for false discovery rate; these examples are only to show functionality of the package).

interactions_df %>% 
    filter(pVal < 0.05)

The results showed the most significant pathway is WP1559 (“TFs Regulate miRNAs related to cardiac hypertrophy”). We can inspect the model results for this pathway directly.

summary(interactions_ls[["WP1559"]]$model)
#> Call:
#> coxph(formula = Surv(time, status) ~ PC1 + male + PC1 * male, 
#>     data = data_df)
#> 
#>   n= 320, number of events= 51 
#> 
#>                 coef exp(coef) se(coef)      z Pr(>|z|)    
#> PC1           0.7534    2.1242   0.1782  4.228 2.36e-05 ***
#> maleTRUE      0.1825    1.2002   0.4268  0.427  0.66903    
#> PC1:maleTRUE -0.5974    0.5502   0.2162 -2.763  0.00573 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>              exp(coef) exp(-coef) lower .95 upper .95
#> PC1             2.1242     0.4708    1.4980    3.0122
#> maleTRUE        1.2002     0.8332    0.5199    2.7705
#> PC1:maleTRUE    0.5502     1.8174    0.3602    0.8406
#> 
#> Concordance= 0.654  (se = 0.047 )
#> Likelihood ratio test= 22.51  on 3 df,   p=5e-05
#> Wald test            = 30.03  on 3 df,   p=1e-06
#> Score (logrank) test = 33.09  on 3 df,   p=3e-07

The results showed that although sex is not significantly associated with survival outcome, the association of pathway gene expression (PC1) with survival is highly dependent on sex of the samples.

Survival curves by sex interaction

For visualization, we first match the PC expression to survival outcome.

# Bind the Pheno Data to WP1559's First PC
kidneyWP1559_df <- inner_join(
  kirpPheno_df,
  getPathPCLs(kirpRNAseq_aespcOut, "WP1559")$PCs,
  by = c("Sample" = "sampleID")
)

Now we divide subjects according to sex and high or low PC-expression, and add a group indicator for the four groups.

# Add Grouping Feature
kidneySurvWP1559grouped_df <- 
  kidneyWP1559_df %>% 
  rename(PC1 = V1) %>% 
  # add strength indicator
  mutate(direction = ifelse(PC1 > median(PC1), "High", "Low")) %>% 
  # group by interaction of sex and strength on PC
  mutate(group = paste0(direction, ifelse(male, "_Male", "_Female"))) %>% 
  # recode time in years
  mutate(time = time / 365.25) %>% 
  # remove summarized columns
  select(-male, -PC1, -direction)

Now we can plot survival curves for the four groups.

# Fit the survival model
fit <- survfit(
  Surv(time, status) ~ group,
  data = kidneySurvWP1559grouped_df
)

ggsurvplot(
  fit, conf.int = FALSE, 
  xlab = "Time in Years",
  ylim = c(0.4, 1),
  xlim = c(0, 10)
)

These Kaplan-Meier curves show that while high or low pathway activities were not significantly associated with survival in male subjects (green and purple curves, respectively), female subjects with high pathway activities (red) had significantly worse survival outcomes than those with low pathway activities (blue).