satijalab / seurat

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

PrepSCTIntegration Error in scale.data[anchor.features, ] : subscript out of bounds #8216

Closed LukeJ89 closed 10 months ago

LukeJ89 commented 10 months ago

Hi,

I'm fairly new to R and Seurat, so I hope this will be an easy fix. I've a list of 7 Seurat objects from 10X spatial datasets, and ran SCTransform on them. The list is "Prep_Healthy". Running this code:

Prep_Healthy<- SCTransform(Prep_Healthy, vst.flavor = "v2", assay="Spatial") features <- SelectIntegrationFeatures(object.list = Prep_healthy, nfeatures = 3000) Prep_healthy <- PrepSCTIntegration(object.list = Prep_healthy, anchor.features = features)

But running this I get the error code "| 0 % ~calculating Error in scale.data[anchor.features, ] : subscript out of bounds" If I set the nfeatures lower, then it works, the highest is nfeatures = 690 after which the error occurs. So it seems that the list for the anchor.features is too long, not sure if it's a memory thing but I've tried it on 2 computers.

I also tried merging and splitting the list as in #3085 and also tried to use intersect in case some names were missing across the datasets as in #6707. I'm at a bit of a loss, any help would be greatly appreciated.

" Healthy_Samples <- list(HEALTHY.Male.s1, HEALTHY.Male.s2, HEALTHY.Female1.s1, HEALTHY.Female1.s2.clean, HEALTHY.Female1.s3, HEALTHY.Female2.s1, HEALTHY.Female2.s2)

Filter for counts greater than 200

i <- 1 while(i <= length(Healthy_Samples)){ filtered_data <- st_filter_by_genes(st.data = Healthy_Samples[[i]],x = 200) Healthy_Samples[[i]] <- filtered_data i <- i+1 }

st_SCT <- function(data.list){ for (i in 1:length(data.list)) { data.list[[i]] <- SCTransform(data.list[[i]], vst.flavor = "v2", verbose = T,assay="Spatial")} data.list[[i]] <- UpdateSeuratObject(data.list[[i]]) return(data.list) }

Prep_healthy <- st_SCT(Healthy_Samples)

features <- SelectIntegrationFeatures(object.list = Prep_healthy, nfeatures = 3000) Prep_healthy <- PrepSCTIntegration(object.list = Prep_healthy, anchor.features = features)

saketkc commented 10 months ago
features <- SelectIntegrationFeatures(object.list = Prep_healthy, nfeatures = 3000)

expects that object.list to be list. Assuming this was an issue while pasting, can you paste your entire code here along with sessionInfo()

LukeJ89 commented 10 months ago
st_combine <- function(data.list,ndim,res,seed=100){
  set.seed(seed)
  for (i in 1:length(data.list)) {
    data.list[[i]] <- SCTransform(data.list[[i]], vst.flavor = "v2", verbose = T,assay="Spatial")
  }
  features <- SelectIntegrationFeatures(object.list =  data.list, nfeatures = 3000)
  data.list <- PrepSCTIntegration(object.list = data.list, anchor.features = features)
  data.anchors <- FindIntegrationAnchors(object.list = data.list, normalization.method = "SCT", 
                                         anchor.features = features)
  data.combined.sct <- IntegrateData(anchorset = data.anchors, normalization.method = "SCT")
  data.combined.sct <- RunPCA(data.combined.sct, verbose = FALSE)
  data.combined.sct <- FindNeighbors(data.combined.sct, dims = 1:ndim)
  data.combined.sct <- FindClusters(data.combined.sct, verbose = FALSE,resolution = res)
  data.combined.sct <- RunUMAP(data.combined.sct, dims = 1:ndim)
  return(data.combined.sct)
}

Prep_healthy <- c(HEALTHY.Male.s1, HEALTHY.Male.s2,
                        HEALTHY.Female1.s1, HEALTHY.Female1.s2.clean,
                        HEALTHY.Female1.s3, HEALTHY.Female2.s1, HEALTHY.Female2.s2)

#Filter for counts greater than 200
i <- 1
while(i <= length(Prep_healthy)){
  filtered_data <- st_filter_by_genes(st.data = Prep_healthy[[i]],x = 200)
  Prep_healthy[[i]] <- filtered_data
  i <- i+1
}

Prep_healthy <- st_combine(Prep_healthy, ndim = 20, res = 0.6)

Each name in the "Prep_healthy" list is a Seurat object

  |                                                  | 0 % ~calculating  Error in scale.data[anchor.features, ] : subscript out of bounds
In addition: There were 50 or more warnings (use warnings() to see the first 50)

> sessionInfo()
R version 4.2.3 (2023-03-15 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22635)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.utf8  LC_CTYPE=English_United Kingdom.utf8   
[3] LC_MONETARY=English_United Kingdom.utf8 LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.utf8    

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

other attached packages:
 [1] harmony_1.1.0      Rcpp_1.0.11        ggsci_3.0.0        scales_1.3.0       pheatmap_1.0.12   
 [6] RColorBrewer_1.1-3 cowplot_1.1.2      Seurat_5.0.1.9001  SeuratObject_5.0.1 sp_2.1-1          
[11] lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2       
[16] readr_2.1.4        tidyr_1.3.0        tibble_3.2.1       ggplot2_3.4.4      tidyverse_2.0.0   

