⑤ Single-cell learning-cellchat intergroup communication difference analysis

catalogs

1, data loading

2, create CellChat objects

① LS group data

② NL group data

3, build the cellchat file separately

① Construct LS group cellchat data

② Constructing NL group cellchat data

4, perform cellchat analysis

① Cellchat analysis of the LS group

② Cellchat analysis of NL groups

5, merge and group cellchat objects

①Data loading

②Merge cellchat objects

analyze

1. Overall comparison: differences in the number and intensity of communications

2. Differences in communication at the level of cell subpopulations

2.1 Network diagram of cellular communication differences

2.2 Heatmap of cellular communication differences

3. Communication differences at the level of signaling pathways

3.1 Bar graph of differences in enriched signaling pathways between groups

3.2 Horizontal heat map of outgoing signaling pathways

3.3 Horizontal heat map of afferent signaling pathways

3.4 Heat map of overall signaling pathway levels

4. Differences in horizontal communication between ligand-receptor pairs

4.1 Bubble diagram of probability differences between total paired receptor pairs

4.2 Distinguishing between upper and lower marshaled receptor pairs

5. Visualization of differences at the level of individual/specific signaling pathways

5.1 Network Diagram

5.2 Heat map

④-2 Single Cell Learning - cellchat Single Data Code Supplement (Communication Network) - CSDN Blogs

1, data loading

# Package loading




#CellChatIntergroup Difference Analysis of Cellular Communication#




rm(list=ls())




library(CellChat)




library(patchwork)




library(ggplot2)




library(Seurat)




library(ggalluvial)# Plotting the Sankey map




library(expm)




library(sna)




library(NMF)




library(ComplexHeatmap)




options(stringsAsFactors = FALSE)## input data is not automatically converted to factors (to prevent data formatting errors)

Data loading

#DataLoad ################



load("data_humanSkin.Rdata")# Data loading: here's thecount datadigital




data.input = data_humanSkin$data#Normalized gene expression matrix and cell grouping information files required



#dat <- as.data.frame(data.input)# You can view the expression matrix [1] 17328  7563




meta = data_humanSkin$meta




data.input[1:6,1:3] #Expressioncount




head(meta);table(meta$condition) #with normal(NL) and disequations(LS)



unique(meta$labels) #check cell subpopulation label type

---

2, create CellChat objects

① LS group data

# LS group data



cell.use1 = rownames(meta)[meta$condition == 'LS'] #Name of the cell from which the LS was extracted




data.input1 = data.input[, cell.use1]# Extract the LS expression matrix



meta1 = meta[cell.use1, ]#Extract LS cell information



identical(rownames(meta1),colnames(data.input1)) #Check that matrix column names and grouped file line names are identical



unique(meta1$labels) #Check cell subpopulation label types

② NL group data

# NL group data



cell.use2 = rownames(meta)[meta$condition == 'NL'] #Name of cell from which NL was extracted




data.input2 = data.input[, cell.use2]# Extract the NL expression matrix



meta2 = meta[cell.use2, ]#Extract NL cell information



identical(rownames(meta2),colnames(data.input2)) #Check that matrix column names and grouped file line names are identical



unique(meta2$labels) #Check cell subpopulation label types



#Data checking: cell line checking in different subgroups



setdiff(unique(meta1$labels),unique(meta2$labels))#character(0) Cell classifications are all consistent

3, build the cellchat file separately

① Construct LS group cellchat data

# Build cellchat files separately



## Build LS group cellchat data



cellchatLS <- createCellChat(object = data.input1, # Support for normalized expression matrices, Seurat objects, and SingleCellExperiment objects



                             meta = meta1, #meta file



                             group.by = 'labels') #cell classification columns in meta

② Constructing NL group cellchat data

## Build NL group cellchat data



cellchatNL <- createCellChat(object = data.input2, # Support for normalized expression matrices, Seurat objects, and SingleCellExperiment objects



                             meta = meta2, #meta file



                             group.by = 'labels') #cell classification columns in meta



