cellchat-(1)单数据集的细胞通讯分析

library(tidyverse)
library(Seurat)
library(CellChat)
pbmc3k <- read_rds('data/pbmc3k_anno.rds')

# 加个label
pbmc3k$label_cc <- Idents(pbmc3k)

# 创建cellchat对象，注意需要用Normalized data, 这里用cellchat提供的`normalizeData`函数
cellchat <- createCellChat(normalizeData(pbmc3k[['RNA']]@counts), 
                           meta = pbmc3k[[]],
                           group.by = 'label_cc')

showDatabaseCategory(CellChatDB.human) #替换后缀的物种名即可查看其他物种的分类情况

glimpse(CellChatDB.human$interaction)

#以下是终端输出内容，为了方便浏览放到了代码框中
Rows: 1,939
Columns: 11
$ interaction_name   <chr> "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2", "TGFB3_TGFBR1_TGFBR2", "…
$ pathway_name       <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb"…
$ ligand             <chr> "TGFB1", "TGFB2", "TGFB3", "TGFB1", "TGFB1", "TGFB2", "TGFB2", "TGFB3"…
$ receptor           <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGFbR2", "ACVR1C_TGFbR2…
$ agonist            <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb …
$ antagonist         <chr> "TGFb antagonist", "TGFb antagonist", "TGFb antagonist", "TGFb antagon…
$ co_A_receptor      <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""…
$ co_I_receptor      <chr> "TGFb inhibition receptor", "TGFb inhibition receptor", "TGFb inhibiti…
$ evidence           <chr> "KEGG: hsa04350", "KEGG: hsa04350", "KEGG: hsa04350", "PMID: 27449815"…
$ annotation         <chr> "Secreted Signaling", "Secreted Signaling", "Secreted Signaling", "Sec…
$ interaction_name_2 <chr> "TGFB1 - (TGFBR1+TGFBR2)", "TGFB2 - (TGFBR1+TGFBR2)", "TGFB3 - (TGFBR1…

cellchat@DB <- CellChatDB.human
#也可以只选择其中某种类型的相互作用,如：
# cellchat@DB <- subsetDB(CellChatDB.human, search = 'Secreted Signaling')

# 首先取一下要计算的基因子集，因为CellChatDB.human中包含4万多基因，
# 而我们的pbmc3k中只有1万多基因，这样可以节省接下来的计算时间
cellchat <- subsetData(cellchat)
future::plan("multisession", workers = 4) #设置使用的电脑核心数，并行计算节省时间
# 注意，windows系统需将'multisession'改为'multiprocess'

cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

# trim参数用来过滤掉只在少量细胞表达的基因，
# 0.1表示如果某基因在一种细胞类型中表达比例少于10%，
# 则该基因在这种细胞中平均表达量归0，也就是该配体/受体不在这种细胞中表达
# population.size设置为TRUE表示需要考虑每种细胞类型中细胞数量的多少，对于分选的细胞需要设置为FALSE
cellchat <- computeCommunProb(cellchat, trim = 0.1, population.size = T) 
# min.cells用来过滤掉细胞数目较少的细胞类型
cellchat <- filterCommunication(cellchat, min.cells = 10)

# 在信号通路水平推断相互作用
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

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"
)

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])
}