smorabit / hdWGCNA

High dimensional weighted gene co-expression network analysis
https://smorabit.github.io/hdWGCNA/
Other
316 stars 31 forks source link

Error in Matrix(from, sparse = TRUE, doDiag = FALSE) : invalid 'data' from MetaspotsByGroups() #173

Closed weiyouzi321 closed 5 months ago

weiyouzi321 commented 6 months ago

Describe the bug Error in Matrix(from, sparse = TRUE, doDiag = FALSE) : invalid 'data'

exciting work!i have tried the example 10X data (sagittal brain section) it worked well and showed well reproducibility,so package installation and R enviroment may be approriate;but when trying hdWGCNA on my own project,i encounter error in MetaspotsByGroups()

Steps to reproduce i used another public spatial transcriptomic data,Stereoseq of a half brain section of adult mouse from MOSTA(Cell,2022),and here is the link https://db.cngb.org/stomics/mosta/download/ the seurat object was bulit manually so i suspect the error may be caused by the procedure of the seurat object. i will be appreciate any suggestion from you!

# Put code in this box
seurat_obj <- MetaspotsByGroups(
  seurat_obj,
  group.by = c("orig.ident"),
  ident.group = "orig.ident",
  assay = 'Spatial',
  slot = 'counts',min_spots = 30
)

