SCTransform gives 'BLAS/LAPACK routine 'DLASCLS'" error

[spadavec]

I'm trying to load multiple scRNA datasets, merge them, and run SCTransform, but I'm getting the following error:

warning: solve(): system seems singular; attempting approx solution
Error in qpois_reg(regressor_data, y, 1e-09, 100, 1.0001, TRUE) : 
  BLAS/LAPACK routine 'DLASCLS' gave error code -4

If I try to follow up my SCTransform w/ the PrepSCTIntegration step with the following command

obj_l <- PrepSCTIntegration(object.list = obj_l, anchor.features = prostate.features, verbose = FALSE)

I get the following error:

Error: The following assays have not been processed with SCTransform:
 object: 1 - assay: RNA
 object: 2 - assay: RNA
 object: 3 - assay: RNA
 object: 4 - assay: RNA

My full code is :

enza.data <- Read10X(data.dir = '/data/Enza/')
enza_h.data <- Read10X(data.dir = '/data/HypoxEnza/')
hyp.data <- Read10X(data.dir = '/data/LNCAP/Hypox/')
control.data <- Read10X(data.dir = '/data/Control/')

enza <- CreateSeuratObject(counts = enza.data, project = "enza", min.cells = 3)
enzaH <- CreateSeuratObject(counts = enza_h.data, project = "enzaH", min.cells = 3)
Hyp <- CreateSeuratObject(counts = hyp.data, project='Hyp', min.cells= 3)
Control <- CreateSeuratObject(counts = control.data, project='Control', min.cells=3)

obj_l = merge(object1= enza, object2= enzaH, object3=Hyp, objec4=Control, add.cell.id1='enza', add.cell.id2='enzaH', add.cell.id3='Hyp', add.cell.id4='Control')

for (i in 1:length(obj_l)) {
  obj_l[[i]] <- SCTransform(obj_l[[i]], verbose = FALSE)
}

[ChristophH]

I have not seen this error before. Can you share the data and the exact steps that led to this problem?

In the meantime, make sure to install the latest sctransform version (remotes::install_github("ChristophH/sctransform@develop")). You can also try method = 'glmGamPoi' when calling SCTransform

[ahmedmahfouz]

I am getting the same error on a visium dataset. Tried also method = 'glmGamPoi' with no success. This only happens to one out of 4 datasets that I am working with.

Happy to share the data if necessary.

What I ran so far:
library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)

DBA_2J <- Load10X_Spatial(data.dir = paste0(data.dir,samples[3],"/outs"))
DBA_2J <- SCTransform(DBA_2J, assay = "Spatial", verbose = TRUE)

Output:
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 12362 by 1519
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 1519 cells
| | 0%
Error in qpois_reg(regressor_data, y, 1e-09, 100, 1.0001, TRUE): BLAS/LAPACK routine 'DLASCL' gave error code -4

Session info:
R version 3.6.1 (2019-07-05)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 8

Matrix products: default
BLAS/LAPACK: /home/amahfouz/miniconda3/envs/rnotebook/lib/R/lib/libRblas.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] dplyr_0.8.0.1 patchwork_1.1.0 ggplot2_3.3.2 Seurat_3.2.2

loaded via a namespace (and not attached):
[1] nlme_3.1-150 matrixStats_0.57.0 bit64_4.0.5
[4] RcppAnnoy_0.0.17 RColorBrewer_1.1-2 httr_1.4.0
[7] repr_0.19.2 sctransform_0.3.1 tools_3.6.1
Here is my session info:
R version 3.6.1 (2019-07-05)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 8

Matrix products: default
BLAS/LAPACK: /home/amahfouz/miniconda3/envs/rnotebook/lib/R/lib/libRblas.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] dplyr_0.8.0.1 patchwork_1.1.0 ggplot2_3.3.2 Seurat_3.2.2

