PancCanNet - Workflow 1

General instructions (read before running the code snippets)

In this workflow, we will run identify affected pathways in two pancreatic cancer subtypes, and visualize the data on the pathway models.

• The script contains several code snippets which should be run one after the other.
• Make sure all the required packages are installed beforehand (BiocManager::install(…)).
• Make sure you have Cytoscape installed (version 3.8.0+) and running before you start running the script.

Setup

Loading libraries

options(connectionObserver = NULL)
# uncomment if you need to install the libraries
#BiocManager::install(c("GO.db","DO.db","readr","here","dplyr","tibble","clusterProfiler","org.Hs.eg.db","tidyverse","rWikiPathways","RCy3","RColorBrewer","EnhancedVolcano","reshape","data.table","svglite","VennDiagram","BridgeDbR","rstudioapi"))


library(readr)
library(here)
library(dplyr)
library(tibble)
library(clusterProfiler)
library(org.Hs.eg.db)
library(tidyverse)
library(rWikiPathways)
library(RCy3)
library(RColorBrewer)
library(EnhancedVolcano)
library(reshape)
library(data.table)
library(svglite)
library(VennDiagram)
library(BridgeDbR)
library(rstudioapi)

setwd(dirname(getActiveDocumentContext()$path))

Differential expression data visualization

We will use a publicly available dataset, which identified two different subtypes of pancreatic cancer - classical and basal (GSE71729).

First, let’s import the data and use a volcano plot to visualize the result of the differential gene expression analysis result, and use a Venn diagram to study how many differentially expressed genes are shared between the subtypes.

Interpretation

• Q1: Are there more up- or down-regulated genes? What cutoffs have been used for pvalue and log2FC (see line 63)?
• Q2: What are the most significantly changed genes? Do a quick Pubmed search to see if they have a known link to pancreatic cancer.
• Q3: How many differentially expressed genes are shared between the two subtypes?

Pathway enrichment analysis

We will perform pathway enrichment with the gene sets of all pathway models in WikiPathways (human only).

gmt <- rWikiPathways::downloadPathwayArchive(organism = "Homo sapiens",format = "gmt")
wp2gene <- readPathwayGMT(gmt)
wpid2gene <- wp2gene %>% dplyr::select(wpid,gene) #TERM2GENE
wpid2name <- wp2gene %>% dplyr::select(wpid,name) #TERM2NAME

bkgd.genes <- unique(data.panc[,c(1,2)])

Interpretation

• Q4: How many pathway models are in the collection? (wpid2name)

The clusterProfiler R-package is used to perform overrepresentation analysis (ORA). The function can be easily replaced to use other enrichment methods (GSEA / rSEA / etc). We will run the analysis separately for basal and classical subtype.

ewp.basal <- clusterProfiler::enricher(
  deg.basal$GeneId,
  universe = bkgd.genes$GeneId,
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.02,
  TERM2GENE = wpid2gene,
  TERM2NAME = wpid2name)
ewp.basal.res <- as.data.frame(ewp.basal) 

# number of genes measured in pathways
length(ewp.basal@universe)
# number of DEG in pathways
length(deg.basal$GeneId[deg.basal$GeneId %in% unique(wp2gene$gene)])

num.pathways.basal <- dim(ewp.basal.res)[1]

# export enrichment result
png('output1/basal_barplot.png', width = 1200, height=1000)
ggplot(ewp.basal[1:num.pathways.basal], aes(x=reorder(Description, -pvalue), y=Count)) +
  geom_bar(stat ="identity", fill="#BA8CD7") +
  coord_flip() +
  labs(x="", y="Basal DEG gene count", fill="") +
  theme(axis.text=element_text(size=25)) + 
  theme(legend.position="none")
dev.off()
write.table(ewp.basal.res, file="output1/basal_enrich_res.txt", sep="\t", quote=FALSE, row.names = FALSE)

Interpretation

• Q5: How many pathways are altered in the basal subtype and how do they link to pancreatic cancer (expected or unexpected)?

ewp.classical <- clusterProfiler::enricher(
  deg.classical$GeneId,
  universe = bkgd.genes$GeneId,
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.2,
  TERM2GENE = wpid2gene,
  TERM2NAME = wpid2name)
ewp.classical.res <- as.data.frame(ewp.classical) 