save(cellchatLS, cellchatNL, file = "")# Staging grouped cellchat data

4, perform cellchat analysis

① Cellchat analysis of the LS group

#ReloadData for cellchat analysis #################################################



rm(list=ls())




load("")




Cellchat analysis of #LS groups



cellchat <- cellchatLS



cellchat <- setIdent(cellchat,  = 'labels') # Set the labels to the default order of display




levels(cellchat@idents) #View celltype and factor order




table(cellchat@idents) # of cells per celltype




cellchat@DB <- ## Setting up the ligand receptor database (CellChatDB).




cellchat <- subsetData(cellchat)## A subset of the expression matrix of the signaling genes Assignment to cellchat@




cellchat <- identifyOverExpressedGenes(cellchat)## Identify overexpressed signaling genes associated with each cell subpopulation




cellchat <- identifyOverExpressedInteractions(cellchat)## Recognize overexpressed gene ligand-receptor interactions




cellchat <- projectData(cellchat, )#Mapping gene expression data to PPI networks (skippable)




cellchat <- computeCommunProb(cellchat,  = TRUE) #Calculating cellular communication probabilities




cellchat <- filterCommunication(cellchat,  = 10)## Cellular communication filtering




cellchat <- computeCommunProbPathway(cellchat)# Calculate the probability of communication at the level of the signaling pathway




cellchat <- aggregateNet(cellchat)# Calculate the number and probability strength of inter-cell pair communication




cellchat <- netAnalysis_computeCentrality(cellchat,#Computing network centrality weights: identifying roles/activities of each cell type in signaling pathways




                                           = "netP")




cellchatLS <- cellchat



saveRDS(cellchatLS, "")

② Cellchat analysis of the NL group

Cellchat analysis of #cellchatNL groups



rm(list=ls())



load("")



cellchat <- cellchatNL



cellchat <- setIdent(cellchat,  = 'labels') # Set the labels to the default order of display




levels(cellchat@idents) #View celltype and factor order




table(cellchat@idents) # of cells per celltype




cellchat@DB <- ## Setting up the ligand receptor database (CellChatDB).




cellchat <- subsetData(cellchat)## A subset of the expression matrix of the signaling genes Assigned to cellchat@




cellchat <- identifyOverExpressedGenes(cellchat)## Identify overexpressed signaling genes associated with each cell subpopulation




cellchat <- identifyOverExpressedInteractions(cellchat)## Recognize overexpressed gene ligand-receptor interactions




cellchat <- projectData(cellchat, )#Mapping gene expression data to PPI networks (skippable)




cellchat <- computeCommunProb(cellchat,  = TRUE) #Calculating cellular communication probabilities




cellchat <- filterCommunication(cellchat,  = 10)## Cellular communication filtering




cellchat <- computeCommunProbPathway(cellchat)# Calculate the probability of communication at the level of the signaling pathway




cellchat <- aggregateNet(cellchat)# Calculate the number and probability strength of inter-cell pair communication




cellchat <- netAnalysis_computeCentrality(cellchat,#Computing network centrality weights: identifying roles/activities of each cell type in signaling pathways




                                           = "netP")



cellchatNL <- cellchat



saveRDS(cellchatNL, "")




#######################################################

---

5, merge and group cellchat objects

①Data loading

# Two datasets of skin samples from patients with atopic dermatitis: non-lesional (NL, normal) and lesional (LS, diseased)




rm(list=ls())



cellchatNL <- readRDS("")



cellchatLS <- readRDS("")




levels(cellchatNL@idents)




levels(cellchatLS@idents)




identical(levels(cellchatNL@idents),levels(cellchatLS@idents))#Data checking for consistency [1] TRUE

②Merge cellchat objects

#Merge cellchat objects.




object.list <- list(NL = cellchatNL,



                    LS = cellchatLS) # control group (NL) first, comparison group (LS) second, note order



cellchat <- mergeCellChat(object.list, add.names = names(object.list))



