常用的細(xì)胞通訊軟件：
- CellphoneDB：是公開的人工校正的敛熬，儲(chǔ)存受體、配體以及兩種相互作用的數(shù)據(jù)庫启搂。此外，還考慮了結(jié)構(gòu)組成，能夠描述異構(gòu)復(fù)合物。（配體-受體+多聚體）
- iTALK：通過平均表達(dá)量方式伐蒂，篩選高表達(dá)的胚體和受體，根據(jù)結(jié)果作圈圖肛鹏。（配體-受體）
- CellChat：CellChat將細(xì)胞的基因表達(dá)數(shù)據(jù)作為輸入逸邦，并結(jié)合配體受體及其輔助因子的相互作用來模擬細(xì)胞間通訊。（配體-受體+多聚體+輔因子）
- NicheNet // NicheNet多樣本分析 // 目標(biāo)基因的配體和靶基因活性預(yù)測(cè)：通過將相互作用細(xì)胞的表達(dá)數(shù)據(jù)與信號(hào)和基因調(diào)控網(wǎng)絡(luò)的先驗(yàn)知識(shí)相結(jié)合來預(yù)測(cè)相互作用細(xì)胞之間的配體-靶標(biāo)聯(lián)系的方法在扰。( 配體-受體+信號(hào)通路）
附：NicheNet使用的常見問題匯總其它細(xì)胞互作軟件還包括
Celltalker缕减，SingleCellSignalR，scTensor和SoptSC（這幾個(gè)也是基于配體-受體相互作用）
CellChat通過綜合信號(hào)配體芒珠、受體及其輔因子基因的表達(dá)只與它們之間互作的先驗(yàn)知識(shí)對(duì)細(xì)胞通訊概率建模桥狡。在推斷出細(xì)胞間通訊網(wǎng)絡(luò)后，CellChat提供了進(jìn)一步數(shù)據(jù)探索、分析和可視化的功能裹芝。
文章鏈接：https://www.nature.com/articles/s41467-021-21246-9
視頻鏈接：https://www.youtube.com/watch?v=kc45au1RhNs
CellChat工作流程圖：
在進(jìn)行細(xì)胞交互分析的時(shí)候部逮，不同分組的樣本盡量不要一起進(jìn)行分析，想要一起分析的時(shí)候需要保證不同分組間的細(xì)胞種類一致嫂易。同一分組的不同生物學(xué)重復(fù)可以一起分析（很多文獻(xiàn)這樣做）兄朋。
一、單樣本分析（可以是同一組的多個(gè)生物學(xué)重復(fù)一起分析）
1. 安裝：
devtools::install_github("sqjin/CellChat")
2. 數(shù)據(jù)準(zhǔn)備：
Sys.setenv(RETICULATE_PYTHON="/usr/bin/python3")
# ??要根據(jù)自己python3的路徑來設(shè)置怜械，可以在終端使用which python3來查看
library(Seurat)
library(SeuratData)
library(tidyverse)
library(CellChat)
library(NMF)
library(ggalluvial)
library(patchwork)
library(ggplot2)
library(svglite)
options(stringsAsFactors = FALSE)
pbmc3k.final <- readRDS("pbmc.rds")
pbmc3k.final@commands$FindClusters # 也可以看一看作者的其他命令颅和，Seurat是記錄其分析過程的。
#以思維導(dǎo)圖方式查看pbm3k.final的結(jié)構(gòu)
library(mindr)
(out <- capture.output(str(pbmc3k.final)))
out2 <- paste(out, collapse="\n")
mm(gsub("\\.\\.@","# ",gsub("\\.\\. ","#",out2)),type ="text",root= "Seurat")
3. 創(chuàng)建一個(gè)Cell Chat對(duì)象
從Seurat對(duì)象直接創(chuàng)建：
??：構(gòu)建Cell Chat對(duì)象時(shí)缕允，輸入的是log后的數(shù)據(jù)峡扩。
cellchat <- createCellChat(object=pbmc3k.final,group.by = "cell_type")
#cellchat <- createCellChat(pbmc3k.final@assays$RNA@data, meta = pbmc3k.final@meta.data, group.by = "cell_type")
cellchat
summary(cellchat)
str(cellchat)
levels(cellchat@idents)
#cellchat <- setIdent(cellchat, ident.use = "cell_type")
groupSize <- as.numeric(table(cellchat@idents))
#查看每個(gè)cluster有多少個(gè)細(xì)胞，后面畫圖的時(shí)候需要用到這個(gè)值
groupSize
# [1] 711 480 472 344 279 162 144 32 14
4. 導(dǎo)入配體受體數(shù)據(jù)庫
CellChatDB <- CellChatDB.human
#導(dǎo)入小鼠是CellChatDB <- CellChatDB.mouse
str(CellChatDB) #查看數(shù)據(jù)庫信息
#包含interaction灼芭、complex有额、cofactor和geneInfo這4個(gè)dataframe
colnames(CellChatDB$interaction)
CellChatDB$interaction[1:4,1:4]
head(CellChatDB$cofactor)
head(CellChatDB$complex)
head(CellChatDB$geneInfo)
#dev.new() #下一步不出圖的時(shí)候運(yùn)行
showDatabaseCategory(CellChatDB)
在CellChat中，我們還可以先擇特定的信息描述細(xì)胞間的相互作用彼绷，可以理解為從特定的側(cè)面來刻畫細(xì)胞間相互作用巍佑，比用一個(gè)大的配體庫又精細(xì)了許多。
unique(CellChatDB$interaction$annotation)#查看可以選擇的側(cè)面寄悯，也就是上圖左中的三種
#選擇"Secreted Signaling"進(jìn)行后續(xù)細(xì)胞互作分析
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling")
# use all CellChatDB for cell-cell communication analysis
# CellChatDB.use <- CellChatDB # simply use the default CellChatDB
# set the used database in the object
cellchat@DB <- CellChatDB.use # set the used database in the object
5. 預(yù)處理
對(duì)表達(dá)數(shù)據(jù)進(jìn)行預(yù)處理萤衰，用于細(xì)胞間的通信分析。首先在一個(gè)細(xì)胞組中識(shí)別過表達(dá)的配體或受體猜旬，然后將基因表達(dá)數(shù)據(jù)投射到蛋白-蛋白相互作用(PPI)網(wǎng)絡(luò)上脆栋。如果配體或受體過表達(dá)，則識(shí)別過表達(dá)配體和受體之間的相互作用洒擦。
## 在矩陣的所有的基因中提取signaling gene椿争，結(jié)果保存在data.signaling(13714個(gè)基因，過濾完只有270個(gè)）
cellchat <- subsetData(cellchat)
future::plan("multiprocess", workers = 4)
#相當(dāng)于Seurat的FindMarkers熟嫩，找每個(gè)細(xì)胞群中高表達(dá)的配體受體
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat) #Identify over-expressed ligand-receptor interactions (pairs) within the used CellChatDB
#上一步運(yùn)行的結(jié)果儲(chǔ)存在cellchat@LR$LRsig
cellchat <- projectData(cellchat, PPI.human)
#找到配體受體關(guān)系后秦踪，projectData將配體受體對(duì)的表達(dá)值投射到PPI上，來對(duì)@data.signaling中的表達(dá)值進(jìn)行校正掸茅。結(jié)果保存在@data.project
6. 推斷細(xì)胞通訊網(wǎng)絡(luò)
通過為每個(gè)相互作用分配一個(gè)概率值并進(jìn)行置換檢驗(yàn)來推斷生物意義上的細(xì)胞-細(xì)胞通信椅邓。
-
推斷配體-受體水平細(xì)胞通訊網(wǎng)絡(luò)（結(jié)果儲(chǔ)存在@net下面，有一個(gè)概率值和對(duì)應(yīng)的pval）
??這一步也是CellChat比CellPhoneDB多的一步
通過計(jì)算與每個(gè)信號(hào)通路相關(guān)的所有配體-受體相互作用的通信概率來推斷信號(hào)通路水平上的通信概率昧狮。
#根據(jù)表達(dá)值推測(cè)細(xì)胞互作的概率（cellphonedb是用平均表達(dá)值代表互作強(qiáng)度）景馁。
cellchat <- computeCommunProb(cellchat, raw.use = FALSE, population.size = TRUE) #如果不想用上一步PPI矯正的結(jié)果，raw.use = TRUE即可逗鸣。
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
cellchat <- filterCommunication(cellchat, min.cells = 10)
df.net <- subsetCommunication(cellchat)
write.csv(df.net, "net_lr.csv")
net_lr（配體-受體水平細(xì)胞通訊網(wǎng)絡(luò)）
- 推斷信號(hào)通路水平的細(xì)胞通訊網(wǎng)絡(luò)（結(jié)果儲(chǔ)存在@netP下面合住，有一個(gè)概率值和對(duì)應(yīng)的pval）
我們可以通過計(jì)算鏈路的數(shù)量或匯總通信概率來計(jì)算細(xì)胞間的聚合通信網(wǎng)絡(luò)绰精。
cellchat <- computeCommunProbPathway(cellchat)
df.netp <- subsetCommunication(cellchat, slot.name = "netP")
write.csv(df.netp, "net_pathway.csv")
net_pathway（信號(hào)通路水平的細(xì)胞通訊網(wǎng)絡(luò)）
至此，統(tǒng)計(jì)分析部分已經(jīng)完成聊疲。
7. 細(xì)胞互作關(guān)系展示
7.1 所有細(xì)胞群總體觀：細(xì)胞互作數(shù)量與強(qiáng)度統(tǒng)計(jì)分析：
#統(tǒng)計(jì)細(xì)胞和細(xì)胞之間通信的數(shù)量（有多少個(gè)配體-受體對(duì)）和強(qiáng)度（概率）
cellchat <- aggregateNet(cellchat)
#計(jì)算每種細(xì)胞各有多少個(gè)
groupSize <- as.numeric(table(cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
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")
# save as TIL/net_number_strength.pdf
左圖：外周各種顏色圓圈的大小表示細(xì)胞的數(shù)量茬底，圈越大，細(xì)胞數(shù)越多获洲。發(fā)出箭頭的細(xì)胞表達(dá)配體阱表，箭頭指向的細(xì)胞表達(dá)受體。配體-受體對(duì)越多贡珊，線越粗最爬。右圖：互作的概率/強(qiáng)度值（強(qiáng)度就是概率值相加）
檢查每種細(xì)胞發(fā)出的信號(hào)
mat <- cellchat@net$count
par(mfrow = c(3,3), xpd=TRUE)
for (i in 1:nrow(mat)) {
# i = 1
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, arrow.width = 0.2,
arrow.size = 0.1, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}
# save as TIL/net_number_individual.pdf
## 運(yùn)行上述代碼出現(xiàn)報(bào)錯(cuò) Error in plot.new() : figure margins too large
# par("mar")
## [1] 5.1 4.1 4.1 2.1
# par(mar=c(1,1,1,1))
# 重新運(yùn)行上面的代碼
每個(gè)細(xì)胞如何跟別的細(xì)胞互作（number of interaction圖）
mat <- cellchat@net$weight
par(mfrow = c(3,3), xpd=T)
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, arrow.width = 0.2,
arrow.size = 0.1, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}
# save as TIL/net_strength_individual.pdf
每個(gè)細(xì)胞如何跟別的細(xì)胞互作（互作的強(qiáng)度/概率圖）
mat <- cellchat@net$weight
par(mfrow = c(2,3), xpd=T)
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, arrow.width = 0.2,
arrow.size = 0.1, edge.weight.max = max(mat), title.name = rownames(mat)[i])
}
# save as TIL/net_strength_individual.pdf
7.2 單個(gè)信號(hào)通路或配體-受體介導(dǎo)的細(xì)胞互作可視化（層次圖舆蝴、網(wǎng)絡(luò)圖隘蝎、和弦圖、熱圖）
??注：層次圖智嚷、網(wǎng)絡(luò)圖寒随、和弦圖糠悯、熱圖只是不同的展示方法，展示的內(nèi)容和代表的意思一模一樣
比如在前面的功能富集分析或case control的比較中找到了一些信號(hào)通路差異妻往，就可以進(jìn)一步在細(xì)胞水平上驗(yàn)證互艾。
cellchat@netP$pathways #查看都有哪些信號(hào)通路
# [1] "TGFb" "NRG" "PDGF" "CCL" "CXCL" "MIF" "IL2" "IL6"
# [9] "IL10" "IL1" "CSF" "IL16" "IFN-II" "LT" "LIGHT" "FASLG"
# [17] "TRAIL" "BAFF" "CD40" "VISFATIN" "COMPLEMENT" "PARs" "FLT3" "ANNEXIN"
# [25] "GAS" "GRN" "GALECTIN" "BTLA" "BAG"
# 選擇其中一個(gè)信號(hào)通路，比如說TGFb
pathways.show <- c("TGFb")
層次圖（Hierarchy plot）
levels(cellchat@idents) # show all celltype
# [1] "Naive CD4 T" "Memory CD4 T" "CD14+ Mono" "B" "CD8 T"
# [6] "FCGR3A+ Mono" "NK" "DC" "Platelet"
vertex.receiver = c(1,2,4,6) # define a numeric vector （淋系細(xì)胞）giving the index of the celltype as targets
#par(mar=c(5.1,4.1,4.1,2.1))
netVisual_aggregate(cellchat, signaling = pathways.show, vertex.receiver = vertex.receiver)
# save as TIL/CXCL_hierarchy.pdf
在層次圖中讯泣，實(shí)體圓和空心圓分別表示源和目標(biāo)纫普。圓的大小與每個(gè)細(xì)胞組的細(xì)胞數(shù)成比例。線越粗好渠，互作信號(hào)越強(qiáng)昨稼。左圖中間的target是我們選定的靶細(xì)胞。右圖是選中的靶細(xì)胞之外的另外一組放在中間看互作拳锚。
網(wǎng)絡(luò)圖（Circle plot）
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle")
# save as TIL/CXCL_circle.pdf
和弦圖（Chord diagram）
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "chord")
# save as TIL/CXCL_chord.pdf
熱圖（Heatmap）
par(mfrow=c(1,1))
netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
# save as TIL/CXCL_heatmap.pdf
縱軸是發(fā)出信號(hào)的細(xì)胞假栓，橫軸是接收信號(hào)的細(xì)胞，熱圖顏色深淺代表信號(hào)強(qiáng)度霍掺。上側(cè)和右側(cè)的柱子是縱軸和橫軸強(qiáng)度的累積
7.3 配體-受體層級(jí)的可視化（計(jì)算各個(gè)ligand-receptor pair對(duì)信號(hào)通路的貢獻(xiàn)）
#計(jì)算配體受體對(duì)選定信號(hào)通路的貢獻(xiàn)值（在這里就是查看哪條信號(hào)通路對(duì)TGFb貢獻(xiàn)最大）
netAnalysis_contribution(cellchat, signaling = pathways.show)
pairLR.TGFb <- extractEnrichedLR(cellchat, signaling = pathways.show, geneLR.return = FALSE) #提取對(duì)TGFb有貢獻(xiàn)的所有配體受體
# save as TIL/CXCL_LR_contribution.pdf
層次圖（ Hierarchy plot）
#提取對(duì)這個(gè)通路貢獻(xiàn)最大的配體受體對(duì)來展示（也可以選擇其他的配體受體對(duì)）
LR.show <- pairLR.TGFb[1,]
vertex.receiver = c(1,2,4,6) # a numeric vector
netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, vertex.receiver = vertex.receiver)
# save as TIL/CXCL_hierarchy2.pdf
網(wǎng)絡(luò)圖（Circle plot）
netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
# save as TIL/CXCL_circle2.pdf
和弦圖（Chord diagram）
netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "chord")
# save as TIL/CXCL_chord2.pdf
Error: not enough space for cells at track index '1'.
7.4 自動(dòng)（批量）保存每個(gè)信號(hào)通路的互作結(jié)果
# Access all the signaling pathways showing significant communications將所有信號(hào)通路找出來
pathways.show.all <- cellchat@netP$pathways
# check the order of cell identity to set suitable vertex.receiver
levels(cellchat@idents)
vertex.receiver = c(1,2,4,6) #不畫層次圖就不需要這一步
dir.create("all_pathways_com_circle") #創(chuàng)建文件夾保存批量畫圖結(jié)果
setwd("all_pathways_com_circle")
for (i in 1:length(pathways.show.all)) {
# Visualize communication network associated with both signaling pathway and individual L-R pairs
netVisual(cellchat, signaling = pathways.show.all[i], out.format = c("pdf"),
vertex.receiver = vertex.receiver, layout = "circle") #繪制網(wǎng)絡(luò)圖
# Compute and visualize the contribution of each ligand-receptor pair to the overall signaling pathway
gg <- netAnalysis_contribution(cellchat, signaling = pathways.show.all[i])
ggsave(filename=paste0(pathways.show.all[i], "_L-R_contribution.pdf"),
plot=gg, width = 5, height = 2.5, units = 'in', dpi = 300)
}
setwd("../")
7.5多個(gè)配體-受體介導(dǎo)的細(xì)胞互作關(guān)系可視化
氣泡圖（全部配體受體）
levels(cellchat@idents)
# show all the significant interactions (L-R pairs)
#需要指定受體細(xì)胞和配體細(xì)胞
p = netVisual_bubble(cellchat, sources.use = c(3,5,7,8,9),
targets.use = c(1,2,4,6), remove.isolate = FALSE)
ggsave("Mye_Lymph_bubble.pdf", p, width = 8, height = 12) #髓系對(duì)淋巴的調(diào)節(jié)
# save as TIL/Mye_Lymph_bubble.pdf
展現(xiàn)細(xì)胞之間配體受體關(guān)系最直觀的圖
氣泡圖（指定信號(hào)通路或配體-受體）
#比如制定CCL和CXCL這兩個(gè)信號(hào)通路
netVisual_bubble(cellchat, sources.use = c(3,5,7,8,9), targets.use = c(1,2,4,6),
signaling = c("CCL","CXCL"), remove.isolate = FALSE)
氣泡圖（指定信號(hào)通路或配體-受體并指定細(xì)胞）
# show all the significant interactions (L-R pairs) based on user's input
pairLR.use <- extractEnrichedLR(cellchat, signaling = c("CCL","CXCL","TGFb"))
netVisual_bubble(cellchat, sources.use = c(3,6,8), targets.use = c(1,4,5),
pairLR.use = pairLR.use, remove.isolate = TRUE)
CCL匾荆，CXCL和TGFb信號(hào)通路參與單核細(xì)胞和樹突狀細(xì)胞對(duì)T細(xì)胞的調(diào)控作用情況
參與某條信號(hào)通路（如TGFb）的所有基因在細(xì)胞群中的表達(dá)情況展示（小提琴圖和氣泡圖）
## Plot the signaling gene expression distribution
p = plotGeneExpression(cellchat, signaling = "TGFb")
ggsave("TGFb_GeneExpression_vln.pdf", p, width = 8, height = 8)
p = plotGeneExpression(cellchat, signaling = "TGFb", type = "dot")
ggsave("TGFb_GeneExpression_dot.pdf", p, width = 8, height = 6)
8. 通訊網(wǎng)絡(luò)系統(tǒng)分析（可信度有待考量）
通訊網(wǎng)絡(luò)系統(tǒng)分析使用了三種算法：社會(huì)網(wǎng)絡(luò)分析、NMF分析和流行學(xué)習(xí)與分類
??：不同的算法算出來的結(jié)果可能會(huì)相互矛盾抗楔，需要結(jié)合生物學(xué)知識(shí)加以判斷
8.1 社會(huì)網(wǎng)絡(luò)分析（通訊網(wǎng)絡(luò)中的角色識(shí)別）
計(jì)算網(wǎng)絡(luò)中心性權(quán)重
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
通過計(jì)算每個(gè)細(xì)胞群的網(wǎng)絡(luò)中心性指標(biāo)棋凳，識(shí)別每類細(xì)胞在信號(hào)通路中的角色/作用C（發(fā)送者拦坠、接收者连躏、調(diào)解者和影響者）
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show,
width = 15, height = 6, font.size = 10)
# # save as TIL/SNA_CXCL_signalingRole.pdf
識(shí)別細(xì)胞的信號(hào)流模式
ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing", font.size = 5)
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming", font.size = 5)
ht1 + ht2
# save as TIL/SNA_SignalingPattern.pdf
上圖橫軸是細(xì)胞類型，縱軸是pathway贞滨。左圖是各個(gè)細(xì)胞類型中各個(gè)通路發(fā)出信號(hào)的強(qiáng)度入热，由圖是各個(gè)細(xì)胞類型中各個(gè)通路接受信號(hào)的強(qiáng)度（上面和右側(cè)的條柱和熱圖不對(duì)應(yīng)拍棕？）
8.2 非負(fù)矩陣分解（NMF）識(shí)別細(xì)胞的通訊模式（這里是一個(gè)比較標(biāo)準(zhǔn)的NMF的應(yīng)用方式）
- 信號(hào)輸出細(xì)胞的模式識(shí)別
#計(jì)算分解成幾個(gè)因子(pattern)比較合適（這一步運(yùn)行比較慢 。在使用NMF對(duì)細(xì)胞進(jìn)行亞群細(xì)分時(shí)勺良，如果不測(cè)試的話绰播，最好選擇比細(xì)胞類型多一點(diǎn)的值）
selectK(cellchat, pattern = "outgoing")
# save as TIL/NMF_outgoing_selectK.pdf
熱圖
nPatterns = 2 # 挑選曲線中第一個(gè)出現(xiàn)下降的點(diǎn)（從2就開始下降了）
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = nPatterns,
width = 5, height = 9, font.size = 6)
# save as TIL/NMF_outgoing_comPattern_heatmap.pdf
按selectK算出來的pattern值分為了2個(gè)大的pattern，左圖縱軸是細(xì)胞類型尚困，算法自動(dòng)將DC蠢箩、Memory CD4T、FCGR3A+ Mono和Platelet分為了一個(gè)pattern（和Pattern 2對(duì)應(yīng)）事甜，剩下的五種細(xì)胞分為了一個(gè)pattern（和Pattern 1對(duì)應(yīng)）谬泌。右圖縱軸是信號(hào)通路，也是被分為了兩個(gè)大的pattern逻谦，上面一個(gè)部分是在Pattern 1中比較活躍的通路掌实，下面一個(gè)部分是在Pattern 2中比較活躍的通路
沖擊圖/河流圖（river plot）
netAnalysis_river(cellchat, pattern = "outgoing")
# save as TIL/NMF_outgoing_comPattern_river.pdf
氣泡圖
netAnalysis_dot(cellchat, pattern = "outgoing")
# save as TIL/NMF_outgoing_comPattern_dotplot.pdf
- 信號(hào)輸入細(xì)胞的模式識(shí)別
selectK(cellchat, pattern = "incoming")
# save as TIL/NMF_incoming_selectK.pdf
熱圖
nPatterns = 2
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = nPatterns,
width = 5, height = 9, font.size = 6)
# save as TIL/NMF_incoming_comPattern_heatmap.pdf
沖擊圖
netAnalysis_river(cellchat, pattern = "incoming")
# save as TIL/NMF_incoming_comPattern_river.pdf
氣泡圖
netAnalysis_dot(cellchat, pattern = "incoming")
# save as TIL/NMF_incoming_comPattern_dotplot.pdf
8.3 信號(hào)網(wǎng)絡(luò)的流行學(xué)習(xí)與分類
把共同起作用的信號(hào)通路歸納在一起，分為基于功能的歸納和基于拓?fù)浣Y(jié)構(gòu)的歸納
- 基于功能相似性的流行學(xué)習(xí)與分類
cellchat <- computeNetSimilarity(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
cellchat <- netClustering(cellchat, type = "functional")
#Error in do_one(nmeth) : NA/NaN/Inf in foreign function call (arg 1)
p = netVisual_embedding(cellchat, type = "functional", label.size = 3.5)
ggsave("Manifold_functional_cluster.pdf", p, width = 8, height = 6)
#netVisual_embeddingZoomIn(cellchat, type = "functional", nCol = 2)
- 基于拓?fù)湎嗨菩缘牧餍袑W(xué)習(xí)與分類
cellchat <- computeNetSimilarity(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
cellchat <- netClustering(cellchat, type = "structural")
p = netVisual_embedding(cellchat, type = "structural", label.size = 3.5)
ggsave("Manifold_structural_cluster.pdf", p, width = 8, height = 6)
## The end
saveRDS(cellchat, file = "cellchat.rds")
二邦马、不同分組之間的配對(duì)分析
??：配對(duì)分析必須保證細(xì)胞類型是一樣的贱鼻，才可以進(jìn)行配對(duì)。如果 兩個(gè)樣本的細(xì)胞類型不一樣又想進(jìn)行配對(duì)分析時(shí)滋将，可以用subset把兩個(gè)樣本的細(xì)胞類型取成一致的邻悬。
1. 數(shù)據(jù)準(zhǔn)備，分別創(chuàng)建CellChat對(duì)象
Sys.setenv(RETICULATE_PYTHON="/usr/bin/python3")
library(Seurat)
library(tidyverse)
library(CellChat)
library(NMF)
library(ggalluvial)
library(patchwork)
rm(list = ls())
options(stringsAsFactors = FALSE)
## 創(chuàng)建cellchat對(duì)象
### 提取數(shù)據(jù)子集
scRNA <- readRDS("~/project/Integrate/scRNA.classified.rds")
scRNA$celltype <- scRNA$SingleR
table(scRNA$celltype)
Idents(scRNA) <- 'celltype'
scRNA <- subset(scRNA, idents = c('B cells','CD4+ T cells','CD8+ T cells','Dendritic cells','Monocytes','NK cells'))
scRNA$celltype <- as.factor(as.character(scRNA$celltype))
table(scRNA$orig.ident)
Idents(scRNA) <- 'orig.ident'
sco.til <- subset(scRNA, idents = c('HNC01TIL', 'HNC10TIL', 'HNC20TIL'))
sco.pbmc <- subset(scRNA, idents = c('HNC01PBMC', 'HNC10PBMC', 'HNC20PBMC'))
### 創(chuàng)建cellchat對(duì)象
cco.til <- createCellChat(sco.til@assays$RNA@data, meta = sco.til@meta.data, group.by = "celltype")
cco.pbmc <- createCellChat(sco.pbmc@assays$RNA@data, meta = sco.pbmc@meta.data, group.by = "celltype")
save(cco.til, cco.pbmc, file = "cco.rda")
2. 細(xì)胞通訊網(wǎng)絡(luò)分析
- 2.1 數(shù)據(jù)準(zhǔn)備和路徑切換
dir.create("./Compare")
setwd("./Compare")
# load("../cco.rda")
# cco.pbmc <- setIdent(cco.pbmc, ident.use = "celltype")
# cco.til <- setIdent(cco.til, ident.use = "celltype")
- 2.2 分析樣本cco.pbmc的細(xì)胞通訊網(wǎng)絡(luò)
??：cellchat不太穩(wěn)定耕渴，identifyOverExpressedGenes常出錯(cuò)（不出現(xiàn)進(jìn)度條就是出錯(cuò)了）拘悦，重啟Rstudio后再運(yùn)算。
cellchat <- cco.pbmc
cellchat@DB <- CellChatDB.human
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
#cellchat <- projectData(cellchat, PPI.human)
cellchat <- computeCommunProb(cellchat, raw.use = TRUE, population.size = TRUE)
#cellchat <- filterCommunication(cellchat, min.cells = 5)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
#cellchat <- computeNetSimilarity(cellchat, type = "functional")
#cellchat <- netEmbedding(cellchat, type = "functional")
#cellchat <- netClustering(cellchat, type = "functional")
#cellchat <- computeNetSimilarity(cellchat, type = "structural")
#cellchat <- netEmbedding(cellchat, type = "structural")
#cellchat <- netClustering(cellchat, type = "structural")
cco.pbmc <- cellchat
saveRDS(cco.pbmc, "cco.pbmc.rds")
- 2.3 分析樣本cco.til的細(xì)胞通訊網(wǎng)絡(luò)
cellchat <- cco.til
cellchat@DB <- CellChatDB.human
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
#cellchat <- projectData(cellchat, PPI.human)
cellchat <- computeCommunProb(cellchat, raw.use = TRUE, population.size = TRUE)
#cellchat <- filterCommunication(cellchat, min.cells = 5)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
#cellchat <- computeNetSimilarity(cellchat, type = "functional")
#cellchat <- netEmbedding(cellchat, type = "functional")
#cellchat <- netClustering(cellchat, type = "functional")
#cellchat <- computeNetSimilarity(cellchat, type = "structural")
#cellchat <- netEmbedding(cellchat, type = "structural")
#cellchat <- netClustering(cellchat, type = "structural")
cco.til <- cellchat
saveRDS(cco.til, "cco.til.rds")
- 2.4 合并cellchat對(duì)象
cco.list <- list(pbmc=cco.pbmc, til=cco.til)
cellchat <- mergeCellChat(cco.list, add.names = names(cco.list), cell.prefix = TRUE)
3. 可視化
3.1 所有細(xì)胞群總體觀：通訊數(shù)量與強(qiáng)度對(duì)比
gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "count")
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "weight")
p <- gg1 + gg2
ggsave("Overview_number_strength.pdf", p, width = 6, height = 4)
左圖展示通訊數(shù)量之間的差異橱脸，右圖展示通訊強(qiáng)度之間的差異础米。本例中信號(hào)通路強(qiáng)度weight值過低，導(dǎo)致顯示時(shí)均為0（實(shí)際上有數(shù)值的添诉，只是過小屁桑，顯示為0）
數(shù)量與強(qiáng)度差異網(wǎng)絡(luò)圖
par(mfrow = c(1,2))
netVisual_diffInteraction(cellchat, weight.scale = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")
# save as Diff_number_strength_net.pdf
紅色是case相對(duì)于control上調(diào)的，藍(lán)色是下調(diào)的栏赴。
數(shù)量與強(qiáng)度差異熱圖
par(mfrow = c(1,1))
h1 <- netVisual_heatmap(cellchat)
h2 <- netVisual_heatmap(cellchat, measure = "weight")
h1+h2
# save as Diff_number_strength_heatmap.pdf
case和control對(duì)比蘑斧，紅色是上調(diào)，藍(lán)色是下調(diào)须眷。
細(xì)胞互作數(shù)量對(duì)比網(wǎng)絡(luò)圖
par(mfrow = c(1,2))
weight.max <- getMaxWeight(cco.list, attribute = c("idents","count"))
for (i in 1:length(cco.list)) {
netVisual_circle(cco.list[[i]]@net$count, weight.scale = T, label.edge= F,
edge.weight.max = weight.max[2], edge.width.max = 12,
title.name = paste0("Number of interactions - ", names(cco.list)[i]))
}
# save as Counts_Compare_net.pdf
左圖是control竖瘾，右圖是case，可以直接對(duì)比數(shù)量變化花颗。
3.2 指定細(xì)胞互作數(shù)量對(duì)比網(wǎng)絡(luò)圖
par(mfrow = c(1,2))
s.cell <- c("CD4+ T cells", "CD8+ T cells", "Monocytes")
count1 <- cco.list[[1]]@net$count[s.cell, s.cell]
count2 <- cco.list[[2]]@net$count[s.cell, s.cell]
weight.max <- max(max(count1), max(count2))
netVisual_circle(count1, weight.scale = T, label.edge= T, edge.weight.max = weight.max, edge.width.max = 12,
title.name = paste0("Number of interactions-", names(cco.list)[1]))
netVisual_circle(count2, weight.scale = T, label.edge= T, edge.weight.max = weight.max, edge.width.max = 12,
title.name = paste0("Number of interactions-", names(cco.list)[2]))
# save as Counts_Compare_select.pdf 10*6.5
3.3 保守和特異性信號(hào)通路的識(shí)別與可視化
## 通路信號(hào)強(qiáng)度對(duì)比分析
gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
p <- gg1 + gg2
ggsave("Compare_pathway_strengh.pdf", p, width = 10, height = 6)
左圖最下面5個(gè)信號(hào)通路是case組獨(dú)有的
3.4 流行學(xué)習(xí)識(shí)別差異信號(hào)通路
這里function的圖出不來捕传，只有structural的圖可以出來
cellchat <- computeNetSimilarityPairwise(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
cellchat <- netClustering(cellchat, type = "functional")
#netVisual_embeddingPairwise(cellchat, type = "functional", label.size = 3.5)
#netVisual_embeddingPairwiseZoomIn(cellchat, type = "functional", nCol = 2)
cellchat <- computeNetSimilarityPairwise(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
cellchat <- netClustering(cellchat, type = "structural")
#netVisual_embeddingPairwise(cellchat, type = "structural", label.size = 3.5)
#netVisual_embeddingPairwiseZoomIn(cellchat, type = "structural", nCol = 2)
p <- rankSimilarity(cellchat, type = "structural") + ggtitle("Structural similarity of pathway")
ggsave("Pathway_Similarity.pdf", p, width = 8, height = 5)
saveRDS(cellchat, "cellchat.rds")
case和control之間信號(hào)通路差異相差程度排行，在這張圖中扩劝，ICAM相差最大庸论，其次是SELPLG职辅。??與上面那張圖的結(jié)果不相符??♀?，個(gè)人更傾向于相信上一張圖的結(jié)果
3.5 細(xì)胞信號(hào)模式對(duì)比
library(ComplexHeatmap)
總體信號(hào)模式對(duì)比
pathway.union <- union(cco.list[[1]]@netP$pathways, cco.list[[2]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(cco.list[[1]], pattern = "all", signaling = pathway.union,
title = names(cco.list)[1], width = 8, height = 10)
ht2 = netAnalysis_signalingRole_heatmap(cco.list[[2]], pattern = "all", signaling = pathway.union,
title = names(cco.list)[2], width = 8, height = 10)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))
# save as Compare_signal_pattern_all.pdf 10*6
左圖下面幾個(gè)信號(hào)通路在pbmc組中沒有聂示，和3.3中的圖相符
輸出信號(hào)模式對(duì)比
pathway.union <- union(cco.list[[1]]@netP$pathways, cco.list[[2]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(cco.list[[1]], pattern = "outgoing", signaling = pathway.union,
title = names(cco.list)[1], width = 8, height = 10)
ht2 = netAnalysis_signalingRole_heatmap(cco.list[[2]], pattern = "outgoing", signaling = pathway.union,
title = names(cco.list)[2], width = 8, height = 10)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))
# save as Compare_signal_pattern_outgoing.pdf 10*6
輸入信號(hào)模式對(duì)比
pathway.union <- union(cco.list[[1]]@netP$pathways, cco.list[[2]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(cco.list[[1]], pattern = "incoming", signaling = pathway.union,
title = names(cco.list)[1], width = 8, height = 10)
ht2 = netAnalysis_signalingRole_heatmap(cco.list[[2]], pattern = "incoming", signaling = pathway.union,
title = names(cco.list)[2], width = 8, height = 10)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))
# save as Compare_signal_pattern_incoming.pdf 10*6
3.6 特定信號(hào)通路的對(duì)比
網(wǎng)絡(luò)圖
pathways.show <- c("IL16")
weight.max <- getMaxWeight(cco.list, slot.name = c("netP"), attribute = pathways.show)
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(cco.list)) {
netVisual_aggregate(cco.list[[i]], signaling = pathways.show, layout = "circle",
edge.weight.max = weight.max[1], edge.width.max = 10,
signaling.name = paste(pathways.show, names(cco.list)[i]))
}
# save as Compare_IL16_net.pdf 10*6.5
熱圖
par(mfrow = c(1,2), xpd=TRUE)
ht <- list()
for (i in 1:length(cco.list)) {
ht[[i]] <- netVisual_heatmap(cco.list[[i]], signaling = pathways.show, color.heatmap = "Reds",
title.name = paste(pathways.show, "signaling ",names(cco.list)[i]))
}
ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))
# save as Compare_IL16_heatmap.pdf 12*6.5
和弦圖
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(cco.list)) {
netVisual_aggregate(cco.list[[i]], signaling = pathways.show, layout = "chord", pt.title = 3, title.space = 0.05,
vertex.label.cex = 0.6, signaling.name = paste(pathways.show, names(cco.list)[i]))
}
# save as Compare_IL16_chord.pdf 10*6.5
3.7 配體-受體對(duì)比分析
氣泡圖展示所有配體受體對(duì)的差異
levels(cellchat@idents$joint)
p <- netVisual_bubble(cellchat, sources.use = c(4,5), targets.use = c(1,2,3,6), comparison = c(1, 2), angle.x = 45)
ggsave("Compare_LR_bubble.pdf", p, width = 12, height = 8)
case和control配對(duì)的圖域携，文章中常見。cellphoneDB找到的結(jié)果鱼喉，也可以用這種方式呈現(xiàn)秀鞭。
氣泡圖展示上調(diào)或下調(diào)的配體受體對(duì)
p1 <- netVisual_bubble(cellchat, sources.use = c(4,5), targets.use = c(1,2,3,6), comparison = c(1, 2),
max.dataset = 2, title.name = "Increased signaling in TIL", angle.x = 45, remove.isolate = T)
p2 <- netVisual_bubble(cellchat, sources.use = c(4,5), targets.use = c(1,2,3,6), comparison = c(1, 2),
max.dataset = 1, title.name = "Decreased signaling in TIL", angle.x = 45, remove.isolate = T)
pc <- p1 + p2
ggsave("Compare_LR_regulated.pdf", pc, width = 12, height = 5.5)
上調(diào)下調(diào)分開展示
和弦圖
par(mfrow = c(1, 2), xpd=TRUE)
for (i in 1:length(cco.list)) {
netVisual_chord_gene(cco.list[[i]], sources.use = c(4,5), targets.use = c(1,2,3,6), signaling = "MHC-I",
lab.cex = 0.6, legend.pos.x = 10, legend.pos.y = 20,
title.name = paste0("Signaling from Treg - ", names(cco.list)[i]))
}
# save as Compare_LR_chord.pdf 10*6.5
小思考：但組間細(xì)胞類型占比變化似乎對(duì)CellChat組間分析結(jié)果影響較大，尤其是細(xì)胞互作強(qiáng)度扛禽。
如處理組和對(duì)照組相比气筋，巨噬細(xì)胞顯著減少，單核細(xì)胞顯著增多的情況下旋圆，蛋殼圖中巨噬細(xì)胞和巨噬細(xì)胞的互作強(qiáng)度（注意 不是數(shù)量）就會(huì)顯著下降宠默，單核和單核的互作強(qiáng)度就會(huì)顯著上升。
由于互作強(qiáng)度是根據(jù)表達(dá)的受體配體數(shù)決定的灵巧，那么細(xì)胞占比增多搀矫，強(qiáng)度就會(huì)增加，細(xì)胞占比下降刻肄，強(qiáng)度就下降瓤球？這樣似乎不是很科學(xué)
