Running into Error: std::bad_alloc while using FindMultiModalNeighbors and DSB

[smk5g5]

Hi,

Here is my code. I keep running into Error: std::bad_alloc and am not able to resolve it.

library(dsb)
library(ComplexHeatmap)
library(ggplot2)
library(pheatmap)
library(scCustomize)
library(Seurat)

# scp sak4832@transfer.rc.hms.harvard.edu:/n/scratch/users/s/sak4832/Cite_seq_latest/majority_doublet_predictions.R .
# apply DSB using full background for each batch separately.  

# neg_adt = a list of protein data indexed by batch (background / empty droplets)
# pos_adt = a list of protein data indexed by batch (cell containing droplets)

# modify to match isotype control names in protein data matrix
# isotypes = c("Mouse IgG2bkIsotype_PROT", "MouseIgG1kappaisotype_PROT",
#             "MouseIgG2akappaisotype_PROT", "RatIgG2bkIsotype_PROT")

myoutdir = paste0('/n/scratch/users/s/sak4832/Cite_seq_latest/DSB_analysis_multisample')

create_dir <- function(inp_dir) {
  if (!dir.exists(inp_dir)){
    dir.create(inp_dir)
  } else {
    print("Dir already exists!")
  }
}

create_dir(myoutdir)

isotype.controls = c("Mouse_IgG1","Mouse_IgG2b","Mouse_IgG2a","Rat_IgG2b","Rat_IgG2a")

# raw = Seurat::Read10X(raw_data)
# cells = Seurat::Read10X(filtered_data)

c = list()
neg_adt = list()
c_rna = list()
cellmd_list = list()
# prot_list = list()
# rna_list = list()
# doublet_removed_srt = readRDS(doublet_removed_seurat)
# selected_barcodes = intersect(rownames(cellmd),Cells(doublet_removed_srt))
# selected_barcodes_df = as.data.frame(selected_barcodes)


mypat = paste0('GBM',1:10)
# c('GBM1','GBM2','GBM3','GBM4','GBM5')


doublet_removed_srt_list <- lapply(mypat, function(x) {
mysrt = readRDS(paste0('/n/scratch/users/s/sak4832/Cite_seq_latest/Seurat_per_sample/',x,'/',x,'.doublet_removed.majority_doublet_predictions.250724.rds'))
return(mysrt)
})

names(doublet_removed_srt_list) = mypat

doublet_removed_cells = lapply(doublet_removed_srt_list,Cells)

for(i in mypat){
	 raw = Seurat::Read10X(paste0('/n/data1/mgh/neuro/petti/lab/data/human/gbm_citeseq/SAMPLES/',i,'/outs/raw_feature_bc_matrix/'))
	cells = Seurat::Read10X(paste0('/n/data1/mgh/neuro/petti/lab/data/human/gbm_citeseq/SAMPLES/',i,'/outs/filtered_feature_bc_matrix/'))
	stained_cells = colnames(cells$`Gene Expression`)
	background = setdiff(colnames(raw$`Gene Expression`), stained_cells)
	prot = raw$`Antibody Capture`
	rna = raw$`Gene Expression`
	mynew_colnames_prot = paste0(i,'_',colnames(rna))
	mynew_colnames_rna = paste0(i,'_',colnames(rna))
	mtgene = grep(pattern = "^MT-", rownames(rna), value = TRUE) # used below
	md = data.frame(
	  rna.size = log10(Matrix::colSums(rna)), 
	  prot.size = log10(Matrix::colSums(prot)), 
	  n.gene = Matrix::colSums(rna > 0), 
	  mt.prop = Matrix::colSums(rna[mtgene, ]) / Matrix::colSums(rna)
	)

	md$drop.class = ifelse(rownames(md) %in% stained_cells, 'cell', 'background')

	md = md[md$rna.size > 0 & md$prot.size > 0, ]

	background_drops = rownames(
	  md[ md$prot.size > 1.5 & 
	      md$prot.size < 3 & 
	      md$rna.size < 2.5, ]
	  ) 
	background.adt.mtx = as.matrix(prot[ , background_drops])

	cellmd = md[md$drop.class == 'cell', ]

	rna.mult = (3*mad(cellmd$rna.size))
	prot.mult = (3*mad(cellmd$prot.size))
	rna.lower = median(cellmd$rna.size) - rna.mult
	rna.upper = median(cellmd$rna.size) + rna.mult
	prot.lower = median(cellmd$prot.size) - prot.mult
	prot.upper = median(cellmd$prot.size) + prot.mult

	qc_cells = rownames(
	  cellmd[cellmd$prot.size > prot.lower & 
	         cellmd$prot.size < prot.upper & 
	         cellmd$rna.size > rna.lower & 
	         cellmd$rna.size < rna.upper & 
	         cellmd$mt.prop < 0.14, ]
	  )

	myselected_cells = intersect(qc_cells,doublet_removed_cells[[i]])

	cell.adt.raw = as.matrix(prot[ , myselected_cells])
	cell.rna.raw = rna[ ,myselected_cells]
	cellmd = cellmd[myselected_cells, ]

	pm = sort(apply(cell.adt.raw, 1, max))
	pm2 = apply(background.adt.mtx, 1, max)

	c[[i]] = cell.adt.raw
	neg_adt[[i]] = background.adt.mtx
	c_rna[[i]] = cell.rna.raw
	cellmd_list[[i]] = cellmd
}