# number of genes measured in pathways
length(ewp.classical@universe)
# number of DEG in pathways
length(deg.classical$GeneId[deg.classical$GeneId %in% unique(wp2gene$gene)])

num.pathways.classical <- dim(ewp.classical.res)[1]

# export enrichment result
png('output1/classical_barplot.png', width = 1200, height=1000)
ggplot(ewp.classical[1:num.pathways.classical], aes(x=reorder(Description, -pvalue), y=Count)) +
  geom_bar(stat ="identity", fill="#BA8CD7") +
  coord_flip() +
  labs(x="", y="Classical DEG gene count", fill="") +
  theme(axis.text=element_text(size=25)) + 
  theme(legend.position="none")
dev.off()
write.table(ewp.classical.res, file="output1/classical_enrich_res.txt", sep="\t", quote=FALSE, row.names = FALSE)

Interpretation

• Q6: How many pathways are altered in the classical subtype and how do they link to pancreatic cancer (expected or unexpected)?

Venn diagrams

venn.diagram(x = list(ewp.basal.res$ID, ewp.classical.res$ID),
  category.names = c("Basal" , "Classical"),
  filename = 'output1/venn_pathways.png',
  output=TRUE,
  col=c("#440154ff", '#21908dff'),
  fill = c(alpha("#440154ff",0.3), alpha('#21908dff',0.3)),
  cex = 1.5,
)

Interpretation

• Q7: How many pathways are altered in both subtypes?

Pathway visualization

The pathways can then be visualized with the gene expression data as shown with the “Vitamin D receptor pathway” (WP2877) pathway from WikiPathways. The pathway was altered in both subtypes.

You can easily replace the WikiPathways identifier in line 209 with another pathway of interest and visualize the data on that pathway.

RCy3::cytoscapePing()
RCy3::installApp(c("wikipathways","CyTargetLinker"))

RCy3::commandsRun('wikipathways import-as-pathway id=WP2877') 
toggleGraphicsDetails()
RCy3::mapTableColumn("Ensembl", "Homo sapiens", "Ensembl", "Entrez Gene")
loadTableData(data.panc, data.key.column = "GeneId", table.key.column = "Entrez Gene")

RCy3::installApp("enhancedGraphics")
RCy3::copyVisualStyle("WikiPathways", "my_style_heatmap")

RCy3::setNodeCustomHeatMapChart(c("B_logFC","C_logFC"), slot = 2, style.name = "my_style_heatmap", colors = c("#CC3300","#FFFFFF","#6699FF","#CCCCCC"))

RCy3::setVisualStyle("my_style_heatmap")

# Saving output
#png.file <- file.path("figures/PathwayVisualization.png")
#exportImage(png.file,'PNG', zoom = 500)
#cys.file <- file.path(getwd(),"output/PathwayVisualization.cys")
#saveSession(cys.file) 
#comment following line if you want to manipulate the visualization in Cytoscape
#RCy3::closeSession(save.before.closing = F)

Interpretation

• Q8: Do you see any differences between the two subtypes in the vitamin D receptor pathway? (left is basal and right is classical)

Pathway overlap visualization

There is often crosstalk and overlap between pathways enriched in gene expression analyses. The following step visualizes the overlap between the enriched pathways in a pathway-gene network.

The genes not present in any pathway are included in the visualization but can be filtered in a follow-up step if preferred.

pwy <- unique(ewp.basal.res[,c(1,2)])
colnames(pwy) <- c("id","label")
pwy$type <- 'pathway'
edges <- wpid2gene[wpid2gene$wpid %in% pwy$id,]
colnames(edges) <- c("source", "target")
genes <- unique(deg.basal)
colnames(genes) <- c("id","label")
genes <- transform(genes, id = as.character(id))
genes$type <- 'gene'
edges <- unique(edges[edges$target %in% genes$id,])
genes <- genes[genes$id %in% edges$target,]
nodes <- dplyr::bind_rows(genes, pwy)
rownames(nodes) <- NULL
createNetworkFromDataFrames(nodes=nodes,edges=edges,title="Pathway-Gene-Associations", collection="PathwayGeneCrosstalk")
loadTableData(data.panc, data.key.column = "GeneId", table.key.column = "id")
# Visual style
RCy3::copyVisualStyle("default","wp.vis")
RCy3::setNodeLabelMapping("label", style.name="wp.vis")
RCy3::lockNodeDimensions(TRUE, style.name="wp.vis")
RCy3::setNodeShapeMapping('type', c('gene','pathway'), c("ellipse","hexagon"), style.name="wp.vis")
RCy3::setNodeSizeMapping('type', c('gene','pathway'), c(40,25), mapping.type = "d", style.name = "wp.vis")
data.values<-c(-1,0,1) 
node.colors <- c(rev(brewer.pal(length(data.values), "RdBu")))
setNodeColorMapping("B_logFC", data.values, node.colors, default.color = "#99FF99", style.name = "wp.vis")
RCy3::setVisualStyle("wp.vis")
RCy3::toggleGraphicsDetails()

