inst/doc/ActivePathways-vignette.R

In ActivePathways: Integrative Pathway Enrichment Analysis of Multivariate Omics Data

## ----echo=FALSE---------------------------------------------------------------
knitr::opts_chunk$set(warning=FALSE, 
                      message=FALSE, 
                      width=500)
options(max.print=35)
library("ggplot2")
library("data.table")

## -----------------------------------------------------------------------------
scores <- read.table(
system.file('extdata', 'Adenocarcinoma_scores_subset.tsv', package = 'ActivePathways'), 
header = TRUE, sep = '\t', row.names = 'Gene')
scores <- as.matrix(scores)
scores

## -----------------------------------------------------------------------------
scores[is.na(scores)] <- 1

## -----------------------------------------------------------------------------
library(ActivePathways)
gmt_file <- system.file('extdata', 'hsapiens_REAC_subset.gmt', package = 'ActivePathways')
ActivePathways(scores, gmt_file)

## -----------------------------------------------------------------------------
fname_data_matrix <- system.file('extdata', 
		'Differential_expression_rna_protein.tsv',
		package = 'ActivePathways')

pvals_FCs <- read.table(fname_data_matrix, header = TRUE, sep = '\t')

example_genes <- c('ACTN4','PIK3R4','PPIL1','NELFE','LUZP1','ITGB2')
pvals_FCs[pvals_FCs$gene %in% example_genes,]

## -----------------------------------------------------------------------------
pval_matrix <- data.frame(
		row.names = pvals_FCs$gene, 
		rna = pvals_FCs$rna_pval, 
		protein = pvals_FCs$protein_pval)

pval_matrix <- as.matrix(pval_matrix)
pval_matrix[is.na(pval_matrix)] <- 1

pval_matrix[example_genes,]

## -----------------------------------------------------------------------------
dir_matrix <- data.frame(
		row.names = pvals_FCs$gene, 
		rna = pvals_FCs$rna_log2fc, 
		protein = pvals_FCs$protein_log2fc)

dir_matrix <- as.matrix(dir_matrix)
dir_matrix <- sign(dir_matrix)
dir_matrix[is.na(dir_matrix)] <- 0

dir_matrix[example_genes,]

## -----------------------------------------------------------------------------
constraints_vector <- c(1,1)
constraints_vector

# constraints_vector <- c(1,-1)

## -----------------------------------------------------------------------------
directional_merged_pvals <- merge_p_values(pval_matrix, 
		method = "DPM", dir_matrix, constraints_vector)

merged_pvals <- merge_p_values(pval_matrix, method = "Brown")

sort(merged_pvals)[1:5]
sort(directional_merged_pvals)[1:5]

## -----------------------------------------------------------------------------
pvals_FCs[pvals_FCs$gene == "PIK3R4",]

pval_matrix["PIK3R4",]
dir_matrix["PIK3R4",]

merged_pvals["PIK3R4"]
directional_merged_pvals["PIK3R4"]

## -----------------------------------------------------------------------------
lineplot_df <- data.frame(original = -log10(merged_pvals),
			  modified = -log10(directional_merged_pvals))

ggplot(lineplot_df) +
	geom_point(size = 2.4, shape = 19,
		aes(original, modified,
		    color = ifelse(original <= -log10(0.05),"gray",
                                    ifelse(modified > -log10(0.05),"#1F449C","#F05039")))) +
	labs(title = "",
		 x ="Merged -log10(P)",
		 y = "Directional Merged -log10(P)") + 
            geom_hline(yintercept = 1.301, linetype = "dashed",
		       col = 'black', size = 0.5) +
            geom_vline(xintercept = 1.301, linetype = "dashed",
		       col = "black", size = 0.5) + 
            geom_abline(size = 0.5, slope = 1,intercept = 0) +
	    scale_color_identity()

## -----------------------------------------------------------------------------
constraints_vector <- c(-1,1)

constraints_vector <- c(1,-1)

## -----------------------------------------------------------------------------
constraints_vector <- c(0,0)

## -----------------------------------------------------------------------------
constraints_vector <- c(1,-1)

## -----------------------------------------------------------------------------
constraints_vector <- c(1,1,-1)

## -----------------------------------------------------------------------------
fname_GMT2 <- system.file("extdata", "hsapiens_REAC_subset2.gmt", 
		package = "ActivePathways")

## -----------------------------------------------------------------------------
enriched_pathways <- ActivePathways(
		pval_matrix, gmt = fname_GMT2, cytoscape_file_tag = "Original_")

constraints_vector <- c(1,1)
constraints_vector
dir_matrix[example_genes,]

enriched_pathways_directional <- ActivePathways(
		pval_matrix, gmt = fname_GMT2, cytoscape_file_tag = "Directional_", merge_method = "DPM",
		scores_direction = dir_matrix, constraints_vector = constraints_vector)

## -----------------------------------------------------------------------------
pathways_lost_in_directional_integration = 
		setdiff(enriched_pathways$term_id, enriched_pathways_directional$term_id)

pathways_lost_in_directional_integration
enriched_pathways[enriched_pathways$term_id %in% pathways_lost_in_directional_integration,] 

## -----------------------------------------------------------------------------
wnt_pathway_id <- "REAC:R-HSA-3858494"

enriched_pathway_genes <- unlist(
		enriched_pathways[enriched_pathways$term_id == wnt_pathway_id,]$overlap)
enriched_pathway_genes

