library(Seurat)library(scCustomize)library(harmony)library(ggpubr)# 第一部分：基因表达矩阵构建与热图绘制（基于目标基因的组织表达水平）# 读取组织特异性新增差异基因（DEGs）数据add_gene <- read.csv("./results/tissue based add DEGs.csv")# 合并top10基因与新增差异基因，得到最终用于分析的基因集合（去重）gene_used <- union(top10, add_gene$gene)# 计算目标基因在各组织中的平均表达量# eos_0.8是单细胞RNA测序数据对象，提取RNA assay中的目标基因表达矩阵，按组织分组求均值markers_mean_exp <-   aggregate(    as.matrix(t(eos_0.8@assays$RNA@data[gene_used, ])),  # 转置表达矩阵使基因为列、细胞为行    list(Tissue = eos_0.8@meta.data[ ,"Tissue"]),        # 按 Tissue 列分组    mean)                                                 # 计算每组（组织）的平均表达量row.names(markers_mean_exp) <- markers_mean_exp$Tissue  # 将 Tissue 列设为行名markers_mean_exp$Tissue <- c()                           # 删除原 Tissue 列，保留纯表达矩阵# 定义z-score标准化函数（消除基因表达量量级差异，便于跨基因比较）zscore <- function(x){  return( (x - mean(x))/sd(x) )  # (原始值 - 均值) / 标准差}# 对表达矩阵按列（基因）进行z-score标准化，然后转置使基因为行、组织为列markers_plot_matrix <- apply(markers_mean_exp, 2, zscore) %>% t()# 确定每个基因表达量最高的组织（用于后续基因排序）max.avg <- apply(markers_plot_matrix, 1, which.max)# 按「基因高表达组织」分组，对基因进行排序（同组织高表达基因聚类）gene_order <- c()for(i in 1:ncol(markers_plot_matrix)){  # 遍历每个组织（列）  if(sum(max.avg == i) > 1){            # 若多个基因在该组织高表达    # 按该组织的表达量降序排序    temp <- data.frame(      gene = names(sort(markers_plot_matrix[names(max.avg)[max.avg == i],i], decreasing = T)),       Tissue = colnames(markers_plot_matrix)[i], stringsAsFactors = F)  }else if(sum(max.avg == i) == 1){      # 若仅1个基因在该组织高表达    temp <- data.frame(gene = names(max.avg)[max.avg == i],                        Tissue = colnames(markers_plot_matrix)[i])  }else{    temp <- c()  # 无基因在该组织高表达时为空  }  gene_order <- rbind(gene_order, temp)  # 合并所有基因的排序信息}# 定义热图颜色映射（8级渐变色，基于标准化后的表达量范围）color_used <- circlize::colorRamp2(seq(min(markers_plot_matrix),                                        max(markers_plot_matrix),                                        length = 8),                                    RColorBrewer::brewer.pal(9,'BuPu')[1:8])# 按排序后的基因顺序整理绘图数据plotdata <- markers_plot_matrix[gene_order$gene, ]# 绘制基因表达热图（不聚类行和列，保留手动排序结果）p <- Heatmap(plotdata, cluster_rows = FALSE, cluster_columns = FALSE,             show_row_names = T, column_names_gp = gpar(fontsize = 10),  # 列名（组织）字体大小             row_names_gp = gpar(fontsize = 6),  # 行名（基因）字体大小             col = color_used, name = "Exp")     # 颜色映射及图例名称# 第二部分：GO功能富集分析与通路热图绘制（对应Figure 1h）# 构建各组织差异基因列表（DEGs_BM、DEGs_Blood等为各组织预计算的差异基因数据框）csv_list <- list(BM = DEGs_BM, Blood = DEGs_Blood, eWAT = DEGs_eWAT,                 BAT = DEGs_BAT, scWAT = DEGs_scWAT, SI = DEGs_SI,                 Colon = DEGs_Colon, Lung = DEGs_Lung)# 加载功能富集分析所需包（clusterProfiler用于富集分析，org.Mm.eg.db为小鼠基因注释数据库）library(clusterProfiler)library(org.Mm.eg.db)# 准备背景基因集（所有组织差异基因的集合，转换为ENTREZID格式）all_genes <- unique(unlist(lapply(csv_list, rownames)))  # 提取所有组织的差异基因符号（SYMBOL）universe_entrezid <- bitr(all_genes, fromType = "SYMBOL",  # 符号转ENTREZID                          toType = c("ENTREZID"),                           OrgDb = org.Mm.eg.db)universe_entrezid <- unique(universe_entrezid$ENTREZID)  # 去重背景基因# 筛选各组织的top200差异基因（过滤线粒体基因、核糖体基因等，调整后P值<0.05）gene_list <- list()for(tissue in names(csv_list)){  # 遍历每个组织  top_gene <-    csv_list[[tissue]] %>%    filter(!grepl("^mt-|^Rps|^Rpl|^ENS|Rik", Symbol),  # 排除线粒体、核糖体等无关基因           adj.P.Val < 0.05, Grp != "others") %>%       # 差异显著且非"others"分组    slice_max(order_by = abs(logFC), n = 200) %>%      # 按logFC绝对值取top200    pull(Symbol)  # 提取基因符号  gene.df <- bitr(top_gene,  # 符号转ENTREZID                  fromType = "SYMBOL",                  toType = c("ENTREZID"),                  OrgDb = org.Mm.eg.db)  gene_list[[tissue]] <- unique(gene.df$ENTREZID)  # 去重后存入列表}# 批量进行GO Biological Process（BP）功能富集分析cluster_go <- compareCluster(gene_list, fun = "enrichGO",                             OrgDb = org.Mm.eg.db, ont = "BP",  # 小鼠数据库，BP本体                             universe = universe_entrezid,      # 背景基因集                             pvalueCutoff = 0.05,              # P值阈值                             qvalueCutoff = 0.05,              # 校正后P值阈值                             readable = TRUE, pool = TRUE)     # 输出基因符号，合并富集结果# 提取富集分析结果并转换为数据框BP_result<-as.data.frame(cluster_go@compareClusterResult)# （可选）从文件读取预计算的富集结果（避免重复分析）BP_result<-read.csv("./results/BP_compareCluster_results.csv",row.names = 1)# 读取手动筛选的关键通路列表（根据研究重点选择的BP通路）pathway<-read.csv("./results/chosen BP pathway new.csv")pathway_used <- pathway %>%  pull(Description) %>%  unique()  # 提取通路描述并去重# 整理通路富集热图数据（按组织×通路整理校正后P值的-log10转换值）heatmap_df <- BP_result %>%   filter(Description %in% pathway_used) %>%  # 仅保留筛选后的通路  dplyr::select(Cluster, Description, p.adjust) %>%  # 选择组织、通路、校正后P值  mutate(logp = -log10(p.adjust)) %>%  # 转换为-log10(p.adjust)（值越大富集越显著）  dplyr::select(-p.adjust) %>%  # 删除原P值列  pivot_wider(names_from = Cluster, values_from = logp, values_fill = 0) %>%  # 宽格式转换  tibble::column_to_rownames(var = "Description")  # 通路设为行名# 定义组织顺序（按研究设计的逻辑顺序排列）Tissue_order <- c("BM","Blood","eWAT","BAT","scWAT","SI","Colon","Lung")heatmap_df <- heatmap_df[ ,Tissue_order]  # 按顺序调整列# 对通路富集矩阵按行（通路）进行z-score标准化markers_plot_matrix <- apply(heatmap_df, 1, zscore) %>% t()# 确定每个通路富集最显著的组织（用于通路排序）max.avg <- apply(markers_plot_matrix, 1, which.max)# 按「通路富集最显著组织」分组，对通路进行排序（同组织富集通路聚类）pathway_order <- c()for(i in 1:ncol(markers_plot_matrix)){  # 遍历每个组织（列）  if(sum(max.avg == i) > 1){            # 若多个通路在该组织富集显著    # 按该组织的富集分数降序排序    temp <- data.frame(      pathway = names(sort(markers_plot_matrix[names(max.avg)[max.avg == i],i], decreasing = T)),       tissue = colnames(markers_plot_matrix)[i], stringsAsFactors = F)  }else if(sum(max.avg == i) == 1){      # 若仅1个通路在该组织富集显著    temp <- data.frame(pathway = names(max.avg)[max.avg == i],                        tissue = colnames(markers_plot_matrix)[i])  }else{    temp <- c()  # 无通路在该组织富集显著时为空  }  pathway_order <- rbind(pathway_order, temp)  # 合并所有通路的排序信息}pathway_order <- pathway_order[nrow(pathway_order):1,]  # 反转顺序（可选，调整绘图布局）pathways_labeled <- pathway_used  # 通路标签# 定义通路热图颜色映射（同基因表达热图）color_used <- circlize::colorRamp2(seq(min(markers_plot_matrix),                                        max(markers_plot_matrix),                                        length = 8),                                    RColorBrewer::brewer.pal(9,'BuPu')[1:8])# 按排序后的通路顺序整理绘图数据plotdata_BP<-as.data.frame(markers_plot_matrix[pathway_order$pathway,])# 绘制通路富集热图（不聚类行和列，保留手动排序结果）p <- Heatmap(as.matrix(plotdata_BP),             cluster_rows = FALSE,              cluster_columns = FALSE,             show_row_names = T,             column_names_gp = gpar(fontsize = 10),  # 组织列名字体大小             row_names_gp = gpar(fontsize = 9),       # 通路行名字体大小             row_names_max_width = max_text_width(rownames(plotdata_BP)),  # 适配长通路名称             col = color_used,             show_heatmap_legend = FALSE)  # 隐藏图例（可根据需要调整为TRUE）