cellchat

### perform cellchat
p_load(Seurat, tidyverse, CellChat, qs)
options(future.globals.maxSize = 1 * 1024^4)
sce <- qread('sce.qs', nthreads = 40)
data.input <- GetAssayData(sce, layer = "data")
identity <- data.frame(group = sce$cellType, row.names = names(sce$cellType))
unique(identity$group)
cellchat <- createCellChat(data.input)
cellchat <- addMeta(cellchat, meta = identity, meta.name = "labels")
cellchat <- setIdent(cellchat, ident.use = "labels")
# cellchat@DB <- CellChatDB.mouse #.mouse or .human
cellchat@DB <- subsetDB(CellChatDB.mouse, search = 'Secreted Signaling')
future::plan("multisession", workers = 25)
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
cellchat <- projectData(cellchat, PPI.mouse) #.mouse or .human
cellchat <- computeCommunProb(cellchat, type = "triMean")
cellchat <- filterCommunication(cellchat, min.cells = 5)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
qsave(cellchat, 'cellchat.qs', nthreads = 40)

### visualization
p_load(tidyverse, CellChat)
cellchat <- readRDS("cellchat.rds")

groupSize <- as.numeric(table(cellchat@idents))
netVisual_circle(cellchat@net$count, vertex.weight = groupSize, weight.scale = T, 
                 label.edge= F, title.name = "Number of interactions")
netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, 
                 label.edge= F, title.name = "Interaction weights/strength")

mat <- cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat)) {
  mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
  mat2[i, ] <- mat[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, 
                   edge.weight.max = max(mat), title.name = rownames(mat)[i])
}

ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing", 
                                         height = 36)
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming", 
                                         height = 36, color.heatmap = 'Reds')
pdf('cellchat_LRpair.pdf', height=16, width=12)
ht1 + ht2
dev.off()

# show all the significant interactions (L-R pairs) from some cell groups 
# (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(cellchat, sources.use = c(6), targets.use = c(1:13))

# show all the interactions sending from
netVisual_chord_gene(cellchat, sources.use = c(6), targets.use = c(1:3))
# show all the interactions received by 
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(8), lab.cex = 0.5)

netVisual_chord_gene(cellchat, sources.use = c(1:12), targets.use = c(13), 
                     signaling = c('MHC-I'),legend.pos.x = 8)

gg1 <- netAnalysis_signalingRole_scatter(cellchat)
gg2 <- netAnalysis_signalingRole_scatter(cellchat, signaling = c("CXCL", "CCL"))
gg1 + gg2

pathways.show <- c("MHC-I")
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "chord")
netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
netAnalysis_contribution(cellchat, signaling = pathways.show)


plotGeneExpression(cellchat, signaling = pathways.show, enriched.only = FALSE)

p_load(NMF, ggalluvial)
selectK(cellchat, pattern = "outgoing")
nPatterns = 6
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", 
                                          k = nPatterns)
netAnalysis_river(cellchat, pattern = "outgoing")
netAnalysis_dot(cellchat, pattern = "outgoing")

# Identify signaling groups based on their functional similarity
cellchat <- computeNetSimilarity(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
# Manifold learning of the signaling networks for a single dataset
cellchat <- netClustering(cellchat, type = "functional")
# Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(cellchat, type = "functional", label.size = 3.5)

# Identify signaling groups based on structure similarity
cellchat <- computeNetSimilarity(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
# Manifold learning of the signaling networks for a single dataset
cellchat <- netClustering(cellchat, type = "structural")
# Classification learning of the signaling networks for a single dataset
# Visualization in 2D-space
netVisual_embedding(cellchat, type = "structural", label.size = 3.5)

### for multiple cellchat object
p_load(tidyverse, CellChat, qs)
options(future.globals.maxSize = 1 * 1024^4)
ctrl <- qread('./ctrl.qs', nthreads = 40)
treat <- qread('./pa.qs', nthreads = 40)
object_list <- list(ctrl = ctrl, treat = treat)
cellchat <- mergeCellChat(object_list, add.names = names(object_list), 
                          merge.data = T, cell.prefix = T)

# 1.overview_interactions
gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2))
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), 
                           measure = "weight")
