Last updated: 2023-10-23

Checks: 7 0

Knit directory: muse/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20200712) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 9b9327a. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    r_packages_4.3.0/
    Ignored:    r_packages_4.3.1/

Untracked files:
    Untracked:  analysis/cell_ranger.Rmd
    Untracked:  analysis/complex_heatmap.Rmd
    Untracked:  analysis/sleuth.Rmd
    Untracked:  analysis/tss_xgboost.Rmd
    Untracked:  code/multiz100way/
    Untracked:  data/HG00702_SH089_CHSTrio.chr1.vcf.gz
    Untracked:  data/HG00702_SH089_CHSTrio.chr1.vcf.gz.tbi
    Untracked:  data/ncrna_NONCODE[v3.0].fasta.tar.gz
    Untracked:  data/ncrna_noncode_v3.fa
    Untracked:  data/netmhciipan.out.gz
    Untracked:  data/test
    Untracked:  export/davetang039sblog.WordPress.2023-06-30.xml
    Untracked:  export/output/
    Untracked:  women.json

Unstaged changes:
    Modified:   analysis/graph.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/gsva.Rmd) and HTML (docs/gsva.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 9b9327a Dave Tang 2023-10-23 No progress bar
html b0af093 Dave Tang 2023-10-23 Build site.
Rmd aeb90e1 Dave Tang 2023-10-23 GSVA

Following the vignette.

Gene set variation analysis (GSVA) is a particular type of gene set enrichment method that works on single samples and enables pathway-centric analyses of molecular data by performing a conceptually simple but powerful change in the functional unit of analysis, from genes to gene sets. The GSVA package provides the implementation of four single-sample gene set enrichment methods, concretely zscore, plage, ssGSEA and its own called GSVA. While this methodology was initially developed for gene expression data, it can be applied to other types of molecular profiling data. In this vignette we illustrate how to use the GSVA package with bulk microarray and RNA-seq expression data.

Package

Install GSVA. (Dependencies are listed in the Imports section in the DESCRIPTION file.)

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

if (!require("GSVA", quietly = TRUE))
  BiocManager::install("GSVA")

Load package.

library(GSVA)

Quick start

Generate example expression matrix.

p <- 10000
n <- 30
set.seed(1984)
X <- matrix(
  rnorm(p*n),
  nrow=p,
  dimnames=list(paste0("g", 1:p), paste0("s", 1:n))
)
X[1:5, 1:5]
           s1         s2         s3          s4           s5
g1  0.4092032  1.4676435  0.3515056  1.53512312 -1.279009469
g2 -0.3230250 -1.8501416 -0.9198650  1.40036448  0.086613315
g3  0.6358523  1.6084120  1.6380322  0.23799146  0.216628121
g4 -1.8461288 -0.2928844  0.4651573 -0.09766558 -0.009887299
g5  0.9536474 -0.4816006  0.1807824  1.03141311  0.206414282

Generate 100 gene sets that are contain from 10 to up to 100 genes sampled from 1:p.

set.seed(1984)
gs <- as.list(sample(10:100, size=100, replace=TRUE))

gs <- lapply(gs, function(n, p){
  paste0("g", sample(1:p, size=n, replace=FALSE))
}, p)
names(gs) <- paste0("gs", 1:length(gs))

sapply(gs, length)
  gs1   gs2   gs3   gs4   gs5   gs6   gs7   gs8   gs9  gs10  gs11  gs12  gs13 
   49    29    67    90    94    87    41    26    86    77    97    90    45 
 gs14  gs15  gs16  gs17  gs18  gs19  gs20  gs21  gs22  gs23  gs24  gs25  gs26 
   47    54    83    11    75    95    99    94    89    93    50    49    87 
 gs27  gs28  gs29  gs30  gs31  gs32  gs33  gs34  gs35  gs36  gs37  gs38  gs39 
   36    61    84    99    58    30    63    29    35    29    69    41    46 
 gs40  gs41  gs42  gs43  gs44  gs45  gs46  gs47  gs48  gs49  gs50  gs51  gs52 
   38    17    48    72    15    81   100    93    37    99    89    43    36 
 gs53  gs54  gs55  gs56  gs57  gs58  gs59  gs60  gs61  gs62  gs63  gs64  gs65 
   84    83    40    72    90    86    37    23    69    96    20    93    36 
 gs66  gs67  gs68  gs69  gs70  gs71  gs72  gs73  gs74  gs75  gs76  gs77  gs78 
   21    46    76    71    57    48    25    73    26    46    29    53    69 
 gs79  gs80  gs81  gs82  gs83  gs84  gs85  gs86  gs87  gs88  gs89  gs90  gs91 
   69    42    76    30    16    49    35    12    83    99    88    66    10 
 gs92  gs93  gs94  gs95  gs96  gs97  gs98  gs99 gs100 
   51    82    73    97    59    59    42    10    64 

Calculate GSVA enrichment scores.

The first argument to the gsva() function is the gene expression data matrix and the second the collection of gene sets. The gsva() function can take the input expression data and gene sets using different specialized containers that facilitate the access and manipulation of molecular and phenotype data, as well as their associated metadata. Another advanced features include the use of on-disk and parallel backends to enable, respectively, using GSVA on large molecular data sets and speed up computing time.

es_gsva <- gsva(X, gs, verbose=FALSE)
dim(es_gsva)
[1] 100  30

Median enrichment scores.

apply(es_gsva, 2, median)
          s1           s2           s3           s4           s5           s6 
 0.009061514 -0.008458424 -0.005713628 -0.021937531  0.006417182  0.019422486 
          s7           s8           s9          s10          s11          s12 
 0.010140618  0.005812097  0.006495584  0.008887644  0.024577619 -0.031697634 
         s13          s14          s15          s16          s17          s18 
 0.001642502  0.008919786 -0.022622470  0.027695420 -0.015799537 -0.011108686 
         s19          s20          s21          s22          s23          s24 
-0.013956442 -0.015493300 -0.004809844 -0.014081494  0.026845336  0.023895676 
         s25          s26          s27          s28          s29          s30 
 0.006358240  0.010642450 -0.012690144 -0.005999451 -0.005058572 -0.012403422 

ssgsea (Barbie et al. 2009). Single sample GSEA (ssGSEA) is a non-parametric method that calculates a gene set enrichment score per sample as the normalized difference in empirical cumulative distribution functions (CDFs) of gene expression ranks inside and outside the gene set. By default, the implementation in the GSVA package follows the last step described in (Barbie et al. 2009, online methods, pg. 2) by which pathway scores are normalized, dividing them by the range of calculated values. This normalization step may be switched off using the argument ssgsea.norm in the call to the gsva() function; see below.

es_ssgsea <- gsva(X, gs, method = "ssgsea", verbose=FALSE)
[1] "Calculating ranks..."
[1] "Calculating absolute values from ranks..."
[1] "Normalizing..."
apply(es_ssgsea, 2, median)
       s1        s2        s3        s4        s5        s6        s7        s8 
0.1056332 0.1105148 0.1149790 0.1045303 0.1192256 0.1274702 0.1192217 0.1287768 
       s9       s10       s11       s12       s13       s14       s15       s16 
0.1316429 0.1257247 0.1293241 0.1083643 0.1222038 0.1012195 0.1018693 0.1174057 
      s17       s18       s19       s20       s21       s22       s23       s24 
0.1103184 0.1112410 0.1164373 0.1126160 0.1212393 0.1067666 0.1167117 0.1418775 
      s25       s26       s27       s28       s29       s30 
0.1221909 0.1235449 0.1112199 0.1012381 0.1254514 0.1210392 

sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [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: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] GSVA_1.48.3         BiocManager_1.30.22 workflowr_1.7.1    

loaded via a namespace (and not attached):
 [1] blob_1.2.4                  Biostrings_2.68.1          
 [3] bitops_1.0-7                fastmap_1.1.1              
 [5] SingleCellExperiment_1.22.0 RCurl_1.98-1.12            
 [7] promises_1.2.1              rsvd_1.0.5                 
 [9] XML_3.99-0.14               digest_0.6.33              
[11] lifecycle_1.0.3             processx_3.8.2             
[13] KEGGREST_1.40.1             RSQLite_2.3.1              
[15] magrittr_2.0.3              compiler_4.3.1             
[17] rlang_1.1.1                 sass_0.4.7                 
[19] tools_4.3.1                 utf8_1.2.3                 
[21] yaml_2.3.7                  knitr_1.44                 
[23] S4Arrays_1.0.6              bit_4.0.5                  
[25] DelayedArray_0.26.7         abind_1.4-5                
[27] BiocParallel_1.34.2         HDF5Array_1.28.1           
[29] BiocGenerics_0.46.0         grid_4.3.1                 
[31] stats4_4.3.1                fansi_1.0.5                
[33] git2r_0.32.0                beachmat_2.16.0            
[35] xtable_1.8-4                Rhdf5lib_1.22.1            
[37] SummarizedExperiment_1.30.2 cli_3.6.1                  
[39] rmarkdown_2.25              crayon_1.5.2               
[41] rstudioapi_0.15.0           httr_1.4.7                 
[43] DelayedMatrixStats_1.22.6   DBI_1.1.3                  
[45] cachem_1.0.8                rhdf5_2.44.0               
[47] stringr_1.5.0               zlibbioc_1.46.0            
[49] parallel_4.3.1              AnnotationDbi_1.62.2       
[51] XVector_0.40.0              matrixStats_1.0.0          
[53] vctrs_0.6.4                 Matrix_1.5-4.1             
[55] jsonlite_1.8.7              BiocSingular_1.16.0        
[57] callr_3.7.3                 IRanges_2.34.1             
[59] S4Vectors_0.38.2            bit64_4.0.5                
[61] irlba_2.3.5.1               GSEABase_1.62.0            
[63] jquerylib_0.1.4             annotate_1.78.0            
[65] glue_1.6.2                  codetools_0.2-19           
[67] ps_1.7.5                    stringi_1.7.12             
[69] later_1.3.1                 GenomeInfoDb_1.36.4        
[71] GenomicRanges_1.52.1        ScaledMatrix_1.8.1         
[73] tibble_3.2.1                pillar_1.9.0               
[75] htmltools_0.5.6.1           rhdf5filters_1.12.1        
[77] graph_1.78.0                GenomeInfoDbData_1.2.10    
[79] R6_2.5.1                    sparseMatrixStats_1.12.2   
[81] rprojroot_2.0.3             evaluate_0.22              
[83] Biobase_2.60.0              lattice_0.21-8             
[85] png_0.1-8                   memoise_2.0.1              
[87] httpuv_1.6.11               bslib_0.5.1                
[89] Rcpp_1.0.11                 whisker_0.4.1              
[91] xfun_0.40                   fs_1.6.3                   
[93] MatrixGenerics_1.12.3       getPass_0.2-2              
[95] pkgconfig_2.0.3