for (i in mypat) {
# colnames(c[[i]]) = rename_columns(mymatrix=c[[i]],mysample=i)
mynewcolnames_c = paste0(i,'_',colnames(c[[i]]))
colnames(c[[i]]) = mynewcolnames_c

mynewcolnames_negadt = paste0(i,'_',colnames(neg_adt[[i]]))
colnames(neg_adt[[i]]) = mynewcolnames_negadt

mynewcolnames_crna = paste0(i,'_',colnames(c_rna[[i]]))
colnames(c_rna[[i]]) = mynewcolnames_crna

myrownames = paste0(i,'_',rownames(cellmd_list[[i]]))
rownames(cellmd_list[[i]]) = myrownames
}
# cells.dsb.norm = DSBNormalizeProtein(
#   cell_protein_matrix = cell.adt.raw, 
#   empty_drop_matrix = background.adt.mtx, 
#   denoise.counts = TRUE,scale.factor = 'standardize',
#   use.isotype.control = TRUE, return.stats = TRUE,isotype.control.name.vec = isotype.controls
#   )

# apply(c[['GBM1']], 1, quantile)
# apply(c[['GBM1']], 1, mean)

# sel_genes = Reduce(intersect,genes_by_sample)
# names(abc_meanexp) = mypat

myprots_selected = lapply(mypat,function(x){
	pmat_grt0 <- apply(c[[x]], 1, max)
	pmat_grt0_sub = pmat_grt0[pmat_grt0>0]
	pmat_background_grt0 <- apply(neg_adt[[x]], 1, max)
	pmat_background_grt0_sub = pmat_background_grt0[pmat_background_grt0>0]

	selected_abs = intersect(names(pmat_grt0_sub),names(pmat_background_grt0_sub))
	return(selected_abs)
	})

names(myprots_selected) = mypat

selected_proteins = Reduce(intersect,myprots_selected)

proteins_to_remove = unique(unlist(lapply(mypat,function(x){
	setdiff(rownames(c[[x]]),selected_proteins)
	})))
# abc_quantile = lapply(mypat,function(x){
# 	apply(c[[x]], 1, quantile)
# 	})

# names(abc_quantile) = mypat


# pm = sort(apply(c[[i]], 1, max))
# pm2 = apply(neg_adt[[i]], 1, max)
# head(pm)

# Proteins without raw data in stained cells 
# (below, the maximum UMI count for one protein was 4, 
# an order of magnitude below even the isotype controls) 
# can be removed from both matrices prior to normalization.
# In many cases, removing proteins is not necessary, but we
#  recommend checking your raw data.
# # flter 
# pm = sort(apply(cell.adt.raw, 1, max))
# pm2 = apply(background.adt.mtx, 1, max)
# head(pm)
# prot pmax
# CD34_TotalSeqB 4
# CD80_TotalSeqB 60
# CD274_TotalSeqB 75
# IgG2b_control_TotalSeqB 90

