if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("roastgsa")
This vignette explains broadly the main functions for applying roastgsa
in RNA-seq data. A more exhaustive example to explore the roastgsa
functionality is presented in the “roastgsa vignette (main)”. All the analyses
explained in the main vignette can be reproduced for RNA-seq data, after
undertaking the steps covered here in the section
“Data normalization and filtering”.
library( roastgsa )
We consider the first dataset available in the tcga
compendium
from the GSEABenchmarkeR
package [1], which consists of a RNA-seq study
with 19 tumor Bladder Urothelial Carcinoma samples and 19 adjacent
healthy tissues.
#library(GSEABenchmarkeR)
#tcga <- loadEData("tcga", nr.datasets=1,cache = TRUE)
#ysel <- assays(tcga[[1]])$expr
#fd <- rowData(tcga[[1]])
#pd <- colData(tcga[[1]])
data(fd.tcga)
data(pd.tcga)
data(expr.tcga)
fd <- fd.tcga
ysel <- expr.tcga
pd <- pd.tcga
N <- ncol(ysel)
head(pd)
## DataFrame with 6 rows and 4 columns
## sample type GROUP
## <character> <factor> <numeric>
## TCGA-K4-A3WV-01A-11R-A22U-07 TCGA-K4-A3WV-01A-11R.. BLCA 1
## TCGA-BT-A20W-01A-21R-A14Y-07 TCGA-BT-A20W-01A-21R.. BLCA 1
## TCGA-K4-A5RI-01A-11R-A28M-07 TCGA-K4-A5RI-01A-11R.. BLCA 1
## TCGA-BT-A20N-01A-11R-A14Y-07 TCGA-BT-A20N-01A-11R.. BLCA 1
## TCGA-BL-A13J-01A-11R-A277-07 TCGA-BL-A13J-01A-11R.. BLCA 1
## TCGA-BT-A20U-01A-11R-A14Y-07 TCGA-BT-A20U-01A-11R.. BLCA 1
## BLOCK
## <character>
## TCGA-K4-A3WV-01A-11R-A22U-07 TCGA-K4-A3WV
## TCGA-BT-A20W-01A-21R-A14Y-07 TCGA-BT-A20W
## TCGA-K4-A5RI-01A-11R-A28M-07 TCGA-K4-A5RI
## TCGA-BT-A20N-01A-11R-A14Y-07 TCGA-BT-A20N
## TCGA-BL-A13J-01A-11R-A277-07 TCGA-BL-A13J
## TCGA-BT-A20U-01A-11R-A14Y-07 TCGA-BT-A20U
cnames <- c("BLOCK","GROUP")
covar <- data.frame(pd[,cnames,drop=FALSE])
covar$GROUP <- as.factor(covar$GROUP)
colnames(covar) <- cnames
print(table(covar$GROUP))
##
## 0 1
## 19 19
To apply roastgsa
, the expression data should be approximately normally
distributed, at least in their univariate form. Depending on the user’s
preferred method for differential expression analysis, counts transformation
methods such as rlog
or vst
(DESeq2
) [2], zscoreDGE
(edgeR
) [3] or
voom
(limma
) [4], can be applied. In the paper we explored the type I
and type II errors when applying the rlog
or vst
transformation followed by
roastgsa
, showing both good control of type I errors and decent true
discovery rates. In the example presented here we transform the expression
data with vst
function from DESeq2
R package
library(DESeq2)
dds1 <- DESeqDataSetFromMatrix(countData=ysel,colData=pd,
design= ~ BLOCK + GROUP)
dds1 <- estimateSizeFactors(dds1)
ynorm <- assays(vst(dds1))[[1]]
colnames(ynorm) <- rownames(covar) <- paste0("s",1:ncol(ynorm))
Another key step before using the roastgsa methods for enrichment analysis
is to filter out low expressed genes, where coverage might be a
limitation for detecting true differentially expressed genes. For the
TCGA data considered here, the default filter employed by the authors when
loading the data was to exclude genes with cpm < 2 in more than half of
the samples. A short discussion about the relationship between gene coverage
and statistical power for the roastgsa
approach is available in our article
presenting the roastgsa
package.
threshLR <- 10
dim(ysel)
## [1] 3621 38
min(apply(ysel,1,mean))
## [1] 88.26316
We consider a classic repository of general biological functions for battery gene set analysis such as broad hallmarks [5]. The gene sets for human are saved within the roastgsa package and can be loaded by
data(hallmarks.hs)
head(names(hallmarks.hs))
## [1] "HALLMARK_TNFA_SIGNALING_VIA_NFKB" "HALLMARK_HYPOXIA"
## [3] "HALLMARK_CHOLESTEROL_HOMEOSTASIS" "HALLMARK_MITOTIC_SPINDLE"
## [5] "HALLMARK_WNT_BETA_CATENIN_SIGNALING" "HALLMARK_TGF_BETA_SIGNALING"
In this case, hallmarks.hs
contains gene symbols whereas the row
names for ynorm
are entrez identifiers. We can set the row names to
symbols, which in this case presents a one-to-one relationship
rownames(ynorm) <-fd[rownames(ynorm),1]
Other gene set databases that could be applied to these data for battery
testing are presented in the roastgsa
vignette (gene set collections).
The comparison of interest can be specified by a numeric vector with length matching the number of columns in the design.
form <- as.formula(paste0("~ ", paste0(cnames, collapse = "+")))
design <- model.matrix(form , data = covar)
terms <- colnames(design)
contrast <- rep(0, length(terms))
contrast[length(colnames(design))] <- 1
Below, there is the standard roastgsa
instruction (under competitive
testing) for maxmean
and mean
statistics.
fit.maxmean <- roastgsa(ynorm, form = form, covar = covar,
contrast = contrast, index = hallmarks.hs, nrot = 500,
mccores = 1, set.statistic = "maxmean",
self.contained = FALSE, executation.info = FALSE)
f1 <- fit.maxmean$res
rownames(f1) <- gsub("HALLMARK_","",rownames(f1))
head(f1)
## total_genes measured_genes est nes
## E2F_TARGETS 200 194 1.4694371 4.060285
## G2M_CHECKPOINT 200 188 1.1752646 3.989106
## MYOGENESIS 200 155 -0.8304752 -2.827713
## MYC_TARGETS_V2 58 58 1.0304122 2.705374
## UNFOLDED_PROTEIN_RESPONSE 113 104 0.3749170 2.697455
## MYC_TARGETS_V1 200 192 0.7727287 2.534479
## pval adj.pval
## E2F_TARGETS 0.001996008 0.03326680
## G2M_CHECKPOINT 0.001996008 0.03326680
## MYOGENESIS 0.001996008 0.03326680
## MYC_TARGETS_V2 0.013972056 0.09980040
## UNFOLDED_PROTEIN_RESPONSE 0.005988024 0.05988024
## MYC_TARGETS_V1 0.013972056 0.09980040
fit.mean <- roastgsa(ynorm, form = form, covar = covar,
contrast = contrast, index = hallmarks.hs, nrot = 500,
mccores = 1, set.statistic = "mean",
self.contained = FALSE, executation.info = FALSE)
f2 <- fit.mean$res
rownames(f2) <- gsub("HALLMARK_","",rownames(f2))
head(f2)
## total_genes measured_genes est nes
## E2F_TARGETS 200 194 1.1896796 2.863615
## G2M_CHECKPOINT 200 188 0.9287256 2.732674
## UNFOLDED_PROTEIN_RESPONSE 113 104 0.4270303 2.476257
## MYOGENESIS 200 155 -0.7076989 -2.361425
## UV_RESPONSE_DN 144 134 -0.5941011 -2.346723
## MYC_TARGETS_V2 58 58 0.9175999 2.233710
## pval adj.pval
## E2F_TARGETS 0.001996008 0.0499002
## G2M_CHECKPOINT 0.005988024 0.0748503
## UNFOLDED_PROTEIN_RESPONSE 0.005988024 0.0748503
## MYOGENESIS 0.001996008 0.0499002
## UV_RESPONSE_DN 0.009980040 0.0998004
## MYC_TARGETS_V2 0.033932136 0.1885119
Several graphics can be obtained to complement the table results in
f1
and f2
. Here we only show the heatmaps that summarize the
expression patterns obtained for all tested hallmarks. Full description
and usage of all graphical options available in the roastgsa
package
are considered in the roastgsa
vignette for arrays data and the
roastgsa
manual
hm1 <- heatmaprgsa_hm(fit.maxmean, ynorm, intvar = "GROUP", whplot = 1:50,
toplot = TRUE, pathwaylevel = TRUE, mycol = c("orange","green",
"white"), sample2zero = FALSE)
hm2 <- heatmaprgsa_hm(fit.mean, ynorm, intvar = "GROUP", whplot = 1:50,
toplot = TRUE, pathwaylevel = TRUE, mycol = c("orange","green",
"white"), sample2zero = FALSE)
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB 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
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DESeq2_1.46.0 SummarizedExperiment_1.36.0
## [3] Biobase_2.66.0 MatrixGenerics_1.18.0
## [5] matrixStats_1.4.1 GenomicRanges_1.58.0
## [7] GenomeInfoDb_1.42.0 IRanges_2.40.0
## [9] S4Vectors_0.44.0 BiocGenerics_0.52.0
## [11] roastgsa_1.4.0 knitr_1.48
## [13] BiocStyle_2.34.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.48 bslib_0.8.0
## [4] ggplot2_3.5.1 caTools_1.18.3 lattice_0.22-6
## [7] vctrs_0.6.5 tools_4.4.1 bitops_1.0-9
## [10] generics_0.1.3 parallel_4.4.1 tibble_3.2.1
## [13] fansi_1.0.6 highr_0.11 pkgconfig_2.0.3
## [16] Matrix_1.7-1 KernSmooth_2.23-24 RColorBrewer_1.1-3
## [19] lifecycle_1.0.4 GenomeInfoDbData_1.2.13 farver_2.1.2
## [22] compiler_4.4.1 gplots_3.2.0 tinytex_0.53
## [25] statmod_1.5.0 munsell_0.5.1 codetools_0.2-20
## [28] htmltools_0.5.8.1 sass_0.4.9 yaml_2.3.10
## [31] pillar_1.9.0 crayon_1.5.3 jquerylib_0.1.4
## [34] BiocParallel_1.40.0 cachem_1.1.0 DelayedArray_0.32.0
## [37] limma_3.62.0 magick_2.8.5 abind_1.4-8
## [40] gtools_3.9.5 locfit_1.5-9.10 tidyselect_1.2.1
## [43] digest_0.6.37 dplyr_1.1.4 bookdown_0.41
## [46] labeling_0.4.3 fastmap_1.2.0 grid_4.4.1
## [49] colorspace_2.1-1 cli_3.6.3 SparseArray_1.6.0
## [52] magrittr_2.0.3 S4Arrays_1.6.0 utf8_1.2.4
## [55] withr_3.0.2 scales_1.3.0 UCSC.utils_1.2.0
## [58] rmarkdown_2.28 XVector_0.46.0 httr_1.4.7
## [61] evaluate_1.0.1 rlang_1.1.4 Rcpp_1.0.13
## [64] glue_1.8.0 BiocManager_1.30.25 jsonlite_1.8.9
## [67] R6_2.5.1 zlibbioc_1.52.0
[1] Geistlinger L, Csaba G, Santarelli M, Schiffer L, Ramos M, Zimmer R, Waldron L (2019). GSEABenchmarkeR: Reproducible GSEA Benchmarking. R package version 1.6.0, https://github.com/waldronlab/GSEABenchmarkeR.
[2] Love MI, Huber W, Anders S (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550. doi:10.1186/s13059-014-0550-8.
[3] Robinson MD, McCarthy DJ, Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1), 139-140. doi:10.1093/bioinformatics/btp616.
[4] M. E. Ritchie, B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, and G. K. Smyth. limma powers differential expression analyses for RNAsequencing and microarray studies. Nucleic acids research, 43(7):e47, 2015.
[5] A. Liberzon, C. Birger, H. Thorvaldsdottir, M. Ghandi, J. P. Mesirov, and P. Tamayo. The Molecular Signatures Database Hallmark Gene Set Collection. Cell Systems, 1(6):417-425, 2015.