search

satijalab / seurat

R toolkit for single cell genomics
http://www.satijalab.org/seurat
Other
2.27k stars 910 forks source link

Problem with Leverage Score Calculation for Integration of Single Cell Seurat Objects #7690

Open SSridhar234 opened 1 year ago

SSridhar234 commented 1 year ago

I am trying to use the SketchData function with the Leverage Score method, and the computation will not go through. All of my Seurat objects seem consistent with each other, and the matrices don't seem to have any major issues. Since everything seems okay, it perplexes me why I am having this issue.

This is the code that I am currently running (the error occurs on the last line):

### This script takes as input the Seurat Objects from each CODEX sample in the atlas and creates a combined CODEX Atlas Object
library(Seurat)
options(Seurat.object.assay.version = 'v5')
library(tidyverse)
library(patchwork)
library(readr)
library(dplyr) 
library(ggrepel)
library(ggplot2)
library(RColorBrewer)
options(future.globals.maxSize = 1e9)
options(Seurat.object.assay.version = "v5")
source("/mnt/isilon/some_lab/myname/CODEX_Functions.R")
#Import RDS Files
#load in every individual R file for each sample 
Sample1 <- readRDS("/mnt/isilon/some_lab/myname/Sample1.RDS")
Sample2 <- readRDS("/mnt/isilon/some_lab/myname/Sample2.RDS")       
#Work with Combined Data
ob.list <- c()
ob.list[[1]] <- Sample1
ob.list[[2]] <- Sample2
combined <- merge(x = ob.list[[1]], y = c(ob.list[[2]]), add.cell.ids = c("Sample1", "Sample2"))
combined <- JoinLayers(combined)
#important_features  <- c("AB", "CD", "EF", "GH", "IJ", "KL", "MN", "OP", "QR", "ST", "UV",  "WX")

important_features = rownames(combined)
VariableFeatures(combined) = important_features
combined <- NormalizeData(object = combined, normalization.method = "CLR", margin = 1)
combined <- split(combined, f = combined$orig.ident)
combined <- SketchData(object = combined, ncells = 50000, method = "LeverageScore", sketched.assay = "sketch")

The error that I am getting is the following:

Calcuating Leverage Score
Running LeverageScore for layer data.pHGG_Autopsy
sampling 5000 cells
Performing QR decomposition
'as(<dgeMatrix>, "dgCMatrix")' is deprecated.
Use 'as(., "CsparseMatrix")' instead.
See help("Deprecated") and help("Matrix-deprecated").
Error in backsolve(r = R, x = diag(x = ncol(x = R))) : 
  singular matrix in 'backsolve'. First zero in diagonal [52]

This is the SessionInfo():

R version 4.2.3 (2023-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 9.1 (Plow)

Matrix products: default
BLAS:   /cm/shared/apps_chop/R/4.2.3/lib64/R/lib/libRblas.so
LAPACK: /cm/shared/apps_chop/R/4.2.3/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] RColorBrewer_1.1-3      ggrepel_0.9.3           patchwork_1.1.2         lubridate_1.9.2        
 [5] forcats_1.0.0           stringr_1.5.0           dplyr_1.1.2             purrr_1.0.1            
 [9] readr_2.1.4             tidyr_1.3.0             tibble_3.2.1            ggplot2_3.4.2          
[13] tidyverse_2.0.0         Seurat_4.9.9.9049       SeuratObject_4.9.9.9086 sp_2.0-0               

