SpaCET.deconvolution.matched.scRNAseq hanging at specific cell type

[ahwanpandey]

SpaCET/R/extensions.R

Line 496 in 0199d93

markers <- pbmcapply::pbmclapply(setdiff(cellTypes_level_1,cellType), run_limma, mc.cores=coreNo)

I think this step is hanging at reference generation at "Endothelial Cells" when I use my custom reference. I am not sure why! Please help.

I have attached the relevant objects in the Dropbox link below:
https://www.dropbox.com/scl/fo/erizsk0pighy3cbv32q7p/AENZqTpFtZA_sNpchJsvs0s?rlkey=4jyd8bimmcudrgseee1hyfu4l&dl=0

And the code I ran is:

SpaCET_obj <- SpaCET.deconvolution.matched.scRNAseq(
      SpaCET_obj, 
      sc_counts=sc_counts, 
      sc_annotation=sc_annotation, 
      sc_lineageTree=sc_lineageTree, 
      sc_includeMalignant = TRUE,
      coreNo = 8,
      sc_nCellEachLineage = 100,
      sc_downsampling = TRUE # make sure to downsample tha data beforehand!
    )

Thanks!
Ahwan

[ahwanpandey]

OK I removed "Endothelial Cells" and now its stuck at "Fibroblasts". Also removed that and then its stuck at "Macrophages". It seems to be getting stuck at the third Cell Type every time. I also tried with 1 thread, same issue.

Here is my sessionInfo

R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Rocky Linux 9.6 (Blue Onyx)

Matrix products: default
BLAS/LAPACK: FlexiBLAS OPENBLAS-OPENMP;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C               LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8     LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8    LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       

time zone: Australia/Melbourne
tzcode source: system (glibc)

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

other attached packages:
[1] limma_3.62.2       readxl_1.4.5       dplyr_1.1.4        data.table_1.17.2  qs_0.27.3          Seurat_5.4.0       SeuratObject_5.3.0 sp_2.2-0           SpaCET_1.4.0      

loaded via a namespace (and not attached):
  [1] deldir_2.0-4           pbapply_1.7-2          gridExtra_2.3          rlang_1.1.6            magrittr_2.0.3         RcppAnnoy_0.0.22       matrixStats_1.5.0      ggridges_0.5.6         compiler_4.4.2        
 [10] spatstat.geom_3.3-6    png_0.1-8              vctrs_0.6.5            reshape2_1.4.4         stringr_1.5.1          pkgconfig_2.0.3        fastmap_1.2.0          promises_1.3.2         purrr_1.0.4           
 [19] jsonlite_2.0.0         goftest_1.2-3          later_1.4.2            spatstat.utils_3.1-3   irlba_2.3.5.1          parallel_4.4.2         cluster_2.1.6          R6_2.6.1               ica_1.0-3             
 [28] stringi_1.8.7          RColorBrewer_1.1-3     spatstat.data_3.1-6    reticulate_1.46.0      parallelly_1.44.0      spatstat.univar_3.1-3  cellranger_1.1.0       lmtest_0.9-40          scattermore_1.2       
 [37] Rcpp_1.0.14            tensor_1.5             future.apply_1.11.3    zoo_1.8-14             sctransform_0.4.2      httpuv_1.6.16          Matrix_1.7-1           splines_4.4.2          igraph_2.1.4          
 [46] tidyselect_1.2.1       rstudioapi_0.17.1      dichromat_2.0-0.1      abind_1.4-8            stringfish_0.16.0      spatstat.random_3.3-3  codetools_0.2-20       miniUI_0.1.2           spatstat.explore_3.4-2
 [55] listenv_0.9.1          lattice_0.22-6         tibble_3.2.1           plyr_1.8.9             shiny_1.10.0           S7_0.2.0               ROCR_1.0-11            Rtsne_0.17             future_1.49.0         
 [64] fastDummies_1.7.5      survival_3.7-0         RcppParallel_5.1.10    polyclip_1.10-7        fitdistrplus_1.2-2     pillar_1.10.2          KernSmooth_2.23-24     plotly_4.10.4          generics_0.1.4        
 [73] RcppHNSW_0.6.0         ggplot2_4.0.2          scales_1.4.0           RApiSerialize_0.1.4    globals_0.18.0         xtable_1.8-4           glue_1.8.0             lazyeval_0.2.2         tools_4.4.2           
 [82] RSpectra_0.16-2        RANN_2.6.2             dotCall64_1.2          cowplot_1.1.3          grid_4.4.2             tidyr_1.3.1            nlme_3.1-166           patchwork_1.3.2        cli_3.6.5             
 [91] spatstat.sparse_3.1-0  spam_2.11-1            viridisLite_0.4.2      uwot_0.2.3             gtable_0.3.6           digest_0.6.37          progressr_0.15.1       ggrepel_0.9.6          htmlwidgets_1.6.4     
