Pathway analysis using SBGNview gene set and visualization
Xiaoxi Dong, Weijun Luo
29 October, 2019
1 Install R packages
1.1 Install SBGNview
Install SBGNview through Bioconductor
BiocManager::install(c("SBGNview"))1.2 Install package ‘gage’ for pathway enrichment analysis
Install gage through Bioconductor
BiocManager::install(c("gage"))2 IFNg dataset
In this example, we analyze a RNA-SEQ dataset of wild type mice and IFNg KO mice. It contains normalized RNA-seq gene expression data described in this study: Greer, Renee L., Xiaoxi Dong, et al. “Akkermansia muciniphila mediates negative effects of IFNg on glucose metabolism.” Nature communications 2016.
2.1 Load RNA-seq data
The RNA abundance data was quantile normalized and log2 transformed, stored in a “SummarizedExperiment” object
library(SBGNview)
library(SummarizedExperiment)
library(SBGNview.data)
data("IFNg")
count.data <- assays(IFNg)$counts
wt <- colnames(IFNg)[IFNg$group == "wt"]
ko <- colnames(IFNg)[IFNg$group == "ko"]
head(count.data)2.2 Pathway enrichment analysis: run SBGNview gene set
2.2.1 Load gene list: mouse ensemble IDs.
ensemble.to.pathway <- getMolList(
database = "pathwayCommons"
,mol.list.ID.type = "ENSEMBL"
,org = "mmu"
,cpd.or.gene = "gene"
,output.pathway.name = TRUE
)
head(ensemble.to.pathway[[2]])2.2.2 Run gage
gage is an R package for pathway enrichment analysis.
library(gage)
degs <- gage(exprs = count.data
,gsets = ensemble.to.pathway
,ref = which(colnames(count.data) %in% wt)
,samp = which(colnames(count.data) %in% ko)
,compare = "as.group"
)
head(degs$greater)[,c(1,4:5)]
head(degs$less)
down.pathways <- degs$less
head(down.pathways)2.3 Visualize fold change on pathway maps
2.3.1 Generate fold change data for visualization.
The abundance table was log2 transformed. Here we calculate the fold change of KO group v.s. WT group.
mean.wt <- apply(count.data[,wt] ,1 ,"mean")
head(mean.wt)
mean.ko <- apply(count.data[,ko],1,"mean")
head(mean.ko)
ensemble.to.koVsWt <- mean.ko - mean.wt
head(ensemble.to.koVsWt)2.3.2 Visualize fold change data
down.pathways <- do.call(rbind,strsplit(row.names(down.pathways),"::"))[,1]
head(down.pathways)
sbgnview.obj <- SBGNview(
gene.data = ensemble.to.koVsWt,
gene.id.type = "ENSEMBL",
input.sbgn = down.pathways[1],
output.file = "ifn.sbgnview.less",
show.pathway.name = TRUE,
max.gene.value = 2,
min.gene.value = -2,
mid.gene.value = 0,
node.sum = "mean",
label.spliting.string = c("-","_"," "),
output.format = c("png"),
inhibition.edge.end.shift = 1,
font.size = 2,
logic.node.font.scale = 6,
font.size.scale.compartment = 1.5,
org = "mmu",
text.length.factor.macromolecule = 0.8,
text.length.factor.compartment = 2,
text.length.factor.complex = 3,
if.scale.compartment.font.size = TRUE,
node.width.adjust.factor.compartment = 0.058
)
sbgnview.obj
Figure 2.1: SBGNview
2.3.3 Visualize RNA abundance data: SummarizedExperiment object as input.
data("cancer.ds")
sbgnview.obj <- SBGNview(
gene.data = cancer.ds,
gene.id.type = "ENTREZID",
input.sbgn = "R-HSA-877300",
output.file = "demo.SummarizedExperiment",
show.pathway.name = TRUE,
max.gene.value = 1,
min.gene.value = -1,
mid.gene.value = 0,
node.sum = "mean",
label.spliting.string = c("-","_"," "),
output.format = c("png"),
inhibition.edge.end.shift = 1,
font.size = 2,
logic.node.font.scale = 6,
font.size.scale.compartment = 1.5,
org = "hsa",
text.length.factor.macromolecule = 0.8,
text.length.factor.compartment = 2,
text.length.factor.complex = 3,
if.scale.compartment.font.size = TRUE,
node.width.adjust.factor.compartment = 0.058
)
sbgnview.obj3 Session Info
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] gage_2.36.0 SummarizedExperiment_1.16.0
## [3] DelayedArray_0.12.0 BiocParallel_1.20.0
## [5] matrixStats_0.55.0 GenomicRanges_1.38.0
## [7] GenomeInfoDb_1.22.0 SBGNview_1.0.0
## [9] SBGNview.data_0.99.18 pathview_1.26.0
## [11] org.Hs.eg.db_3.10.0 AnnotationDbi_1.48.0
## [13] IRanges_2.20.0 S4Vectors_0.24.0
## [15] Biobase_2.46.0 BiocGenerics_0.32.0
## [17] knitr_1.25
##
## loaded via a namespace (and not attached):
## [1] KEGGREST_1.26.0 xfun_0.10 lattice_0.20-38
## [4] vctrs_0.2.0 htmltools_0.4.0 yaml_2.2.0
## [7] blob_1.2.0 XML_3.98-1.20 rlang_0.4.1
## [10] pillar_1.4.2 DBI_1.0.0 Rgraphviz_2.30.0
## [13] bit64_0.9-7 GenomeInfoDbData_1.2.2 stringr_1.4.0
## [16] zlibbioc_1.32.0 Biostrings_2.54.0 memoise_1.1.0
## [19] evaluate_0.14 rsvg_1.3 gbRd_0.4-11
## [22] highr_0.8 Rcpp_1.0.2 backports_1.1.5
## [25] graph_1.64.0 XVector_0.26.0 bit_1.1-14
## [28] png_0.1-7 digest_0.6.22 stringi_1.4.3
## [31] bookdown_0.14 grid_3.6.1 bibtex_0.4.2
## [34] Rdpack_0.11-0 tools_3.6.1 bitops_1.0-6
## [37] magrittr_1.5 RCurl_1.95-4.12 tibble_2.1.3
## [40] RSQLite_2.1.2 crayon_1.3.4 pkgconfig_2.0.3
## [43] zeallot_0.1.0 Matrix_1.2-17 xml2_1.2.2
## [46] KEGGgraph_1.46.0 rmarkdown_1.16 httr_1.4.1
## [49] R6_2.4.0 igraph_1.2.4.1 compiler_3.6.1