醋醋百科网

# 抑制加载包时的提示信息，避免输出冗余suppressMessages({    library(plyr)          # 数据处理工具包    library(CellChat)      # 单细胞通讯分析核心包    library(patchwork)     # 图形组合工具    library(Seurat)        # 单细胞数据分析基础包    library(SingleCellExperiment)  # 单细胞数据结构包    library(ggalluvial)    # 桑基图绘制工具    library(repr)          # 图形尺寸调整工具    library(ggplot2)       # 基础可视化包    library(gplots)        # 热图绘制（heatmap.2函数）    library(RColorBrewer)  # 颜色 palette 生成    library(tidyr)         # 数据重塑（长表/宽表转换）    library(stringr)       # 字符串处理    library(ComplexHeatmap)# 复杂热图绘制})# 全局设置：字符串不自动转换为因子型，避免数据处理错误options(stringsAsFactors = FALSE)# 自定义运算符：判断元素是否"不在"向量中（与%in%相反）'%!in%' <- function(x,y)!('%in%'(x,y))# 定义数据输出路径path_data="/cellchat/output_files/"# 加载预处理后的CellChat对象（分别对应肥胖、瘦、减重三种状态）# 数据已过滤：排除淋巴B细胞，仅保留特定细胞类型，经过严格子采样obese=readRDS(paste0(path_data,"cellchat_obese_no_lymphatic_bcells_kit_scott_only_ct_fine_stringent_subsampled.rds"))lean=readRDS(paste0(path_data,"cellchat_lean_no_lymphatic_bcells_kit_scott_only_ct_fine_stringent_subsampled.rds"))wl=readRDS(paste0(path_data,"cellchat_weightloss_no_lymphatic_bcells_kit_scott_only_ct_fine_stringent_subsampled.rds"))# 减重组（wl）缺少Mu5细胞集群，需与肥胖组保持细胞类型一致（通过liftCellChat函数对齐）group.new = levels(obese@idents)  # 以肥胖组的细胞类型为参考wl <- liftCellChat(wl, group.new)  # 对齐减重组细胞类型lean <- liftCellChat(lean, group.new)  # 对齐瘦组细胞类型# 合并不同状态的CellChat对象（用于组间比较）object.list <- list(Obese = obese, Lean = lean)obese_lean <- mergeCellChat(object.list, add.names = names(object.list))  # 肥胖vs瘦object.list <- list(Obese = obese, Weightloss = wl)obese_wl <- mergeCellChat(object.list, add.names = names(object.list))  # 肥胖vs减重object.list <- list(Obese = obese, Lean = lean, Weightloss = wl)cellchat_all <- mergeCellChat(object.list, add.names = names(object.list))  # 三者合并# 计算各状态下细胞通讯网络的中心性指标（基于推断的通讯网络netP）# 中心性指标反映细胞在通讯中的"角色"（如发送者、接收者等）obese <- netAnalysis_computeCentrality(obese, slot.name = "netP")lean <- netAnalysis_computeCentrality(lean, slot.name = "netP")wl <- netAnalysis_computeCentrality(wl, slot.name = "netP")# --------------------------# 提取并整理中心性指标数据# --------------------------# 获取所有组中出现的通讯通路（避免通路遗漏）pathways.show.all <- unique(c(obese@netP$pathways,wl@netP$pathways,lean@netP$pathways))# 定义中心性指标及其对应中文含义（用于结果标注）measure = c("outdeg","indeg","flowbet","info")  # 算法层面的指标名称measure.name = c("Sender","Receiver","Mediator","Influencer")  # 生物学解读名称# 提取肥胖组的中心性数据并标准化mat2=NULL  # 临时存储数据框for(j in pathways.show.all){  # 遍历所有通路    centr=obese@netP$centr[j]  # 获取该通路的中心性结果    for(i in 1:length(centr)) {        centr0 <- centr[[i]]        if(!is.null(centr0)){  # 跳过空结果            # 将列表格式的中心性指标转换为矩阵            mat <- matrix(unlist(centr0), ncol = length(centr0), byrow = FALSE)            mat <- t(mat)  # 转置矩阵（行：指标，列：细胞类型）            rownames(mat) <- names(centr0); colnames(mat) <- names(centr0$outdeg)
            # 筛选需要的指标并替换为生物学名称            if (!is.null(measure)) {              mat <- mat[measure,]              if (!is.null(measure.name)) {                rownames(mat) <- measure.name              }            }            # 按行标准化（每行最大值为1），消除不同指标的量纲差异            mat <- sweep(mat, 1L, apply(mat, 1, max), '/', check.margin = FALSE)            mat=as.data.frame(mat)  # 转换为数据框            mat$Pathway=j  # 添加通路标签            mat$Role=rownames(mat)  # 添加角色标签（Sender等）
            # 合并到结果数据框            if (is.null(mat2)){                mat2=data.frame(mat)            } else {                mat2=rbind(mat2,mat)            }        }        }}mat2$Condition="Obese"  # 添加组别标签mat_obese=mat2  # 保存肥胖组结果# 重复上述步骤，提取瘦组的中心性数据mat2=NULLfor(j in pathways.show.all){    centr=lean@netP$centr[j]    for(i in 1:length(centr)) {        centr0 <- centr[[i]]        if(!is.null(centr0)){            mat <- matrix(unlist(centr0), ncol = length(centr0), byrow = FALSE)            mat <- t(mat)            rownames(mat) <- names(centr0); colnames(mat) <- names(centr0$outdeg)            if (!is.null(measure)) {              mat <- mat[measure,]              if (!is.null(measure.name)) {                rownames(mat) <- measure.name              }            }            mat <- sweep(mat, 1L, apply(mat, 1, max), '/', check.margin = FALSE)            mat=as.data.frame(mat)            mat$Pathway=j            mat$Role=rownames(mat)            if (is.null(mat2)){                mat2=data.frame(mat)            } else {                mat2=rbind(mat2,mat)            }        }        }}mat2$Condition="Lean"mat_lean=mat2# 重复上述步骤，提取减重组的中心性数据mat2=NULLfor(j in pathways.show.all){    centr=wl@netP$centr[j]    for(i in 1:length(centr)) {        centr0 <- centr[[i]]        if(!is.null(centr0)){            mat <- matrix(unlist(centr0), ncol = length(centr0), byrow = FALSE)            mat <- t(mat)            rownames(mat) <- names(centr0); colnames(mat) <- names(centr0$outdeg)            if (!is.null(measure)) {              mat <- mat[measure,]              if (!is.null(measure.name)) {                rownames(mat) <- measure.name              }            }            mat <- sweep(mat, 1L, apply(mat, 1, max), '/', check.margin = FALSE)            mat=as.data.frame(mat)            mat$Pathway=j            mat$Role=rownames(mat)            if (is.null(mat2)){                mat2=data.frame(mat)            } else {                mat2=rbind(mat2,mat)            }        }        }}mat2$Condition="Weightloss"mat_wl=mat2# 合并三组中心性数据并输出为CSVmat_final=rbind(mat_obese,mat_lean,mat_wl)head(mat_final)  # 预览前几行write.csv(mat_final, paste0(path_data,"cellchat_ct_fine_no_lymphatic_bcells_kit_stringent_subsampled_scott_only_centrality_all.csv"), row.names=FALSE)# --------------------------# 提取通讯网络的详细连接信息# --------------------------# 提取肥胖组的通讯连接数据（如发送细胞、接收细胞、通路、强度等）df.net <- subsetCommunication(obese)unique(df.net$pathway_name)  # 查看该组涉及的通路write.csv(df.net, paste0(path_data,"cellchat_ct_fine_no_lymphatic_bcells_kit_stringent_subsampled_scott_only_net_obese.csv"), row.names=FALSE)# 提取瘦组的通讯连接数据df.net <- subsetCommunication(lean)write.csv(df.net, paste0(path_data,"cellchat_ct_fine_no_lymphatic_bcells_kit_stringent_subsampled_scott_only_net_lean.csv"), row.names=FALSE)# 提取减重组的通讯连接数据df.net <- subsetCommunication(wl)write.csv(df.net, paste0(path_data,"cellchat_ct_fine_no_lymphatic_bcells_kit_stringent_subsampled_scott_only_net_wl.csv"), row.names=FALSE)# --------------------------# 组间通讯差异分析（通路水平）# --------------------------options(repr.plot.width = 6, repr.plot.height = 6)  # 设置图形尺寸# 初始化存储差异分析结果的列表（贡献度和P值）heatmap.df <- list(contrib = data.frame(matrix(ncol = 0, nrow = 0)), pval = data.frame(matrix(ncol = 0, nrow = 0)))cellchat=cellchat_all  # 使用之前合并的所有组数据genotypes <- names(cellchat@net)  # 获取组别名称（Obese, Lean, Weightloss）# 比较肥胖组（1）与瘦组（2）、肥胖组（1）与减重组（3）for (i in c(2:3)){    # 计算两组间的通讯差异（以肥胖组为对照）    gg <- rankNet(cellchat, mode = "comparison", comparison=c(1,i)            , title= genotypes[i]            , cutoff.pvalue = 0.05  # 初步筛选P值阈值            , stacked = T, show.raw = T, do.stat = TRUE, return.data = FALSE)
    # 提取处理组（Lean/Weightloss）和对照组（Obese）的贡献度及P值    temp <- data.frame(        gg$data[gg$data$group == levels(gg$data$group)[1],c("contribution", "pvalues")],  # 处理组的贡献度和P值        control.contribution=gg$data[gg$data$group == levels(gg$data$group)[2],c("contribution")],  # 对照组的贡献度        row.names = rownames(gg$data)[gg$data$group == levels(gg$data$group)[2]]  # 行名设为通路    )
    # 计算处理组/对照组的贡献度比值（反映相对变化）    temp2 <- data.frame((temp[,1]/temp[,3]), row.names = rownames(temp))    colnames(temp2) <- genotypes[i]  # 列名设为组别（Lean/Weightloss）
    # 合并贡献度结果    if(sum(dim(heatmap.df$contrib) == c(0,0))){        heatmap.df$contrib <- temp2    } else{        heatmap.df$contrib <- merge(heatmap.df$contrib, temp2, all = TRUE, by = "row.names")        rownames(heatmap.df$contrib) <- heatmap.df$contrib$Row.names        heatmap.df$contrib <- heatmap.df$contrib[,c(-1)]  # 移除合并产生的Row.names列    }
    # 对P值进行多重检验校正（Bonferroni法）    temp2 <- data.frame(p.adjust(temp[,2], method="bonferroni"), row.names = rownames(temp))    colnames(temp2) <- genotypes[i]
    # 合并校正后的P值结果    if(sum(dim(heatmap.df$pval) == c(0,0))){        heatmap.df$pval <- temp2    } else {        heatmap.df$pval <- merge(heatmap.df$pval, temp2, all = TRUE, by = "row.names")        rownames(heatmap.df$pval) <- heatmap.df$pval$Row.names        heatmap.df$pval <- heatmap.df$pval[,c(-1)]    }}# 对贡献度比值取log2（将倍数变化转换为对称分布）heatmap.df$contrib <- log2(heatmap.df$contrib)# 确保数据格式正确（数值型）heatmap.df$pval <- data.frame(sapply(heatmap.df$pval, as.numeric), row.names=rownames(heatmap.df$pval))heatmap.df$contrib <- data.frame(sapply(heatmap.df$contrib, as.numeric), row.names=rownames(heatmap.df$contrib))# 将宽表转换为长表（便于后续可视化或统计）# 贡献度长表（包含通路、组别、log2倍数变化）temp <- heatmap.df$contribtemp$Pathway <- rownames(temp)heatmap.df.contrib.wide <- gather(temp, Group, log2FC_contrib, Lean:Weightloss, factor_key=TRUE)print(dim(heatmap.df.contrib.wide))head(heatmap.df.contrib.wide)# P值长表（包含通路、组别、校正后P值）temp <- heatmap.df$pvaltemp$Pathway <- rownames(temp)heatmap.df.pval.wide <- gather(temp, Group, adjusted_pval,Lean:Weightloss, factor_key=TRUE)print(dim(heatmap.df.pval.wide))head(heatmap.df.pval.wide)# 合并贡献度和P值的长表heatmap.df.wide <- merge(heatmap.df.contrib.wide, heatmap.df.pval.wide, by=c("Pathway","Group"), all.x = TRUE, all.y = TRUE)head(heatmap.df.wide)# 查看特定通路（如NRG）的结果heatmap.df.wide[which(heatmap.df.wide$Pathway=="NRG"),]# --------------------------# 绘制通路水平的差异热图# --------------------------# 确定颜色映射的范围（排除无限值和NA）MAX.CONTRIB <- max(abs(heatmap.df$contrib[!(sapply(heatmap.df$contrib, is.infinite) | is.na(heatmap.df$contrib))]))MAX.CONTRIB <- round(MAX.CONTRIB)  # 取整便于设置颜色梯度INF.VALUE <- MAX.CONTRIB+1  # 用该值替代无限值（避免绘图错误）# 过滤弱差异（log2FC绝对值<0.2视为无差异，P值设为1）heatmap.df.wide$adjusted_pval[heatmap.df.wide$log2FC_contrib<0.2]=1heatmap.df.wide$adjusted_pval[heatmap.df.wide$log2FC_contrib<0.2]=1  # 重复过滤确保生效options(repr.plot.width = 20, repr.plot.height = 5)  # 调整热图尺寸# 处理热图数据（替换无限值、NA，标记显著性）heatmap.df2 <- heatmap.df# 替换无限值（正无限→INF.VALUE，负无限→-INF.VALUE）heatmap.df2$contrib[sapply(heatmap.df2$contrib, is.infinite) & heatmap.df2$contrib > 0] <- INF.VALUEheatmap.df2$contrib[sapply(heatmap.df2$contrib, is.infinite) & heatmap.df2$contrib < 0] <- -INF.VALUE# 弱差异的P值设为1（不显著），NA的P值设为2（后续标记为“-”）heatmap.df2$pval[(abs(heatmap.df2$contrib) < 0.2)] <- 1heatmap.df2$pval[is.na(heatmap.df2$pval)] <- 2heatmap.df2$pval <- heatmap.df2$pval[(apply(is.na(heatmap.df2$contrib), 1, sum) < 2),]heatmap.df2$contrib <- heatmap.df2$contrib[(apply(is.na(heatmap.df2$contrib), 1, sum) < 2),]