# cell.adt.raw  = cell.adt.raw[!rownames(cell.adt.raw) == 'CD34_TotalSeqB', ]
# background.adt.mtx = background.adt.mtx[!rownames(background.adt.mtx) == 'CD34_TotalSeqB', ]

# cell.adt.raw  = cell.adt.raw[!rownames(cell.adt.raw) == 'CD34_TotalSeqB', ]
# background.adt.mtx = background.adt.mtx[!rownames(background.adt.mtx) == 'CD34_TotalSeqB', ]


c_modified = lapply(mypat,function(x){
	return(c[[x]][!rownames(c[[x]])==proteins_to_remove,])
	})

names(c_modified) = mypat

neg_adt_modified = lapply(mypat,function(x){
	return(neg_adt[[x]][!rownames(neg_adt[[x]])==proteins_to_remove, ])
	})

names(neg_adt_modified) = mypat

lapply(c_modified,function(x){ print(x[1:4,1:5])})

# return.stats = TRUE



dsb_norm = list()
for (i in mypat) {
  dsb_norm[[i]] = 
    DSBNormalizeProtein(cell_protein_matrix = c_modified[[i]], 
                        empty_drop_matrix = neg_adt_modified[[i]], 
                        denoise.counts = TRUE, 
                        use.isotype.control = TRUE, quantile.clipping = TRUE,
                        isotype.control.name.vec = isotype.controls,return.stats = TRUE)
}


lapply(dsb_norm,function(x){ dim(x$technical_stats)})


dsb_norm_matrixlist = lapply(dsb_norm,function(x){ 
	return(x$dsb_normalized_matrix)
	})

 # denoise.counts = TRUE, 
 #  use.isotype.control = TRUE, 

dsb_multi = do.call(cbind, dsb_norm_matrixlist)

c_rna_multi = do.call(cbind,c_rna)

c_modified_adt_multi = do.call(cbind,c_modified)

myrownames_cellmdmulti = unname(unlist(lapply(cellmd_list,rownames)))
cellmd_multi = do.call(rbind,cellmd_list)
rownames(cellmd_multi) = myrownames_cellmdmulti
# min.cells=10, min.features=100
# Seurat workflow 
library(Seurat)

# # integrating with Seurat
# stopifnot(isTRUE(all.equal(rownames(cellmd), colnames(cell.adt.raw))))
# stopifnot(isTRUE(all.equal(rownames(cellmd), colnames(cell.rna.raw))))

# create Seurat object note: min.cells is a gene filter, not a cell filter
seurat_obj = Seurat::CreateSeuratObject(counts = c_rna_multi, 
                               meta.data = cellmd_multi,
                               assay = "RNA", 
                               min.cells=10, min.features=100,project='CITE_seq_DSBmulti')


seurat_obj_adt2 = Seurat::CreateSeuratObject(counts = c_rna_multi, 
                               meta.data = cellmd_multi,
                               assay = "RNA", 
                               min.cells=10, min.features=100,project='CITE_seq_DSBmulti')


adt_assay <- CreateAssay5Object(counts = c_modified_adt_multi[,Cells(seurat_obj_adt2)])

# add this assay to the previously created Seurat object
seurat_obj_adt2[["ADT"]] <- adt_assay

# Validate that the object now contains multiple assays
# Assays(cbmc)

# Now we cluster cells based on dsb normalized protein levels. 
# Similar to workflow used in our paper Kotliarov et al. 2020
#  we don’t cluster based on principal components from ADT, 
#  instead directly using the normalized values.

seurat_obj_adt2_filt <- subset(seurat_obj_adt2, subset = nCount_ADT > 100)