[100] farver_2.1.2           htmltools_0.5.8.1      lifecycle_1.0.4        httr_1.4.7             statmod_1.5.0          mime_0.13              MASS_7.3-61           

[ahwanpandey]

I also changed the order to have "B cells" third and "Endothelial" first, now its stuck at B cells.. So weird!

[beibeiru]

Hi, @ahwanpandey ,

It is very weird!

The following code works really well on both my Mac and HPC.

library(SpaCET)

visiumPath <- file.path(system.file(package = "SpaCET"), "extdata/Visium_BC")
SpaCET_obj <- create.SpaCET.object.10X(dataPath = visiumPath, platform = "Visium", organism = "human")

load("my_data.RData")
SpaCET_obj <- SpaCET.deconvolution.matched.scRNAseq(
      SpaCET_obj, 
       sc_counts=sc_counts, 
       sc_annotation=sc_annotation, 
       sc_lineageTree=sc_lineageTree, 
       sc_includeMalignant = TRUE,
       coreNo = 8,
       sc_nCellEachLineage = 100,
       sc_downsampling = TRUE # make sure to downsample tha data beforehand!
     )

I noticed you opened an issue a few months ago.
Just wondering if you ever got it to run successfully?

Best,
Beibei

[ahwanpandey]

Hey @beibeiru ,

So I asked google gemini by providing the code and the issue, and it solved it for me! Please see the attached doc, and the updated code is below. Also it says "rowMeans" is much faster than the "apply" step and has Incorporated that.

google_gemini_fix.docx