## -----------------------------------------------------------------------------
pathway_gene_pvals = pval_matrix[enriched_pathway_genes,]
pathway_gene_directions = dir_matrix[enriched_pathway_genes,]

directional_conflict_genes = names(which(
		pathway_gene_directions[,1] != pathway_gene_directions[,2] &
		pathway_gene_directions[,1] != 0 & pathway_gene_directions[,2] != 0))

pathway_gene_pvals[directional_conflict_genes,]
pathway_gene_directions[directional_conflict_genes,]
length(directional_conflict_genes)

## -----------------------------------------------------------------------------
nrow(ActivePathways(scores, gmt_file, significant = 0.05))
nrow(ActivePathways(scores, gmt_file, significant = 0.1))

## -----------------------------------------------------------------------------
gmt <- read.GMT(gmt_file)
names(gmt[[1]])

# Pretty-print the GMT
gmt[1:3]

# Look at the genes annotated to the first term
gmt[[1]]$genes

# Get the full name of Reactome pathway 2424491
gmt$`REAC:2424491`$name

## -----------------------------------------------------------------------------
gmt <- Filter(function(term) length(term$genes) >= 10, gmt)
gmt <- Filter(function(term) length(term$genes) <= 500, gmt)

## -----------------------------------------------------------------------------
ActivePathways(scores, gmt)

## -----------------------------------------------------------------------------
ActivePathways(scores, gmt_file, geneset_filter = c(10, 500))

## ----eval=FALSE---------------------------------------------------------------
#  write.GMT(gmt, 'hsapiens_REAC_subset_filtered.gmt')

## -----------------------------------------------------------------------------
background <- makeBackground(gmt)
background <- background[background != 'TP53']
ActivePathways(scores, gmt_file, background = background)

## -----------------------------------------------------------------------------
sort(merge_p_values(scores, 'Fisher'))[1:5]
sort(merge_p_values(scores, 'Brown'))[1:5]

sort(merge_p_values(scores, 'Stouffer'))[1:5]
sort(merge_p_values(scores, 'Strube'))[1:5]

## -----------------------------------------------------------------------------
scores2 <- cbind(scores[, 'CDS'], merge_p_values(scores[, c('X3UTR', 'X5UTR', 'promCore')], 'Brown'))
colnames(scores2) <- c('CDS', 'non_coding')
scores[c(2179, 1760),]
scores2[c(2179, 1760),]

ActivePathways(scores, gmt_file)
ActivePathways(scores2, gmt_file)

## -----------------------------------------------------------------------------

nrow(ActivePathways(scores, gmt_file))
nrow(ActivePathways(scores, gmt_file, cutoff = 0.01))


## -----------------------------------------------------------------------------

nrow(ActivePathways(scores, gmt_file))
nrow(ActivePathways(scores, gmt_file, correction_method = 'none'))


## -----------------------------------------------------------------------------
res <- ActivePathways(scores, gmt_file)
res

## -----------------------------------------------------------------------------
res$overlap[1:3]

## -----------------------------------------------------------------------------
unlist(res[res$term_id == "REAC:422475","evidence"])

## ----eval = FALSE-------------------------------------------------------------
#  result_file <- paste('ActivePathways_results.csv', sep = '/')
#  export_as_CSV (res, result_file) # remove comment to run
#  read.csv(result_file, stringsAsFactors = F)[1:3,]

## ----eval=FALSE---------------------------------------------------------------
#  result_file <- paste('ActivePathways_results2.txt', sep = '/')
#  data.table::fwrite(res, result_file, sep = '\t', sep2 = c('', ',', ''))
#  cat(paste(readLines(result_file)[1:2], collapse = '\n'))

## ----eval=FALSE---------------------------------------------------------------
#  res <- ActivePathways(scores, gmt_file, cytoscape_file_tag = "enrichmentMap__")

## -----------------------------------------------------------------------------
files <- c(system.file('extdata', 'enrichmentMap__pathways.txt', package='ActivePathways'),
           system.file('extdata', 'enrichmentMap__subgroups.txt', package='ActivePathways'),
           system.file('extdata', 'enrichmentMap__pathways.gmt', package='ActivePathways'),
           system.file('extdata', 'enrichmentMap__legend.pdf', package='ActivePathways'))

## ----eval=FALSE---------------------------------------------------------------
#  gmt_file <- system.file('extdata', 'hsapiens_REAC_subset.gmt', package = 'ActivePathways')
#  scores_file <- system.file('extdata', 'Adenocarcinoma_scores_subset.tsv', package = 'ActivePathways')
#  
#  scores <- read.table(scores_file, header = TRUE, sep = '\t', row.names = 'Gene')
#  scores <- as.matrix(scores)
#  scores[is.na(scores)] <- 1
#  
#  res <- ActivePathways(scores, gmt_file, cytoscape_file_tag = "enrichmentMap__")

## -----------------------------------------------------------------------------
cat(paste(readLines(files[1])[1:5], collapse='\n'))
cat(paste(readLines(files[2])[1:5], collapse='\n'))
cat(paste(readLines(files[3])[18:19], collapse='\n'))

## -----------------------------------------------------------------------------
res <- ActivePathways(scores, gmt_file, cytoscape_file_tag = "enrichmentMap__", color_palette = "Pastel1")

## -----------------------------------------------------------------------------
res <- ActivePathways(scores, gmt_file, cytoscape_file_tag = "enrichmentMap__", custom_colors = c("violet","green","orange","red"))