In pouryany/pathwayDisturbance: Pathway enrichment analysis of genomic using Causal Disturbance Analysis

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

Overview

This document contains the walk-through/tutorial for the Causal Disturbance Analysis (CADIA) as described in the corresponding publication by Naderi and Mostafavi (References to be provided). CADIA is an enrichment analysis tool for interpreting gene perturbations by contrasting them with underlying graphs of annotated pathways. This program takes an input list of differentially expressed gene IDs and a gene universe and produces p-values that indicate pathway enrichments.

Dependencies and Installation guide

CADIA depends on the following packages: KEGGgraph, RBGL, graph, stringr, dplyr, magrittr. Make sure the packages are installed before using CADIA. To install the packages you can use the following code chunk.

dependencies  <- c("KEGGgraph", "RBGL", "graph", "stringr", "dplyr", "magrittr")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Please note the parameters of BiocManager::install() and modify accordingly

for (i in dependencies) {
    if(!require(i))
        BiocManager::install(i,update = F)
}

Package Installation and Preparation

You can install CADIA using devtools library. To install devtool run the following commnad:

install.packages("devtools")

After having devtools installed, you can install CADIA from its github repository

library(devtools)
devtools::install_github("pouryany/CADIA")

File guides

The github page of CADIA packages contains associated the code and documentation. The folder data-raw and its subfolders contains relevant data and code accompanying the package and the original manuscripts. In particular, KEGGPATHS contains the raw unprocessed XML files of the original pathways. The TestCodes folder contains the test cases and codes for the accompaniying publications. Additional instructions and guides regarding the specific test cases and codes are provided in the GuideMe.R document inthe data-raw folder.

Enrichment Analysis with CADIA

The current version of CADIA works by taking two inputs of differentially expressed genes (DEG) and the gene universe. The output of CADIA is a list of KEGG pathways along with a group of p-values that describe the association of the pathway with a the DEG.

This document contains a test case of CADIA using an ovarian cancer dataset by Bowtell and colleagues, with 60 High-grade serous ovarian cancer and 30 Low malignant potential tumors. This data is available from NCBI-GEO portal through accession code GSE12172 [@ang08]. The datasets platform us affymetrix HG-U133b.

Finding differentially expressed genes

NCBI-GEO portal provides and R-code generation tool for differential expression analysis of its datasets. The below code shows the automatically generated script for differential expression analysis of the GSE12172 dataset. In general, it uses the limma package to contrast the normalized expressions of two predefined subsets of the samples. For running the next code chunk you would need to have these packages installed: Biobase, GEOquery, and limma. If you already have a limma topTable object you may skip to the next section where we use a predefined table of this analysis.

# The following lines are based on the code generated from NCBI GEO portal
# Version info: R 3.2.3, Biobase 2.30.0, GEOquery 2.40.0, limma 3.26.8
# R scripts generated  Tue Feb 20 21:18:14 EST 2018

################################################################
#   Differential expression analysis with limma
library(Biobase)
library(GEOquery)
library(limma)
library(KEGGgraph)


# load series and platform data from GEO
# 30 ovarian LMP vs 60 ovarian cancer
gset <- getGEO("GSE12172", GSEMatrix =TRUE, AnnotGPL=TRUE)
if (length(gset) > 1) idx <- grep("GPL570", attr(gset, "names")) else idx <- 1
gset <- gset[[idx]]

# make proper column names to match toptable
fvarLabels(gset) <- make.names(fvarLabels(gset))

# group names for all samples, 1 For cancer and 0 for non-cancer
gsms <- paste0("11111010110111001100111111110101001010101111101111",
               "1101111010011110011110000011111010010111")
sml <- c()
for (i in 1:nchar(gsms)) { sml[i] <- substr(gsms,i,i) }

