前言
在前面的文章中，我们已经分别介绍了单样本 CellChat 分析以及多样本 CellChat 的比较分析思路。然而，在真正进入多样本比较时，研究者往往会遇到一个不可回避的现实问题：CellChat 在进行样本间对比时，默认要求各个样本具有完全一致的细胞类型组成。这一前提在理想化的数据中或许成立，但在真实的生物学研究中却几乎难以满足。
以肿瘤研究为例，肿瘤组织中通常富含多种免疫细胞，如 CD4⁺ T 细胞、NK 细胞、浆细胞等，而在对应的正常组织中，这些细胞类型可能完全不存在或仅占极低比例。细胞组成的这种显著差异，直接导致两个样本在 CellChat 框架下无法进行直接的细胞通讯比较——并不是因为通讯机制本身不可比，而是因为参与通讯的“细胞角色”本身就不一致。
如果在这种情况下强行进行比较，就如同拿两张结构完全不同的地图进行对照：一张多了道路和城市，另一张却缺少关键节点，最终得到的差异很可能并非真实的生物学变化，而只是由结构不一致所引入的分析偏差。正是为了解决这一核心问题，CellChat 提供了liftCellChat方法，通过先统一不同数据集的细胞类型框架，使所有样本在同一套“细胞语言体系”下进行分析，从而保证后续观察到的通讯差异，真正反映的是生物学状态的变化，而不是细胞组成差异带来的技术性假象。
一.数据准备
基于CellChat v2.20
依旧是我们的七例甲状腺肿瘤
library(Seurat) library(CellChat) # 读取之前的数据 scRNA <- readRDS("D:/datat/新建文件夹/GSE191288_RAW/scRNA.rds") T_sce <- readRDS("D:/datat/新建文件夹/GSE191288_RAW/T_sce.rds") B_sce <- readRDS("D:/datat/新建文件夹/GSE191288_RAW/B_sce.rds") # 子集注释回写母对象 phe <- scRNA@meta.data phe_T <- T_sce@meta.data # 如有因子，先转换为字符 phe$celltype <- as.character(phe$celltype) phe_T$celltype.1 <- as.character(phe_T$celltype.1) # 找到匹配的行 common_T <- intersect(rownames(phe), rownames(phe_T)) # 回写细胞类型 phe[common_T, "celltype"] <- phe_T[common_T, "celltype.1"] scRNA@meta.data <- phe table(phe$orig.ident,phe$celltype) # B Cells CD 4 CD 8 Endothelial Epithelial Cells Fibroblasts Cells Myeloid Cells NK NKT γδ T # NT 12 55 79 1122 909 878 51 58 1 28 # T1L 60 139 111 663 135 1390 67 25 4 30 # T1R 300 637 363 222 66 425 192 22 16 82 # T2L 48 834 570 191 164 406 150 261 20 191 # T2R 559 1005 398 191 31 166 114 100 11 67 # T3L 11 260 160 60 1506 336 269 40 13 29 # T3R 5 175 332 159 3113 188 164 62 6 54 ## NT为癌旁组织，其余为甲状腺组织（临床分型为PTC） ## 本次流程使用癌旁组织（NT）和 甲状腺组织（T1L）这次我们挑选NT（癌旁）和T2L（甲状腺癌）做对比，NT样本的B细胞（12个），NKT（1个），我们可以将NT样本中的B Cells和NKT删除，这样NT样本中就只有CD 4 ；CD 8； Endothelial ；Epithelial Cells ；Fibroblasts Cells； Myeloid Cells ； NK ； γδ T，8种细胞类型，而T2L则相比较于NT多了两种细胞类型（B Cells ，NKT）
实际操作（提取矩阵和meta）
## 拆分样本 obj.list <- SplitObject(scRNA, split.by = "orig.ident") ## 删除NT样本中的B Cells 和NKT obj.list[["NT"]] <- subset( obj.list[["NT"]], subset = !(celltype %in% c("NKT", "B Cells")) ) ## 查看两个样本的细胞类型和数量 table(obj.list[["NT"]]$celltype) # # CD 4 CD 8 Endothelial Epithelial Cells Fibroblasts Cells Myeloid Cells NK # 55 79 1122 909 878 51 58 # γδ T # 28 table(obj.list[["T2L"]]$celltype) # # B Cells CD 4 CD 8 Endothelial Epithelial Cells Fibroblasts Cells Myeloid Cells # 48 834 570 191 164 406 150 # NK NKT γδ T # 261 20 191 ## 现在NT样本中有8种细胞类型，而T2L种则有10种细胞类型（注意！！ T2L中的10中细胞类型必须包含 NT样本中的8种细胞类型） # 根据研究情况进行细胞排序（T2L） celltype_order <- c( "CD 4", "CD 8", "γδ T", "NK", "NKT", "B Cells", "Myeloid Cells", "Fibroblasts Cells", "Epithelial Cells", "Endothelial" ) ## 细胞的基因表达数据，由于每个样本只删除几个细胞，所以不需要再做标准化，直接提取就可以 data.input_T2L <- GetAssayData(obj.list[["T2L"]], slot = 'data') # normalized data matrix meta_T2L <- obj.list[["T2L"]]@meta.data[,c("orig.ident","celltype")] # 对 meta 和 data.input 进行排序 identical(rownames(meta_T2L), colnames(data.input_T2L)) meta_T2L$celltype <- factor(meta_T2L$celltype ,levels = celltype_order) ordered_indices <- order(meta_T2L$celltype) meta_T2L <- meta_T2L[ordered_indices, ] data.input_T2L <- data.input_T2L[, ordered_indices] identical(rownames(meta_T2L),colnames(data.input_T2L)) # 根据研究情况进行细胞排序（NT） celltype_order <- c( "CD 4", "CD 8", "γδ T", "NK", "Myeloid Cells", "Fibroblasts Cells", "Epithelial Cells", "Endothelial" ) data.input_NT <- GetAssayData(obj.list[["NT"]], slot = 'data') # normalized data matrix ## metadata文件 meta_NT <- obj.list[["NT"]]@meta.data[,c("orig.ident","celltype")] # 对 meta 和 data.input 进行排序 identical(rownames(meta_NT), colnames(data.input_NT)) meta_NT$celltype <- factor(meta_NT$celltype ,levels = celltype_order) ordered_indices <- order(meta_NT$celltype) meta_NT <- meta_NT[ordered_indices, ] data.input_NT <- data.input_NT[, ordered_indices] identical(rownames(meta_NT),colnames(data.input_NT))创建cellchat
# 加载数据库 CellChatDB <- CellChatDB.human #人源样本 # 查看数据的分类情况与 # 具体内容showDatabaseCategory(CellChatDB) dplyr::glimpse(CellChatDB$interaction) #筛选数据库，只保留“Secreted Signaling” CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling") ## NT样品 cellchat_NT <- createCellChat(object = data.input_NT, meta = meta_NT, group.by = "celltype") # 添加细胞信息 cellchat_NT <- addMeta(cellchat_NT, meta = meta_NT) cellchat_NT <- setIdent(cellchat_NT, ident.use = "celltype") levels(cellchat_NT@idents) groupSize_NT <- as.numeric(table(cellchat_NT@idents)) cellchat_NT@DB <- CellChatDB.use cellchat_NT <- subsetData(cellchat_NT) # 只保留 配体/受体基因的表达矩阵 # 类似于找高变基因 cellchat_NT <- identifyOverExpressedGenes(cellchat_NT) cellchat_NT <- identifyOverExpressedInteractions(cellchat_NT) cellchat_NT <- computeCommunProb(cellchat_NT, type = "triMean") cellchat_NT <- filterCommunication(cellchat_NT, min.cells = 10) cellchat_NT <- computeCommunProbPathway(cellchat_NT) cellchat_NT <- aggregateNet(cellchat_NT) cellchat_NT <- netAnalysis_computeCentrality(cellchat_NT, slot.name = "netP") ## T2L样品 cellchat_T2L <- createCellChat(object = data.input_T2L, meta = meta_T2L, group.by = "celltype") # 添加细胞信息 cellchat_T2L <- addMeta(cellchat_T2L, meta = meta_T2L) cellchat_T2L <- setIdent(cellchat_T2L, ident.use = "celltype") levels(cellchat_T2L@idents) groupSize_T2L <- as.numeric(table(cellchat_T2L@idents)) cellchat_T2L@DB <- CellChatDB.use cellchat_T2L <- subsetData(cellchat_T2L) # 只保留 配体/受体基因的表达矩阵 # 类似于找高变基因 cellchat_T2L <- identifyOverExpressedGenes(cellchat_T2L) cellchat_T2L <- identifyOverExpressedInteractions(cellchat_T2L) cellchat_T2L <- computeCommunProb(cellchat_T2L, type = "triMean") cellchat_T2L <- filterCommunication(cellchat_T2L, min.cells = 10) cellchat_T2L <- computeCommunProbPathway(cellchat_T2L) cellchat_T2L <- aggregateNet(cellchat_T2L) cellchat_T2L <-netAnalysis_computeCentrality(cellchat_T2L, slot.name = "netP")合并样本重点！！
## 更新为当前版本兼容的对象结构,如果是自己从头跑且期间没有更新过cellchat包则不用跑 ## 如果是加载比别人的数据或者中途更新了包或者R版本则需要updateCellChat()更新数据类型 cellchat_NT <- updateCellChat(cellchat_NT) cellchat_T2L <- updateCellChat(cellchat_T2L) # 添加细胞类型前 table(cellchat_NT@idents) CD 4 CD 8 γδ T NK Myeloid Cells Fibroblasts Cells Epithelial Cells 55 79 28 58 51 878 909 Endothelial 1122 ## 使用liftCellChat()给NT样本添加和N2L一样细胞类型 group.new = levels(cellchat_T2L@idents) cellchat_NT <- liftCellChat(cellchat_NT, group.new) # 查看添加后（细胞数量都是0） table(cellchat_NT@idents) CD 4 CD 8 γδ T NK NKT B Cells Myeloid Cells 55 79 28 58 0 0 51 Fibroblasts Cells Epithelial Cells Endothelial 878 909 1122 ## 然后合并样本 object.list <- list(NT = cellchat_NT, T2L = cellchat_T2L) cellchat <- mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE)可视化
Circle plot
## 查看两个样本有什么通路 lapply(object.list, function(x) x@netP$pathways) $NT [1] "MIF" "CypA" "VEGF" "ANNEXIN" "MK" "CXCL" "GALECTIN" "PTN" "TNF" "SPP1" "VISFATIN" [12] "CCL" "ANGPTL" "PARs" "PDGF" "GRN" "OSM" "IFN-II" "EGF" "PLAU" "TGFb" "GAS" [23] "ANGPT" "SEMA3" "CX3C" "EDN" "IL16" "TWEAK" "CALCR" "FGF" "BMP" "BAG" "PROS" [34] "KIT" "HGF" $T2L [1] "MIF" "CypA" "MK" "PARs" "GALECTIN" "VISFATIN" "TNF" "VEGF" "CCL" "TGFb" [11] "CXCL" "IFN-II" "ANNEXIN" "PTN" "EGF" "GRN" "ANGPT" "OSM" "GAS" "ANGPTL" [21] "LIGHT" "TWEAK" "COMPLEMENT" "PDGF" "PROS" "IL16" "CD137" "CD70" "SLIT" "CX3C" [31] "BAFF" "TRAIL" "FASLG" "SEMA3" "PLAU" "MHC-I" "WNT" "CSF" "NRG" "PERIOSTIN" [41] "APRIL" "BTLA" "BAG" "FGF" ## 必须要挑选两个样本都有的代谢通路 pathways.show <- c("MIF" ) weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) # control the edge weights across different datasets par(mfrow = c(1,2), xpd=TRUE) for (i in 1:length(object.list)) { netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "circle", edge.weight.max = weight.max[1], edge.width.max = 10, signaling.name = paste(pathways.show, names(object.list)[i])) }所有细胞类型围成一个圆； 边表示：MIF 通路下的细胞—细胞通讯； 边粗细代表通讯强度
特点：不区分“谁是信号发送者 / 接收者”
强调的是：网络是否存在； 哪些细胞之间连接更强
Hierarchy plot
pathways.show <- c("MIF") weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) ## 指定接受端，细胞1-5类型为接收端 vertex.receiver = seq(1,5) for (i in 1:length(object.list)) { netVisual_aggregate(object.list[[i]], signaling = pathways.show, vertex.receiver = vertex.receiver, edge.weight.max = weight.max[1], edge.width.max = 10, signaling.name = paste(pathways.show, names(object.list)[i])) }Chord diagram
pathways.show <- c("MIF") par(mfrow = c(1,2), xpd=TRUE) for (i in 1:length(object.list)) { netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "chord", signaling.name = paste(pathways.show, names(object.list)[i])) }重新指定分组
## 跟上一期一样，我们把T/NK细胞合并为一组 group.merged <- c(rep("T/NK Cells", 5), "B Cells", "Myeloid Cells", "Fibroblasts Cells", "Epithelial Cells", "Endothelial") names(group.merged) <- levels(object.list[[1]]@idents) pathways.show <- c("MIF") par(mfrow = c(1,2), xpd=TRUE) for (i in 1:length(object.list)) { netVisual_chord_cell(object.list[[i]], signaling = pathways.show, group = group.merged, title.name = paste0(pathways.show, " signaling network - ", names(object.list)[i])) }于其他的可视化，大家可以参考cellchat的第二篇，双（多）样本对比。
那么到这里cellchat的流程内容就是更新完了，祝大家科研顺利！