空转也可以使用cellchat进行细胞通讯，详细的代码在cellchat官网也有哒
题外话：当然，也有单细胞的细胞cellchat细胞通讯：
制备单张空转片子 并整理空转内部数据
print(getwd())#request 2.libPaths(c( "/home/data/t040413/R/x86_64-pc-linux-gnu-library/4.2","/home/data/t040413/R/yll/usr/local/lib/R/site-library","/home/data/refdir/Rlib/", "/usr/local/lib/R/library"))# SiO2_7, NS_7, SiO2_56, NS_56setwd('/home/data/t040413/silicosis/spatial/')### #############################sio2_7# sio2_7_exp <- Read10X_h5(filename = "./sio7d/2.3.h5_files/filtered_feature_bc_matrix.h5" ) #1815# sio2_7_blank <- as.matrix(read.csv("./blank/sio7_B26.csv",sep=",",header=T));sio2_7_blank #72# sio2_7_valid=setdiff(colnames(sio2_7_exp),sio2_7_blank[,1])# sio2_7_exp_valid=sio2_7_exp[,sio2_7_valid] #1743## sio2_7=CreateSeuratObject(counts = sio2_7_exp_valid, project = 'SiO2_7', assay = 'Spatial',min.cells=3) #18690*1743# sio2_7$stim <- "sio2_7"## sio2_7_img <- Read10X_Image(image.dir="./sio7d/spatial/")# DefaultAssay(object = sio2_7_img) <- 'Spatial'# sio2_7_img <- sio2_7_img[colnames(sio2_7)]# sio2_7[['image']] <- sio2_7_img## SpatialFeaturePlot(sio2_7, features = "nCount_Spatial")# SpatialFeaturePlot(sio2_7, features = c("Grem1")) #画的count值# ############单样本SCT# sio2_7_sct=SCTransform(sio2_7, verbose = FALSE,,assay ="Spatial")# SpatialFeaturePlot(sio2_7, features = c("Grem1")) #画的count值# sio7=sio2_7_sct### ########## #############################sio2_56# sio2_56_exp <- Read10X_h5(filename = "./sio56d/2.3.h5_files/filtered_feature_bc_matrix.h5" )# sio2_56_blank <- as.matrix(read.csv("./blank/sio56_B11.csv",sep=",",header=T));sio2_56_blank #135# sio2_56_valid=setdiff(colnames(sio2_56_exp),sio2_56_blank[,1])# sio2_56_exp_valid=sio2_56_exp[,sio2_56_valid]## sio2_56=CreateSeuratObject(counts = sio2_56_exp_valid, project = 'sio2_56', assay = 'Spatial',min.cells=3)# sio2_56$stim <- "sio2_56"## sio2_56_img <- Read10X_Image(image.dir="./sio56d/spatial/")# DefaultAssay(object = sio2_56_img) <- 'Spatial'# sio2_56_img <- sio2_56_img[colnames(sio2_56)]# sio2_56[['image']] <- sio2_56_img## SpatialFeaturePlot(sio2_56, features = "nCount_Spatial")# SpatialFeaturePlot(sio2_56, features = c("Grem1")) #画的count值# ############单样本SCT# sio2_56_sct=SCTransform(sio2_56, verbose = FALSE,,assay ="Spatial")# SpatialFeaturePlot(sio2_56, features = c("Grem1")) #画的count值# sio56=sio2_56_sct#### ##########3# #############################ns_7# ns_7_exp <- Read10X_h5(filename = "./ns7/2.3.h5_files/filtered_feature_bc_matrix.h5" )# ns_7_blank <- as.matrix(read.csv("./blank/ns7d_A72_.csv",sep=",",header=T));ns_7_blank# ns_7_valid=setdiff(colnames(ns_7_exp),ns_7_blank[,1])# ns_7_exp_valid=ns_7_exp[,ns_7_valid] #1743## ns_7=CreateSeuratObject(counts = ns_7_exp_valid, project = 'ns_7', assay = 'Spatial',min.cells=3) #18690*1743# ns_7$stim <- "ns_7"## ns_7_img <- Read10X_Image(image.dir="./ns7/spatial/")# DefaultAssay(object = ns_7_img) <- 'Spatial'# ns_7_img <- ns_7_img[colnames(ns_7)]# ns_7[['image']] <- ns_7_img## SpatialFeaturePlot(ns_7, features = "nCount_Spatial")# SpatialFeaturePlot(ns_7, features = c("Grem1")) #画的count值# ############单样本SCT# ns_7_sct=SCTransform(ns_7, verbose = FALSE,,assay ="Spatial")# SpatialFeaturePlot(ns_7, features = c("Grem1")) #画的count值# ns7=ns_7_sct### ###### #############################ns_56# ns_56_exp <- Read10X_h5(filename = "./ns56/2.3.h5_files/filtered_feature_bc_matrix.h5" )# ns_56_blank <- as.matrix(read.csv("./blank/ns56d_A73_.csv",sep=",",header=T));ns_56_blank# ns_56_valid=setdiff(colnames(ns_56_exp),ns_56_blank[,1])# ns_56_exp_valid=ns_56_exp[,ns_56_valid] #1743## ns_56=CreateSeuratObject(counts = ns_56_exp_valid, project = 'ns_56', assay = 'Spatial',min.cells=3) #18690*1743# ns_56$stim <- "ns_56"## ns_56_img <- Read10X_Image(image.dir="./ns56/spatial/")# DefaultAssay(object = ns_56_img) <- 'Spatial'# ns_56_img <- ns_56_img[colnames(ns_56)]# ns_56[['image']] <- ns_56_img## SpatialFeaturePlot(ns_56, features = "nCount_Spatial")# SpatialFeaturePlot(ns_56, features = c("Grem1")) #画的count值# ############单样本SCT# ns_56_sct=SCTransform(ns_56, verbose = FALSE,,assay ="Spatial")# SpatialFeaturePlot(ns_56, features = c("Grem1")) #画的count值# ns56=ns_56_sct#library(Seurat)library(dplyr)#d.all=readRDS("~/silicosis/spatial/integrated_slides/integrated_slides.rds")load("/home/data/t040413/silicosis/spatial_transcriptomics/silicosis_ST_harmony_SCT_r0.5.rds")dim(d.all)DefaultAssay(d.all)="Spatial"#visium_slides=SplitObject(object = d.all,split.by = "stim")names(d.all);dim(d.all)d.all@meta.data %>%head()head(colnames(d.all))#给d.all 添加meta信息------adata_obs=read.csv("~/silicosis/spatial/adata_obs.csv")head(adata_obs)mymeta= paste0(d.all@meta.data$orig.ident,"_",colnames(d.all)) %>% gsub("-.*","",.) # %>% head()head(mymeta)tail(mymeta)#掉-及其之后内容adata_obs$col= adata_obs$spot_id %>% gsub("-.*","",.) # %>% head()head(adata_obs)rownames(adata_obs)=adata_obs$coladata_obs=adata_obs[mymeta,]head(adata_obs)identical(mymeta,adata_obs$col)d.all=AddMetaData(d.all,metadata = adata_obs)table(d.all$origional_clusters)d.all$labels=d.all$origional_clusters##从d.all 提取单个样本-----Idents(d.all)=d.all$labelshead(d.all@meta.data)#str(ns7)##自建函数----get_one_slideimage=function(slide,image_name="NS_7"){# slide=eachslidetmpimages=slide@images[[image_name]]slide@images=list()slide@images[[image_name]]=tmpimagesreturn(slide)}ns7= d.all[,d.all@meta.data$orig.ident=='NS_7']ns7=get_one_slideimage(slide = ns7,image_name = "NS_7")ns7=SCTransform(ns7,assay = 'Spatial')ns7@imagessio7=d.all[,d.all@meta.data$orig.ident=='SiO2_7']sio7=get_one_slideimage(slide = sio7,image_name = 'SiO2_7')sio7=SCTransform(sio7,assay = "Spatial")visium_slides=SplitObject(object = d.all,split.by = "stim")visium_slides[[1]]@imagesmyslides=list()for (eachslide in visium_slides ) {#eachslide= visium_slides[[1]]eachslide=get_one_slideimage(slide = eachslide,image_name = unique(eachslide$stim) )eachslide=SCTransform(eachslide,assay = "Spatial")myslides[[unique(eachslide$stim)]]=eachslide}myslides[[1]]@meta.data %>%head()sio7@imagesseu1=ns7seu2=sio7ptm = Sys.time()#library(CellChat)library(patchwork)# show the image and annotated spots# color.use <- scPalette(nlevels(seu1)); names(color.use) <- levels(seu1)# p1 <- Seurat::SpatialDimPlot(seu1, label = F, label.size = 3, cols = color.use)# color.use <- scPalette(nlevels(seu2)); names(color.use) <- levels(seu2)# p2 <- Seurat::SpatialDimPlot(seu2, label = F, label.size = 3, cols = color.use) #+ NoLegend()# p1 + p2
2制作单张空转片子 提取cellchat所需数据
5#开始制备cellchat空转所需数据----5.1##5.1 获取单个片子的表达矩阵---# Prepare input data for CelChat analysisdata.input1 = Seurat::GetAssayData(seu1, slot = "data", assay = "SCT") # normalized data matrixhead(data.input1)[,1:9]data.input2 = Seurat::GetAssayData(seu2, slot = "data", assay = "SCT")#-my_data_inputs=list()for (slide in myslides) {my_data_inputs[[unique(slide$stim)]] = Seurat::GetAssayData(slide, slot = "data", assay = "SCT")}my_data_inputs5.2##5.2获取交集基因-----genes.common <- intersect(rownames(data.input1), rownames(data.input2))length(genes.common)#--genes.common= Reduce(intersect, lapply(myslides, rownames))5.3## 确保不同片子的列名不一样-----colnames(data.input1) <- paste0("NS_7_", colnames(data.input1))colnames(data.input2) <- paste0("SiO2_7_", colnames(data.input2))names(my_data_inputs)class(my_data_inputs)for (eachdata in names(my_data_inputs)) {print(names(eachdata))colnames( my_data_inputs[[eachdata]] ) =paste0(paste0(eachdata,"_"), colnames(my_data_inputs[[eachdata]]) )}my_data_inputs[[1]] %>%head() %>%.[,1:7]# #掉-及其之后内容# adata_obs$col= adata_obs$spot_id %>% gsub("-.*","",.) # %>% head()# head(adata_obs)# #掉-及其之后内容gsub_final_dash=function(string=c("SiO2_7_AAACTTGCAAACGTAT-1_1",'fad-fe')) {string= gsub(pattern = "-.*",replacement = '',x = string)return(string)}gsub_final_dash()for (eachdata in names(my_data_inputs)) {print(names(eachdata))colnames( my_data_inputs[[eachdata]] ) = gsub_final_dash(string = colnames(my_data_inputs[[eachdata]]) )}my_data_inputs[[1]] %>%head() %>%.[,1:7]colnames(my_data_inputs) %>%head()5.4#合并多个空转片子，获取表达矩阵---data.input <- cbind(my_data_inputs[[1]][genes.common, ],my_data_inputs[[2]][genes.common, ],my_data_inputs[[3]][genes.common, ],my_data_inputs[[4]][genes.common, ])data.input %>%head() %>%.[,1:7]stringr::str_split(string = colnames(data.input),pattern = "_[A-Z]",simplify = TRUE)[,1] %>% table()
添加空转meta信息，添加cell2location的细胞丰度信息
6 #添加metadata信息-------# define the meta data# a column named `slices` should be provided for spatial transcriptomics analysis, which is useful for analyzing cell-cell communication by aggregating multiple slices/replicates. Of note, for comparison analysis across different conditions, users still need to create a CellChat object seperately for each condition.meta1 = data.frame(labels = Idents(seu1), slices = "NS_7") # manually create a dataframe consisting of the cell labelsmeta2 = data.frame(labels = Idents(seu2), slices = "SiO2_7")myslides[[1]]@meta.data %>%head()##2 添加细胞丰度信息-------#major_abundance1=read.csv("~/silicosis/spatial_c2l_results/sc_silicosis_20231006/cell2location_map/q05_cell_abundance_w_sf.csv")major_abundance1=read.csv("~/silicosis/spatial_c2l_results/silicosis_mapped_after_kaiti/sc_silicosis_analysis/cell2location_map/q05_cell_abundance_w_sf.csv")head(major_abundance1)selected_abundance=major_abundance1[,c("spot_id",'AM1','AM3','Inmt.fibroblast','Grem1.fibroblast')]rownames(selected_abundance)=gsub(pattern = "-[1-9]",replacement = "",selected_abundance$spot_id)head(selected_abundance)selected_abundance[rowSums(selected_abundance[,-1])==0,]dim(selected_abundance)selected_abundance=selected_abundance[d.all@meta.data$col,]dim(selected_abundance)identical(rownames(selected_abundance),d.all@meta.data$col)# load("~/silicosis/spatial_c2l_results/silicosis_mapped_after_kaiti/barcodes_cellltypes_proportion.Rdata")# head(new_df)# top1_df=new_df %>%group_by(barcodes) %>% slice_max(proportion) %>%as.data.frame()# top1_df %>%head()# table(top1_df$`cell type`)## head(wide_df)# wide_df_new=wide_df %>%as.data.frame()# rownames(wide_df_new)=wide_df_new$barcodes# identical(wide_df_new$barcodes,d.all@meta.data$col)head(selected_abundance)str(selected_abundance)celltype <- apply(selected_abundance[,-1], 1, function(x) colnames(selected_abundance[,-1]) [which.max(x)] );print(length(celltype))selected_abundance$celltype <- celltypehead(selected_abundance)table(selected_abundance$celltype)which.max(c(1,54,3,532,1))meta1 %>%head()# meta <- rbind(meta1, meta2)# meta %>%head()meta= selected_abundance # d.all@meta.data # %>%head()d.all$stim=factor(d.all$stim,levels = unique(d.all$stim))table(d.all$stim)#tail(data.input[,(ncol(data.input)-10):ncol(data.input)])#调整行顺序，使两个数据框行名一致---identical(rownames(meta),colnames(data.input))rownames(meta) <- colnames(data.input)# a factor level should be defined for the `meta$labels` and `meta$slices`#meta$labels <- factor(meta$labels, levels = levels(Idents(d.all) ) )meta$labels=factor(meta$celltype,levels = unique(meta$celltype));head(rownames(meta))#meta$stim= sapply( stringr::str_split(rownames(meta), "_", simplify = TRUE) , `[`, 1:2)meta$stim= stringr:: str_extract(rownames(meta), "^([^_]+_[^_]+)")meta$slices <- factor(meta$stim, levels = unique(d.all$stim) )unique(meta$labels) # check the cell labelshead(meta)
空转空间信息来自于cellranger的下降数据
7#加载空间位置信息----# Spatial locations of spots from full (NOT high/low) resolution images are required. For 10X Visium, this information is in `tissue_positions.csv`.spatial.locs1 = Seurat::GetTissueCoordinates(myslides$SiO2_7, scale = NULL, cols = c("imagerow", "imagecol"))spatial.locs2 = Seurat::GetTissueCoordinates(myslides$NS_7, scale = NULL, cols = c("imagerow", "imagecol"))spatial.locs3 = Seurat::GetTissueCoordinates(myslides$SiO2_56, scale = NULL, cols = c("imagerow", "imagecol"))spatial.locs4 = Seurat::GetTissueCoordinates(myslides$NS_56, scale = NULL, cols = c("imagerow", "imagecol"))spatial.locs <- rbind(spatial.locs1, spatial.locs2,spatial.locs3, spatial.locs4)head(spatial.locs)head(data.input)rownames(spatial.locs) <- colnames(data.input)ns7@images$NS_7@scale.factors8#加载图片缩放因子信息# Scale factors of spatial coordinates# For 10X Visium, the conversion factor of converting spatial coordinates from Pixels to Micrometers can be computed as the ratio of the theoretical spot size (i.e., 65um) over the number of pixels that span the diameter of a theoretical spot size in the full-resolution image (i.e., 'spot_diameter_fullres' in pixels in the 'scalefactors_json.json' file).# Of note, the 'spot_diameter_fullres' factor is different from the `spot` in Seurat object and thus users still need to get the value from the original json file.scalefactors1 = jsonlite::fromJSON(txt = file.path("~/silicosis/spatial/sio7d/spatial/scalefactors_json.json"))scalefactors1myslides$SiO2_7@images$SiO2_7@scale.factorsspot.size = 65 # the theoretical spot size (um) in 10X Visiumconversion.factor1 = spot.size/scalefactors1$spot_diameter_fullresconversion.factor1spatial.factors1 = data.frame(ratio = conversion.factor1, tol = spot.size/2)spatial.factors1######ns7scalefactors2 = jsonlite::fromJSON(txt = file.path("~/silicosis/spatial/ns7/spatial/scalefactors_json.json"))scalefactors2myslides$SiO2_7@images$SiO2_7@scale.factorsspot.size = 65 # the theoretical spot size (um) in 10X Visiumconversion.factor2 = spot.size/scalefactors2$spot_diameter_fullresconversion.factor2spatial.factors2 = data.frame(ratio = conversion.factor2, tol = spot.size/2)spatial.factors2######sio56scalefactors3 = jsonlite::fromJSON(txt = file.path("~/silicosis/spatial/sio56d/spatial/scalefactors_json.json"))scalefactors3myslides$SiO2_7@images$SiO2_7@scale.factorsspot.size = 65 # the theoretical spot size (um) in 10X Visiumconversion.factor3= spot.size/scalefactors3$spot_diameter_fullresconversion.factor3spatial.factors3 = data.frame(ratio = conversion.factor3, tol = spot.size/2)spatial.factors3######ns56scalefactors4 = jsonlite::fromJSON(txt = file.path("~/silicosis/spatial/ns56/spatial/scalefactors_json.json"))scalefactors4myslides$SiO2_7@images$SiO2_7@scale.factorsspot.size = 65 # the theoretical spot size (um) in 10X Visiumconversion.factor4= spot.size/scalefactors4$spot_diameter_fullresconversion.factor4spatial.factors4 = data.frame(ratio = conversion.factor4, tol = spot.size/2)spatial.factors4# scalefactors2 = jsonlite::fromJSON(txt = file.path("/Users/suoqinjin/Library/CloudStorage/OneDrive-Personal/works/CellChat/tutorial/spatial_imaging_data-intestinalA2", 'scalefactors_json.json'))# conversion.factor2 = spot.size/scalefactors2$spot_diameter_fullres# spatial.factors2 = data.frame(ratio = conversion.factor2, tol = spot.size/2)spatial.factors <- rbind(spatial.factors1, spatial.factors2,spatial.factors3, spatial.factors4)rownames(spatial.factors) <- c("SiO2_7", "NS_7","SiO2_56", "NS_56")spatial.factors9#创建cellchat空转对象head(data.input)[,1:6]head(meta)head(spatial.locs)spatial.factors#BiocManager::install("CellChat")getwd().libPaths()dir.create("~/R/mycellchat2").libPaths(c("~/R/mycellchat2",.libPaths()))#devtools::install_github("jinworks/CellChat")library(CellChat)#meta$labels=paste0("Cluster",meta$labels)identical(colnames(data.input),rownames(meta))head(meta);dim(meta)table(meta$labels)cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels",datatype = "spatial", coordinates = spatial.locs, spatial.factors = spatial.factors)#>1010#设置配体受体数据库-----CellChatDB <- CellChatDB.mouse # use CellChatDB.human if running on human datahead(CellChatDB$complex)lapply(CellChatDB,head)head(CellChatDB$interaction)CellChatDB$geneInfotable(CellChatDB$interaction$annotation)# use a subset of CellChatDB for cell-cell communication analysisCellChatDB.use <- CellChatDBno #subsetDB(CellChatDB, search = "Secreted Signaling", key = "annotation") # use Secreted SignalingCellChatDB.use# set the used database in the objectcellchat@DB <- CellChatDB.use11# subset the expression data of signaling genes for saving computation costcellchat <- subsetData(cellchat) # This step is necessary even if using the whole databasefuture::plan("multisession", workers = 20)cellchat <- identifyOverExpressedGenes(cellchat)cellchat <- identifyOverExpressedInteractions(cellchat)execution.time = Sys.time() - ptmprint(as.numeric(execution.time, units = "secs"))print(getwd())12#推断互作网络ptm = Sys.time()trims=0.8cellchat <- computeCommunProb(cellchat, type = "truncatedMean", trim = trims,distance.use = FALSE, interaction.range = 250, scale.distance = NULL,contact.dependent = TRUE, contact.range = 100)#> truncatedMean is used for calculating the average gene expression per cell group.#> [1] ">>> Run CellChat on spatial transcriptomics data without distance values as constraints of the computed communication probability <<< [2023-12-11 09:02:32.712269]"#> [1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2023-12-11 09:05:16.73344]"# Filter out the cell-cell communication if there are only few number of cells in certain cell groupscellchat <- filterCommunication(cellchat, min.cells = 10)#> The cell-cell communication related with the following cell groups are excluded due to the few number of cells: Stromal 4cellchat <- computeCommunProbPathway(cellchat)cellchat <- aggregateNet(cellchat)execution.time = Sys.time() - ptmgsub(pattern = ":",replacement = '_',x = Sys.time())print(as.numeric(execution.time, units = "secs"))print(getwd())mypath=paste0("~/silicosis/spatial/cellchat2_spatial/",gsub(pattern = ":",replacement = '_',x = Sys.time()),'figures')dir.create(mypath)setwd(mypath)print(getwd())save(cellchat,file = paste0(mypath,"/",paste0("trims_",trims),'cellchat.rds') )print('here-====================================================')#load("~/silicosis/spatial/cellchat2_spatial/2024-01-09 18_51_33.849761figures/trims_0.8cellchat.rds")#setwd("~/silicosis/spatial/cellchat2_spatial/2024-01-09 18_51_33.849761figures/")netVisual_heatmap(cellchat, measure = "count", color.heatmap = "Blues")#netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle")# DefaultAssay(d.all)="SCT"# d.all@meta.data %>%head()# Idents(d.all)# brain=d.all# brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)# brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30)# brain <- FindClusters(brain, verbose = FALSE)# brain <- RunUMAP(brain, reduction = "pca", dims = 1:30)# SpatialDimPlot(brain, label = TRUE, label.size = 3)getwd()
画图
# setwd("./cellchat2_spatial/2024-01-09 18_51_33.849761figures/")#输出单细胞类图表------library(CellChat){print(getwd())df.net <- subsetCommunication(cellchat)head(df.net)openxlsx:: write.xlsx(df.net,file='0.Cell-Cell_Communications_At_L-R.xlsx', rowNames=F, colNames=T)df.net <- subsetCommunication(cellchat, slot.name = "netP")openxlsx:: write.xlsx(df.net,file='0.Cell-Cell_Communications_At_Pathway.xlsx', rowNames=F, colNames=T)groupSize <- as.numeric(table(cellchat@idents))## NumberOfInteractionsmat <- cellchat@net$countopenxlsx:: write.xlsx(mat, file='1.NumberOfInteractions.xlsx', rowNames=T, colNames=T)pdf("1.NumberOfInteractions.pdf")netVisual_circle(mat, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")dev.off()pdf("1.NumberOfInteractions_Split.pdf")for (i in 1:nrow(mat)) {mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))mat2[i, ] <- mat[i, ]p = netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat), title.name = rownames(mat)[i])print(p)}dev.off()## InteractionWeightsmat <- cellchat@net$weightopenxlsx:: write.xlsx(mat, file='2.InteractionWeights.xlsx', rowNames=T, colNames=T)pdf("2.InteractionWeights.pdf")netVisual_circle(mat, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")dev.off()pdf("2.InteractionWeights_Split.pdf")for (i in 1:nrow(mat)) {mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))mat2[i, ] <- mat[i, ]p = netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat), title.name = rownames(mat)[i])print(p)}dev.off()## cellchat@netP$pathways: the signaling pathways showing significant communicationspathways = cellchat@netP$pathways## the left portion shows autocrine and paracrine signaling to certain cell groups of interest (i.e, the defined vertex.receiver)## the right portion shows autocrine and paracrine signaling to the remaining cell groups in the datasetvertex.receiver = seq(1,2)pdf("3.Sig_Pathway_Hierarchy_Plot.pdf")for(i in pathways){print(i)p = netVisual_aggregate(cellchat, signaling = i,vertex.receiver = vertex.receiver,layout='hierarchy' #, vertex.label.cex = 0.4)title(main = paste0(i,' signaling'))print(p)}dev.off()pdf("3.Sig_Pathway_Circle_Plot.pdf")for(i in pathways){print(i)p = netVisual_aggregate(cellchat, signaling = i,vertex.receiver = vertex.receiver,layout = "circle")title(main = paste0(i,' signaling'))print(p)}dev.off()pdf("4.Sig_Pathway_L-R_pair_Contribution.pdf")for(i in pathways){print(i)p = netAnalysis_contribution(cellchat, signaling = i, title = paste0(i, " signaling pathway", " Contribution of each L-R pair"))print(p)}dev.off()pdf("4.Sig_Pathway_L-R_pair_bubbleplot.pdf", width=25, height=20)p = netVisual_bubble(cellchat, remove.isolate = FALSE)print(p)dev.off()cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")pdf("5.Signaling_Roles_Of_Cell_Groups_Heatmap.pdf")for(i in pathways){print(i)p = netAnalysis_signalingRole_network(cellchat, signaling = i, width = 8, height = 2.5, font.size = 10)print(p)}dev.off()pdf("4.Sig_Pathway_L-R_pair_bubbleplot.pdf", width=25, height=20)p = netVisual_bubble(cellchat, remove.isolate = FALSE)print(p)dev.off()pdf("5.Signaling_Roles_Of_Cell_Groups_2D.pdf")p = netAnalysis_signalingRole_scatter(cellchat)print(p)dev.off()pdf("5.signals_Contribution_Of_Cell_Groups_Heatmap.pdf", width=10)ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing", font.size = 5)ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming", font.size = 5)print(ht1 + ht2)dev.off()save(cellchat, file = paste0("cellchat_","_.RData"))pdf('6_heatmap_counts_weight.pdf')p1=netVisual_heatmap(cellchat, measure = "count", color.heatmap = 'Blues')p2=netVisual_heatmap(cellchat, measure = "weight", color.heatmap = "Blues")print(p1)print(p2)dev.off()head(cellchat@netP)pdf('7_chord_signalling_cellgroups.pdf')for (i in cellchat@netP$pathways) {p= netVisual_chord_cell(cellchat,signaling = i)print(p)}dev.off()pdf('8_chord_LR_cellgroups.pdf')for (i in cellchat@netP$pathways) {p= netVisual_chord_gene(cellchat,sources.use = 1,targets.use = (3:4),slot.name = "net" )print(p)}for (i in cellchat@netP$pathways) {p= netVisual_chord_gene(cellchat,sources.use = 1,targets.use = (3:4),slot.name = "netP" )print(p)}dev.off()library(NMF)library(ggalluvial)library(CellChat)selectK(cellchat,pattern = "outgoing").libPaths()npatterns=3pdf('9_patterns_outcoming3.pdf' )par(mfrow=c(1,2))cellchat=identifyCommunicationPatterns(cellchat,pattern = "outgoing",k=npatterns,heatmap.show = TRUE,height = 10)dev.off()pdf("9_river_outgoning_pattern3.pdf",width = 9)p=netAnalysis_river(cellchat,pattern = "outgoing")print(p)dev.off()pdf("9_dotplot_outgoing_pattern3.pdf" )p=netAnalysis_dot(cellchat,pattern = 'outgoing')print(p)dev.off()getwd()#selectK(cellchat,pattern = "incoming")pdf('9_patterns_incoming 3.pdf' )par(mfrow=c(1,2))# cellchat=identifyCommunicationPatterns(cellchat,pattern = "incoming",k=npatterns,height = 9)cellchat=identifyCommunicationPatterns(cellchat,pattern = "incoming",k=3,height = 15)dev.off()pdf('9_patterns_incoming 3.pdf' )par(mfrow=c(1,2))# cellchat=identifyCommunicationPatterns(cellchat,pattern = "incoming",k=npatterns,height = 9)cellchat=identifyCommunicationPatterns(cellchat,pattern = "incoming",k=3,height = 20)dev.off()pdf("9_river_incoming_pattern3.pdf",width = 9)p=netAnalysis_river(cellchat,pattern = "incoming")print(p)dev.off()pdf("9_dotplot_incoming_pattern3.pdf" )p=netAnalysis_dot(cellchat,pattern = 'incoming')print(p)dev.off()npatterns=2pdf('patterns_outcoming2.pdf' )par(mfrow=c(1,2))cellchat=identifyCommunicationPatterns(cellchat,pattern = "outgoing",k=npatterns,heatmap.show = TRUE,height = 9)dev.off()pdf("9_river_outgoning_pattern2.pdf",width = 9)p=netAnalysis_river(cellchat,pattern = "outgoing")print(p)dev.off()pdf("9_dotplot_outgoing_pattern2.pdf" )p=netAnalysis_dot(cellchat,pattern = 'outgoing')print(p)dev.off()getwd()#selectK(cellchat,pattern = "incoming")pdf('patterns_incoming2.pdf' )par(mfrow=c(1,2))cellchat=identifyCommunicationPatterns(cellchat,pattern = "incoming",k=npatterns,height = 9)cellchat=identifyCommunicationPatterns(cellchat,pattern = "incoming",k=3,height = 9)dev.off()pdf("9_river_incoming_pattern2.pdf",width = 9)p=netAnalysis_river(cellchat,pattern = "incoming")print(p)dev.off()pdf("9_dotplot_incoming_pattern2.pdf" )p=netAnalysis_dot(cellchat,pattern = 'incoming')print(p)dev.off()}13#可视化---ptm = Sys.time()groupSize <- as.numeric(table(cellchat@idents))par(mfrow = c(1,2), xpd=TRUE)netVisual_circle(cellchat@net$count,vertex.weight = rowSums(cellchat@net$count),weight.scale = T, label.edge= F,title.name = "Number of interactions")netVisual_circle(cellchat@net$weight,vertex.weight = rowSums(cellchat@net$weight),weight.scale = T, label.edge= F,title.name = "Interaction weights/strength")color.heatmap <- colorRampPalette(c("lightpink", "red", "darkred"))(9)netVisual_heatmap(cellchat, measure = "count", color.heatmap = c('white','blue'))netVisual_heatmap(cellchat, measure = "count", color.heatmap = 'Blues')netVisual_heatmap(cellchat, measure = "weight", color.heatmap = "Blues")cellchat@netP$pathwayspathways.show <- c("SPP1")# Circle plotpar(mfrow=c(1,1))
空间cellchat细胞通讯可视化
14#空间可视化# Spatial plotlibrary(ggplot2) # Assuming ggplot2 is used for plotting# Define the slice names to iterate overslice.names <- c("SiO2_7", "SiO2_56", "NS_56", "NS_7")# Initiate a PDF file for outputpdf("cellchat_aggregate_spatial_plots.pdf")# Loop through each slice name and generate the plotfor (slice.name in slice.names) {p <- netVisual_aggregate(cellchat, signaling = pathways.show, slice.use = slice.name,layout = "spatial", edge.width.max = 2, vertex.size.max = 1,alpha.image = 0.2, vertex.label.cex = 0)print(p) # Print the plot to the PDF}# Close the PDF filedev.off()# Setting `vertex.label.cex = 0` to hide the labels on the spatial plotexecution.time = Sys.time() - ptmprint(as.numeric(execution.time, units = "secs"))#> [1] 1.171192# Compute the network centrality scorescellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP") # the slot 'netP' means the inferred intercellular communication network of signaling pathways# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groupspar(mfrow=c(1,1))netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)par(mfrow=c(1,1))netVisual_aggregate(cellchat, signaling = pathways.show,slice.use = c("SiO2_56"), layout = "spatial", edge.width.max = 2, alpha.image = 0.2,vertex.weight = "incoming",vertex.size.max = 6, vertex.label.cex = 0)netAnalysis_contribution(cellchat, signaling = pathways.show)# Take an input of a few genesspatialFeaturePlot(cellchat, features = c("Tgfb1"), slice.use = "SiO2_7",point.size = 0.8, color.heatmap = "Reds", direction = 1)lrs=cellchat@LR$LRsig# Take an input of a ligand-receptor pairspatialFeaturePlot(cellchat, pairLR.use = "TGFB1_TGFBR1_TGFBR2",slice.use = "SiO2_7", point.size = 0.5, do.binary = FALSE,cutoff = 0.05, enriched.only = F, color.heatmap = c('blue','white','red'), direction = 1)# Take an input of a ligand-receptor pair and show expression in binaryspatialFeaturePlot(cellchat, pairLR.use = "BMP2_BMPR1A_ACVR2A",slice.use = "SiO2_7", point.size = 1.5,do.binary = TRUE, cutoff = 0.05, enriched.only = F,color.heatmap = "Reds", direction = 1)# Take an input of a ligand-receptor pair and show expression in binaryspatialFeaturePlot(cellchat, pairLR.use = "BMP2_BMPR1A_ACVR2A", slice.use = "SiO2_7",point.size = 1.5, do.binary = TRUE, cutoff = 0.05,enriched.only = F, color.heatmap = "Reds", direction = 1)