gg1 + gg2

# 2.network_interactions
par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")

# 3.interaction_strength
num.link <- sapply(object_list, 
                   function(x) {rowSums(x@net$count) + 
                       colSums(x@net$count)-diag(x@net$count)})
weight.MinMax <- c(min(num.link), max(num.link)) # control the dot size in the different datasets
gg <- list()
for (i in 1:length(object_list)) {
  gg[[i]] <- netAnalysis_signalingRole_scatter(object_list[[i]], 
                                               title = names(object_list)[i], 
                                               weight.MinMax = weight.MinMax)
}
patchwork::wrap_plots(plots = gg)

# 4.pathway_difference
type = 'functional' 
### type = 'structural'
cellchat <- computeNetSimilarityPairwise(cellchat, type = type)
cellchat <- netEmbedding(cellchat, type = type)
cellchat <- netClustering(cellchat, type = type)
# netVisual_embeddingPairwise(cellchat, type = type, label.size = 3.5)
gg1 <- rankSimilarity(cellchat, type = type)
# gg1 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg1 + gg2

# 5.LR_pair for certain celltype
netVisual_bubble(cellchat, sources.use = c(1:13), 
                 targets.use = c(12),  comparison = c(1, 2), angle.x = 45)

# 6.detailed_chordgram
pathways.show <- c("TGFb") 
par(mfrow = c(1,2), xpd=TRUE)
ht <- list()
for (i in 1:length(object_list)) {
  ht[[i]] <- netVisual_heatmap(object_list[[i]], signaling = pathways.show, 
                               color.heatmap = "Reds",
                               title.name = paste(pathways.show,"_signaling_",names(object_list)[i]))
}
ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))

# 7.detailed_expression_in_certain_pathway
cellchat@meta$datasets = factor(cellchat@meta$datasets, 
                                levels = c("ctrl", "treat")) # set factor level
plotGeneExpression(cellchat, signaling = "CXCL", 
                   split.by = "datasets", colors.ggplot = T)

# 8.Compare outgoing signal patterns
pathway.union <- union(object_list[[1]]@netP$pathways, object_list[[2]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(object_list[[1]], pattern = "outgoing", 
                                        signaling = pathway.union,title = names(object_list)[1],
                                        width=5,height=24,font.size=6, color.heatmap = "Blues")
ht2 = netAnalysis_signalingRole_heatmap(object_list[[1]], pattern = "incoming", 
                                        signaling = pathway.union,title = names(object_list)[1],
                                        width=5,height=24,font.size=6, color.heatmap = "Greens")
ht3 = netAnalysis_signalingRole_heatmap(object_list[[2]], pattern = "outgoing", 
                                        signaling = pathway.union,title = names(object_list)[2],
                                        width=5,height=24,font.size=6, color.heatmap = "Blues")
ht4 = netAnalysis_signalingRole_heatmap(object_list[[2]], pattern = "incoming", 
                                        signaling = pathway.union,title = names(object_list)[2],
                                        width=5,height=24,font.size=6, color.heatmap = "Greens")
pdf('8.compare_signal_patterns', height=12, width=12)
ComplexHeatmap::draw(ht1 + ht3 + ht2 + ht4, ht_gap = unit(0.5, "cm"))
dev.off()

# 9. Comparing communications
gg1 <- netVisual_bubble(cellchat, sources.use = c(1:13), targets.use = c(10),  comparison = c(1, 2), max.dataset = 1, 
                        title.name = paste0("Increased_signaling_in_",names(cellchat@net)[1]), angle.x = 45, remove.isolate = F)
gg2 <- netVisual_bubble(cellchat, sources.use = c(1:13), targets.use = c(10),  comparison = c(1, 2), max.dataset = 2, 
                        title.name = paste0("Increased_signaling_in_",names(cellchat@net)[2]), angle.x = 45, remove.isolate = T)
gg1 + gg2