loaded via a namespace (and not attached):
  [1] spam_2.10-0                 plyr_1.8.9                  igraph_1.5.1               
  [4] lazyeval_0.2.2              splines_4.2.3               RcppHNSW_0.5.0             
  [7] listenv_0.9.0               scattermore_1.2             GenomeInfoDb_1.34.9        
 [10] usethis_2.2.2               digest_0.6.33               htmltools_0.5.7            
 [13] fansi_1.0.5                 magrittr_2.0.3              memoise_2.0.1              
 [16] tensor_1.5                  cluster_2.1.4               ROCR_1.0-11                
 [19] tzdb_0.4.0                  remotes_2.4.2.1             globals_0.16.2             
 [22] matrixStats_1.1.0           timechange_0.2.0            spatstat.sparse_3.0-3      
 [25] prettyunits_1.2.0           colorspace_2.1-0            ggrepel_0.9.4              
 [28] RCurl_1.98-1.13             callr_3.7.3                 crayon_1.5.2               
 [31] jsonlite_1.8.7              progressr_0.14.0            spatstat.data_3.0-3        
 [34] survival_3.5-7              zoo_1.8-12                  glue_1.6.2                 
 [37] polyclip_1.10-6             gtable_0.3.4                zlibbioc_1.44.0            
 [40] XVector_0.38.0              leiden_0.4.3.1              DelayedArray_0.24.0        
 [43] pkgbuild_1.4.2              future.apply_1.11.0         BiocGenerics_0.44.0        
 [46] abind_1.4-5                 spatstat.random_3.2-1       miniUI_0.1.1.1             
 [49] viridisLite_0.4.2           xtable_1.8-4                reticulate_1.34.0          
 [52] dotCall64_1.1-0             stats4_4.2.3                profvis_0.3.8              
 [55] htmlwidgets_1.6.4           httr_1.4.7                  ellipsis_0.3.2             
 [58] ica_1.0-3                   urlchecker_1.0.1            pkgconfig_2.0.3            
 [61] farver_2.1.1                uwot_0.1.16                 deldir_2.0-2               
 [64] utf8_1.2.4                  tidyselect_1.2.0            labeling_0.4.3             
 [67] rlang_1.1.2                 reshape2_1.4.4              later_1.3.1                
 [70] munsell_0.5.0               tools_4.2.3                 cachem_1.0.8               
 [73] cli_3.6.1                   generics_0.1.3              devtools_2.4.5             
 [76] ggridges_0.5.5              fastmap_1.1.1               goftest_1.2-3              
 [79] processx_3.8.2              fs_1.6.3                    fitdistrplus_1.1-11        
 [82] RANN_2.6.1                  sparseMatrixStats_1.10.0    pbapply_1.7-2              
 [85] future_1.33.0               nlme_3.1-163                mime_0.12                  
 [88] compiler_4.2.3              rstudioapi_0.15.0           plotly_4.10.3              
 [91] png_0.1-8                   spatstat.utils_3.0-4        glmGamPoi_1.10.2           
 [94] stringi_1.8.2               ps_1.7.5                    RSpectra_0.16-1            
 [97] lattice_0.22-5              Matrix_1.6-3                vctrs_0.6.4                
[100] pillar_1.9.0                lifecycle_1.0.4             BiocManager_1.30.22        
[103] spatstat.geom_3.2-7         lmtest_0.9-40               RcppAnnoy_0.0.21           
[106] bitops_1.0-7                data.table_1.14.8           irlba_2.3.5.1              
[109] GenomicRanges_1.50.2        httpuv_1.6.12               patchwork_1.1.3            
[112] R6_2.5.1                    promises_1.2.1              KernSmooth_2.23-22         
[115] gridExtra_2.3               IRanges_2.32.0              parallelly_1.36.0          
[118] sessioninfo_1.2.2           codetools_0.2-19            fastDummies_1.7.3          
[121] MASS_7.3-60                 pkgload_1.3.3               SummarizedExperiment_1.28.0
[124] withr_2.5.2                 sctransform_0.4.1           GenomeInfoDbData_1.2.9     
[127] S4Vectors_0.36.2            parallel_4.2.3              hms_1.1.3                  
[130] grid_4.2.3                  DelayedMatrixStats_1.20.0   MatrixGenerics_1.10.0      
[133] Rtsne_0.16                  spatstat.explore_3.2-5      Biobase_2.58.0             
[136] shiny_1.8.0   
saketkc commented 10 months ago

Ahh, I see. This is related to #8203. My suggestion for now would be to create a "RNA" assay, by copying the "Spatial" assay and doing everything on the "RNA" assay.

object[["RNA"]] <- object[["Spatial"]]
LukeJ89 commented 10 months ago

Perfect! That works, thank you so much!

FerrenaAlexander commented 8 months ago

Hi all, thanks for all of your amazing hard work! Just wanted to report something maybe related to this, I was able to get around the following error in the GetResidual function: "Error in .subscript.2ary(x, , j, drop = drop) : subscript out of bounds".

By setting object[["RNA"]] <- object[["Spatial"]] , at two specific times:

  1. before running SCTransform on each indiviudal object
  2. before running PrepSCTIntegration

I was using GetResidual after integration, to plot genes that were not in the scale.data layer after integration. Not sure if the second step was really necessary. Doing the first of these fixed the GetResidual error in integration, but the second fixed only the PrepSCTIntegration error.

Also not sure whether this may be related to #8561.