heatmap.df2$contrib <- heatmap.df2$contrib[(apply(heatmap.df2$pval > 0.01, 1, sum) < 2),]heatmap.df2$pval <- heatmap.df2$pval[(apply(heatmap.df2$pval > 0.01, 1, sum) < 2),]
heatmap.df2$pval[heatmap.df2$pval == 2] <- NA
heatmap.df2$sig <- heatmap.df2$pvalheatmap.df2$sig[] <- ''heatmap.df2$sig[heatmap.df2$pval <= 0.01] <- '*'
heatmap.df2$sig[is.na(heatmap.df2$pval)] <- '-'
heatmap.df2$contrib[is.na(heatmap.df2$contrib)] <- 0
heatmap.df2$sig <- heatmap.df2$sig[,c("Lean","Weightloss")]heatmap.df2$pval <- heatmap.df2$pval[,c("Lean","Weightloss")]heatmap.df2$contrib <- heatmap.df2$contrib[,c("Lean","Weightloss")]
scale.interval.size <- 0.05num_breaks <- (MAX.CONTRIB * 2 / scale.interval.size)+1breaks <- c(-(INF.VALUE+1), seq(from=-MAX.CONTRIB, to=MAX.CONTRIB, length.out=num_breaks), INF.VALUE+1)midpoint <- 0 # the mid of the "real" scale
rampCol2 <- colorRampPalette(c("#6699FF", "white", "#FF6600"))(length(breaks)-1)mypalette <- c(rampCol2)mypalette[1] <- "#5689EF" # just a random extreme color for -infmypalette[num_breaks+1] <- "#EF5600" # just a random extreme color for inf