loaded via a namespace (and not attached):
  [1] Rtsne_0.16             colorspace_2.1-0       deldir_1.0-9           ellipsis_0.3.2        
  [5] ggridges_0.5.4         RcppHNSW_0.4.1         rstudioapi_0.14        spatstat.data_3.0-1   
  [9] leiden_0.4.3           listenv_0.9.0          RSpectra_0.16-1        fansi_1.0.4           
 [13] codetools_0.2-19       splines_4.2.3          polyclip_1.10-4        spam_2.9-1            
 [17] jsonlite_1.8.5         ica_1.0-3              cluster_2.1.4          png_0.1-8             
 [21] uwot_0.1.16            shiny_1.7.4            sctransform_0.3.5      spatstat.sparse_3.0-2 
 [25] BiocManager_1.30.21    compiler_4.2.3         httr_1.4.6             Matrix_1.5-4          
 [29] fastmap_1.1.1          lazyeval_0.2.2         cli_3.6.1              later_1.3.1           
 [33] htmltools_0.5.5        tools_4.2.3            igraph_1.5.0           dotCall64_1.0-2       
 [37] gtable_0.3.3           glue_1.6.2             RANN_2.6.1             reshape2_1.4.4        
 [41] Rcpp_1.0.10            scattermore_1.2        vctrs_0.6.3            spatstat.explore_3.2-1
 [45] nlme_3.1-162           progressr_0.13.0       lmtest_0.9-40          spatstat.random_3.1-5 
 [49] globals_0.16.2         timechange_0.2.0       mime_0.12              miniUI_0.1.1.1        
 [53] lifecycle_1.0.3        irlba_2.3.5.1          renv_0.17.3            goftest_1.2-3         
 [57] future_1.32.0          MASS_7.3-60            zoo_1.8-12             scales_1.2.1          
 [61] hms_1.1.3              promises_1.2.0.1       spatstat.utils_3.0-3   parallel_4.2.3        
 [65] reticulate_1.30        pbapply_1.7-2          gridExtra_2.3          stringi_1.7.12        
 [69] fastDummies_1.6.3      rlang_1.1.1            pkgconfig_2.0.3        matrixStats_1.0.0     
 [73] lattice_0.21-8         ROCR_1.0-11            tensor_1.5             htmlwidgets_1.6.2     
 [77] cowplot_1.1.1          tidyselect_1.2.0       parallelly_1.36.0      RcppAnnoy_0.0.20      
 [81] plyr_1.8.8             magrittr_2.0.3         R6_2.5.1               generics_0.1.3        
 [85] pillar_1.9.0           withr_2.5.0            fitdistrplus_1.1-11    survival_3.5-5        
 [89] abind_1.4-5            future.apply_1.11.0    KernSmooth_2.23-21     utf8_1.2.3            
 [93] spatstat.geom_3.2-1    plotly_4.10.2          tzdb_0.4.0             grid_4.2.3            
 [97] data.table_1.14.8      digest_0.6.32          xtable_1.8-4           httpuv_1.6.11         
[101] munsell_0.5.0          viridisLite_0.4.2
Sandman-1 commented 8 months ago

Hello. I am also experiencing this type of issue. I have a dataset with 30 layers, 3 of which have fewer than 10 K cells and several others that have 10-12 K cells. One of the layers with fewer than 10 K cells returns this same error message. I have resorted to editing the SketchData function to not calculate leverage scores on any layers with fewer than 10 K cells by subsetting the object, as I do not know how to edit the LeverageScore function to return values of 1 if the matrix following QR decomposition is singular. It would be awesome if you guys could edit the function to do so.

levinhein commented 6 months ago

Hello. I'm also experiencing this error. My data has 174 layers.

Any suggestions?

haolangchen commented 5 months ago

I am trying to use the SketchData function with the Leverage Score method, and the computation will not go through. All of my Seurat objects seem consistent with each other, and the matrices don't seem to have any major issues. Since everything seems okay, it perplexes me why I am having this issue.

This is the code that I am currently running (the error occurs on the last line):

### This script takes as input the Seurat Objects from each CODEX sample in the atlas and creates a combined CODEX Atlas Object
library(Seurat)
options(Seurat.object.assay.version = 'v5')
library(tidyverse)
library(patchwork)
library(readr)
library(dplyr) 
library(ggrepel)
library(ggplot2)
library(RColorBrewer)
options(future.globals.maxSize = 1e9)
options(Seurat.object.assay.version = "v5")
source("/mnt/isilon/some_lab/myname/CODEX_Functions.R")
#Import RDS Files
#load in every individual R file for each sample 
Sample1 <- readRDS("/mnt/isilon/some_lab/myname/Sample1.RDS")
Sample2 <- readRDS("/mnt/isilon/some_lab/myname/Sample2.RDS")       
#Work with Combined Data
ob.list <- c()
ob.list[[1]] <- Sample1
ob.list[[2]] <- Sample2
combined <- merge(x = ob.list[[1]], y = c(ob.list[[2]]), add.cell.ids = c("Sample1", "Sample2"))
combined <- JoinLayers(combined)
#important_features  <- c("AB", "CD", "EF", "GH", "IJ", "KL", "MN", "OP", "QR", "ST", "UV",  "WX")

important_features = rownames(combined)
VariableFeatures(combined) = important_features
combined <- NormalizeData(object = combined, normalization.method = "CLR", margin = 1)
combined <- split(combined, f = combined$orig.ident)
combined <- SketchData(object = combined, ncells = 50000, method = "LeverageScore", sketched.assay = "sketch")