# seurat_obj = Seurat::CreateSeuratObject(counts = c_rna_multi, 
#                                meta.data = cellmd_multi,
#                                assay = "RNA", 
#                                min.cells=10, min.features=100,project='CITE_seq_DSBmulti')

seurat_obj_sub = subset(seurat_obj,cells=Cells(seurat_obj_adt2_filt))

dsb_multi_sub = dsb_multi[,Cells(seurat_obj_sub)]
# add dsb normalized matrix "cell.adt.dsb" to the "CITE" data (not counts!) slot
seurat_obj_sub[["CITE"]] = Seurat::CreateAssayObject(data = dsb_multi_sub)


# define proteins to use in clustering (non-isptype controls)
prots = rownames(seurat_obj_sub@assays$CITE@data)[1:42]

# cluster and run umap 
seurat_obj_sub = Seurat::FindNeighbors(object = seurat_obj_sub, dims = NULL,assay = 'CITE', 
                          features = prots, k.param = 30, 
                          verbose = FALSE)


# direct graph clustering 
seurat_obj_sub = Seurat::FindClusters(object = seurat_obj_sub, resolution = c(0.1,0.3,0.5, 0.7, 0.9,1,1.5,2), 
                         algorithm = 3, 
                         graph.name = 'CITE_snn', 
                         verbose = T)
# umap (optional)
# s = Seurat::RunUMAP(object = s, assay = "CITE", features = prots,
#                     seed.use = 1990, min.dist = 0.2, n.neighbors = 30,
#                     verbose = FALSE)

# make results dataframe 
d_seurat = cbind(seurat_obj_sub@meta.data, 
          as.data.frame(t(seurat_obj_sub@assays$CITE@data))
          # s@reductions$umap@cell.embeddings)
          )


library(magrittr)
# calculate the median protein expression separately for each cluster 
adt_plot = d_seurat %>% 
  dplyr::group_by(CITE_snn_res.0.5) %>% 
  dplyr::summarize_at(.vars = prots, .funs = median) %>% 
  tibble::remove_rownames() %>% 
  tibble::column_to_rownames("CITE_snn_res.0.5") 
# plot a heatmap of the average dsb normalized values for each cluster
pdf('adt_plot_pheatmap_multi.pdf',width=30, height=40)
print(pheatmap::pheatmap(t(adt_plot), 
                   color = viridis::viridis(25, option = "B"), 
                   fontsize_row = 8, border_color = NA))
dev.off()

DefaultAssay(seurat_obj_sub) = "RNA"
seurat_obj_sub = NormalizeData(seurat_obj_sub, verbose = FALSE) %>% 
  FindVariableFeatures(selection.method = 'vst', verbose = FALSE) %>% 
  ScaleData(verbose = FALSE) %>%
  RunPCA(verbose = FALSE,irlba = FALSE)

# set up dsb values to use in WNN analysis (do not normalize with CLR, use dsb normalized values)
DefaultAssay(seurat_obj_sub) = "CITE"
VariableFeatures(seurat_obj_sub) = prots


seurat_obj_sub = seurat_obj_sub %>% 
  ScaleData(verbose = TRUE) %>% 
  RunPCA(reduction.name = 'apca', verbose = TRUE)

seurat_obj_sub2 = FindMultiModalNeighbors(
  seurat_obj_sub, reduction.list = list("pca", "apca"), 
  dims.list = list(1:10, 1:10), 
  modality.weight.name = "RNA.weight", weighted.nn.name = "weighted.nn",
  snn.graph.name = "wsnn",
  verbose = TRUE
)

I have tried everything that has been suggested here #6473
or here #7265 but nothing seems to work.

I am unable to resolve this issue and would appreciate any help in this matter.

Here is my session info :

sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

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

time zone: Etc/UTC
tzcode source: system (glibc)

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

other attached packages:
[1] scCustomize_2.1.2     Seurat_5.0.0          SeuratObject_5.0.2
[4] sp_2.1-4              pheatmap_1.0.12       ggplot2_3.5.1.9000
[7] ComplexHeatmap_2.18.0 dsb_1.0.4

Kindly let me know how I can resolve it.