**R session info**
R version 4.2.0 (2022-04-22 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=Chinese (Simplified)_China.utf8  LC_CTYPE=Chinese (Simplified)_China.utf8   
[3] LC_MONETARY=Chinese (Simplified)_China.utf8 LC_NUMERIC=C                               
[5] LC_TIME=Chinese (Simplified)_China.utf8    

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

other attached packages:
[1] hdWGCNA_0.2.24        harmony_1.1.0         Rcpp_1.0.10           WGCNA_1.72-1         
[5] fastcluster_1.2.3     dynamicTreeCut_1.63-1 SeuratObject_4.1.3    Seurat_4.3.0         

loaded via a namespace (and not attached):
  [1] utf8_1.2.4              spatstat.explore_3.2-1  reticulate_1.31        
  [4] tidyselect_1.2.0        RSQLite_2.3.1           AnnotationDbi_1.62.2   
  [7] htmlwidgets_1.6.2       grid_4.2.0              Rtsne_0.16             
 [10] devtools_2.4.5          munsell_0.5.0           codetools_0.2-19       
 [13] ica_1.0-3               preprocessCore_1.62.1   future_1.33.0          
 [16] miniUI_0.1.1.1          withr_2.5.2             tester_0.1.7           
 [19] spatstat.random_3.1-5   colorspace_2.1-0        progressr_0.14.0       
 [22] Biobase_2.62.0          knitr_1.45              rstudioapi_0.15.0      
 [25] stats4_4.2.0            ROCR_1.0-11             tensor_1.5             
 [28] listenv_0.9.0           labeling_0.4.3          GenomeInfoDbData_1.2.11
 [31] polyclip_1.10-4         farver_2.1.1            bit64_4.0.5            
 [34] parallelly_1.36.0       vctrs_0.6.1             generics_0.1.3         
 [37] xfun_0.41               R6_2.5.1                doParallel_1.0.17      
 [40] GenomeInfoDb_1.38.1     bitops_1.0-7            spatstat.utils_3.0-3   
 [43] cachem_1.0.7            promises_1.2.0.1        scales_1.2.1           
 [46] nnet_7.3-19             gtable_0.3.4            globals_0.16.2         
 [49] processx_3.8.2          goftest_1.2-3           rlang_1.1.2            
 [52] splines_4.2.0           lazyeval_0.2.2          impute_1.74.1          
 [55] spatstat.geom_3.2-4     checkmate_2.3.0         BiocManager_1.30.22    
 [58] reshape2_1.4.4          abind_1.4-5             backports_1.4.1        
 [61] httpuv_1.6.9            Hmisc_5.1-1             tools_4.2.0            
 [64] usethis_2.2.2           ggplot2_3.4.4           ellipsis_0.3.2         
 [67] RColorBrewer_1.1-3      proxy_0.4-27            BiocGenerics_0.48.1    
 [70] sessioninfo_1.2.2       ggridges_0.5.4          plyr_1.8.9             
 [73] base64enc_0.1-3         zlibbioc_1.48.0         purrr_1.0.1            
 [76] RCurl_1.98-1.13         ps_1.7.5                prettyunits_1.2.0      
 [79] rpart_4.1.19            deldir_1.0-9            pbapply_1.7-2          
 [82] cowplot_1.1.1           urlchecker_1.0.1        S4Vectors_0.40.1       
 [85] zoo_1.8-12              ggrepel_0.9.4           cluster_2.1.4          
 [88] fs_1.6.3                magrittr_2.0.3          data.table_1.14.8      
 [91] scattermore_1.2         lmtest_0.9-40           RANN_2.6.1             
 [94] fitdistrplus_1.1-11     matrixStats_1.1.0       pkgload_1.3.3          
 [97] patchwork_1.1.3         mime_0.12               evaluate_0.23          
[100] xtable_1.8-4            IRanges_2.36.0          gridExtra_2.3          
[103] compiler_4.2.0          tibble_3.2.1            KernSmooth_2.23-22     
[106] crayon_1.5.2            htmltools_0.5.5         later_1.3.0            
[109] Formula_1.2-5           tidyr_1.3.0             DBI_1.1.3              
[112] corrplot_0.92           MASS_7.3-60             Matrix_1.5-3           
[115] cli_3.4.0               parallel_4.2.0          igraph_1.5.1           
[118] pkgconfig_2.0.3         foreign_0.8-85          sp_2.1-1               
[121] plotly_4.10.3           spatstat.sparse_3.0-2   foreach_1.5.2          
[124] XVector_0.42.0          stringr_1.5.0           callr_3.7.3            
[127] digest_0.6.31           sctransform_0.4.1       RcppAnnoy_0.0.21       
[130] spatstat.data_3.0-3     Biostrings_2.68.1       rmarkdown_2.25         
[133] leiden_0.4.3            htmlTable_2.4.2         uwot_0.1.16            
[136] shiny_1.7.5.1           lifecycle_1.0.4         nlme_3.1-162           
[139] jsonlite_1.8.7          viridisLite_0.4.2       fansi_1.0.5            
[142] pillar_1.9.0            lattice_0.21-8          KEGGREST_1.40.0        
[145] fastmap_1.1.1           httr_1.4.7              pkgbuild_1.4.2         
[148] survival_3.5-5          GO.db_3.17.0            glue_1.6.2             
[151] remotes_2.4.2.1         png_0.1-8               iterators_1.0.14       
[154] bit_4.0.5               stringi_1.7.12          profvis_0.3.8          
[157] blob_1.2.4              memoise_2.0.1           dplyr_1.1.3            
[160] irlba_2.3.5.1           future.apply_1.11.0    

Screenshots 1700055543022

traceback 7: stop("invalid 'data'") 6: Matrix(from, sparse = TRUE, doDiag = FALSE) 5: asMethod(object) 4: as(agg_X, "sparseMatrix") 3: (function (cur_seurat, mode = "sum", assay = "Spatial", slot = "counts") { X <- GetAssayData(cur_seurat, slot = "counts") if (sum(unlist(lapply(Images(cur_seurat), function(x) { nrow(cur_seurat@images[[x]]@coordinates) != 0 }))) != 1) { stop("More than one sample present in grouping. Please specify a metadata column with group.by indicating different ST samples.") } row_range <- range(cur_seurat$row) col_range <- range(cur_seurat$col) col_bounds <- col_range[1]:col_range[2] if (length(col_bounds) >= 4) { col_bounds <- col_bounds[which(1:length(col_bounds)%%4 == 0)] } else if (length(col_bounds) == 3) { warning("The selected grouping is spatially constrained and might fail the metaspot aggregation step. You may wish to change the group.by parameter to form larger groups.") col_bounds <- col_bounds[2] } else if (length(col_bounds) >= 1) { ... 2: mapply(ConstructMetaspots, cur_seurat = seurat_list, MoreArgs = list(mode = mode, assay = assay, slot = slot)) 1: MetaspotsByGroups(seurat_obj, group.by = c("orig.ident"), ident.group = "orig.ident", assay = "Spatial", slot = "counts", min_spots = 30)

smorabit commented 6 months ago

Thank you for taking the time to write your issue. Can you please provide the code that you used to create the Seurat object?

weiyouzi321 commented 6 months ago

thank you for your attention! there two ways to build my own Seurat object;the object with which i encounter to error was built by first way,i will try whether the object built by second way can sovle this error first way is to use the transform code to trans .h5ad to rds and then use readRDS to get Seurat object

library(SeuratDisk)
library(Seurat)
library(ggplot2)
library(dplyr)
library(rjson)
library(argparser)
args <- arg_parser("Converting h5ad file(.h5ad) to RDS.")
args <- add_argument(args, "--infile", help = "input .h5ad file")
args <- add_argument(args, "--outfile", help = "output RDS file")
argv <- parse_args(args)

infile <- argv$infile
outfile <- argv$outfile

print("Converting h5ad to h5seurat file.")
Convert(infile, dest = "h5seurat",assay = "Spatial", overwrite = TRUE)
h5file<-paste(paste(unlist(strsplit(infile,"h5ad",fixed=TRUE)),collapse='h5ad'),"h5seurat",sep="")
print("Loading h5seurat file.")
Stdata <- LoadH5Seurat(h5file)

if(!is.null(Stdata@reductions$ignore)){
  Stdata@reductions$ignore <- NULL
}
if(!is.null(Stdata@reductions$spatial)){
  Stdata@reductions$spatial <- NULL
}
if(!is.null(Stdata@misc$raw_counts)){
  if(!is.null(Stdata@misc$raw_genename)){
    Stdata@misc$raw_counts@Dimnames[[1]]  = Stdata@misc$raw_genename
    Stdata@misc$raw_genename<-NULL
  }
  if(!is.null(Stdata@misc$raw_cellname)){
    Stdata@misc$raw_counts@Dimnames[[2]]  = Stdata@misc$raw_cellname
    Stdata@misc$raw_cellname<-NULL
  }
  print("Adding misc raw_counts to Spatial assay as counts.")
  Stdata@assays$Spatial@counts=Stdata@misc$raw_counts
  if(dim(Stdata[['Spatial']]@scale.data)[1]>0){
    print("Adding misc raw_counts to Spatial assay as data.")
    Stdata@assays$Spatial@data=Stdata@misc$raw_counts
  }
  Stdata@misc$raw_counts<-NULL
}
if(!is.null(Stdata@misc$sct_data)){
  print("Creating SCT assay.")
  Stdata@assays$SCT=Stdata@assays$Spatial
  Stdata@assays$SCT@data=Stdata@misc$sct_data
  Stdata@misc$sct_data<-NULL
  Stdata@assays$Spatial@scale.data=matrix(,nrow=10,ncol=0)

  f_len<-length(Stdata@assays$SCT@meta.features)
  if(f_len>0){
    Stdata[['SCT']]@meta.features <- Stdata[['SCT']]@meta.features[,-(1:f_len)]
  }
  if(!is.null(Stdata@misc$sct_genename)){
    print("Adding gene dimnames to SCT data.")
    Stdata@assays$SCT@data@Dimnames[[1]]  = Stdata@misc$sct_genename
  }
  if(!is.null(Stdata@misc$sct_cellname)){
    print("Adding cell dimnames to SCT data.")
    Stdata@assays$SCT@data@Dimnames[[2]]  = Stdata@misc$sct_cellname
  }
  #Stdata@assays$SCT@data@Dimnames=Stdata@assays$Spatial@data@Dimnames

  Stdata@assays$SCT@counts=Stdata@misc$sct_counts
  Stdata@misc$sct_counts<-NULL
  if(!is.null(Stdata@misc$sct_genename)){
    print("Adding gene dimnames to SCT counts.")
    Stdata@assays$SCT@counts@Dimnames[[1]]  = Stdata@misc$sct_genename
  }
  if(!is.null(Stdata@misc$sct_cellname)){
    print("Adding cell dimnames to SCT counts.")
    Stdata@assays$SCT@counts@Dimnames[[2]]  = Stdata@misc$sct_cellname
  }
  #Stdata@assays$SCT@counts@Dimnames=Stdata@assays$Spatial@counts@Dimnames
  DefaultAssay(Stdata)<-"SCT"
  Stdata@misc$sct_genename<-NULL
  Stdata@misc$sct_cellname<-NULL
}
## add image
cell_coords=unique(Stdata@meta.data[, c('x', 'y')])
cell_coords['cell']=row.names(cell_coords)
cell_coords$x <- cell_coords$x - min(cell_coords$x) + 1
cell_coords$y <- cell_coords$y - min(cell_coords$y) + 1
# object of images$slice1@image, all illustrated as 1 since no concrete pic
tissue_lowres_image <- matrix(1, max(cell_coords$y), max(cell_coords$x))

# object of images$slice1@coordinates, concrete coordinate of X and Y
tissue_positions_list <- data.frame(row.names = cell_coords$cell,
                                    tissue = 1,
                                    row = cell_coords$y, col = cell_coords$x,
                                    imagerow = cell_coords$y, imagecol = cell_coords$x)
##@images$slice1@scale.factors
scalefactors_json <- toJSON(list(fiducial_diameter_fullres = 1,tissue_hires_scalef = 1,tissue_lowres_scalef = 1))

# generate object @images$slice1
generate_BGI_spatial <- function(image, scale.factors, tissue.positions, filter.matrix = TRUE){
  if (filter.matrix) {
    tissue.positions <- tissue.positions[which(tissue.positions$tissue == 1), , drop = FALSE]
  }
  unnormalized.radius <- scale.factors$fiducial_diameter_fullres * scale.factors$tissue_lowres_scalef
  spot.radius <- unnormalized.radius / max(dim(x = image))
  return(new(Class = 'VisiumV1',
             image = image,
             scale.factors = scalefactors(spot = scale.factors$tissue_hires_scalef,
                                          fiducial = scale.factors$fiducial_diameter_fullres,
                                          hires = scale.factors$tissue_hires_scalef,
                                          lowres = scale.factors$tissue_lowres_scalef),
             coordinates = tissue.positions,
             spot.radius = spot.radius))
}

BGI_spatial <- generate_BGI_spatial(image = tissue_lowres_image,
                                    scale.factors = fromJSON(scalefactors_json),
                                    tissue.positions = tissue_positions_list)
# can be thought of as a background of spatial
#' import image into seurat object
Stdata@images[['slice1']] <-BGI_spatial
Stdata@images$slice1@key<-"slice1_"
Stdata@images$slice1@assay<-"Spatial"
# conversion done, save
saveRDS(Stdata,outfile)
print("Finished RDS.")

another way is to read [.h5ad] by scanpy and built Seurat object from original matrix

library("reticulate")
  unloadNamespace("ggplot2") 
 library("scater")
    #library("Seurat")
  Sys.setenv(RETICULATE_PYTHON="D:\\Anaconda3\\python.exe")
  sc<-import("scanpy")
  pd<-import("pandas")
  adata <- sc$read_h5ad("D:\\HDMouse\\2023bin100\\cellbin\\Mouse_brain_cell_bin.h5ad")
  adata <- sc$read_h5ad("E:\\Mouse_embroys\\data\\E16.5_E1S3_cell_bin_whole_brain.h5ad")

  counts1 <- t(adata$layers["counts"])
  colnames(counts1) <- adata$obs_names$to_list()
  rownames(counts1) <- adata$var_names$to_list()

  counts1 <- Matrix::Matrix(as.matrix(counts1), sparse = T)

  seurat <- CreateSeuratObject(counts1)
  seurat <- AddMetaData(seurat,adata$obsm['spatial'])
  seurat <- AddMetaData(seurat,adata$uns['sim anno_colors'])
smorabit commented 5 months ago

Based on this information it seems like this issue is outside of the scope of hdWGCNA since it has to do with the way that you have chosen to construct your Seurat object rather than to do with the behavior of hdWGCNA. Personally I have not had success with packages that try to convert betweeen scanpy/R formats like SeuratDisk for example.