# Saving output
svg.file <- file.path(getwd(), "output1/PathwayCrosstalk.svg")
exportImage(svg.file,'SVG')
png.file <- file.path(getwd(), "output1/PathwayCrosstalk.png")
exportImage(png.file,'PNG', zoom = 500)
cys.file <- file.path(getwd(), "output1/PathwayCrosstalk.cys")
saveSession(cys.file) 
#comment following line if you want to manipulate the visualization in Cytoscape
#RCy3::closeSession(save.before.closing = F)

Interpretation

• Q9: Do the affected pathways share any genes?

Drug target information

Next, we will add information about known drug-target interactions for the genes in the affected pathways using information from DrugBank using the CyTargetLinker app.

We will show this for the classical subtype.

RCy3::cytoscapePing()
installApp('CyTargetLinker') 

pwy <- unique(ewp.classical.res[,c(1,2)])
colnames(pwy) <- c("id","label")
pwy$type <- 'pathway'
edges <- wpid2gene[wpid2gene$wpid %in% pwy$id,]
colnames(edges) <- c("source", "target")
genes <- unique(deg.basal)
colnames(genes) <- c("id","label")
genes <- transform(genes, id = as.character(id))
genes$type <- 'gene'
edges <- unique(edges[edges$target %in% genes$id,])
genes <- genes[genes$id %in% edges$target,]
nodes <- dplyr::bind_rows(genes, pwy)
rownames(nodes) <- NULL
createNetworkFromDataFrames(nodes=nodes,edges=edges,title="Pathway-Gene-Associations", collection="PathwayGeneCrosstalk")
loadTableData(data.panc, data.key.column = "GeneId", table.key.column = "id")
# Visual style
RCy3::copyVisualStyle("default","wp.vis")
RCy3::setNodeLabelMapping("label", style.name="wp.vis")
RCy3::lockNodeDimensions(TRUE, style.name="wp.vis")
RCy3::setNodeShapeMapping('type', c('gene','pathway'), c("ellipse","hexagon"), style.name="wp.vis")
RCy3::setNodeSizeMapping('type', c('gene','pathway'), c(40,25), mapping.type = "d", style.name = "wp.vis")
data.values<-c(-1,0,1) 
node.colors <- c(rev(brewer.pal(length(data.values), "RdBu")))
setNodeColorMapping("C_logFC", data.values, node.colors, default.color = "#99FF99", style.name = "wp.vis")
RCy3::setVisualStyle("wp.vis")
RCy3::toggleGraphicsDetails()

drugbank <- file.path(getwd(), "data/drugbank-5.1.0.xgmml")

# run CyTargetLinker
commandsRun(paste0('cytargetlinker extend idAttribute="id" linkSetFiles="', drugbank, '"'))
commandsRun('cytargetlinker applyLayout network="current"')
RCy3::setVisualStyle("wp.vis")

#let's change the visualization of the drugs in the network using the ByPass option
selected <- RCy3::selectNodes(nodes="drug", by.col = "CTL.Type")
RCy3::setNodeShapeBypass(node.names = selected$nodes, new.shapes = "Triangle")
RCy3::setNodeColorBypass(node.names = selected$nodes, "#FFFFCE")
RCy3::setNodeBorderColorBypass(node.names = selected$nodes, "#000000")
RCy3::setNodeBorderWidthBypass(node.names = selected$nodes, 4)
RCy3::clearSelection()
RCy3::toggleGraphicsDetails()

Interpretation

• Q10: Are there proteins targeted by many known drugs? Can you find supporting literature regarding these genes or drugs in the context of pancreatic cancer?