generateRef <- function(
    sc.matrix.data = sc.matrix.data,
    sc.matrix.anno = sc.matrix.anno,
    sc.matrix.tree = sc.matrix.tree,
    coreNo = coreNo
)
{
  # 1. LOAD THE LIBRARY AT THE TOP
  library(limma)
  
  # 2. PREVENT CORES JAMMING (Turns off Blas/MKL multi-threading)
  Sys.setenv(MKL_NUM_THREADS = 1)
  Sys.setenv(OPENBLAS_NUM_THREADS = 1)
  Sys.setenv(OMP_NUM_THREADS = 1)
  if (requireNamespace("RhpcBLASctl", quietly = TRUE)) {
    RhpcBLASctl::blas_set_num_threads(1)
    RhpcBLASctl::omp_set_num_threads(1)
  }

  sc.matrix.data.norm <- sweep(sc.matrix.data, 2, Matrix::colSums(sc.matrix.data), "/") * 1e5
  sc.matrix.data.log2 <- log2(sc.matrix.data.norm + 1)

  cellTypes_level_1 <- names(sc.matrix.tree)
  cellTypes_level_1_toBeSplit <- names(sc.matrix.tree)[sapply(sc.matrix.tree, function(x) length(x) != 1)]

  refProfiles <- data.frame(row.names = rownames(sc.matrix.data.norm))
  sigGenes <- list()

  # Define the parallel function OUTSIDE the loop so it compiles cleanly once
  run_limma <- function(cellTypeOther, cellType, sc.matrix.anno, sc.matrix.tree, sc.matrix.data.log2)
  {
    # DO NOT put library(limma) here anymore
    TT <- as.numeric(sc.matrix.anno[, 2] %in% sc.matrix.tree[[cellType]])
    WT <- as.numeric(sc.matrix.anno[, 2] %in% sc.matrix.tree[[cellTypeOther]])
    design <- cbind(TT, WT)
    fit <- lmFit(sc.matrix.data.log2, design)
    cont.matrix <- makeContrasts(TTvsWT = TT - WT, levels = design)
    fit2 <- contrasts.fit(fit, cont.matrix)
    fit2 <- eBayes(fit2)
    res <- topTable(fit2, coef = 1, number = nrow(sc.matrix.data.log2))
    res <- res[order(res[, "t"], decreasing = TRUE), ]
    res <- res[1:500, ]

    rownames(res)[res[, "logFC"] > 0.25 & res[, "adj.P.Val"] < 0.01]
  }

  for (cellType in cellTypes_level_1)
  {
    message(paste0("  Processing: ", cellType))

    # SPEED UP: Replaced slow apply(..., 1, mean) with blazing fast rowMeans
    subset_matrix <- sc.matrix.data.norm[, sc.matrix.anno[, 2] %in% sc.matrix.tree[[cellType]], drop = FALSE]
    refProfiles[[cellType]] <- Matrix::rowMeans(subset_matrix)

    if (length(cellTypes_level_1) > 1) 
    {
      # Explicitly passing variables down into the parallel landscape prevents pbmclapply hanging
      markers <- pbmcapply::pbmclapply(
        setdiff(cellTypes_level_1, cellType), 
        run_limma, 
        cellType = cellType,
        sc.matrix.anno = sc.matrix.anno,
        sc.matrix.tree = sc.matrix.tree,
        sc.matrix.data.log2 = sc.matrix.data.log2,
        mc.cores = coreNo
      )
      
      temp <- table(unlist(markers))
      sigGenes[[cellType]] <- names(temp)[temp >= length(cellTypes_level_1) - 1]
    }

    if (cellType %in% cellTypes_level_1_toBeSplit)
    {
      cellTypeSubs <- sc.matrix.tree[[cellType]]

      for (cellTypeSub in cellTypeSubs)
      {
        # SPEED UP: Replaced slow apply with rowMeans
        sub_subset_matrix <- sc.matrix.data.norm[, sc.matrix.anno[, 2] %in% cellTypeSub, drop = FALSE]
        refProfiles[[cellTypeSub]] <- Matrix::rowMeans(sub_subset_matrix)

        TT <- as.numeric(sc.matrix.anno[, 2] %in% cellTypeSub)
        WT <- as.numeric(sc.matrix.anno[, 2] %in% setdiff(sc.matrix.tree[[cellType]], cellTypeSub))
        design <- cbind(TT, WT)
        fit <- lmFit(sc.matrix.data.log2, design)
        cont.matrix <- makeContrasts(TTvsWT = TT - WT, levels = design)
        fit2 <- contrasts.fit(fit, cont.matrix)
        fit2 <- eBayes(fit2)
        res <- topTable(fit2, coef = 1, number = nrow(sc.matrix.data.log2))
        res <- res[order(res[, "t"], decreasing = TRUE), ]
        res <- res[1:500, ]

        markers <- rownames(res)[res[, "logFC"] > 0.25 & res[, "adj.P.Val"] < 0.01]
        sigGenes[[cellTypeSub]] <- markers
      }
    }
  }

  Ref <- list(refProfiles = refProfiles, sigGenes = sigGenes, lineageTree = sc.matrix.tree)
  return(Ref)
}

I updated your code with the custom code above and it ran to completion. Hopefully you could see what the issue was?

In the word doc you will see that in the end it asks to compare what the different BLAS version was, mine is as follows in my HPC.

> extSoftVersion()["BLAS"]
                       BLAS
"FlexiBLAS OPENBLAS-OPENMP"

Best,
Ahwan