diff --git a/DESCRIPTION b/DESCRIPTION
index 0af09d4..fa2d664 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -7,7 +7,7 @@ Description: Cell-cell Interactions at Single-Cell Resolution
 License: GPL-3
 Encoding: UTF-8
 LazyData: yes
-RoxygenNote: 7.1.2
+RoxygenNote: 7.2.3
 VignetteBuilder: knitr
 Depends:
     R (>= 4.0)
diff --git a/R/RunCellToCell.R b/R/RunCellToCell.R
index b293b7a..29b1478 100644
--- a/R/RunCellToCell.R
+++ b/R/RunCellToCell.R
@@ -38,7 +38,7 @@ RunCellToCell <- function(sys.small,
       colnames(subunit.list[[s]]) <- colnames(temp)
       rownames(subunit.list[[s]]) <- rownames(ground.truth$source.subunits)
       non.na.indices <- !is.na(ground.truth$source.subunits[,s]) #Identify rows in the s-th column of the ground truth which are not NA
-      subunit.list[[s]][non.na.indices,] <- as.matrix(temp@assays[[assay]]@data[ground.truth$source.subunits[non.na.indices,s],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+      subunit.list[[s]][non.na.indices,] <- as.matrix(getSeuratAssay(temp,assay,"data")[ground.truth$source.subunits[non.na.indices,s],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
     }
     lig.list[[i]] <- Reduce('*',subunit.list)
     rm(subunit.list)
@@ -54,7 +54,7 @@ RunCellToCell <- function(sys.small,
       colnames(subunit.list[[t]]) <- colnames(temp)
       rownames(subunit.list[[t]]) <- rownames(ground.truth$target.subunits)
       non.na.indices <- !is.na(ground.truth$target.subunits[,t]) #Identify rows in the t-th column of the ground truth which are not NA
-      subunit.list[[t]][non.na.indices,] <- as.matrix(temp@assays[[assay]]@data[ground.truth$target.subunits[non.na.indices,t],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+      subunit.list[[t]][non.na.indices,] <- as.matrix(getSeuratAssay(temp,assay,"data")[ground.truth$target.subunits[non.na.indices,t],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
     }
     rec.list[[i]] <- Reduce('*',subunit.list)
     rm(subunit.list)
@@ -108,6 +108,13 @@ RunCellToCell <- function(sys.small,
   
   #Use this matrix to create a Seurat object:
   demo <- Seurat::CreateSeuratObject(counts = as.matrix(scc),assay = 'CellToCell')
+  # JC: Seurat V5 will not create data slot automatically, the following step is to manually add this slot
+  if(SeuratObject::Version(demo) >= 5){
+    demo <- NormalizeData(demo,assay = "CellToCell")  # Seura Object need to be >= 5.0.1
+    demo@assays$CellToCell@layers$data <- demo@assays$CellToCell@layers$counts # Seura Object need to be>= 5.0.1
+    
+  }
+  
   
   # Gather and assemble metadata based on "ident" slot
   sending.cell.idents.2 <- do.call(c,sending.cell.idents)
@@ -139,7 +146,8 @@ RunCellToCell <- function(sys.small,
   # Compile
   meta.data.to.add.also <- cbind(sending.metadata,receiving.metadata,joint.metadata)
   rownames(meta.data.to.add.also) <- paste(sending.barcodes,receiving.barcodes,sep='—')
-  # Add additional metadata
+  # Add additional metadata 
+  ### COMMENT 2024-02-03 MSBR: Is this compatible with Seurat v5?###
   demo <- Seurat::AddMetaData(demo,metadata = as.data.frame(meta.data.to.add.also))
   }
   # Define initial identity
@@ -152,7 +160,7 @@ RunCellToCell <- function(sys.small,
   else{
     output_list <- vector(mode = "list",length=2)
     names(output_list) <- c("CellToCellMatrix","metadata")
-    output_list[["CellToCellMatrix"]] <- demo[["CellToCell"]]@counts
+    output_list[["CellToCellMatrix"]] <- getSeuratAssay(demo,"CellToCell","counts")
     output_list[["metadata"]] <- demo@meta.data
     return(output_list)
   }
diff --git a/R/RunCellToCellSpatial.R b/R/RunCellToCellSpatial.R
index 25fedc7..9230ddd 100644
--- a/R/RunCellToCellSpatial.R
+++ b/R/RunCellToCellSpatial.R
@@ -26,11 +26,11 @@ RunCellToCellSpatial <- function(sys.small,
 
     subunit.list <- list() # Builds sending (ligand) data for any number of ligand subunits
     for (s in 1:ncol(ground.truth$source.subunits)){ #For each subunit column...
-      subunit.list[[s]] <- matrix(data = 1,nrow = nrow(ground.truth$source.subunits),ncol = ncol(sys.small@assays[[assay]]@data[,edgelist$from])) #initialize a mechanism x barcode matrix of all NAs
-      colnames(subunit.list[[s]]) <- colnames(sys.small@assays[[assay]]@data[,edgelist$from])
+      subunit.list[[s]] <- matrix(data = 1,nrow = nrow(ground.truth$source.subunits),ncol = ncol(getSeuratAssay(sys.small,assay,"data")[,edgelist$from])) #initialize a mechanism x barcode matrix of all NAs
+      colnames(subunit.list[[s]]) <- colnames(getSeuratAssay(sys.small,assay,"data")[,edgelist$from])
       rownames(subunit.list[[s]]) <- rownames(ground.truth$source.subunits)
       non.na.indices <- !is.na(ground.truth$source.subunits[,s]) #Identify rows in the s-th column of the ground truth which are not NA
-      subunit.list[[s]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$source.subunits[non.na.indices,s],edgelist$from])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+      subunit.list[[s]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$source.subunits[non.na.indices,s],edgelist$from])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
     }
     lig.data <- Reduce('*',subunit.list)
     rm(subunit.list)
@@ -41,11 +41,11 @@ RunCellToCellSpatial <- function(sys.small,
     
     subunit.list <- list() # Builds receiving (receptor) data for any number of receptor subunits
     for (t in 1:ncol(ground.truth$target.subunits)){
-      subunit.list[[t]] <- matrix(data = 1,nrow = nrow(ground.truth$target.subunits),ncol = ncol(sys.small@assays[[assay]]@data[,edgelist$to])) #initialize a mechanism x barcode matrix of all NAs
-      colnames(subunit.list[[t]]) <- colnames(sys.small@assays[[assay]]@data[,edgelist$to])
+      subunit.list[[t]] <- matrix(data = 1,nrow = nrow(ground.truth$target.subunits),ncol = ncol(getSeuratAssay(sys.small,assay,"data")[,edgelist$to])) #initialize a mechanism x barcode matrix of all NAs
+      colnames(subunit.list[[t]]) <- colnames(getSeuratAssay(sys.small,assay,"data")[,edgelist$to])
       rownames(subunit.list[[t]]) <- rownames(ground.truth$target.subunits)
       non.na.indices <- !is.na(ground.truth$target.subunits[,t]) #Identify rows in the t-th column of the ground truth which are not NA
-      subunit.list[[t]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$target.subunits[non.na.indices,t],edgelist$to])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+      subunit.list[[t]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$target.subunits[non.na.indices,t],edgelist$to])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
     }
     rec.data <- Reduce('*',subunit.list)
     rm(subunit.list)
@@ -61,6 +61,13 @@ RunCellToCellSpatial <- function(sys.small,
   # Use this matrix to create a Seurat object:
   demo <- Seurat::CreateSeuratObject(counts = as.matrix(scc),assay = 'CellToCellSpatial')
 
+  # JC: Seurat V5 will not create data slot automatically, the following step is to manually add this slot
+  if(SeuratObject::Version(demo) >= 5){
+    demo <- NormalizeData(demo,assay = "CellToCellSpatial")  # Seura Object need to be >= 5.0.1
+    demo@assays$CellToCellSpatial@layers$data <- demo@assays$CellToCellSpatial@layers$counts # Seura Object need to be >= 5.0.1
+    
+  }
+  
   # Add key metadata
 
   meta.data.to.add <- data.frame(SendingType = sending.cell.idents,
@@ -73,7 +80,7 @@ RunCellToCellSpatial <- function(sys.small,
 
   #Add metadata to the Seurat object
   demo <- Seurat::AddMetaData(demo,metadata = meta.data.to.add)
-
+  
   # Gather and assemble additional metadata
   if (!is.null(meta.data.to.map)){
     # Identify sending and receiving barcodes
@@ -107,7 +114,7 @@ RunCellToCellSpatial <- function(sys.small,
   else{
     output_list <- vector(mode = "list",length=2)
     names(output_list) <- c("CellToCellSpatialMatrix","metadata")
-    output_list[["CellToCellSpatialMatrix"]] <- demo[["CellToCellSpatial"]]@counts
+    output_list[["CellToCellSpatialMatrix"]] <- getSeuratAssay(demo,"CellToCellSpatial","counts")
     output_list[["metadata"]] <- demo@meta.data
     return(output_list)
   }
diff --git a/R/RunCellToNeighborhood.R b/R/RunCellToNeighborhood.R
index f5db948..9779572 100644
--- a/R/RunCellToNeighborhood.R
+++ b/R/RunCellToNeighborhood.R
@@ -25,11 +25,11 @@ RunCellToNeighborhood <- function(sys.small,
   
   subunit.list <- list() # Builds sending (ligand) data for any number of ligand subunits
   for (s in 1:ncol(ground.truth$source.subunits)){ #For each subunit column...
-    subunit.list[[s]] <- matrix(data = 1,nrow = nrow(ground.truth$source.subunits),ncol = ncol(sys.small@assays[[assay]]@data[,edgelist$from])) #initialize a mechanism x barcode matrix of all NAs
-    colnames(subunit.list[[s]]) <- colnames(sys.small@assays[[assay]]@data[,edgelist$from])
+    subunit.list[[s]] <- matrix(data = 1,nrow = nrow(ground.truth$source.subunits),ncol = ncol(getSeuratAssay(sys.small,assay,"data")[,edgelist$from])) #initialize a mechanism x barcode matrix of all NAs
+    colnames(subunit.list[[s]]) <- colnames(getSeuratAssay(sys.small,assay,"data")[,edgelist$from])
     rownames(subunit.list[[s]]) <- rownames(ground.truth$source.subunits)
     non.na.indices <- !is.na(ground.truth$source.subunits[,s]) #Identify rows in the s-th column of the ground truth which are not NA
-    subunit.list[[s]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$source.subunits[non.na.indices,s],edgelist$from])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+    subunit.list[[s]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$source.subunits[non.na.indices,s],edgelist$from])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
   }
   lig.data <- Reduce('*',subunit.list)
   rm(subunit.list)
@@ -40,11 +40,11 @@ RunCellToNeighborhood <- function(sys.small,
   
   subunit.list <- list() # Builds receiving (receptor) data for any number of receptor subunits
   for (t in 1:ncol(ground.truth$target.subunits)){
-    subunit.list[[t]] <- matrix(data = 1,nrow = nrow(ground.truth$target.subunits),ncol = ncol(sys.small@assays[[assay]]@data[,edgelist$to])) #initialize a mechanism x barcode matrix of all NAs
-    colnames(subunit.list[[t]]) <- colnames(sys.small@assays[[assay]]@data[,edgelist$to])
+    subunit.list[[t]] <- matrix(data = 1,nrow = nrow(ground.truth$target.subunits),ncol = ncol(getSeuratAssay(sys.small,assay,"data")[,edgelist$to])) #initialize a mechanism x barcode matrix of all NAs
+    colnames(subunit.list[[t]]) <- colnames(getSeuratAssay(sys.small,assay,"data")[,edgelist$to])
     rownames(subunit.list[[t]]) <- rownames(ground.truth$target.subunits)
     non.na.indices <- !is.na(ground.truth$target.subunits[,t]) #Identify rows in the t-th column of the ground truth which are not NA
-    subunit.list[[t]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$target.subunits[non.na.indices,t],edgelist$to])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+    subunit.list[[t]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$target.subunits[non.na.indices,t],edgelist$to])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
   }
   rec.data <- Reduce('*',subunit.list)
   rm(subunit.list)
@@ -67,7 +67,13 @@ RunCellToNeighborhood <- function(sys.small,
 
   # Use this matrix to create a Seurat object:
   demo <- Seurat::CreateSeuratObject(counts = as.matrix(scc),assay = 'CellToNeighborhood')
-
+  # JC: Seurat V5 will not create data slot automatically, the following step is to manually add this slot
+  if(SeuratObject::Version(demo) >= 5){
+    demo <- NormalizeData(demo,assay = "CellToNeighborhood")  # Seura Object need to be >= 5.0.1
+    demo@assays$CellToNeighborhood@layers$data <- demo@assays$CellToNeighborhood@layers$counts # Seura Object need to be >= 5.0.1
+    
+  }
+  
   # Add metadata based on ident slot
   demo <- Seurat::AddMetaData(demo,metadata = barcodes,col.name = 'SendingCell')
   # bug fix: add the Neighborhood - prefix
@@ -107,7 +113,7 @@ RunCellToNeighborhood <- function(sys.small,
   else{
     output_list <- vector(mode = "list",length=2)
     names(output_list) <- c("CellToNeighborhoodMatrix","metadata")
-    output_list[["CellToNeighborhoodMatrix"]] <- demo[["CellToNeighborhood"]]@counts
+    output_list[["CellToNeighborhoodMatrix"]] <- getSeuratAssay(demo,"CellToNeighborhood","counts")
     output_list[["metadata"]] <- demo@meta.data
     return(output_list)
   }
diff --git a/R/RunCellToSystem.R b/R/RunCellToSystem.R
index e5f1e04..3ca87c8 100644
--- a/R/RunCellToSystem.R
+++ b/R/RunCellToSystem.R
@@ -37,7 +37,7 @@ RunCellToSystem <- function(sys.small,
     colnames(subunit.list[[t]]) <- colnames(sys.small)
     rownames(subunit.list[[t]]) <- rownames(ground.truth$target.subunits)
     non.na.indices <- !is.na(ground.truth$target.subunits[,t]) #Identify rows in the s-th column of the ground truth which are not NA
-    subunit.list[[t]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$target.subunits[non.na.indices,t],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+    subunit.list[[t]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$target.subunits[non.na.indices,t],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
   }
   rec.map <- Reduce('*',subunit.list)
   rm(subunit.list)
@@ -65,7 +65,7 @@ RunCellToSystem <- function(sys.small,
     colnames(subunit.list[[s]]) <- colnames(sys.small)
     rownames(subunit.list[[s]]) <- rownames(ground.truth$source.subunits)
     non.na.indices <- !is.na(ground.truth$source.subunits[,s]) #Identify rows in the s-th column of the ground truth which are not NA
-    subunit.list[[s]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$source.subunits[non.na.indices,s],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+    subunit.list[[s]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$source.subunits[non.na.indices,s],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
   }
   lig.map <- Reduce('*',subunit.list)
   rm(subunit.list)
@@ -82,6 +82,12 @@ RunCellToSystem <- function(sys.small,
   
   # Use this matrix to create a Seurat object:
   demo <- Seurat::CreateSeuratObject(counts = as.matrix(sc.connectome),assay = 'CellToSystem')
+  # JC: Seurat V5 will not create data slot automatically, the following step is to manually add this slot
+  if(SeuratObject::Version(demo) >= 5){
+    demo <- NormalizeData(demo,assay = "CellToSystem")  # Seura Object need to be >= 5.0.1
+    demo@assays$CellToSystem@layers$data <- demo@assays$CellToSystem@layers$counts # Seura Object need to be >= 5.0.1
+    
+  }
   
   # Add metadata to the Seurat object
   meta.data.to.add <- data.frame(as.character(colnames(lig.map)))
@@ -126,7 +132,7 @@ RunCellToSystem <- function(sys.small,
   else{
     output_list <- vector(mode = "list",length=2)
     names(output_list) <- c("CellToSystemMatrix","metadata")
-    output_list[["CellToSystemMatrix"]] <- demo[["CellToSystem"]]@counts
+    output_list[["CellToSystemMatrix"]] <- getSeuratAssay(demo,"CellToSystem","counts")
     output_list[["metadata"]] <- demo@meta.data
     return(output_list)
   }
diff --git a/R/RunNeighborhoodToCell.R b/R/RunNeighborhoodToCell.R
index a768288..183198b 100644
--- a/R/RunNeighborhoodToCell.R
+++ b/R/RunNeighborhoodToCell.R
@@ -25,11 +25,11 @@ RunNeighborhoodToCell <- function(sys.small,
   
   subunit.list <- list() # Builds sending (ligand) data for any number of ligand subunits
   for (s in 1:ncol(ground.truth$source.subunits)){ #For each subunit column...
-    subunit.list[[s]] <- matrix(data = 1,nrow = nrow(ground.truth$source.subunits),ncol = ncol(sys.small@assays[[assay]]@data[,edgelist$from])) #initialize a mechanism x barcode matrix of all NAs
-    colnames(subunit.list[[s]]) <- colnames(sys.small@assays[[assay]]@data[,edgelist$from])
+    subunit.list[[s]] <- matrix(data = 1,nrow = nrow(ground.truth$source.subunits),ncol = ncol(getSeuratAssay(sys.small,assay,"data")[,edgelist$from])) #initialize a mechanism x barcode matrix of all NAs
+    colnames(subunit.list[[s]]) <- colnames(getSeuratAssay(sys.small,assay,"data")[,edgelist$from])
     rownames(subunit.list[[s]]) <- rownames(ground.truth$source.subunits)
     non.na.indices <- !is.na(ground.truth$source.subunits[,s]) #Identify rows in the s-th column of the ground truth which are not NA
-    subunit.list[[s]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$source.subunits[non.na.indices,s],edgelist$from])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+    subunit.list[[s]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$source.subunits[non.na.indices,s],edgelist$from])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
   }
   lig.data <- Reduce('*',subunit.list)
   rm(subunit.list)
@@ -40,11 +40,11 @@ RunNeighborhoodToCell <- function(sys.small,
   
   subunit.list <- list() # Builds receiving (receptor) data for any number of receptor subunits
   for (t in 1:ncol(ground.truth$target.subunits)){
-    subunit.list[[t]] <- matrix(data = 1,nrow = nrow(ground.truth$target.subunits),ncol = ncol(sys.small@assays[[assay]]@data[,edgelist$to])) #initialize a mechanism x barcode matrix of all NAs
-    colnames(subunit.list[[t]]) <- colnames(sys.small@assays[[assay]]@data[,edgelist$to])
+    subunit.list[[t]] <- matrix(data = 1,nrow = nrow(ground.truth$target.subunits),ncol = ncol(getSeuratAssay(sys.small,assay,"data")[,edgelist$to])) #initialize a mechanism x barcode matrix of all NAs
+    colnames(subunit.list[[t]]) <- colnames(getSeuratAssay(sys.small,assay,"data")[,edgelist$to])
     rownames(subunit.list[[t]]) <- rownames(ground.truth$target.subunits)
     non.na.indices <- !is.na(ground.truth$target.subunits[,t]) #Identify rows in the t-th column of the ground truth which are not NA
-    subunit.list[[t]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$target.subunits[non.na.indices,t],edgelist$to])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+    subunit.list[[t]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$target.subunits[non.na.indices,t],edgelist$to])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
   }
   rec.data <- Reduce('*',subunit.list)
   rm(subunit.list)
@@ -67,6 +67,12 @@ RunNeighborhoodToCell <- function(sys.small,
   
   # Use this matrix to create a Seurat object:
   demo <- Seurat::CreateSeuratObject(counts = as.matrix(scc),assay = 'NeighborhoodToCell')
+  # JC: Seurat V5 will not create data slot automatically, the following step is to manually add this slot
+  if(SeuratObject::Version(demo) >= 5){
+    demo <- NormalizeData(demo,assay = "NeighborhoodToCell")  # Seura Object need to be >= 5.0.1
+    demo@assays$NeighborhoodToCell@layers$data <- demo@assays$NeighborhoodToCell@layers$counts # Seura Object need to be >= 5.0.1
+    
+  }
   
   # Add metadata based on ident slot
   demo <- Seurat::AddMetaData(demo,metadata = barcodes,col.name = 'ReceivingCell')
@@ -107,7 +113,7 @@ RunNeighborhoodToCell <- function(sys.small,
   else{
     output_list <- vector(mode = "list",length=2)
     names(output_list) <- c("NeighborhoodToCellMatrix","metadata")
-    output_list[["NeighborhoodToCellMatrix"]] <- demo[["NeighborhoodToCell"]]@counts
+    output_list[["NeighborhoodToCellMatrix"]] <- getSeuratAssay(demo,"NeighborhoodToCell","counts")
     output_list[["metadata"]] <- demo@meta.data
     return(output_list)
   }
diff --git a/R/RunSystemToCell.R b/R/RunSystemToCell.R
index 6b08145..e69abab 100644
--- a/R/RunSystemToCell.R
+++ b/R/RunSystemToCell.R
@@ -37,7 +37,7 @@ RunSystemToCell <- function(sys.small,
     colnames(subunit.list[[s]]) <- colnames(sys.small)
     rownames(subunit.list[[s]]) <- rownames(ground.truth$source.subunits)
     non.na.indices <- !is.na(ground.truth$source.subunits[,s]) #Identify rows in the s-th column of the ground truth which are not NA
-    subunit.list[[s]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$source.subunits[non.na.indices,s],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+    subunit.list[[s]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$source.subunits[non.na.indices,s],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
   }
   lig.map <- Reduce('*',subunit.list)
   rm(subunit.list)
@@ -64,7 +64,7 @@ RunSystemToCell <- function(sys.small,
     colnames(subunit.list[[t]]) <- colnames(sys.small)
     rownames(subunit.list[[t]]) <- rownames(ground.truth$target.subunits)
     non.na.indices <- !is.na(ground.truth$target.subunits[,t]) #Identify rows in the s-th column of the ground truth which are not NA
-    subunit.list[[t]][non.na.indices,] <- as.matrix(sys.small@assays[[assay]]@data[ground.truth$target.subunits[non.na.indices,t],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
+    subunit.list[[t]][non.na.indices,] <- as.matrix(getSeuratAssay(sys.small,assay,"data")[ground.truth$target.subunits[non.na.indices,t],])   #For every row in the initialized matrix corresponding to the indices of the ground.truth which are not NA, replace with the rows from the Seurat object corresponding to the genes in the ground.truth at those indices
   }
   rec.map <- Reduce('*',subunit.list)
   rm(subunit.list)
@@ -80,6 +80,12 @@ RunSystemToCell <- function(sys.small,
   
   #Use this matrix to create a Seurat object:
   demo <- Seurat::CreateSeuratObject(counts = as.matrix(sc.connectome),assay = 'SystemToCell')
+  # JC: Seurat V5 will not create data slot automatically, the following step is to manually add this slot
+  if(SeuratObject::Version(demo) >= 5){
+    demo <- NormalizeData(demo,assay = "SystemToCell")  # Seurat Object need to be >= 5.0.1
+    demo@assays$SystemToCell@layers$data <- demo@assays$SystemToCell@layers$counts # Seurat Object need to be >= 5.0.1
+    
+  }
   
   # Add metadata to the Seurat object
   meta.data.to.add <- data.frame(as.character(colnames(rec.map)))
@@ -126,7 +132,7 @@ RunSystemToCell <- function(sys.small,
   else{
     output_list <- vector(mode = "list",length=2)
     names(output_list) <- c("SystemToCellMatrix","metadata")
-    output_list[["SystemToCellMatrix"]] <- demo[["SystemToCell"]]@counts
+    output_list[["SystemToCellMatrix"]] <- getSeuratAssay(demo,"SystemToCell","counts")
     output_list[["metadata"]] <- demo@meta.data
     return(output_list)
   }
diff --git a/R/utils.R b/R/utils.R
index 2c36af4..11b3b02 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,3 +1,15 @@
+#' Wrapper function for GetAssayData based on the version of the SeuratObject
+#'
+#' @param object input Seurat object 
+#' @param assay which assay to use in object
+#' @param slot "data","counts",or "scale.data"
+#'
+#' @return the specified assay data
+#' @export
+getSeuratAssay <- function(object,assay,slot){
+  if(SeuratObject::Version(object) >= 5) return(Seurat::GetAssayData(object,assay=assay,layer=slot))
+  else return(Seurat::GetAssayData(object,assay=assay,slot=slot))
+}
 
 #' Seurat input preprocessing
 #'
@@ -16,6 +28,7 @@ prepSeurat <- function(object,assay,min.cells.per.ident,min.cells.per.gene){
   # Stash object
   sys.small <- object
   
+  
   # Limit object to cell populations larger than requested minimum
   if (!is.null(min.cells.per.ident)){
     message(paste("\n",'Subsetting to populations with greater than',min.cells.per.ident,'cells'))
@@ -28,9 +41,9 @@ prepSeurat <- function(object,assay,min.cells.per.ident,min.cells.per.gene){
     message(paste("\n",'Subsetting to genes expressed in greater than',min.cells.per.gene,'cells'))
     
     # jc: return an error if the count matrix in the given assay is empty
-    if(!length(sys.small@assays[[assay]]@counts)) stop("Unable to subset: the count matrix in the given assay is empty.")
+    if(!length(getSeuratAssay(sys.small,assay,"counts"))) stop("Unable to subset: the count matrix in the given assay is empty.")
     
-    cells.per.gene <- data.frame(non.zero.cells = Matrix::rowSums(sys.small@assays[[assay]]@counts>0))
+    cells.per.gene <- data.frame(non.zero.cells = Matrix::rowSums(getSeuratAssay(sys.small,assay,"counts")>0))
     GOI <- subset(cells.per.gene,non.zero.cells > min.cells.per.gene)
     sys.small <- sys.small[rownames(GOI),]
   }
diff --git a/man/RunCellToSystem.Rd b/man/RunCellToSystem.Rd
index 37a105e..4a5e6fa 100644
--- a/man/RunCellToSystem.Rd
+++ b/man/RunCellToSystem.Rd
@@ -20,7 +20,7 @@ RunCellToSystem(
 
 \item{assay}{The assay to run the CellToSystem transformation on. Defaults to "RNA."}
 
-\item{blend}{Choice of linear operator to combine edges. Defaults to "mean", also accepts "sum"}
+\item{blend}{Choice of linear operator to combine edges. Defaults to "mean", also accepts "sum","mean.adj"}
 
 \item{meta.data.to.map}{A character vector of metadata names present in the original object which will be carried to the NICHES objects}
 
diff --git a/man/RunNICHES.Rd b/man/RunNICHES.Rd
index fe81f5c..613dcbc 100644
--- a/man/RunNICHES.Rd
+++ b/man/RunNICHES.Rd
@@ -41,7 +41,7 @@ RunNICHES(object, ...)
   species,
   min.cells.per.ident = NULL,
   min.cells.per.gene = NULL,
-  meta.data.to.map = NULL,
+  meta.data.to.map = colnames(object@meta.data),
   position.x = NULL,
   position.y = NULL,
   cell_types = NULL,
diff --git a/man/RunSystemToCell.Rd b/man/RunSystemToCell.Rd
index b0ea582..97c73e8 100644
--- a/man/RunSystemToCell.Rd
+++ b/man/RunSystemToCell.Rd
@@ -20,7 +20,7 @@ RunSystemToCell(
 
 \item{assay}{The assay to run the SystemToCell transformation on. Defaults to "RNA."}
 
-\item{blend}{Choice of linear operator to combine edges. Defaults to "mean", also accepts "sum"}
+\item{blend}{Choice of linear operator to combine edges. Defaults to "mean", also accepts "sum", "mean.adj"}
 
 \item{meta.data.to.map}{A character vector of metadata names present in the original object which will be carried to the NICHES objects}
 
diff --git a/man/getSeuratAssay.Rd b/man/getSeuratAssay.Rd
new file mode 100644
index 0000000..7b00437
--- /dev/null
+++ b/man/getSeuratAssay.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{getSeuratAssay}
+\alias{getSeuratAssay}
+\title{Wrapper function for GetAssayData based on the version of the SeuratObject}
+\usage{
+getSeuratAssay(object, assay, slot)
+}
+\arguments{
+\item{object}{input Seurat object}
+
+\item{assay}{which assay to use in object}
+
+\item{slot}{"data","counts",or "scale.data"}
+}
+\value{
+the specified assay data
+}
+\description{
+Wrapper function for GetAssayData based on the version of the SeuratObject
+}
diff --git a/omnipathr-log/omnipathr-20231109-2228.log b/omnipathr-log/omnipathr-20231109-2228.log
new file mode 100644
index 0000000..cee732b
--- /dev/null
+++ b/omnipathr-log/omnipathr-20231109-2228.log
@@ -0,0 +1,2 @@
+[2023-11-09 22:28:09] [INFO]    [OmnipathR] Initialized cache: `/home/rstudio/.cache/OmnipathR`.
+[2023-11-09 22:28:09] [INFO]    [OmnipathR] Welcome to OmnipathR!
diff --git a/omnipathr-log/omnipathr-20231113-0320.log b/omnipathr-log/omnipathr-20231113-0320.log
new file mode 100644
index 0000000..1713278
--- /dev/null
+++ b/omnipathr-log/omnipathr-20231113-0320.log
@@ -0,0 +1,2 @@
+[2023-11-13 03:20:24] [INFO]    [OmnipathR] Initialized cache: `/home/rstudio/.cache/OmnipathR`.
+[2023-11-13 03:20:24] [INFO]    [OmnipathR] Welcome to OmnipathR!
diff --git a/tests/knn_search_test.R b/tests/knn_search_test.R
new file mode 100644
index 0000000..817c17c
--- /dev/null
+++ b/tests/knn_search_test.R
@@ -0,0 +1,42 @@
+# JC: Test the scaling of different knn methods
+# current method: brute-force
+# target method: kd-tree search
+
+library(Seurat)
+library(Matrix)
+devtools::load_all()
+## first simulate a giant ST dataset
+n_cells <- 400
+n_genes <- 20000
+# nb parameters
+mu <- 1
+size <- 20
+# assuming a gaussian 
+set.seed(42)
+
+count_data <- as(matrix(rnbinom(n = (n_cells*n_genes),mu=mu,size=size),nrow = n_genes,ncol = n_cells),"sparseMatrix")
+genes <- read.csv("/data/test_data/test_data/filtered_gene_bc_matrices/hg19/genes.tsv",sep = '\t',header = F)
+rownames(count_data) <- unique(genes$V2)[1:n_genes]
+
+st_data <- CreateSeuratObject(counts = count_data)
+
+st_data@meta.data$x <- c(1:200)
+st_data@meta.data$y <- c(1:200)
+st_data@meta.data$cell_types <- c(rep("cp1",n_cells/2),rep("cp2",n_cells/2))
+st_data <- NormalizeData(st_data)
+
+NICHES_output <- RunNICHES(object = st_data,
+                           LR.database = "fantom5",
+                           species = "human",
+                           assay = "RNA",
+                           position.x = 'x',
+                           position.y = 'y',
+                           k = 4, 
+                           cell_types = "cell_types",
+                           min.cells.per.ident = 0,
+                           min.cells.per.gene = NULL,
+                           meta.data.to.map = NULL,
+                           CellToCell = T,CellToSystem = F,SystemToCell = F,
+                           CellToCellSpatial = F,CellToNeighborhood = F,NeighborhoodToCell = F)
+
+
diff --git a/vignettes/01 NICHES Spatial.Rmd b/vignettes/01 NICHES Spatial.Rmd
index 6568ce7..0ecaede 100644
--- a/vignettes/01 NICHES Spatial.Rmd	
+++ b/vignettes/01 NICHES Spatial.Rmd	
@@ -38,6 +38,7 @@ Next, we load the data, perform basic pre-processing, and cluster the data so th
 ```{r message=F, warning=F}
 InstallData("stxBrain")
 brain <- LoadData("stxBrain", type = "anterior1")
+brain <- UpdateSeuratObject(brain) # JC: need to update seurat obj
 # Normalization 
 brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)
 SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))
@@ -89,7 +90,7 @@ NICHES outputs a list of objects. Each object contains a certain style of cell-s
 
 ```{r message=F, warning=F}
 niche <- NICHES_output[['NeighborhoodToCell']]
-Idents(niche) <- niche[['ReceivingType']]
+Idents(niche) <- niche@meta.data$ReceivingType
 
 # Scale and visualize
 niche <- ScaleData(niche)
diff --git a/vignettes/02 NICHES Single.Rmd b/vignettes/02 NICHES Single.Rmd
index 63893b4..34c5159 100644
--- a/vignettes/02 NICHES Single.Rmd	
+++ b/vignettes/02 NICHES Single.Rmd	
@@ -32,7 +32,9 @@ Here, we use the 'pbmc3k' dataset from SeuratData. This is a commonly-used datas
 ```{r message=F, warning=F}
 InstallData("pbmc3k")
 data("pbmc3k")
+pbmc3k <- UpdateSeuratObject(pbmc3k) # JC: need to update pbmc3k
 Idents(pbmc3k) <- pbmc3k$seurat_annotations
+pbmc3k <- NormalizeData(pbmc3k) # JC: to be added
 pbmc3k <- ScaleData(pbmc3k)
 pbmc3k <- FindVariableFeatures(pbmc3k)
 pbmc3k <- RunPCA(pbmc3k)
@@ -58,7 +60,9 @@ scc <- RunNICHES(sub,
                  species = 'human',
                  LR.database = 'omnipath',
                  cell_types = 'seurat_annotations',
-                 CellToCell = T)
+                 CellToCell = T,
+                 SystemToCell = T
+                  )
 ```
 
 ## Visualize Cell-Cell Signaling Relationships using UMAP
@@ -86,6 +90,7 @@ scc.imputed <- RunNICHES(imputed,
                  cell_types = 'seurat_annotations',
                  CellToCell = T)
 demo.2 <- scc.imputed$CellToCell
+
 demo.2 <- ScaleData(demo.2)
 demo.2 <- RunPCA(demo.2,features = rownames(demo.2))
 ElbowPlot(demo.2,ndims=50)