heatmap.df3 <- heatmap.df2heatmap.df4 <- heatmap.df2
heatmap.df3$contrib <- heatmap.df2$contrib[(apply(heatmap.df2$contrib > 0 & heatmap.df2$sig != '', 1, sum) >=2 & apply(heatmap.df2$contrib < 0 & heatmap.df2$sig != '', 1, sum) == 0)                                           |(apply(heatmap.df2$contrib < 0 & heatmap.df2$sig != '', 1, sum) >=2 & apply(heatmap.df2$contrib > 0 & heatmap.df2$sig != '', 1, sum) == 0) ,]heatmap.df3$sig <- heatmap.df2$sig[(apply(heatmap.df2$contrib > 0 & heatmap.df2$sig != '', 1, sum) >=2 & apply(heatmap.df2$contrib < 0 & heatmap.df2$sig != '', 1, sum) == 0)                                           |(apply(heatmap.df2$contrib < 0 & heatmap.df2$sig != '', 1, sum) >=2 & apply(heatmap.df2$contrib > 0 & heatmap.df2$sig != '', 1, sum) == 0) ,]heatmap.df3$sig <- heatmap.df3$sig[order(apply(heatmap.df3$contrib, 1, mean)),]heatmap.df3$contrib <- heatmap.df3$contrib[order(apply(heatmap.df3$contrib, 1, mean)),]
heatmap.df4$contrib <- heatmap.df2$contrib[!((apply(heatmap.df2$contrib > 0 & heatmap.df2$sig != '', 1, sum) >=2 & apply(heatmap.df2$contrib < 0 & heatmap.df2$sig != '', 1, sum) == 0)                                           |(apply(heatmap.df2$contrib < 0 & heatmap.df2$sig != '', 1, sum) >=2 & apply(heatmap.df2$contrib > 0 & heatmap.df2$sig != '', 1, sum) == 0)) ,]heatmap.df4$sig <- heatmap.df2$sig[!((apply(heatmap.df2$contrib > 0 & heatmap.df2$sig != '', 1, sum) >=2 & apply(heatmap.df2$contrib < 0 & heatmap.df2$sig != '', 1, sum) == 0)                                           |(apply(heatmap.df2$contrib < 0 & heatmap.df2$sig != '', 1, sum) >=2 & apply(heatmap.df2$contrib > 0 & heatmap.df2$sig != '', 1, sum) == 0)) ,]
heatmap.df4$contrib <- heatmap.df4$contrib[with(data.frame(heatmap.df4$sig == '*'), order(-Lean,-Weightloss)),]heatmap.df4$sig <- heatmap.df4$sig[with(data.frame(heatmap.df4$sig == '*'), order(-Lean,-Weightloss)),]
heatmap.df4$change <- heatmap.df4$contribheatmap.df4$change[] <- 0heatmap.df4$change[heatmap.df4$contrib > 0 & heatmap.df4$sig != ''] <- 1heatmap.df4$change[heatmap.df4$contrib < 0 & heatmap.df4$sig != ''] <- -1heatmap.df4$sig <- heatmap.df4$sig[order(apply(heatmap.df4$change, 1, sum), decreasing = TRUE),]heatmap.df4$contrib <- heatmap.df4$contrib[order(apply(heatmap.df4$change, 1, sum), decreasing = TRUE),]
heatmap.df4$contrib <- heatmap.df4$contrib[order(apply(heatmap.df4$sig == '*', 1, sum)),]heatmap.df4$sig <- heatmap.df4$sig[order(apply(heatmap.df4$sig == '*', 1, sum)),]
heatmap.2(t(heatmap.df3$contrib)        , col=mypalette        , Rowv=NULL          , Colv=NULL        , dendrogram="none"        , na.rm = TRUE        , breaks=breaks, density.info="none", trace="none"        , symm=F,symkey=F,symbreaks=T, scale="none"        , margins = c(33,35)        , lhei=c(1,5), lwid=c(1,7)        , cellnote = t(heatmap.df3$sig)        ,notecex=3.0        ,notecol="black"        ,cexRow=3        ,cexCol=1.5       )
heatmap.2(t(heatmap.df4$contrib)        , col=mypalette        , Rowv=NULL          , Colv=NULL        , dendrogram="none"        , na.rm = TRUE        , breaks=breaks, density.info="none", trace="none"        , symm=F,symkey=F,symbreaks=T, scale="none"        , margins = c(33,35)        , lhei=c(1,5), lwid=c(1,7)        , cellnote = t(heatmap.df4$sig)        ,notecex=3.0        ,notecol="black"        ,cexRow=3        ,cexCol=1.5       )
lean=heatmap.final.corr$pval[,"Lean",drop=FALSE]colnames(lean)=c("pval")lean$contrib=heatmap.final.corr$contrib$Leanlean$Condition="Lean"lean$Group=rownames(lean)rownames(lean)=NULL
head(lean)
wl=heatmap.final.corr$pval[,"Weightloss",drop=FALSE]colnames(wl)=c("pval")wl$contrib=heatmap.final.corr$contrib$Weightlosswl$Condition="Weightloss"wl$Group=rownames(wl)rownames(wl)=NULL
head(wl)
merged.df=rbind(lean,wl)
merged.df[c('Pathway', 'Pair')] <- str_split_fixed(merged.df$Group, '__', 2)
merged.df[c('Sender', 'Receiver')] <- str_split_fixed(merged.df$Pair, '\\|', 2)
rownames(merged.df)=NULL
head(merged.df)
write.csv(merged.df,paste0(path_data,"rankNet_vs_obese_ct_fine_pairwise_scott_only_subsampled_stringent_no_lymphatic_bcells_kit.csv")) # use this data for Sankey plots