gabrielodom/Bioc2019pathwayPCA
A workshop on the pathwayPCA package

Package index
Search the gabrielodom/Bioc2019pathwayPCA package
Vignettes
Functions
10
Source code
13
Man pages
9
Home
/
GitHub
/
gabrielodom/Bioc2019pathwayPCA
/

In gabrielodom/Bioc2019pathwayPCA: A workshop on the pathwayPCA package 

knitr::opts_chunk$set(
  collapse = TRUE,
  cache = FALSE,
  comment = "#>"
)

library(knitr)
library(papeR)

figDir_char <- system.file(package = "Bioc2019pathwayPCA", "vignette")

Integrative Pathway Analysis with pathwayPCA

Instructors

Workshop Description

With the advance in high-throughput technology for molecular assays, multi-omics datasets have become increasingly available. In this workshop, we will demonstrate using the pathwayPCA package to perform integrative pathway-based analyses of multi-omics datasets. In particular, we will demonstrate through three case studies the capabilities of pathwayPCA to

There is a built version of this workshop available at RPubs.^[http://rpubs.com/gabrielodom/Bioc2019_pathwayPCA] The source is available on GitHub.^[https://github.com/gabrielodom/Bioc2019pathwayPCA]

Pre-requisites

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.^[https://gabrielodom.github.io/pathwayPCA/index.html#installing-the-package]

R/Bioconductor packages used

To install the workshop dependencies, assuming you have already installed Bioconductor:^[https://www.bioconductor.org/install]

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

We also strongly recommend the RStudio^[https://www.rstudio.com/products/rstudio/download/] integrated development environment as a user-friendly graphical wrapper for R. We require R version 3.6 or later.

Package descriptions

This workshop package, Bioc2019pathwayPCA, depends on two packages:

The tidyverse package suite is a collection of utility packages for data science: ggplot2, dplyr, tidyr, readr, purrr, tibble, stringr, and forcats. We will make use of some data constructs and ideas from these packages, but we do not expect users to be intimately familiar with them.

Additionally, because we will use survival data in our examples, we need survival^[https://CRAN.R-project.org/package=survival] and survminer^[https://cran.r-project.org/package=survminer] (for drawing survival curves with ggplot2).

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

(OPTIONAL) Development version

Because we are currently in the development phase for version 2 of this package, you can install the package from GitHub. In order to install a package from GitHub, you will need the devtools::^[https://github.com/r-lib/devtools] package and either Rtools^[https://cran.r-project.org/bin/windows/Rtools/] (for Windows) or Xcode^[https://developer.apple.com/xcode/] (for Mac). Then you can install the development version of the pathwayPCA package from GitHub via:

devtools::install_github("gabrielodom/pathwayPCA")

Time outline

50m total

| Activity | Time | |----------------------------------------------------------------------|------| | Introduction to pathwayPCA | 5m | | Case Study 1: Estimating sample-specific pathway activities | 15m | | Case Study 2: Pathway analysis for multi-omics datasets | 10m | | Case Study 3: Analyzing experiments with covariates and interactions | 10m | | Summary and Conclusion | 5m | | Questions and Comments | 5m |

Workshop Goals and Objectives

Learning goals

Learning objectives



Introduction to pathwayPCA

pathwayPCA is an integrative analysis tool that implements PCA-based pathway analysis approaches.^[ Chen et al. (2008) ] ^[ Chen et al. (2010) ] ^[ Chen (2011) ]

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?

Caveats for using pathwayPCA:

Main functions

include_graphics(file.path(figDir_char, "pathwayPCA_workflow.png"))

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,^[https://doi.org/10.18129/B9.bioc.SummarizedExperiment] 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

Assay data

We define our assay data format in R as follows:

Use the two utility functions, TransposeAssay() or SE2Tidy(), to wrangle your assay data into appropriate form for further analysis. Here is an example of "tidying" the first assay from a SummarizedExperiment object:

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

airway_df <- SE2Tidy(airway, whichAssay = 1)

Pathway collections

Based on expert knowledge of the problem at hand, curate a set of biological pathways $1, \ldots, K$, where pathway $k$ contains $p_k$ -ome measurements. These pathway collections are often stored in a .gmt file, a text file with each row corresponding to one pathway. Each row contains an ID (column TERMS), an optional description (column description), and the genes in the pathway (all subsequent columns). This file format is independently defined by the Broad Institute.^[http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats] Pathway collections in .gmt form can be downloaded from the MSigDB database.^[http://software.broadinstitute.org/gsea/msigdb/collections.jsp]

For WikiPathways,^[https://www.wikipathways.org/index.php/WikiPathways] a collection of community-defined biological pathways, one can download monthly data releases in .gmt format using the dowloadPathwayArchive() function in the rWikiPathways package from Bioconductor.^[https://bioconductor.org/packages/release/bioc/html/rWikiPathways.html] For example, the following commands^[https://bioconductor.org/packages/release/bioc/vignettes/rWikiPathways/inst/doc/Overview.html#4_give_me_more] download the latest release of the human pathways from WikiPathways to your current directory (we do not execute this code):

# 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

print("wikipathways-20190110-gmt-Homo_sapiens.gmt")

pathwayPCA includes the June 2018 Wikipathways collection for homo sapiens, which can be loaded using the read_gmt function: 

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

This pathway collection contains:

  1. pathways: a list of character vectors matching some of the names of the -omes recorded in the assay
  2. TERMS: a character vector of the names of the pathways
  3. description: (OPTIONAL) descriptions or URLs of the pathways
wikipathways_PC

Methods

Pathway-based Adaptive, Elastic-net, Sparse PCA (AESPCA)^[ Chen (2011) ] combines the following methods into a single objective function: Elastic-Net^[ Zou and Hastie, 2005 ], Adaptive Lasso^[ Zou, 2012 ], and Sparse PCA^[ Zou, Hastie, and Tibshirani, 2012 ]. 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):^[ Chen et al. (2008) ] ^[ Chen et al. (2010) ] 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.



Case Study 1: Identifying Significant Pathways in Protein Expressions Associated with Survival Outcome in Ovarian Cancer

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).^[https://proteomics.cancer.gov/programs/cptac] The gene-level, log-ratio normalized protein abundance expression dataset for tumor samples can be obtained from the LinkedOmics database.^[http://linkedomics.org/data_download/TCGA-OV/] We used the dataset “Proteome (PNNL, Gene level)” which was generated by the Pacific Northwest National Laboratory (PNNL).^[https://www.pnnl.gov/] 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

Expression and phenotype data

We can obtain protein expression and phenotype datasets for the TCGA ovarian cancer dataset by downloading the raw files from the LinkedOmics website. However, for ease and to save time, we include cleaned and imputed versions of these datasets within this package (for the cleaning process, please see the clean_multi_omics.R file in inst/scripts).

The ovProteome_df dataset is a data frame with protein expression levels and TCGA sample IDs. The variables (columns) include expression data for 4763 proteins for each of the 83 primary tumour samples.

data("ovProteome_df")
ovProteome_df[, 1:5]

The ovPheno_df data frame contains TCGA sample IDs, overall survival time, and overall survival censoring status.

data("ovPheno_df")
ovPheno_df

Create an OmicsSurv data container

Now that we have these three data components (pathway collection, proteomics, and clinical responses), we create an OmicsSurv data container. Note that when assayData_df and response are supplied from two different files, the user must match and merge these data sets by sample IDs. The assay and response must have matching row names, otherwise the function will error.

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  
)

To see a summary of the Omics data object we just created, simply type the name of the object:

ov_OmicsSurv

Testing pathway association with phenotypes {#test_path_assoc}

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.

Implementation

We estimate pathway significance via the following model: phenotype ~ intercept + PC1. Pathway $p$-values are estimated based on a likelihood ratio test that compares this model to the null model phenotype ~ intercept. Note that when the value supplied to the numReps argument is greater than 0, the AESPCA_pvals() function employs a non-parametric, permutation-based test to assign the statistical significance of this model.

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"
)

This ovarian_aespcOut object contains 3 components: a table of pathway $p$-values, AESPCA-estimated PCs of each sample from each pathway, and the loadings of each protein onto the AESPCs.

names(ovarian_aespcOut)

We haven't yet added a print() method for these analysis outputs, so be careful: use the two accessor functions, getPathpVals() and getPathPCLs() instead. We will show examples of these functions next.

Table of pathway analysis results

For this ovarian cancer dataset, the top 20 most significant pathways identified by AESPCA are:

getPathpVals(ovarian_aespcOut)

kable(
  getPathpVals(ovarian_aespcOut)
)

(OPTIONAL) Column chart of significant pathways

Before constructing a graph of the $p$-values, we extract the top 10 pathways (the default value for numPaths is 20). The score = TRUE returns the negative natural logarithm of the unadjusted $p$-values for each pathway (to enhance graphical display of the top pathways).

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

Now we plot the pathway significance level for the top 20 pathways.

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()

Extract relevant genes from significant pathways {#extr_rel_genes}

Because pathways are defined a priori, typically only a subset of genes within each pathway are relevant to the phenotype and contribute to pathway significance. In AESPCA, these relevant genes are the genes with nonzero loadings in the first PC of AESPCs.

For example, we know that the IL-1 Signaling Pathway is one of cancer's "usual suspects".^[https://doi.org/10.1371/journal.pcbi.1003470] For the "IL-1 signaling pathway" (Wikipathways WP195^[https://www.wikipathways.org/index.php/Pathway:WP195]), we can extract the PCs and their protein Loadings using the getPathPCLs() function:

wp195PCLs_ls <- getPathPCLs(ovarian_aespcOut, "WP195")
wp195PCLs_ls

The proteins with non-zero loadings can be extracted as follows: 

wp195Loadings_df <- 
  wp195PCLs_ls$Loadings %>% 
  filter(PC1 != 0)
wp195Loadings_df

(OPTIONAL) Plot protein loadings for a single pathway

We can also prepare these loadings for graphics: 

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")))

Now we will construct a column chart with ggplot2's geom_col() function.

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)

Alternatively, we can also plot the correlation of each gene with first PC for each gene. These correlations can be computed by using the TidyCorrelation() function in pathwayPCA's vignette “Chapter 5 - Visualizing the Results”, Section 3.3.^[https://gabrielodom.github.io/pathwayPCA/articles/Supplement5-Analyse_Results.html#gene-correlations-with-pcs]

Subject-specific PCA estimates {#subj_spec_scores}

In the study of complex diseases, there is often considerable heterogeneity among different subjects with regard to underlying causes of disease and benefit of particular treatment. Therefore, in addition to identifying disease-relevant pathways for the entire patient group, successful (personalized) treatment regimens will also depend upon knowing if a particular pathway is dysregulated for an individual patient.

To this end, we can also assess subject-specific pathway activity. As we saw earlier, the getPathPCLs() function also returns subject-specific estimates for the individual pathway PCs. We can plot these as follows.

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"
  )

This graph shows there can be considerable heterogeneity in pathway activities between the patients.

Extract analysis dataset for entire pathway

Users are often also interested in examining the actual data used for analysis of the top pathways, especially for the relevant genes with the pathway. To extract this dataset, we can use the SubsetPathwayData() function. These commands extract data for the IL-1 signaling pathway:

wp195Data_df <- SubsetPathwayData(ov_OmicsSurv, "WP195")
wp195Data_df

Gene-specific CoxPH model

We can also perform analysis for individual genes belonging to the pathway: 

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

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

Now we inspect the output from the Cox PH model (we've included the pretty version using the prettify() function from the papeR package^[https://CRAN.R-project.org/package=papeR]):

summary(wp195_mod)

wp195pretty_mod <- prettify(
  summary(wp195_mod)
)

kable(wp195pretty_mod)

(OPTIONAL) Gene-specific survival curves

Additionally, we can estimate Kaplan-Meier survival curves for patients with high or low expression values for individual genes: 

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
)

Finally, we can plot these K-M curves over NFKB1 protein expression.

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")
)



Case study 2: An Integrative Multi-Omics Pathway Analysis of Ovarian Cancer Survival

While copy number alterations are common genomic aberrations in ovarian carcer, recent studies have shown these changes do not necessarily lead to concordant changes in protein expression. In Section 1.5.3 above, we illustrated testing pathway activities in protein expression against survival outcome. In this section, we will additionally test pathway activities in copy number against survival outcome. Moreover, we will perform integrative analysis to identify those survival-associated protein pathways, genes, and samples driven by copy number alterations.

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]

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.

AESPCA pathway analysis for copy-number data

Finally, we can apply the AESPCA method to this copy-number data container. Due to the large sample size, this will take a few moments.

# 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

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

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, ]

kable(
  SigBoth_df[1:10, ]
)

Integrative pathway-based gene detection

Similar to the protein pathway analysis shown in Section 1.5.3, we can also identify relevant genes with nonzero loadings that drives pathway significance in CNV. The "IL-1 signaling pathway" (WP195) is significant in both CNV and protein data, so we look for genes with non-zero loadings in both.

# 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)

The result showed that NFKB1, IKBKB, and other genes are selected by AESPCA when testing IL-1 signaling pathway (WP195) against survival outcome in both CNV and protein pathway analysis.

(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.

include_graphics(file.path(figDir_char, "circosPlot_20190617.png"))

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.



(OPTIONAL) Case Study 3: Analysis of Studies with Complex Designs

pathwayPCA is capable of analyzing studies with multiple experimental factors. In this section, we illustrate using pathwayPCA to test differential association of pathway expression with survival outcome in male and female subjects.

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.^[https://xenabrowser.net/datapages/] (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_0(t) \exp\left[ \beta_1\text{PC}_1 + \beta_2\text{male} + \beta_3(\text{PC}_1 \times \text{male}) \right]. ]

In this model, $h(t)$ is expected hazard at time $t$, $h_0 (t)$ is baseline hazard when all predictors are zero, variable $male$ is an indicator variable for male samples, and $PC_1$ 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

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)

kable(
  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)

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).



Further Reading

For addtional information on pathwayPCA, each of our supplementary vignette chapters contain detailed tutorials on the topics discussed above. These vignettes are:

References

Chen, X., Wang, L., Smith, J.D. and Zhang, B. (2008) Supervised principal component analysis for gene set enrichment of microarray data with continuous or survival outcomes. Bioinformatics, 24, 2474-2481.

Chen, X., Wang, L., Hu, B., Guo, M., Barnard, J. and Zhu, X. (2010) Pathway-based analysis for genome-wide association studies using supervised principal components. Genetic epidemiology, 34, 716-724.

Chen, X. (2011) Adaptive elastic-net sparse principal component analysis for pathway association testing. Statistical applications in genetics and molecular biology, 10.

The R session information for this vignette is:

sessionInfo()


gabrielodom/Bioc2019pathwayPCA documentation built on June 19, 2019, 7:40 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close