Error in fitBeta

[cnk113]

Hello,

I've been trying to use glmGamPoi with SCTransform, but I've gotten an error

warning: solve(): system seems singular; attempting approx solution
Error in fitBeta_fisher_scoring(Y, model_matrix, exp(offset_matrix), dispersions,  :
  BLAS/LAPACK routine 'DLASCL' gave error code -4

I should note the matrix I'm using is pretty big, 142k cells by 30k genes. Any ideas why this might be happening?

[const-ae]

Hi Chang,

indeed this seems to be a bug in my software and should not happen. Based on the error message, I guess that the error occurs here.

I should note the matrix I'm using is pretty big, 142k cells by 30k genes.

This is indeed not a small, but actually shouldn't matter. I have tested the package with datasets of similar size, however it is possible that for one of the 30,000 genes there is an unfortunate combination of counts, means, and overdispersion estimates that lead to the failure.

To properly debug the issue, I would need to reproduce the error. Is the your dataset publicly available? If not, could you maybe, do a bisection on the genes: ie. check if the error occurs for the first 15,000 or second 15,000 genes. Then the first quarter or second quarter and so on? Usually this is an easy and fairly quick way to narrow the error down to a small subset of genes that you could share, so that I can debug the problem.

Best, Constantin

[cnk113]

example.zip I've attached an example where it fails. I think it maybe that I'm fitting on genes that occur in only one cell?

[const-ae]

Thanks this is great.

I ran glmGamPoi on the data that you uploaded and did get an error (though not the one that you reported):

library(glmGamPoi)
sessionInfo()
#> R version 4.0.0 Patched (2020-05-04 r78358)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Mojave 10.14.6
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] glmGamPoi_1.0.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] compiler_4.0.0  magrittr_1.5    tools_4.0.0     htmltools_0.4.0
#>  [5] yaml_2.2.1      Rcpp_1.0.4.6    stringi_1.4.6   rmarkdown_2.1  
#>  [9] highr_0.8       knitr_1.28      stringr_1.4.0   xfun_0.13      
#> [13] digest_0.6.25   rlang_0.4.6     evaluate_0.14

test_dat <- readRDS("~/Downloads/example")
set.seed(1)
mat <- as.matrix(test_dat)
#> Loading required package: Matrix
design <- matrix(rnorm(ncol(mat)), ncol = 1)
fit <- glm_gp(mat, design = design, verbose = TRUE)
#> Calculate Size Factors
#> Make initial dispersion estimate
#> Make initial beta estimate
#> Estimate beta
#> Estimate dispersion
#> Estimate beta again
#> Error in fitBeta_fisher_scoring(Y, model_matrix, exp(offset_matrix), dispersions, : BLAS/LAPACK routine 'DLASCLMOIEDDLASRTDSP' gave error code -4

^{Created on 2020-09-22 by the reprex package (v0.3.0)}

However, I actually made the algorithm more robust since the release in March (I just didn't realize this might be the solution earlier) and with the latest version (1.1.13) it runs fine (at least for me):

library(glmGamPoi)
sessionInfo()
#> R version 4.0.0 Patched (2020-05-04 r78358)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Mojave 10.14.6
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] glmGamPoi_1.1.13
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.4.6                knitr_1.28                 
#>  [3] XVector_0.28.0              magrittr_1.5               
#>  [5] GenomicRanges_1.40.0        BiocGenerics_0.34.0        
#>  [7] zlibbioc_1.34.0             IRanges_2.22.1             
#>  [9] lattice_0.20-41             rlang_0.4.6                
#> [11] stringr_1.4.0               highr_0.8                  
#> [13] GenomeInfoDb_1.24.0         tools_4.0.0                
#> [15] grid_4.0.0                  SummarizedExperiment_1.18.1
#> [17] parallel_4.0.0              Biobase_2.48.0             
#> [19] xfun_0.13                   matrixStats_0.56.0         
#> [21] htmltools_0.4.0             yaml_2.2.1                 
#> [23] digest_0.6.25               Matrix_1.2-18              
#> [25] GenomeInfoDbData_1.2.3      S4Vectors_0.26.0           
#> [27] bitops_1.0-6                RCurl_1.98-1.2             
#> [29] evaluate_0.14               rmarkdown_2.1              
#> [31] DelayedArray_0.14.0         stringi_1.4.6              
#> [33] compiler_4.0.0              stats4_4.0.0

test_dat <- readRDS("~/Downloads/example")
set.seed(1)
mat <- as.matrix(test_dat)
design <- matrix(rnorm(ncol(mat)), ncol = 1)
fit <- glm_gp(mat, design = design, verbose = TRUE)
#> Calculate Size Factors (normed_sum)
#> Make initial dispersion estimate
#> Make initial beta estimate
#> Estimate beta
#> Estimate dispersion
#> Fit dispersion trend
#> Shrink dispersion estimates
#> Estimate beta again

^{Created on 2020-09-22 by the reprex package (v0.3.0)}

Could you maybe try to install the latest version from GitHub (devtools::install_github("const-ae/glmGamPoi")) and see if the error disappears?

If not, could you please provide a short reprex that causes the error for you?

Best, Constantin

[cnk113]

It works with the newest version!

[const-ae]

Oh, great. I am relieved that this solved your problem :)