update

NAMESPACE

@@ -41,6 +41,7 @@ import(ggExtra)
 import(ggalluvial)
 import(ggplot2)
 import(ggrepel)
+import(grDevices)
 import(iterators)
 import(methods)
 import(parallel)

R/methods.R

@@ -806,10 +806,10 @@ get_lr_path <- function(object, celltype_sender, celltype_receiver, ligand, rece
     ggi_tf <- unique(pathways[, c("src", "dest", "src_tf", "dest_tf")])
     st_data <- .get_st_data(object)
     st_meta <- .get_st_meta(object)
-    if (!celltype_sender %in% colnames(object@coef)) {
+    if (!celltype_sender %in% st_meta$celltype) {
         stop("Please provide the correct name of celltype_sender!")
     }
-    if (!celltype_receiver %in% colnames(object@coef)) {
+    if (!celltype_receiver %in% st_meta$celltype) {
         stop("Please provide the correct name of celltype_receiver!")
     }
     if (!ligand %in% rownames(st_data)) {

R/plot.R

@@ -571,13 +571,15 @@ plot_ccdist <- function(object, celltype_sender, celltype_receiver, color = NULL
 #' @param type Set 'sig' to plot significant LR pairs or set 'number' to plot the number of spatial LR interactions.
 #' @param fontsize_number fontsize of the numbers displayed in cells.
 #' @param number_color color of the text.
-#' @import ggplot2 pheatmap
+#' @param color_low For 'number' type, define the color for the lowest value.
+#' @param color_high For 'number' type, define the color for the highest value.
+#' @import ggplot2 pheatmap grDevices
 #' @importFrom ggpubr get_palette
 #' @importFrom reshape2 dcast
 #' @export
 
-plot_cci_lrpairs <- function(object, celltype_sender, celltype_receiver, top_lrpairs = 20, color = NULL,
-    border_color = "black", type = NULL, fontsize_number = 1, number_color = "black") {
+plot_cci_lrpairs <- function(object, celltype_sender, celltype_receiver, top_lrpairs = 20, color = NULL, border_color = "black", type = "sig",
+    fontsize_number = 5, number_color = "black", color_low = NULL, color_high = NULL) {
     # check
     if (!is(object, "SpaTalk")) {
         stop("Invalid class for object: must be 'SpaTalk'!")
@@ -624,25 +626,41 @@ plot_cci_lrpairs <- function(object, celltype_sender, celltype_receiver, top_lrp
     }
     ligand <- unique(lrpair$ligand)
     receptor <- unique(lrpair$receptor)
-    # get lrpair_real
-    lrpair_real <- object@lr_path$lrpairs
-    lrpair_real <- lrpair_real[, c(1, 2)]
-    lrpair_real$score <- 1
-    lrpair_mat <- reshape2::dcast(lrpair_real, formula = ligand ~ receptor, fill = 0, value.var = "score")
-    rownames(lrpair_mat) <- lrpair_mat$ligand
-    lrpair_mat <- lrpair_mat[, -1]
-    lrpair_mat <- lrpair_mat[ligand, receptor]
-    if (!is.data.frame(lrpair_mat)) {
-        stop("Limited number of ligand-receptor interactions!")
-    }
-    plot_res <- matrix("", nrow = length(ligand), ncol = length(receptor))
-    rownames(plot_res) <- ligand
-    colnames(plot_res) <- receptor
-    for (i in 1:nrow(lrpair)) {
-        plot_res[lrpair$ligand[i], lrpair$receptor[i]] <- "*"
-    }
-    pheatmap::pheatmap(lrpair_mat, cluster_cols = F, cluster_rows = F, color = heat_col, border_color = border_color,
-        legend = F, display_numbers = plot_res, fontsize_number = fontsize_number, number_color = number_color)
+    if (type == "sig") {
+        # get lrpair_real
+        lrpair_real <- object@lr_path$lrpairs
+        lrpair_real <- lrpair_real[, c(1, 2)]
+        lrpair_real$score <- 1
+        lrpair_mat <- reshape2::dcast(lrpair_real, formula = ligand ~ receptor, fill = 0, value.var = "score")
+        rownames(lrpair_mat) <- lrpair_mat$ligand
+        lrpair_mat <- lrpair_mat[, -1]
+        lrpair_mat <- lrpair_mat[ligand, receptor]
+        if (!is.data.frame(lrpair_mat)) {
+            stop("Limited number of ligand-receptor interactions!")
+        }
+        plot_res <- matrix("", nrow = length(ligand), ncol = length(receptor))
+        rownames(plot_res) <- ligand
+        colnames(plot_res) <- receptor
+        for (i in 1:nrow(lrpair)) {
+            plot_res[lrpair$ligand[i], lrpair$receptor[i]] <- "*"
+        }
+        pheatmap::pheatmap(lrpair_mat, cluster_cols = F, cluster_rows = F, color = heat_col, border_color = border_color, legend = F, display_numbers = plot_res,
+            fontsize_number = fontsize_number, number_color = number_color, main = "Significantly enriched LRI")
+    } else {
+        if (is.null(color_low)) {
+            color_low <- "orange"
+        }
+        if (is.null(color_high)) {
+            color_high <- "red"
+        }
+        lrpair_real <- lrpair[,c("ligand","receptor","lr_co_exp_num")]
+        lrpair_mat <- reshape2::dcast(lrpair_real, formula = ligand ~ receptor, fill = 0, value.var = "lr_co_exp_num")
+        rownames(lrpair_mat) <- lrpair_mat$ligand
+        lrpair_mat <- lrpair_mat[, -1]
+        heat_color <- grDevices::colorRampPalette(c(color_low, color_high))(max(as.matrix(lrpair_mat))-1)
+        heat_color <- c("white", heat_color)
+        pheatmap::pheatmap(lrpair_mat, cluster_cols = F, cluster_rows = F, border_color = border_color, color = heat_color, main = "Number of spatial LRIs")
+    }
 }
 
 #' @title Plot LR pair

script/cellchat.R

@@ -0,0 +1,37 @@
+#==============================================
+#================== CellChat ==================
+#==============================================
+library(CellChat)
+library(SpaTalk)
+library(tidyverse)
+library(Seurat)
+options(stringsAsFactors = FALSE)
+
+load(paste0(system.file(package = 'SpaTalk'),'/extdata/starmap_data.rda'))  # starmap_data
+load(paste0(system.file(package = 'SpaTalk'),'/extdata/starmap_meta.rda'))  # starmap_meta
+
+
+ndata <- starmap_data
+celltype <- starmap_meta
+rownames(celltype) <- celltype$cell
+all(colnames(ndata) == celltype$cell)
+
+seuratObj <- CreateSeuratObject(counts = ndata)
+seuratObj <- NormalizeData(seuratObj)
+Idents(seuratObj) <- celltype$celltype
+
+data.input <- GetAssayData(seuratObj, assay = "RNA", slot = "data") # normalized data matrix
+labels <- Idents(seuratObj)
+meta <- data.frame(labels = labels, row.names = names(labels)) # create a dataframe of the cell labels
+
+cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")
+CellChatDB <- CellChat::CellChatDB.mouse
+cellchat@DB <- CellChatDB
+
+cellchat <- subsetData(cellchat) 
+future::plan("multisession", workers = 4) # do parallel
+cellchat <- identifyOverExpressedGenes(cellchat)
+cellchat <- identifyOverExpressedInteractions(cellchat)
+cellchat <- computeCommunProb(cellchat)
+cellchat <- filterCommunication(cellchat, min.cells = 10)
+df.net <- res_cellchat[[use_db]][['starmap']][['starmap']] <- df.net %>% dplyr::select(1:6)

script/cellphonedb.sh

@@ -0,0 +1,109 @@
+#-----------------------------------
+#--- R script after
+#-----------------------------------
+# Run before python
+library(SpaTalk)
+library(tidyverse)
+
+load(paste0(system.file(package = 'SpaTalk'),'/extdata/starmap_data.rda'))  # starmap_data
+load(paste0(system.file(package = 'SpaTalk'),'/extdata/starmap_meta.rda'))  # starmap_meta
+
+# orthologs gene
+ndata <- starmap_data
+orthologs_gene <- mapper[rownames(ndata)]
+na_gene <- rownames(ndata)[is.na(orthologs_gene)]
+ndata <- ndata[!is.na(orthologs_gene),]
+rownames(ndata) <- orthologs_gene[!is.na(orthologs_gene)] %>% make.names()
+
+celltype <- starmap_meta
+rownames(celltype) <- celltype$cell
+all(colnames(ndata) == celltype$cell)
+
+# create Seurat object 
+
+seuratObj <- CreateSeuratObject(counts = ndata)
+Idents(seuratObj) <- celltype$celltype
+seuratObj[['celltype']] <- celltype$celltype
+seuratObj <- NormalizeData(seuratObj)
+
+writeMM(seuratObj@assays$RNA@counts, file = 'data/cpdb_input/starmap_counts_mtx/matrix.mtx')
+# save gene and cell names
+write(x = rownames(seuratObj@assays$RNA@counts), file = "data/cpdb_input/starmap_counts_mtx/features.tsv")
+write(x = colnames(seuratObj@assays$RNA@counts), file = "data/cpdb_input/starmap_counts_mtx/barcodes.tsv")
+
+[email protected]$Cell = rownames([email protected])
+df = [email protected][, c('Cell', 'celltype')]
+write.table(df, file ='data/cpdb_input/starmap_meta.tsv', sep = '\t', quote = F, row.names = F)
+
+#===============================================================================================
+
+
+#-----------------------------------
+#--- python script after
+#-----------------------------------
+
+conda activate cpdb
+cd data/cpdb_input/
+cellphonedb method statistical_analysis starmap_meta.tsv starmap_counts_mtx/ --counts-data gene_name
+
+#-----------------------------------
+#--- python script before
+#-----------------------------------
+
+
+#=================================================================================
+
+#-----------------------------------
+#--- R script after
+#-----------------------------------
+
+# standardize cpdb output
+res_cpdb_uniondb <- list()
+# create mapper to convert gene into mouse gene
+mapper <- gene2orthologs$mousegene
+names(mapper) <- gene2orthologs$humangene
+gene2orthologs <- readRDS("~/待处理/Spatalk_Benchmark/verified_database/cpdb/gene2orthologs_human_mouse.rds")
+
+cpdb2df <- function(data) {
+  df_data <- data.frame()
+  Ligand = c()
+  Receptor = c()
+  Ligand.Cluster = c()
+  Receptor.Cluster = c()
+  isReceptor_fst = c()
+  isReceptor_scn = c()
+  MeanLR = c()
+
+  for (i in 1:nrow(data)) {
+    pair <- strsplit(data$`interacting_pair`[i],split = "_")
+    for (j in 13:ncol(data)) {
+      c_pair = strsplit(colnames(data)[j],"\\|")
+      if (data[i,j] != 0.0) {
+        Ligand <- c(Ligand,pair[[1]][1])
+        Receptor <- c(Receptor,pair[[1]][2])
+        Ligand.Cluster <- c(Ligand.Cluster,c_pair[[1]][1])
+        Receptor.Cluster <- c(Receptor.Cluster,c_pair[[1]][2])
+        isReceptor_fst <- c(isReceptor_fst,data$`receptor_a`[i])
+        isReceptor_scn <- c(isReceptor_scn,data$`receptor_b`[i])
+        MeanLR <- c(MeanLR,data[i,j])
+      }
+
+    }
+
+  }
+  df_data <- data.frame(Ligand = Ligand, Receptor = Receptor,
+                        Ligand.Cluster = Ligand.Cluster,Receptor.Cluster = Receptor.Cluster,
+                        isReceptor_fst = isReceptor_fst,
+                        isReceptor_scn = isReceptor_scn,
+                        MeanLR = MeanLR)
+  return(df_data)
+}
+
+#--- starmap
+cpdb_O <- read.table(paste0("cpdb_input/","/out/significant_means.txt"),sep = "\t",header = T,check.names = F,stringsAsFactors = F)
+cpdb_O[is.na(cpdb_O)] <- 0
+cpdb_O <- cpdb2df(cpdb_O)
+
+cpdb_O <- cpdb_O %>% mutate(Ligand = mapper[Ligand], Receptor = mapper[Receptor])
+cpdb_O <- cpdb_O %>% rename(receptor = Ligand,ligand = Receptor, sender = Receptor.Cluster, receiver = Ligand.Cluster)
+