cellchat



dplyr::glimpse(cellchat)#data structure view



## Merged content.'','images','net','netP','meta','idents','','DB',and 'LR'

---

The next step was to start the between-group cellchat analysis

---

1. Overall comparison: differences in number and intensity of communications

#1. Overall comparison: differences in the number and intensity of communications



p1 <- compareInteractions(cellchat,  = F, group = c(1,2))



p2 <- compareInteractions(cellchat,  = F, group = c(1,2), 



                          measure = "weight")



p1 + p2

The picture on the left shows the twodata setbetween the number of communications, and the figure on the right illustrates the difference between the strength of communications.

2. Differences in communication at the level of cellular subpopulations

2.1 Network diagram of cellular communication differences

#2. Differences in communication at the cellular subpopulation level



par(mfrow = c(1,2), xpd = TRUE)



netVisual_diffInteraction(cellchat,  = T)



netVisual_diffInteraction(cellchat,  = T, measure = "weight")

Shows the difference in the number of ligand receptor pairs (left) and the difference in the communication probability (right) in all cell subpopulations between samples; different colors of the outer peripheral circle represent different cell subpopulations, and the size indicates the number of ligand receptor pairs of the subpopulation, and the bigger the circle is, the bigger the ratio of the number of ligand receptor pairs between cells is. The blue line indicates stronger communication in the control group, and the red line indicates stronger communication in the comparison group. The thicker the line, the stronger the degree of variation in communication.

2.2 Heatmap of cellular communication differences

#cellular communication differences heatmap



p3 <- netVisual_heatmap(cellchat)



p4 <- netVisual_heatmap(cellchat, measure = "weight")



p3 + p4




dev.off()

The vertical coordinate is the ligand cell, the horizontal coordinate is the receptor cell, and the heat map color represents the communication probability difference; the darker the color, the stronger the communication. The bars on the top and right are the accumulation of communication probability differences in the vertical and horizontal axes. The left panel shows the difference in the number of ligand-receptor pairs between cells, and the right panel shows the difference in communication probability.

3. Communication differences at the level of signaling pathways

3.1 Bar graph of differences in enriched signaling pathways between groups

#3. Communication differences at the level of signaling pathways



## Sorting signaling pathways based on information flow or number of interactions



p5 <- rankNet(cellchat, mode = "comparison", stacked = T,  = TRUE) #Stacking



p6 <- rankNet(cellchat, mode = "comparison", stacked = F,  = TRUE) #No Stacking



p5 + p6

When the ratio of the sum of the probabilities of the pathways in the comparison (LS) group to the control (NL) group is <0.95 and the difference pval <0.05 (rank sum test), the communication strength of the pathway is significantly increased in the control group (vertical coordinates in red); when the ratio of the sum of the probabilities of the pathways in the comparison (LS) group to the control (NL) group is >1.05 and the difference pval <0.05, the communication strength of the pathway in the comparison group is significantly increased (vertical coordinates in blue); vertical coordinates in black indicate that the pathway did not differ between the two groups. Scale plots are shown on the left, and actual value comparison plots are shown on the right.

3.2 Horizontal heat map of outgoing signaling pathways

# Horizontal heat map of efferent signaling pathways



i = 1




 <- union(object.list[[i]]@netP$pathways,



                       object.list[[i+1]]@netP$pathways)







ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]],



                                        pattern = "outgoing", # Outgoing



                                        signaling = ,



                                        title = names(object.list)[i],



                                        width = 5,



                                        height = 6)



ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]],



                                        pattern = "outgoing", # Outgoing



                                        signaling = ,



                                        title = names(object.list)[i+1],



                                        width = 5,



                                        height = 6)



draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

3.3 Horizontal heat map of afferent signaling pathways

#Heatmap of afferent signaling pathway levels



ht3 = netAnalysis_signalingRole_heatmap(object.list[[i]],



                                        pattern = "incoming", # incoming



                                        signaling = ,



                                        title = names(object.list)[i],



                                        width = 5, height = 6,



                                         = "GnBu")



