`gsva` with `method='ssgsea'` errors when given a literal single sample

[ChristopherEeles]

Hi @rcastelo,

Trying to use the GSVA package to put univariate drug response expression biomarkers into the pathway context for use in molecular tumour boards. Our use case allows a researcher or clinician to upload the expression profile (logTPM or RMA normalized microarray) of a patient and use the data to query PharmacoDB.ca for a list of biomarkers associated with response to a given drug therapy from in vitro high throughput drug screening studies in cancer cell lines.

We want to allow clinicians to take these univariate biomarkers and put them into a pathway context to see if the gene associated with drug response is in a pathway known to be involved in their patient's cancer or cancer more generally. To this end, we are trying to use the ssGSEA method via GSVA to identify relevant gene sets from a list we compiled via msigdbr.

I have found a bug, in the sense that I can reproducibly cause an error in the package by passing a single column matrix into the GSVA method. However, it is not clear to me from the package vignette or a brief review of Barbie _et _al.__ (2009) if the method is actually intended to work with a literal single patient's gene expression profile.

Minimal Reprex

expr <- matrix(rnorm(26), ncol=1, dimnames=list(LETTERS, "sample1"))
gs <- list(set1=c("A", "D", "G"), set2=c("Z", "Y", "X"))
gsea <- gsva(expr, gs, method="ssgsea", min.sz=1, verbose=TRUE)
## Estimating ssGSEA scores for 2 gene sets.
## [1] "Calculating ranks..."
## [1] "Calculating absolute values from ranks..."
## |======================================================================| 100%
## [1] "Normalizing..."
## Error in `rownames<-`(`*tmp*`, value = names(geneSets)) : 
##  attempt to set 'rownames' on an object with no dimensions
## In addition: Warning message:
## In .filterFeatures(expr, method) :
##  26 genes with constant expression values throuhgout the samples.

The actual bug is caused because the matrix gets simplified to a vector before trying to assign rownames on line 885 of gsva.R, inside the ssgsea function, and could be fixed by converting es, which is a vector when the bug occurs, into a single column matrix with something like if (n == 1 && !is.matrix(es)) es <- matrix(es, ncol=1, dimnames=list(names(es), colnames(X)).

It is worth noting that the bug doesn't occur when length(geneSets) == 1, since the condition is already converting es back to being a matrix.

I wanted to confirm that the method is able to handle enrichment analysis using only one patient profile before I open a PR with a fix.

If the method does support the single patient case and you would like a PR please let me know. If it does not and/or you have suggestions for an alternative approach it would be appreciated.

We actually already have a "reference" gene expression population that the user uploads to identify genes of interest, so maybe we could add the mean of that population as a reference population.

R Environment

> sessionInfo()
R version 4.2.0 Patched (2022-04-26 r82266)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04 LTS

Matrix products: default
BLAS:   /opt/R/R-latest/lib/R/lib/libRblas.so
LAPACK: /opt/R/R-latest/lib/R/lib/libRlapack.so

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       

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

other attached packages:
[1] GSVA_1.45.0

loaded via a namespace (and not attached):
 [1] KEGGREST_1.37.0             SummarizedExperiment_1.27.1
 [3] beachmat_2.13.1             BiocSingular_1.13.0        
 [5] HDF5Array_1.25.0            rhdf5_2.41.0               
 [7] lattice_0.20-45             pak_0.3.0.9000             
 [9] vctrs_0.4.1                 stats4_4.2.0               
[11] blob_1.2.3                  XML_3.99-0.9               
[13] rlang_1.0.2                 glue_1.6.2                 
[15] DBI_1.1.2                   BiocParallel_1.31.6        
[17] SingleCellExperiment_1.19.0 BiocGenerics_0.43.0        
[19] bit64_4.0.5                 matrixStats_0.62.0         
[21] GenomeInfoDbData_1.2.8      zlibbioc_1.43.0            
[23] MatrixGenerics_1.9.0        Biostrings_2.65.0          
[25] rsvd_1.0.5                  ScaledMatrix_1.5.0         
[27] memoise_2.0.1               Biobase_2.57.1             
[29] callr_3.7.0                 IRanges_2.31.0             
[31] fastmap_1.1.0               ps_1.7.0                   
[33] GenomeInfoDb_1.33.3         irlba_2.3.5                
[35] parallel_4.2.0              AnnotationDbi_1.59.0       
[37] GSEABase_1.59.0             Rcpp_1.0.8.3               
[39] xtable_1.8-4                cachem_1.0.6               
[41] DelayedArray_0.23.0         S4Vectors_0.35.0           
[43] graph_1.75.0                annotate_1.75.0            
[45] XVector_0.37.0              bit_4.0.4                  
[47] png_0.1-7                   processx_3.5.3             
[49] GenomicRanges_1.49.0        grid_4.2.0                 
[51] cli_3.3.0                   tools_4.2.0                
[53] bitops_1.0-7                rhdf5filters_1.9.0         
[55] RCurl_1.98-1.6              RSQLite_2.2.14             
[57] crayon_1.5.1                Matrix_1.4-1               
[59] DelayedMatrixStats_1.19.0   sparseMatrixStats_1.9.0    
[61] httr_1.4.3                  Rhdf5lib_1.19.2            
[63] R6_2.5.1                    compiler_4.2.0 

Best, Christopher Eeles Software Developer Haibe-Kains Lab | PM-Research University Health Network

[rcastelo]

Dear Christopher,

Thank you very much for your very detailed bug report. The problem was in the normalization step, since right before the cbind call always returns a matrix, and that normalization step was dropping the dimensions when there was only a single column (sample) as in your example. So, just adding drop=FALSE in the associated subsetting operation fixed the problem, which has been pushed to the 1.45.1 devel version and 1.44.1 release version and should become available via BiocManager::install() in the next 24/48 hours. Please let me know if this fixes completely your use case or additional fixes are needed.

[ChristopherEeles]

Hi @rcastelo,

Thanks for the prompt fix. Everything is working as expected now with 1.45.1.

Best, Chris