# log2 transform
ex   <- exprs(gset)
qx   <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
    (qx[6]-qx[1] > 50 && qx[2] > 0) ||
    (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
if (LogC) { ex[which(ex <= 0)] <- NaN
exprs(gset) <- log2(ex) }

# set up the data and proceed with analysis
sml              <- paste("G", sml, sep="")    # set group names
fl               <- as.factor(sml)
gset$description <- fl
design           <- model.matrix(~ description + 0, gset)
colnames(design) <- levels(fl)
fit              <- lmFit(gset, design)
cont.matrix      <- makeContrasts(G1-G0, levels=design)
fit2             <- contrasts.fit(fit, cont.matrix)
fit2             <- eBayes(fit2, 0.01)
tT               <- topTable(fit2, adjust="fdr", sort.by="B", number=Inf)

Preparing differentially expressed genes

Here, we will an topTable object that was produced from an ovarian cancer data set from NCBI GEO (GSE12172) [@ang08]. The instructions for producing the topTable object are provided in the previous sections. The next step for enrichment analysis is to identify a subset of differentially expressed genes. The specific platform of this data has multiple probes for each annotated genes. Also, some of the probes do not correspond to any annotated gene. After appropriate filtering for such instances, we use an adjusted p-value and fold-change criteria to identify the subset of DEG. Note that the end results of this analysis step are two lists tT.all.names and tT.de.names.

library(CADIA)

data(tT)
# Filter out unannotated and duplicated genes
tT.filter    <- tT[!is.na(tT$Gene.ID),]
tT.filter    <- tT.filter[!duplicated(tT.filter$Gene.ID),]

# Select a p-value and fold-change cut-off for differential expressions
tT.deGenes   <- tT.filter[tT.filter$adj.P.Val < 0.05, ]
tT.deGenes   <- tT.deGenes[abs(tT.deGenes$logFC) >1,]

# Parameters needed for CADIA. The name of all genes and DE genes.
tT.all.names <- as.vector(tT.filter$Gene.ID)
tT.de.names  <- as.vector(tT.deGenes$Gene.ID)

Pathway Enrichment Analysis with CADIA

The current version of CADIA only works with ENTREZ IDs. To run the enrichment analysis, simply provide the function causalDisturbance() with minimum of two inputs The first argument is the ENTREZ ID of the DEG. The second argument is the ENTREZ ID of all genes from the experiment. If your data has an output other than ENTREZ IDs refer to the subsection at the end of this sections.

Note that, as described in the manuscript, CADIA uses random sampling for calculating the outputs. To ensure the consistency and reproducibility of your results, always the random seed before running the method.

set.seed(1)
cadia.res <- causalDisturbance(tT.de.names,tT.all.names,iter = 10000)

The following are the additional arguments for causalDisturbance(). iter denotes the number of iterations for random bootstrap sampling (preferrably larger than 2000). alpha is the dampening factor of Source-Sink centrality (described in the manuscript, defaul = 0.1), beta is the relative importance of sink centrality compared to the source centrality (default = 1). statEval denotes whether to use a product or summation for calculating the aggeregate centrality of the perturbations (default = 1, product-based aggregation). fdrMethod is the choice multiple hypothesis correction method (default = "B", see p.adjust function documentaion from stats package for options). verbose denotes whether to generate progress report during calculations.

head(cadia.res)

cadia.res is output table. The column cadia is the FDR corrected causal disturbance, which indicates the statistical significance of a pathway enrichment. The output table is sorted based on this value and can be taken as the enrichment resutls. Additional outputs of the method denote different aspects of the analysis. P_ORA is the p-value of over-representation analysis through hypergeometric test. P_SSC is the topological evidence, p-value of Source-Sink Centrality as described in the original manuscript. The column Causal Disturbance is the combined P_SSC and P_ORA using Fisher's method. The column ORA is the adjusted P_ORA. KEGGID is the associated ID of a pathway in the KEGG Database.

cadia.res is a tibble and, thus, can be manipulated using tibble related methods. For example, to create a table with only pathways with cadia < 0.05.

library(dplyr)
top.res <- cadia.res %>% dplyr::filter(.,cadia < 0.05) %>% select(., Name, KEGGID)
head(top.res)

Translating and preparing inputs for CADIA

As mentioned, the current version of CADIA only works with ENTREZ IDs. If the results of your analysis is in some other format, e.g. symbol and ENSEMBEL, we suggest using the functions of clusterProfiler package for appropriate translations [@yu12]. To do this, you would need appropriate library installations. The code below depicts a procedure for translating between different formats in the hypothetical case where your experimental results are in gene symbol notation.

library(clusterProfiler)
library(org.Hs.eg.db)
library(CADIA)
# the object geneList contains a list of all genes in the universe

# the object deGenes contains a list of differentially expressed genes.


gene.df    <- bitr(geneList, fromType = "SYMBOL",
                   toType = c("ENTREZID","ENSEMBL"),
                   OrgDb = org.Hs.eg.db)
deGenes.df <- bitr(deGenes, fromType = "SYMBOL",
                   toType = c("ENTREZID","ENSEMBL"),
                   OrgDb = org.Hs.eg.db)

set.seed(1)
cadia.res  <- CADIA::causalDisturbance(deGenes.df$ENTREZID,gene.df$ENTREZID,
                                       iter = 5000)

Additional functions in CADIA package

The package CADIA provides additional functionality for pathway enrichment analysis and other applications. This section provides a brief overview of the functions.

Internal data

The graph objects of KEGG pathways that are used in CADIA can be accessed using the internal data of package. The data pathways.collection contains the processed pathway graphs in the graphNEL format.

library(CADIA)
data("pathways.collection")
pathways.collection[[1]]

If you are interested in getting a list of pathways that are analyzed by CADIA you can use the built-in function cadia.paths().

head(CADIA::cadia.paths())

One might be intresed to retrieve the differentially expressed genes associated with the significantly enriched pathways, or a subset of them. The function geneReport() faciliates this operation by returning a list containing the pathways and the list of their DEG (in ENTREZ format) concatenated in the rows. See the following.

reports <- geneReport(tT.de.names,top.res$KEGGID)
rownames(reports) <- NULL
head(reports)

The graph analysis functionality of CADIA is implemented separately for those who wish to utilize them in different lines of research. The function source.sink.centrality() implements the concept of Source/Sink Centrality as describe in the original manuscript. This function returns a matrix whose individual elements corresponding to row i and column j can be interpreted as the influence of node i on node j. One can alternatively calculate the centrality of individual nodes by using the function rowSums()

library(CADIA)
library(graph)
data("pathways.collection")
test.graph    <- pathways.collection[[1]]
test.matrix   <- as(test.graph,"matrix")
ssc.influence <- CADIA::source.sink.centrality(test.matrix,alpha = 0.1, beta =1)
head(rowSums(ssc.influence))

The notion of the aggregate score in the original paper can be used in applications where one is interested in computing a centrality score for a subset of nodes. The following provides a showcase of how to calculate this notion of subgraph centrality using the built-in function pathSampler().

test.nodes <- graph::nodes(test.graph)
set.seed(1)
subs.nodes <- sample(test.nodes,10, replace = F)

set.seed(1)
subs.prob  <- pathSampler(inputGraph = test.graph ,deKID = subs.nodes, 
                          iterationNo = 10000, alpha = 0.1, beta = 1,
                          statEval = 1  )
subs.prob

# sub.prob is a probability, one can turn it into a score as following

-log(subs.prob,base = 10)

References

pouryany/pathwayDisturbance documentation built on Sept. 1, 2020, 12:13 a.m.