The error that I am getting is the following:

Calcuating Leverage Score
Running LeverageScore for layer data.pHGG_Autopsy
sampling 5000 cells
Performing QR decomposition
'as(<dgeMatrix>, "dgCMatrix")' is deprecated.
Use 'as(., "CsparseMatrix")' instead.
See help("Deprecated") and help("Matrix-deprecated").
Error in backsolve(r = R, x = diag(x = ncol(x = R))) : 
  singular matrix in 'backsolve'. First zero in diagonal [52]

This is the SessionInfo():

R version 4.2.3 (2023-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 9.1 (Plow)

Matrix products: default
BLAS:   /cm/shared/apps_chop/R/4.2.3/lib64/R/lib/libRblas.so
LAPACK: /cm/shared/apps_chop/R/4.2.3/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] RColorBrewer_1.1-3      ggrepel_0.9.3           patchwork_1.1.2         lubridate_1.9.2        
 [5] forcats_1.0.0           stringr_1.5.0           dplyr_1.1.2             purrr_1.0.1            
 [9] readr_2.1.4             tidyr_1.3.0             tibble_3.2.1            ggplot2_3.4.2          
[13] tidyverse_2.0.0         Seurat_4.9.9.9049       SeuratObject_4.9.9.9086 sp_2.0-0               

loaded via a namespace (and not attached):
  [1] Rtsne_0.16             colorspace_2.1-0       deldir_1.0-9           ellipsis_0.3.2        
  [5] ggridges_0.5.4         RcppHNSW_0.4.1         rstudioapi_0.14        spatstat.data_3.0-1   
  [9] leiden_0.4.3           listenv_0.9.0          RSpectra_0.16-1        fansi_1.0.4           
 [13] codetools_0.2-19       splines_4.2.3          polyclip_1.10-4        spam_2.9-1            
 [17] jsonlite_1.8.5         ica_1.0-3              cluster_2.1.4          png_0.1-8             
 [21] uwot_0.1.16            shiny_1.7.4            sctransform_0.3.5      spatstat.sparse_3.0-2 
 [25] BiocManager_1.30.21    compiler_4.2.3         httr_1.4.6             Matrix_1.5-4          
 [29] fastmap_1.1.1          lazyeval_0.2.2         cli_3.6.1              later_1.3.1           
 [33] htmltools_0.5.5        tools_4.2.3            igraph_1.5.0           dotCall64_1.0-2       
 [37] gtable_0.3.3           glue_1.6.2             RANN_2.6.1             reshape2_1.4.4        
 [41] Rcpp_1.0.10            scattermore_1.2        vctrs_0.6.3            spatstat.explore_3.2-1
 [45] nlme_3.1-162           progressr_0.13.0       lmtest_0.9-40          spatstat.random_3.1-5 
 [49] globals_0.16.2         timechange_0.2.0       mime_0.12              miniUI_0.1.1.1        
 [53] lifecycle_1.0.3        irlba_2.3.5.1          renv_0.17.3            goftest_1.2-3         
 [57] future_1.32.0          MASS_7.3-60            zoo_1.8-12             scales_1.2.1          
 [61] hms_1.1.3              promises_1.2.0.1       spatstat.utils_3.0-3   parallel_4.2.3        
 [65] reticulate_1.30        pbapply_1.7-2          gridExtra_2.3          stringi_1.7.12        
 [69] fastDummies_1.6.3      rlang_1.1.1            pkgconfig_2.0.3        matrixStats_1.0.0     
 [73] lattice_0.21-8         ROCR_1.0-11            tensor_1.5             htmlwidgets_1.6.2     
 [77] cowplot_1.1.1          tidyselect_1.2.0       parallelly_1.36.0      RcppAnnoy_0.0.20      
 [81] plyr_1.8.8             magrittr_2.0.3         R6_2.5.1               generics_0.1.3        
 [85] pillar_1.9.0           withr_2.5.0            fitdistrplus_1.1-11    survival_3.5-5        
 [89] abind_1.4-5            future.apply_1.11.0    KernSmooth_2.23-21     utf8_1.2.3            
 [93] spatstat.geom_3.2-1    plotly_4.10.2          tzdb_0.4.0             grid_4.2.3            
 [97] data.table_1.14.8      digest_0.6.32          xtable_1.8-4           httpuv_1.6.11         
[101] munsell_0.5.0          viridisLite_0.4.2

hi,you need run FindVariableFeatures(),after finsh, you can run sktchdata