ht4 = netAnalysis_signalingRole_heatmap(object.list[[i+1]],



                                        pattern = "incoming", # incoming



                                        signaling = ,



                                        title = names(object.list)[i+1],



                                        width = 5, height = 6,



                                         = "GnBu")



draw(ht3 + ht4, ht_gap = unit(0.5, "cm"))

3.4 Heat map of overall signaling pathway levels

# Overall signaling pathway level heat map



ht5 = netAnalysis_signalingRole_heatmap(object.list[[i]],



                                        pattern = "all", # Overall



                                        signaling = ,



                                        title = names(object.list)[i],



                                        width = 5, height = 6,



                                         = "OrRd")



ht6 = netAnalysis_signalingRole_heatmap(object.list[[i+1]],



                                        pattern = "all", # Overall



                                        signaling = ,



                                        title = names(object.list)[i+1],



                                        width = 5, height = 6,



                                         = "OrRd")



draw(ht5 + ht6, ht_gap = unit(0.5, "cm"))

The vertical coordinate is the signaling pathway, the horizontal coordinate is the cell subpopulation, and the heatmap color represents the signal intensity; the darker the color, the stronger the communication. The bars on the top and right are the accumulation of the intensity of the vertical and horizontal axes. NL group on the left and LS group on the right.

4. Differences in horizontal communication between ligand-receptor pairs

4.1 Bubble diagram of probability difference between total paired receptor pairs

levels(cellchat@idents$joint) #View cell subpopulations



netVisual_bubble(cellchat,



                 sources.use = 4,



                 targets.use = c(5:12),



                 comparison = c(1, 2),



                  = 45)

4.2 Distinguishing between upper and lower deployed receptor pairs

p7 <- netVisual_bubble(cellchat,



                       sources.use = 4,



                       targets.use = c(5:12),



                       comparison = c(1, 2),



                        = 2,



                        = "Increased signaling in LS",



                        = 45,



                        = T) #Increased is the ligand-receptor pair information with stronger probability of communication for the comparator group



p8 <- netVisual_bubble(cellchat,



                       sources.use = 4,



                       targets.use = c(5:12),



                       comparison = c(1, 2),



                        = 1,



                        = "Decreased signaling in LS",



                        = 45,



                        = T) #Decreased for ligand-receptor pair information with stronger probability of communication for controls



p7 + p8

In the bubble diagram, the horizontal coordinates are cell pairs with color distinguishing samples; the vertical coordinates are ligand receptors. The size of the bubble indicates the p-value; the smaller the p-value the larger the bubble. The color indicates the magnitude of the communication probability.

5. Visualization of differences at the level of individual/specific signaling pathways

5.1network diagram

#Use network diagrams:



 <- c("CXCL") #Select target signaling pathways



 <- getMaxWeight(object.list,



                            = c("netP"),



                           attribute = ) #Control edge weights for different datasets



 



par(mfrow = c(1,2), xpd = TRUE)




for (i in 1:length(object.list)) {



  netVisual_aggregate(object.list[[i]],



                      signaling = ,



                      layout = "circle",



                       = [1],



                       = 10,



                       = paste(, names(object.list)[i]))



}

5.2 Heat map

 <- c("CXCL")



par(mfrow = c(1,2), xpd = TRUE)



ht <- list()




for (i in 1:length(object.list)) {



  ht[[i]] <- netVisual_heatmap(object.list[[i]],



                               signaling = ,



                                = "Reds",



                                = paste(, "signaling ",names(object.list)[i]))



}



ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))

---

Reference:

1：Single Cell Analysis of Cellular Interactions-3: CellChat - simple book ()

2：CellChat Cell Communication Analysis (Next) - Hands-on Code - Comparative Analysis of Multiple Data Sets - Knowledge ()

3：CellChat Intergroup Difference Analysis ()

4：focuslyj/CellChat ()