loaded via a namespace (and not attached):
[1] nlme_3.1-150 matrixStats_0.57.0 bit64_4.0.5
[4] RcppAnnoy_0.0.17 RColorBrewer_1.1-2 httr_1.4.0
[7] repr_0.19.2 sctransform_0.3.1 tools_3.6.1
[10] R6_2.5.0 irlba_2.3.3 rpart_4.1-15
[13] KernSmooth_2.23-15 uwot_0.1.9 mgcv_1.8-33
[16] lazyeval_0.2.2 colorspace_2.0-0 withr_2.3.0
[19] tidyselect_0.2.5 gridExtra_2.3 bit_4.0.4
[22] compiler_3.6.1 hdf5r_1.3.3 plotly_4.9.2.1
[25] scales_1.1.1 lmtest_0.9-38 spatstat.data_1.5-2
[28] ggridges_0.5.2 pbapply_1.4-3 goftest_1.2-2
[31] spatstat_1.64-1 pbdZMQ_0.3-3 stringr_1.4.0
[34] digest_0.6.27 spatstat.utils_1.17-0 base64enc_0.1-3
[37] pkgconfig_2.0.3 htmltools_0.3.6 parallelly_1.21.0
[40] htmlwidgets_1.3 rlang_0.4.9 shiny_1.3.2
[43] zoo_1.8-6 jsonlite_1.7.1 ica_1.0-2
[46] magrittr_2.0.1 Matrix_1.2-18 Rcpp_1.0.5
[49] IRkernel_0.8.15 munsell_0.5.0 abind_1.4-5
[52] reticulate_1.18 lifecycle_0.2.0 stringi_1.5.3
[55] MASS_7.3-53 Rtsne_0.15 plyr_1.8.6
[58] grid_3.6.1 parallel_3.6.1 listenv_0.8.0
[61] promises_1.0.1 ggrepel_0.8.2 crayon_1.3.4
[64] deldir_0.2-3 miniUI_0.1.1.1 lattice_0.20-41
[67] IRdisplay_0.7.0 cowplot_1.1.0 splines_3.6.1
[70] tensor_1.5 pillar_1.4.7 igraph_1.2.6
[73] uuid_0.1-2 future.apply_1.6.0 reshape2_1.4.4
[76] codetools_0.2-16 leiden_0.3.6 glue_1.4.2
[79] evaluate_0.14 data.table_1.12.2 vctrs_0.3.5
[82] png_0.1-7 httpuv_1.5.1 polyclip_1.10-0
[85] gtable_0.3.0 RANN_2.6.1 purrr_0.3.2
[88] tidyr_0.8.3 future_1.20.1 assertthat_0.2.1
[91] rsvd_1.0.3 mime_0.6 xtable_1.8-4
[94] later_0.8.0 survival_2.44-1.1 viridisLite_0.3.0
[97] tibble_3.0.4 cluster_2.0.8 globals_0.14.0
[100] fitdistrplus_1.1-3 ellipsis_0.3.1 ROCR_1.0-11

[ChristophH]

I'd be more than happy to have a look. If you want to share the data privately, use christoph.hafemeister@nyu.edu

[ChristophH]

@ahmedmahfouz Thanks for sharing the data. It turns out that there is a 'cell' with zero counts that is causing the problem. After loading the raw data, I recommend plotting some QC metrics and basic filtering. For example:

VlnPlot(DBA_2J, features = c('nCount_Spatial', 'nFeature_Spatial'))
# get rid of cells with very few counts or extremely high
DBA_2J <- DBA_2J[, DBA_2J$nCount_Spatial >= 4000 & DBA_2J$nCount_Spatial <= 15000 & DBA_2J$nFeature_Spatial >= 700 & DBA_2J$nFeature_Spatial <= 3000]
# proceed with analysis

Of course, adjust these filter based on your data, and don't use strict filters if a signal of interest could be correlated with the QC metrics.

[ahmedmahfouz]

Thanks Christoph, this worked.