1.8 单细胞代谢分析（6 套可复用参考代码）

#paper
# METAFlux: Characterizing metabolism from bulk and single-cell RNA-seq data
# https://www.nature.com/articles/s41467-023-40457-w

#github:https://github.com/KChen-lab/METAFlux?tab=readme-ov-file

devtools::install_github('KChen-lab/METAFlux')
library(METAFlux)
library(Seurat)


load("C:~/Downloads/adj_scRNA.RData")

# 1-scRNA数据代谢通量平衡分析 -------------------------------------------------------

#支持人的代谢分析
#METAFlux只能单个样本进行分析，如果我们的单细胞示例数据是多样本的
#所以有两种办法，第一种是subset单个样本的数据进行分析，或者使用循环

#基本分析步骤：很简单，三个函数，三个步骤
#1、calculate_avg_exp:: 计算celltype  bootstrap 平均表达量。
#为了便于理解参数设置，简单解释一下bootstrap，是一种统计方法，通过从原始数据集中有放回地重复抽样，生成多个新的样本集。
#用于估计原始数据集均值的标准误和置信区间，从而评估均值的稳定性和可靠性。
#calculate_avg_exp中的n_bootstrap就是重复抽样次数，实际使用作者建议100，200.。。。

#2、calculate_reaction_score:: #calculate metabolic reaction scores
#3、compute_sc_flux::计算代谢通量


#我们单独提取一个组的数据进行演示
SD_scRNA <- subset(adj_scRNA, orig.ident=='SD_scRNA')
mean_exp=calculate_avg_exp(myseurat = SD_scRNA, #需要分析的单细胞seurat obj,单细胞obj是normalization 的
                           myident = 'celltype',#celltype分组
                           n_bootstrap=3, #抽样次数
                           seed=1)#随机数

#计算代谢反应分数
scores<-calculate_reaction_score(data=mean_exp)
#计算这个样本中每种celltype占比，每种celltype的数量除以总的细胞数
fr =  round(table(SD_scRNA$celltype)/nrow(SD_scRNA@meta.data),2)
fr = as.vector(fr)

#计算代谢通量
flux=compute_sc_flux(num_cell = 10, #celltype数量,比如这里我的数据中有10种celltype
                     fraction = fr, #celltype比例，和必须是1
                     fluxscore=scores,
                     medium = human_blood)#人来源的细胞使用human_blood，细胞系使用cell_medium

# head(flux)
# V1            V2            V3
# celltype 1 HMR_3905 -1.654499e-08 -2.099210e-07  1.238930e-07
# celltype 1 HMR_3907 -6.997868e-09 -2.390169e-07  1.160639e-07
# celltype 1 HMR_4097  2.634345e-07 -6.731786e-08 -3.903460e-07
# celltype 1 HMR_4099  2.042502e-10  2.426213e-10 -1.625041e-10
# celltype 1 HMR_4108  1.318980e-07 -3.268554e-08 -1.953607e-07
# celltype 1 HMR_4133  1.315366e-07 -3.463232e-08 -1.949853e-07
# dim(flux)
# [1] 132468      3

#分析结果可以看到，有很多。这也就是为什么他最后一步非常耗费时间和内存的原因
#得到的结果维度和分析有关系，(num_cell∗13082+1648)∗number_of_bootstrap
#METAflux提供的代谢modle包含13082个代谢反应， 8378 个代谢物。
#行是celltype和代谢反应的组合，celltype1，2，3.。代表compute_sc_flux步骤celltype的顺序，
#列数和n_bootstrap一样，n_bootstrap设置多少，就有多少列
#正值表示化合物释放，负值表示化合物摄取。绝对值代表通量大小



#代谢model，可以查看代谢反应ID，代谢公式，代谢通路！
data("human_gem")
human_gem <- human_gem

#选择感兴趣的代谢反应，提取数据进行可视化
#比如这里看glucose的代谢反应，我们可以通过nutrient_lookup_files文件查找glucose对应的ID
data("nutrient_lookup_files")
glucose_ID <- nutrient_lookup_files$EQUATION[grepl("glucose", nutrient_lookup_files$EQUATION)]
select_ID <- subset(nutrient_lookup_files, nutrient_lookup_files$EQUATION %in% glucose_ID)
#这样我们就可以得知glucose ID
glucose<-data.frame(glucose=flux[grep("HMR_9034",rownames(flux)),])

#if needed we can apply cubic root normalization to normalize the scores
cbrt <- function(x) {
  sign(x) * abs(x)^(1/3)
}

glucose=cbrt(glucose)

#ggplot做一个箱线图展示
library(ggplot2)
glucose$celltype=rownames(glucose)
long_glucose=reshape2::melt(glucose,id.vars='celltype')
ggplot(long_glucose,aes(y=-value,x=celltype))+
  geom_boxplot()+
  ggtitle("Glucose uptake level for different cell types")+
  xlab("")+ylab("Glucose uptake scores")+ 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))


##还可以计算142个通路的pathway level activity。
##对于给定的途径，活度水平可以计算为与该途径相关的反应的平均通量。
#compute pathway level activity for all samples
#将最后的flux文件里面吧celltype做一下转化对应，不然celltype1，2，3.。不然看起来不方便

flux1 <- rowMeans(flux) %>% as.data.frame()
colnames(flux1) <- "flux"

celltype <- c("Endo","Fib","Musc","Ker","APCs","Other","Tcell","Mast","LY","Mela")#顺序和mean-exp顺序一样
celltype <- paste0(celltype," ")
flux_celltype <- paste("celltype", seq(length(celltype)),"")
for(i in 1:length(celltype)){
  
  rownames(flux1) <-  gsub(flux_celltype[i], celltype[i], rownames(flux1))
  
}

flux1$celltype <- sub(" .*$", "", rownames(flux1))
flux1 <- flux1[!grepl("external_medium", flux1$celltype), ]
flux1$pathway <- rep(human_gem$SUBSYSTEM, length(celltype))

library(tidyr)
library(dplyr)
library(pheatmap)

df_summary <- flux1 %>%
  group_by(celltype, pathway) %>%
  summarize(mean_value = mean(abs(flux)))

df_summary_w <- df_summary  %>%
  pivot_wider(names_from = pathway, values_from = mean_value) %>%
  as.data.frame()


rownames(df_summary_w) <- df_summary_w[,1]
df_summary_w <- df_summary_w[,-1]


#heatmap 
mapal <- colorRampPalette(RColorBrewer::brewer.pal(11,"RdBu"))(256)
pheatmap::pheatmap(t(df_summary_w),cluster_cols = F,color = rev(mapal),scale = "row")



#2-bulk数据代谢通量平衡分析 -------------------------------------------------------
#bulk分析和scRNA一样，只不过不需要计算calculate_avg_exp

"""MEBOCOST: a python-based package to infer metabolite-mediated cell-cell communications in scRNA-seq data.
"""
import os
import setuptools
import configparser

def _make_config(conf_path, workdir=os.getcwd()):
    """
    read config file
    """
    #read config
    cf = configparser.ConfigParser()
    cf.read(conf_path)
    config = cf._sections
    # remove the annotation:
    for firstLevel in config.keys():
        for secondLevel in config[firstLevel]:
            if '#' in config[firstLevel][secondLevel]:
                path = config[firstLevel][secondLevel][:config[firstLevel][secondLevel].index('#')-1].rstrip()
                config[firstLevel][secondLevel] = os.path.join(workdir, path)
            else:
                path = config[firstLevel][secondLevel]
                config[firstLevel][secondLevel] = os.path.join(workdir, path)
    ## re-write
    cf_new = configparser.ConfigParser()
    for firstLevel in config.keys():
        cf_new.add_section(firstLevel)
        for secondLevel in config[firstLevel]:
            cf_new.set(firstLevel, secondLevel, config[firstLevel][secondLevel])
    with open('mebocost.conf', 'w') as f:
        cf.write(f)
    return(config)

## setup
def main():
  setuptools.setup(name="mebocost", 
                  version="1.2.2",
                  description="a python-based method to infer metabolite mediated cell-cell communication",
                  author='Rongbin Zheng, Kaifu Chen',
                  author_email='<redacted-email>',
                  url='https://kaifuchenlab.github.io/',
                  zip_safe=True,
                  package_dir={"": "src"},
                  packages=setuptools.find_packages(where="src"),
                  classifiers=[
                      'Environment::Console',
                      'Operating System:: POSIX',
                      "Topic :: Scientific/Engineering :: Bio-Informatics"],
                  keywords='Metabolism',
                  license='OTHER'
  )
if __name__ == '__main__':
    ## change mebocost.conf to absolute path
    _make_config(conf_path = './src/mebocost.conf')
    ## setup
    main()

# ==============================================================================
# decoupleR 实战 01：数据准备与 PROGENy 通路网络加载
# 来源：KS 教程《跟着 NCB 学 Bulk RNA-seq 数据通路活性推断 | decoupleR 实战》
# 场景：胰腺癌细胞 FOXA2 KO vs WT RNA-seq，推断删除 FOXA2 后哪些信号通路受影响
# 输入：bk_data.rds（含 counts / design / limma_ttop）
# 输出：decoupleR_input.rds（供 02-04 脚本直接读取）
# ==============================================================================

library(decoupleR)   # 通路/TF 活性推断核心包
library(dplyr)
library(tibble)
library(tidyr)
library(ggplot2)
library(pheatmap)
library(ggrepel)

# ---- 1. 读取教程示例数据 ------------------------------------------------------
# bk_data.rds 是一个 list：
#   $counts     基因 x 样本 表达矩阵（带 gene 列）
#   $design     分组信息（样本 x 条件）
#   $limma_ttop limma 差异分析结果（WT vs KO）
# ⚠️ 大数据不入库：请从本机原始目录提供该文件，或换成自己的数据
data <- readRDS("bk_data.rds")

# ---- 2. 表达矩阵 counts：缺失值填 0，gene 列转行名 ----------------------------
counts <- data$counts %>%
  dplyr::mutate_if(~ any(is.na(.x)),
                   ~ dplyr::if_else(is.na(.x), 0, .x)) %>%
  tibble::column_to_rownames(var = "gene") %>%
  as.matrix()

head(counts)

# 分组信息
design <- data$design
head(design)

# ---- 3. 差异分析统计量：只保留 t 值（方向 + 显著性） --------------------------
# 对比水平通路活性推断的输入；limma t 或 DESeq2 stat 均可
deg <- data$limma_ttop %>%
  dplyr::select(ID, t) %>%
  dplyr::filter(!is.na(t)) %>%
  tibble::column_to_rownames(var = "ID") %>%
  as.matrix()

head(deg)

# ---- 4. 加载 PROGENy 信号通路-靶基因网络 --------------------------------------
# PROGENy（Pathway RespOnsive GENes）：人工审编的信号通路及其靶基因，
# weight 为该基因对通路活性的响应强度（可正可负）
# top = 500：每个通路取按 p 值排序的前 500 个最显著响应基因
# 物种换小鼠：get_progeny(organism = 'mouse', top = 500)
net <- decoupleR::get_progeny(organism = 'human',
                              top = 500)

head(net)
dim(net)

# ---- 5. 保存中间对象（供后续脚本直接读取） ------------------------------------
saveRDS(list(counts = counts, design = design, deg = deg, net = net),
        "decoupleR_input.rds")

# ==============================================================================
# decoupleR 实战 02：样本水平通路活性推断（MLM）+ 热图
# 原理：多元线性模型（MLM）——把所有通路的靶基因权重放进同一个模型联合拟合，
#       评估每个样本中各信号通路的活性（与 ULM 单变量逐一拟合不同）
# 输入：decoupleR_input.rds（01 脚本产出）
# 输出：通路活性热图.pdf（矢量 PDF）
# ==============================================================================

library(decoupleR)
library(dplyr)
library(tibble)
library(tidyr)
library(pheatmap)
library(RColorBrewer)

input  <- readRDS("decoupleR_input.rds")
counts <- input$counts
net    <- input$net

# ---- 1. MLM 通路活性推断 ------------------------------------------------------
# .mor = 'weight'：PROGENy 的权重列名（区别于 CollecTRI 的 mor）
# minsize = 5：至少 5 个靶基因的通路才参与推断
sample_acts <- decoupleR::run_mlm(mat = counts,
                                  net = net,
                                  .source = 'source',
                                  .target = 'target',
                                  .mor = 'weight',
                                  minsize = 5)

head(sample_acts)

# ---- 2. 长表转宽矩阵：行 = 样本(condition)，列 = 通路 -------------------------
sample_acts_mat <- sample_acts %>%
  tidyr::pivot_wider(id_cols = 'condition',
                     names_from = 'source',
                     values_from = 'score') %>%
  tibble::column_to_rownames('condition') %>%
  as.matrix()

# ---- 3. 按通路（列）做 z-score 标准化，便于跨通路比较相对活性 -----------------
sample_acts_mat <- scale(sample_acts_mat)

# ---- 4. 配色与断点：RdBu 反转，红 = 高活性，蓝 = 低活性 -----------------------
colors <- rev(RColorBrewer::brewer.pal(n = 11, name = "RdBu"))
colors.use <- grDevices::colorRampPalette(colors = colors)(100)

my_breaks <- c(seq(-2, 0, length.out = ceiling(100 / 2) + 1),
               seq(0.05, 2, length.out = floor(100 / 2)))

# ---- 5. 热图（filename 直接落盘为矢量 PDF） ------------------------------------
pheatmap::pheatmap(mat = sample_acts_mat,
                   color = colors.use,
                   border_color = "white",
                   breaks = my_breaks,
                   cellwidth = 20,
                   cellheight = 20,
                   treeheight_row = 20,
                   treeheight_col = 20,
                   filename = "通路活性热图.pdf")

# ==============================================================================
# decoupleR 实战 03：差异水平通路活性推断（基于 DEGs 的 WT vs KO 对比）+ 条形图
# 输入：差异分析 t 值矩阵 deg——代表 WT vs KO 的整体变化方向与显著性
# 输出：通路活性差异_barplot.pdf（矢量 PDF）
# 注：原教程标题"DEGS方法"应为 DEGs（差异表达基因）输入，此处已更正
# ==============================================================================

library(decoupleR)
library(dplyr)
library(ggplot2)
library(RColorBrewer)

input <- readRDS("decoupleR_input.rds")
deg   <- input$deg
net   <- input$net

# ---- 1. 用 t 值做 MLM：得到每条通路在"对比"水平的活性得分 ---------------------
contrast_acts <- decoupleR::run_mlm(mat = deg,
                                    net = net,
                                    .source = 'source',
                                    .target = 'target',
                                    .mor = 'weight',
                                    minsize = 5)

contrast_acts

# ---- 2. 配色：蓝 = 活性下降，红 = 活性上升 ------------------------------------
colors <- rev(RColorBrewer::brewer.pal(n = 11, name = "RdBu")[c(2, 10)])

# ---- 3. 条形图（本教程展示全部通路，不做 top-n 过滤） -------------------------
p <- ggplot2::ggplot(data = contrast_acts,
                     mapping = ggplot2::aes(x = stats::reorder(source, score),
                                            y = score)) +
  ggplot2::geom_bar(mapping = ggplot2::aes(fill = score),
                    color = "black",
                    stat = "identity") +
  ggplot2::scale_fill_gradient2(low = colors[1],
                                mid = "whitesmoke",
                                high = colors[2],
                                midpoint = 0) +
  ggplot2::theme_minimal() +
  ggplot2::theme(axis.title = element_text(face = "bold", size = 12),
                 axis.text.x = ggplot2::element_text(angle = 45,
                                                     hjust = 1,
                                                     size = 10,
                                                     face = "bold"),
                 axis.text.y = ggplot2::element_text(size = 10,
                                                     face = "bold"),
                 panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank()) +
  ggplot2::xlab("Pathways")

p

# ---- 4. 输出矢量 PDF -----------------------------------------------------------
ggplot2::ggsave("通路活性差异_barplot.pdf", p,
                width = 10, height = 6, device = "pdf")

# ---- 5. 文献对照（NCB Fig1d，HSPCs vs progenitors） ----------------------------
# 条形朝 HSPCs 一侧 = 该通路在 HSPCs 中更活跃；只展示 P < 0.05 的差异通路
# 结果：Hypoxia 在 HSPCs 中更活跃；JAK-STAT、PI3K、MAPK 等增殖/活化信号相对更低

# ==============================================================================
# decoupleR 实战 04：深入解析单条通路的关键基因（以 MAPK 为例）
# 逻辑：把 PROGENy 中该通路的靶基因（带 weight 响应权重）映射到差异分析 t 值上，
#       weight 与 t 同号（一致响应）标红，异号（矛盾）标蓝，其余灰色
# 输入：decoupleR_input.rds
# 输出：通路关键基因散点图_MAPK.pdf（矢量 PDF）
# ==============================================================================

library(dplyr)
library(tibble)
library(ggplot2)
library(ggrepel)
library(RColorBrewer)

input <- readRDS("decoupleR_input.rds")
deg   <- input$deg
net   <- input$net

pathway <- 'MAPK'   # 换成任意感兴趣的通路（如 Hypoxia、JAK-STAT、PI3K）

# ---- 1. 取该通路的靶基因，并准备画图框 ----------------------------------------
df <- net %>%
  dplyr::filter(source == pathway) %>%
  dplyr::arrange(target) %>%
  dplyr::mutate(ID = target,
                color = "3") %>%
  tibble::column_to_rownames('target')

# ---- 2. 只保留差异分析结果里也有的靶基因 --------------------------------------
inter <- sort(dplyr::intersect(rownames(deg), rownames(df)))

df <- df[inter, ]

# deg 为单列矩阵（t 值），按行名对齐后赋值
df['t_value'] <- deg[inter, ]

# ---- 3. 着色逻辑 ---------------------------------------------------------------
# color = 1（红）：weight 与 t 同号 → 靶基因变化支持该通路活性方向（一致响应）
# color = 2（蓝）：weight 与 t 异号 → 与 PROGENy 注释的响应方向矛盾
# color = 3（灰）：weight 或 t 缺失 / 为 0
df <- df %>%
  dplyr::mutate(color = dplyr::if_else(weight > 0 & t_value > 0, '1', color)) %>%
  dplyr::mutate(color = dplyr::if_else(weight > 0 & t_value < 0, '2', color)) %>%
  dplyr::mutate(color = dplyr::if_else(weight < 0 & t_value > 0, '2', color)) %>%
  dplyr::mutate(color = dplyr::if_else(weight < 0 & t_value < 0, '1', color))

colors <- rev(RColorBrewer::brewer.pal(n = 11, name = "RdBu")[c(2, 10)])

# ---- 4. 散点图：x = PROGENy 权重，y = 差异 t 值 --------------------------------
p <- ggplot2::ggplot(data = df,
                     mapping = ggplot2::aes(x = weight,
                                            y = t_value,
                                            color = color)) +
  ggplot2::geom_point(size = 2.5,
                      color = "black") +
  ggplot2::geom_point(size = 1.5) +
  ggplot2::scale_colour_manual(values = c(colors[2], colors[1], "grey")) +
  ggrepel::geom_label_repel(mapping = ggplot2::aes(label = ID)) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none") +
  ggplot2::geom_vline(xintercept = 0, linetype = 'dotted') +
  ggplot2::geom_hline(yintercept = 0, linetype = 'dotted') +
  ggplot2::ggtitle(pathway)

p

# ---- 5. 输出矢量 PDF -----------------------------------------------------------
ggplot2::ggsave("通路关键基因散点图_MAPK.pdf", p,
                width = 8, height = 6, device = "pdf")

#首先下载KEGG分级注释
#构建KEGG代谢通路文件 https://www.kegg.jp/kegg-bin/get_htext?br08901.keg
setwd("~/data_analysis/单细胞代谢分析/")
# library(rjson)
# KEGG_pathway <- rjson::fromJSON(file = "./br08901.json",simplify=F)
# KEGG_pathway <- KEGG_pathway[[2]]#结果在2中
# KEGG_pathway <- KEGG_pathway[[1]]#只取代谢
# KEGG_pathway <- KEGG_pathway[[2]]
# 
# 
# cate1 <- c()
# for (i in 1:length(KEGG_pathway)) {
#   length_names <- length(KEGG_pathway[[i]][[2]])
#   cate1 <- append(cate1, rep(KEGG_pathway[[i]][[1]], length_names))
#   
# }
# 
# cate2 <- c()
# for (i in 1:length(KEGG_pathway)) {
#   
#   cate2 <- append(cate2, KEGG_pathway[[i]][[2]])
#   
# }
# 
# cate3 <- c()
# for (i in 1:length(cate2)) {
#   
#   cate3 <- append(cate3, as.character(cate2[[i]]))
#   
# }
# 
# #去除数字
# library(stringr)
# vec = cate3
# vec <- str_replace(vec, "[0-9]+", "")
# vec <- str_replace(vec, "  ", "")
# 
# KEGG_pathways_Metabolism_anno <- data.frame(cat1  = cate1,
#                                        cat2  = cate3,
#                                        vec  = vec)
# 
# rownames(KEGG_pathways_Metabolism_anno) <- KEGG_pathways_Metabolism_anno$vec

#================================================================================
#================================================================================
#================================================================================
library(stringr)
library(reshape2)
library(scales)
library(scater)
library(pheatmap)
library(ggplot2)
library(dplyr)
library(ggrepel)
library(RColorBrewer)

#1. Loading the data
sce <- load("~/data_analysis/单细胞代谢分析/sce.RData")
# library(fgsea)#读入KEGG gmt
# pathway_file <- "./KEGG_metabolism_nc.gmt"
# pathways <- gmtPathways(pathway_file)
# pathway_names <- names(pathways)

#data1: your celltype or sample group
all_cell_types <- as.vector(sce$celltype)
cell_types <- unique(all_cell_types)
#data2 metabolism pathway gene set
pathway_file <- "./KEGG_metabolism_nc.gmt"
pathways <- gmtPathways(pathway_file)


Pathway_act_Score <- function(data,
                              pathways,
                              assay,
                              filterGene=F,
                              Mean_cut=NULL,
                              percent_cut=NULL,
                              all_cell_types,
                              cell_types
                              ){
  
  
  DefaultAssay(data) <- assay
  norm_tpm <- GetAssayData(data, layer = "data")
  norm_tpm <- as.matrix(norm_tpm)
  norm_tpm <- as.data.frame(norm_tpm)
  
  ###filter genes
  if(filterGene==F){
    
    norm_tpm = norm_tpm
    
  }else{
  
    pctExpr <- function(mat){
      nexpr <- apply(mat, 1, function(x) sum(x > 0))
      nexpr / ncol(mat)
    }
    
    norm_tpm$percent_exp = pctExpr(norm_tpm)
    norm_tpm$mean_exp = rowMeans(norm_tpm[,-ncol(norm_tpm)])
    
    norm_tpm <- norm_tpm[which(norm_tpm$mean_exp >Mean_cut & norm_tpm$percent_exp >percent_cut),]
    norm_tpm <- norm_tpm[,-ncol(norm_tpm)]
    norm_tpm <- norm_tpm[,-ncol(norm_tpm)]
  }
  
  
  ##calculate how many pathways of one gene involved.
  num_of_pathways <- function (pathway_file,overlapgenes){
    pathway_names <- names(pathway_file)
    filter_pathways <- list()
    for (p in pathway_names){
      genes <- pathway_file[[p]]
      common_genes <- intersect(genes,overlapgenes)
      if(length(common_genes>=5)){
        filter_pathways[[p]] <- common_genes
      }
    }
    
    all_genes <- unique(as.vector(unlist(filter_pathways)))
    gene_times <- data.frame(num =rep(0,length(all_genes)),row.names = all_genes)
    for(p in pathway_names){
      for(g in filter_pathways[[p]]){
        gene_times[g,"num"] = gene_times[g,"num"]+1
      }
    }
    gene_times
  } 
  
  
  #some genes occur in multiple pathways.
  gene_pathway_number <- num_of_pathways(pathways,rownames(norm_tpm))
  
  ##Calculate the pathway activities
  #mean ratio of genes in each pathway for each cell type
  pathway_names <- names(pathways)
  mean_expression_shuffle <- matrix(NA,nrow=length(pathway_names),ncol=length(cell_types),dimnames = list(pathway_names,cell_types))
  mean_expression_noshuffle <- matrix(NA,nrow=length(pathway_names),ncol=length(cell_types),dimnames = list(pathway_names,cell_types))
  ###calculate the pvalues using shuffle method
  pvalues_mat <- matrix(NA,nrow=length(pathway_names),ncol=length(cell_types),dimnames = (list(pathway_names, cell_types)))
  
  
  for(p in pathway_names){
    genes <- pathways[[p]]
    genes_comm <- intersect(genes, rownames(norm_tpm))
    if(length(genes_comm) < 5) next
    
    pathway_metabolic_tpm <- norm_tpm[genes_comm, ]
    pathway_metabolic_tpm <- pathway_metabolic_tpm[rowSums(pathway_metabolic_tpm)>0,]
    
    mean_exp_eachCellType <- apply(pathway_metabolic_tpm, 1, function(x)by(x, all_cell_types, mean))
    
    #remove genes which are zeros in any celltype to avoid extreme ratio value
    keep <- colnames(mean_exp_eachCellType)[colAlls(mean_exp_eachCellType>0.001)]
    
    if(length(keep)<3) next
    
    #using the loweset value to replace zeros for avoiding extreme ratio value
    pathway_metabolic_tpm <- pathway_metabolic_tpm[keep,]
    pathway_metabolic_tpm <- t( apply(pathway_metabolic_tpm,1,function(x) {x[x<=0] <- min(x[x>0]);x} ))
    
    
    pathway_number_weight = 1 / gene_pathway_number[keep,]
    #
    mean_exp_eachCellType <- apply(pathway_metabolic_tpm, 1, function(x)by(x, all_cell_types, mean))
    ratio_exp_eachCellType <- t(mean_exp_eachCellType) / colMeans(mean_exp_eachCellType)
    #exclude the extreme ratios
    col_quantile <- apply(ratio_exp_eachCellType,2,function(x) quantile(x,na.rm=T))
    col_q1 <- col_quantile["25%",]
    col_q3 <- col_quantile["75%",]
    col_upper <- col_q3 * 3
    col_lower <- col_q1 / 3
    outliers <- apply(ratio_exp_eachCellType,1,function(x) {any( (x>col_upper)|(x<col_lower) )} )
    
    if(sum(!outliers) < 3) next
    
    keep <- names(outliers)[!outliers]
    pathway_metabolic_tpm <- pathway_metabolic_tpm[keep,]
    pathway_number_weight = 1 / gene_pathway_number[keep,]
    mean_exp_eachCellType <- apply(pathway_metabolic_tpm, 1, function(x)by(x, all_cell_types, mean))
    ratio_exp_eachCellType <- t(mean_exp_eachCellType) / colMeans(mean_exp_eachCellType)
    mean_exp_pathway <- apply(ratio_exp_eachCellType,2, function(x) weighted.mean(x, pathway_number_weight/sum(pathway_number_weight)))
    mean_expression_shuffle[p, ] <-  mean_exp_pathway[cell_types]
    mean_expression_noshuffle[p, ] <-  mean_exp_pathway[cell_types]
    
    ##shuffle 5000 times:  
    ##define the functions 
    group_mean <- function(x){
      sapply(cell_types,function(y) rowMeans(pathway_metabolic_tpm[,shuffle_cell_types_list[[x]]==y,drop=F]))
    }
    column_weigth_mean <- function(x){
      apply(ratio_exp_eachCellType_list[[x]],2, function(y) weighted.mean(y, weight_values))
    }
    #####  
    times <- 1:5000
    weight_values <- pathway_number_weight/sum(pathway_number_weight)
    shuffle_cell_types_list <- lapply(times,function(x) sample(all_cell_types)) 
    names(shuffle_cell_types_list) <- times
    mean_exp_eachCellType_list <- lapply(times,function(x) group_mean(x))
    ratio_exp_eachCellType_list <- lapply(times,function(x) mean_exp_eachCellType_list[[x]] / rowMeans(mean_exp_eachCellType_list[[x]]))
    mean_exp_pathway_list <- lapply(times,function(x) column_weigth_mean(x))
    
    shuffle_results <- matrix(unlist(mean_exp_pathway_list),ncol=length(cell_types),byrow = T) 
    rownames(shuffle_results) <- times
    colnames(shuffle_results) <- cell_types
    for(c in cell_types){
      if(is.na(mean_expression_shuffle[p,c])) next
      if(mean_expression_shuffle[p,c]>1){
        pval <- sum(shuffle_results[,c] > mean_expression_shuffle[p,c]) / 5000 
      }else if(mean_expression_shuffle[p,c]<1){
        pval <- sum(shuffle_results[,c] < mean_expression_shuffle[p,c]) / 5000
      }
      if(pval>0.01) mean_expression_shuffle[p, c] <- NA  ### NA is  blank in heatmap
      pvalues_mat[p,c] <- pval
    }
  }
  
  
  analysis_res <- list(mean_expression_shuffle,
                       mean_expression_noshuffle,
                       pvalues_mat)
  
  names(analysis_res) <- c("mean_expression_shuffle",
                           "mean_expression_noshuffle",
                           "pvalues_mat")
  return(analysis_res)
  
}


#================================================================================
#analysis------------------------------------------------------------------------
sce$newcelltype <- paste0(sce$orig.ident,"_", sce$celltype)
#data1: your celltype or sample group
all_cell_types <- as.vector(sce$newcelltype)
cell_types <- unique(all_cell_types)
#data2 metabolism pathway gene set
pathway_file <- "./KEGG_metabolism_nc.gmt"
pathways <- gmtPathways(pathway_file)


metabolism_activaty <- Pathway_act_Score(sce,
                                         pathways=pathways,
                                         assay = "RNA",
                                         filterGene=T,
                                         Mean_cut = 0.001,
                                         percent_cut =0.1,
                                         all_cell_types = all_cell_types,
                                         cell_types = cell_types)


#plot1----------------------------------------------------------------------------
#通路活性数据就是mean_expression_shuffle
pt <- metabolism_activaty[[1]]
#去除所有列都是NA的数据，为什么会有NA，原因是我们前面筛选表达矩阵了，那么
#有些通路的基因在我们的数据中不存在，所以评分就会是NA呀！我们去除这样的数据
all_NA <- rowAlls(is.na(pt))
pt <- pt[!all_NA,]


#然后就可以进行可视化了，首先使用ggplot进行作图
dat <- pt
#以下内容是原作者提供的数据排序，因为我们要添加注释，所以排序这里对我们没用
#不过感兴趣的还是可以学习以下这个代码，应用在其他的内容中
# sort_row <- c()
# sort_column <- c()
# 
# for(i in colnames(dat)){
#   select_row <- which(rowMaxs(dat,na.rm = T) == dat[,i])
#   tmp <- rownames(dat)[select_row][order(dat[select_row,i],decreasing = T)]
#   sort_row <- c(sort_row,tmp)
# }
# sort_column <- apply(dat[sort_row,],2,function(x) order(x)[nrow(dat)])
# sort_column <- names(sort_column)
# dat[is.na(dat)] <- 1
# dat <- as.data.frame(dat)

#读入KEGG annotation文件，挑选代谢通路的注释
KEGG_anno <- read.csv("kegg_annotation.csv", header = T)
rownames(KEGG_anno) <- KEGG_anno$directory3
Anno_pathway <- KEGG_anno[rownames(dat),]
dat$pathway_anno <- Anno_pathway$directory2


#宽数据转化为长数据
dat$pathway <- rownames(dat)
melt_dat = melt(dat,
                id.vars = c("pathway","pathway_anno"),
                measure.vars = 1:(length(dat)-2),
                variable.name = c('celltype'),
                value.name = 'score')
#按照分组排个序
melt_dat <- melt_dat %>%arrange(pathway_anno, desc(score))
melt_dat$celltype <- factor(melt_dat$celltype, levels = c("HC_UEC","HC_CEC","HC_EC","HC_SMC","HC_SF","HC_MAC","HC_LY",
                                                          "AEH_UEC","AEH_CEC","AEH_EC","AEH_SMC","AEH_SF","AEH_MAC","AEH_LY"))
melt_dat$pathway <- factor(melt_dat$pathway, levels = unique(melt_dat$pathway))

p = ggplot(melt_dat,aes(x=celltype,y=pathway, fill=score))+
  geom_tile(color='black')+
  theme(panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title = element_blank(),
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_text(angle = 90,hjust=1,vjust=0.5,color = "black", size=8),
        axis.text.y = element_text(size=8,vjust=1,color = "black"),
        legend.margin = margin(-0.2,-0.2,0,0,'cm'))+
  theme(panel.border = element_rect(fill=NA,color="black", linewidth=0.5, linetype="solid"))+
  scale_y_discrete(position = 'right',expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0))+
  scale_fill_gradientn(colors=colorRampPalette(c("#1AA3FF","white","#FF6B67"))(100),
                       guide = guide_colorbar(ticks.colour = "black",
                                              frame.colour = "black"),
                       name = "Pathway activity score",
                       values = c(0,0.15,0.3,0.5,1))

p


#做一个注释：
Anno_row <- rownames(dat) %>% as.data.frame() %>%
  mutate(group=dat$pathway_anno) %>%
  mutate(p="")
colnames(Anno_row)[1] <- 'pathway'
Anno_row$pathway <- factor(Anno_row$pathway, levels = Anno_row$pathway)
Anno_rows <- ggplot(Anno_row, aes(p,pathway,fill=group))+
  geom_tile() + 
  scale_y_discrete(position="right") +
  theme_minimal()+xlab(NULL) + ylab(NULL) +
  theme(axis.text = element_blank())+
  theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))+
  scale_fill_manual(values = c('#855C75', '#D9AF6B', '#AF6458', '#736F4C', '#526A83', '#625377', '#68855C','#9C9C5E', '#8C785D', '#467378'))+
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0))+
  guides(fill=guide_legend(title="Metabolic pathway"))


library(patchwork)
plot <- Anno_rows+p+plot_layout(ncol  = 2, widths  = c(0.2,5),guides = 'collect')
plot & theme(plot.margin = margin(0.5,0.5,0.5,0.5))



#plot2----------------------------------------------------------------------------
#pheatmap

dat1 = dat[,1:14]
dat1 = dat1[,c("HC_UEC","HC_CEC","HC_EC","HC_SMC","HC_SF","HC_MAC","HC_LY","AEH_UEC","AEH_CEC","AEH_EC","AEH_SMC","AEH_SF","AEH_MAC","AEH_LY")]

mybreaks <- c(
  seq(0.5, 1, length.out=50),
  seq(1.1, max(dat1),length.out=50)
) 
color <- colorRampPalette(c("#1AA3FF","white","#FF6B67"))(100)
pheatmap(dat1,
         cluster_cols = F,cluster_rows = T,
         color=color, breaks = mybreaks, border_color = 'black')

#plot3----------------------------------------------------------------------------
#box plot

#boxplot show the distribution of pathway activity
scRNA_dat <- as.data.frame(metabolism_activaty[[2]])
scRNA_dat$X <- NULL

scRNA_df <- melt(scRNA_dat)
scRNA_df <- scRNA_df[!is.na(scRNA_df$value),]
g <- ggplot(scRNA_df,aes(x=variable,y=value,fill=variable)) +
  # scale_y_continuous(limits=c(0,2),breaks=0:2,labels=0:2)+
  geom_violin(trim=F,size=0.2,show.legend = F,width=1.0) + labs(y=NULL,x=NULL) + 
  stat_summary(fun.y = median,geom="point",size=1,color="blue")+
  scale_fill_brewer(palette="Set5") +
  theme_classic() + 
  theme(legend.position="none",
        axis.text.x=element_text(colour="black", size = 6,angle=45,hjust=1,vjust=1),
        axis.text.y=element_text(colour="black", size = 6),
        axis.line=element_line(size=0.2,color="black"),
        axis.ticks = element_line(colour = "black",size=0.2),
        panel.border = element_blank(), panel.background = element_blank(),
        axis.ticks.length= unit(.5, "mm"))+
  geom_hline(yintercept = 1)+
  labs(y="Pathway activity score")


##小提琴图
library(ggpubr)
CEC <- scRNA_dat[,c("HC_CEC","AEH_CEC")]
CEC <- na.omit(CEC)
CEC$X <- NULL
CEC <- melt(CEC)
colnames(CEC)
# [1] "variable" "value"  
gg1 <- ggplot(data=CEC, aes(x = variable, y = value,
                    fill=variable)) +
  geom_boxplot(alpha =0.5,size=0.5,outlier.shape = NA)+
  scale_fill_manual(values=c("#E29827","#922927"))+
  stat_compare_means(method = "t.test",paired = F, 
                     comparisons=list(c("HC_CEC", "AEH_CEC")))+
  geom_jitter(alpha = 0.3,size=3, shape=21)+
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))+
  theme_bw() + 
  theme(panel.grid =element_blank(),
        axis.text = element_text(size = 10,colour = "black"),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = 'none')+
  labs(y='Pathway activity score')


gg2 <- ggplot(data=CEC, aes(x = variable, y = value,
                            fill=variable)) +
  geom_violin(alpha =0.5,size=0.5,outlier.shape = NA)+
  geom_boxplot(alpha =0.5,size=0.5,outlier.shape = NA, width=0.2, fill="white",color="black")+
  scale_fill_manual(values=c("#E29827","#922927"))+
  stat_compare_means(method = "t.test",paired = F, 
                     comparisons=list(c("HC_CEC", "AEH_CEC")))+
  geom_jitter(alpha = 0.3,size=3, shape=21)+
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))+
  theme_bw() + 
  theme(panel.grid =element_blank(),
        axis.text = element_text(size = 10,colour = "black"),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = 'none')+
  labs(y='Pathway activity score')


gg1+gg2


#######################################################################################
#######################################################################################
#######################################################################################

#筛选
norm_tpm <- GetAssayData(sce, layer = "data")
norm_tpm <- as.matrix(norm_tpm)
norm_tpm <- as.data.frame(norm_tpm)
pctExpr <- function(mat){
  nexpr <- apply(mat, 1, function(x) sum(x > 0))
  nexpr / ncol(mat)
}
norm_tpm$percent_exp = pctExpr(norm_tpm)
norm_tpm$mean_exp = rowMeans(norm_tpm[,-ncol(norm_tpm)])

norm_tpm <- norm_tpm[which(norm_tpm$mean_exp >0.001 & norm_tpm$percent_exp >0.1),]
norm_tpm <- norm_tpm[,-ncol(norm_tpm)]
norm_tpm <- norm_tpm[,-ncol(norm_tpm)]

sce_sec <- sce[rownames(norm_tpm),colnames(norm_tpm)]
dim(sce_sec)
# [1]  8354 15881
dim(sce)
# [1] 27001 15881

save(sce_sec, file = 'sce_sec.RData')

library(reshape2)
library(pheatmap)
library(RColorBrewer)
library(ggtree)
library(aplot)
library(dplyr)
library(ggplot2)
library(ggheatmap)



ks_gg_heatmap <- function(mat,#矩阵，类似于行是基因，列是分组样本
                          scaleData,#是否需要scale data，不需要“none", 按照行"row", 按照列“column"
                          cluster_method,#聚类方法，参考hclust
                          cluster_row=T,#是否聚类，默认T
                          cluster_col=T,#是否聚类，默认T
                          cols,#热图颜色
                          show_cluster_rows=T,#是否展示聚类树，默认T
                          show_cluster_cols=T,#是否展示聚类树，默认T
                          text_size,#x,y文字大小
                          legend_title,#热图legend标题
                          Discrete_data=F,#输入数据类型，是连续性还是离散型
                          rowname_pos="right",#rowname的位置，默认右侧
                          shape="S",#热图形状，默认为S，S表示矩形，C表示⚪
                          myorder=F,#是否自行设定行列顺序，默认F
                          level_row=NULL,#myorder=T时，行的排序
                          level_col=NULL,#myorder=T时，列的排序
                          anno_col=NULL,#列注释
                          anno_row=NULL,#行注释
                          anno_col_color=NULL,#列注释颜色
                          anno_row_color=NULL,#行注释颜色
                          anno_col_title=NULL,#列注释标题
                          anno_row_title=NULL,#行注释标题
                          label_color=F,#是否标记行名特殊颜色，默认F
                          label1=NULL,#label_color=T时，标记的第一组行名
                          label2=NULL,#label_color=T时，标记的第2组行名
                          label3=NULL,#label_color=T时，标记的第3组行名
                          values
                          
){
  
  scale_rows = function(x){
    m = apply(x, 1, mean, na.rm = T)
    s = apply(x, 1, sd, na.rm = T)
    return((x - m) / s)}
  
  df = switch(scaleData, 
              none = mat, 
              row = scale_rows(mat), 
              column = t(scale_rows(t(mat))))  
  
  
  if(cluster_row==T){
    
    row_clust = hclust(dist(df), method = cluster_method)
    roword = row_clust$order
    
    if(show_cluster_rows==T){
      tree <- ape::as.phylo(row_clust)
      row_tree <- ggtree(tree)+
        theme(legend.position = "none")
    }else{
      
      row_tree=NULL
      
    }
    
    
  }else{
    
    roword=seq(1:nrow(df))
    
  }
  
  if(cluster_col==T){
    
    row_clust = hclust(dist(t(df)), method = cluster_method)
    colord = row_clust$order
    
    if(show_cluster_cols==T){
      tree <- ape::as.phylo(row_clust)
      col_tree <- ggtree(tree,layout = "dendrogram")+
        theme(legend.position = "none")
    }else{
      
      col_tree=NULL
      
    }
    
    
  }else{
    
    colord = seq(1:length(df))
    
  }
  
  
  
  dat=df[roword,colord]
  
  dat <- cbind(rownames(dat),dat)
  colnames(dat)[1] <- 'gene'
  
  melt_dat = melt(dat,
                  id.vars = c("gene"),
                  measure.vars = 2:length(dat),
                  variable.name = c('group'),
                  value.name = 'Exp')
  
  
  if(myorder==F){
    
    melt_dat$gene <- factor(melt_dat$gene, levels = rownames(dat))
    melt_dat$group <- factor(melt_dat$group, levels = colnames(dat))
    
  }else{
    
    melt_dat$gene <- factor(melt_dat$gene, levels = level_row)
    melt_dat$group <- factor(melt_dat$group, levels = level_col)
    
  }
  
  
  if(label_color == T){
    
    melt_dat$color <-   ifelse(melt_dat$gene %in% label1, "#D52126",
                               ifelse(melt_dat$gene %in% label2, "#117733",
                                      ifelse(melt_dat$gene %in% label3,"#E68316","black")))
    
    
  }else{
    
    melt_dat$color <- 'black'
  }
  
  
  
  
  if(Discrete_data==F){
    if(shape=="S"){
      p = ggplot(melt_dat,aes(x=group,y=gene, fill=Exp))+
        geom_tile(color='black')+
        theme(panel.background = element_blank(),
              panel.border = element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              axis.title = element_blank(),
              axis.line = element_blank(),
              axis.ticks = element_blank(),
              axis.text.x = element_text(angle = 90,hjust=1,vjust=0.5,color = "black", size=text_size),
              axis.text.y = element_text(size=text_size,vjust=1,
                                         colour =melt_dat$color),
              legend.margin = margin(-0.2,-0.2,0,0,'cm'),
              legend.key.height = unit(0.5,'cm'),
              legend.key.width = unit(0.3,'cm'))+
        scale_y_discrete(position = rowname_pos,expand = c(0,0)) + 
        scale_x_discrete(expand = c(0,0))+
        scale_fill_gradientn(colors=cols,
                             guide = guide_colorbar(ticks.colour = "black",
                                                    frame.colour = "black"),
                             name = legend_title,
                             values = values)
      
    }
    
    if(shape=="C"){
      p = ggplot(melt_dat,aes(x=group,y=gene))+
        geom_point(aes(size=abs(Exp), color=Exp))+
        theme(panel.background = element_blank(),
              panel.border = element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              axis.title = element_blank(),
              axis.line = element_blank(),
              axis.ticks = element_blank(),
              axis.text.x = element_text(angle = 90,hjust=1,vjust=0.5,color = "black", size=text_size),
              axis.text.y = element_text(size=text_size,hjust=1,
                                         colour =melt_dat$color),
              legend.margin = margin(-0.2,-0.2,0,0,'cm'),
              legend.key.height = unit(0.5,'cm'),
              legend.key.width = unit(0.3,'cm'))+
        scale_y_discrete(position = rowname_pos, expand = c(0,0)) + 
        scale_x_discrete(expand = c(0,0))+
        scale_color_gradientn(colors=cols,
                              guide = guide_colorbar(ticks.colour = "black",
                                                     frame.colour = "black"),
                              name = legend_title,
                              values = values)
    }
    
    
    
  }else{
    
    melt_dat$Exp <- as.factor(melt_dat$Exp)
    if(shape=="S"){
      p = ggplot(melt_dat,aes(x=group,y=gene,fill=Exp))+
        geom_tile(color='black')+
        theme(panel.background = element_blank(),
              panel.border = element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              axis.title = element_blank(),
              axis.line = element_blank(),
              axis.ticks = element_blank(),
              axis.text.x = element_text(angle = 90,hjust=1,vjust=0.5,color = "black", size=text_size),
              axis.text.y = element_text(size=text_size,hjust=1,
                                         colour =melt_dat$color),
              legend.margin = margin(-0.2,-0.2,0,0,'cm'),
              legend.key.height = unit(0.5,'cm'),
              legend.key.width = unit(0.3,'cm'))+
        scale_y_discrete(position = rowname_pos,expand = c(0,0)) + 
        scale_x_discrete(expand = c(0,0))+
        scale_fill_manual(values =cols,
                          name = legend_title)
      
    }
    
    if(shape=="C"){
      p = ggplot(melt_dat,aes(x=group,y=gene,color=Exp))+
        geom_point(aes(size=5))+
        theme(panel.background = element_blank(),
              panel.border = element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              axis.title = element_blank(),
              axis.line = element_blank(),
              axis.ticks = element_blank(),
              axis.text.x = element_text(angle = 90,hjust=1,vjust=0.5,color = "black", size=text_size),
              axis.text.y = element_text(size=text_size,vjust=1,
                                         colour =melt_dat$color),
              legend.margin = margin(-0.2,-0.2,0,0,'cm'),
              legend.key.height = unit(0.5,'cm'),
              legend.key.width = unit(0.3,'cm'))+
        scale_y_discrete(position = rowname_pos,expand = c(0,0)) + 
        scale_x_discrete(expand = c(0,0))+
        scale_color_manual(values =cols,
                           name = legend_title)
      
    }
    
    
  }
  # 
  # # plot =p%>%insert_top(col_tree,height = 0.1)
  # # plot = ggplotify::as.ggplot(plot)
  # plot =row_tree+p+plot_layout(ncol = 2, widths  = c(0.5,4),guides = 'collect')
  # plot = plot & theme(plot.margin = margin(0.8,0.8,0.8,0.8))
  
  if(!is.null(anno_col)){
    
    Anno_col <- colnames(mat) %>% as.data.frame() %>% 
      mutate(group=anno_col) %>%
      mutate(p="") 
    Anno_cols <- ggplot(Anno_col, aes(.,y=p,fill=group))+
      geom_tile() + 
      scale_y_discrete(position="right") +
      theme_minimal()+xlab(NULL) + ylab(NULL) +
      theme(axis.text = element_blank())+
      theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))+
      scale_fill_manual(values = anno_col_color)+
      scale_y_discrete(expand = c(0,0)) + 
      scale_x_discrete(expand = c(0,0))+
      guides(fill=guide_legend(title=anno_col_title))
    
  }
  
  
  
  
  if(!is.null(anno_row)){
    Anno_row <- rownames(mat) %>% as.data.frame() %>%
      mutate(group=anno_row) %>%
      mutate(p="")
    Anno_rows <- ggplot(Anno_row, aes(p,.,fill=group))+
      geom_tile() + 
      scale_y_discrete(position="right") +
      theme_minimal()+xlab(NULL) + ylab(NULL) +
      theme(axis.text = element_blank())+
      theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))+
      scale_fill_manual(values = anno_row_color)+
      scale_y_discrete(expand = c(0,0)) + 
      scale_x_discrete(expand = c(0,0))+
      guides(fill=guide_legend(title=anno_row_title))
    
  }
  
  
  # 
  # 
  # if(show_cluster_cols==F & show_cluster_cols==F){
  #   
  #   p <- Anno_cols+p+plot_layout(nrow = 2, heights = c(0.15,5),guides = 'collect')
  #   p <- Anno_rows+p+plot_layout(ncol  = 2, widths  = c(0.15,5),guides = 'collect')
  #   p = p & theme(plot.margin = margin(0.5,0.5,0.5,0.5))
  #   
  # }else{
  
  #layout
  if(show_cluster_rows==F&show_cluster_cols==T){
    if(!is.null(anno_col) &!is.null(anno_row)){
      
      p <-  p%>%insert_left(Anno_rows,width = 0.05)%>%
        insert_top(Anno_cols,height = 0.02)%>%
        insert_top(col_tree,height = 0.1)
    }
    
    if(!is.null(anno_col) & is.null(anno_row)){
      
      p <-  p%>%insert_top(Anno_cols,height = 0.02)%>%
        insert_top(col_tree,height = 0.1)
    }
    
    
    if(is.null(anno_col) &!is.null(anno_row)){
      
      p <-  p%>%insert_left(Anno_rows,width = 0.05)%>%
        insert_top(col_tree,height = 0.1)
    }
    
    
    if(is.null(anno_col) &is.null(anno_row)){
      
      p <-  p%>%insert_top(col_tree,height = 0.1)
    }
    
    
    
  }
  
  
  
  if(show_cluster_rows==T&show_cluster_cols==F){
    if(!is.null(anno_col) &!is.null(anno_row)){
      p <-  p%>%insert_left(Anno_rows,width = 0.05)%>%
        insert_left(row_tree,width = 0.1)%>%
        insert_top(Anno_cols,height = 0.02)
    }
    
    if(!is.null(anno_col) & is.null(anno_row)){
      p <-  p%>%insert_left(row_tree,width = 0.1)%>%
        insert_top(Anno_cols,height = 0.02)
    }
    
    
    if(is.null(anno_col) &!is.null(anno_row)){
      
      p <-  p%>%insert_left(Anno_rows,width = 0.05)%>%
        insert_left(row_tree,width = 0.1)
    }
    
    
    if(is.null(anno_col) &is.null(anno_row)){
      
      p <-  p%>%insert_left(row_tree,width = 0.1)
    }
    
    
    
  }
  
  
  if(show_cluster_rows==F&show_cluster_cols==F){
    
    if(!is.null(anno_col) &!is.null(anno_row)){
      p <-  p%>%insert_left(Anno_rows,width = 0.05)%>%
        insert_top(Anno_cols,height = 0.02)
    }
    
    if(!is.null(anno_col) & is.null(anno_row)){
      p <-  p%>%insert_top(Anno_cols,height = 0.02)
    }
    
    
    if(is.null(anno_col) &!is.null(anno_row)){
      
      p <-  p%>%insert_left(Anno_rows,width = 0.05)
    }
    
    
    if(is.null(anno_col) &is.null(anno_row)){
      
      p <-  p
    }
    
    
  }
  
  
  
  if(show_cluster_rows==T&show_cluster_cols==T){
    if(!is.null(anno_col) &!is.null(anno_row)){
      
      p <-  p%>%insert_left(Anno_rows,width = 0.05)%>%
        insert_left(row_tree,width = 0.1)%>%
        insert_top(Anno_cols,height = 0.02)%>%
        insert_top(col_tree,height = 0.1)
    }
    
    if(!is.null(anno_col) & is.null(anno_row)){
      
      p <-  p%>%insert_left(row_tree,width = 0.1)%>%
        insert_top(Anno_cols,height = 0.02)%>%
        insert_top(col_tree,height = 0.1)
    }
    
    
    if(is.null(anno_col) &!is.null(anno_row)){
      p <-  p%>%insert_left(Anno_rows,width = 0.05)%>%
        insert_left(row_tree,width = 0.1)%>%
        insert_top(col_tree,height = 0.1)
    }
    
    
    if(is.null(anno_col) &is.null(anno_row)){
      
      p <-  p%>%insert_left(row_tree,width = 0.1)%>%
        insert_top(col_tree,height = 0.1)
    }
    
    
  }
  
  
  
  return(p)
  
}


ks_gg_heatmap(dat[,1:14], scaleData="none",
              # cluster_method='ward.D2',
              cluster_row = F, 
              cluster_col = F,
              show_cluster_rows = F,
              show_cluster_cols =F,
              anno_col = c(rep("HC",7),rep("AEH",7)),
              text_size = 8,
              legend_title = "Pathway activity score",
              cols = colorRampPalette(c("#1AA3FF","white","#FF6B67"))(100),
              rowname_pos = "right",
              shape="S",
              anno_col_color  = c("#ABDEB6", "#60C1BF"),
              anno_row = Anno_pathway$directory2,
              anno_row_color = dittoColors(),
              anno_col_title = "group",
              anno_row_title = "Metabolic pathway",
              label_color = F,
              myorder = T,
              level_col =  c("HC_UEC","HC_CEC","HC_EC","HC_SMC","HC_SF","HC_MAC","HC_LY",
                             "AEH_UEC","AEH_CEC","AEH_EC","AEH_SMC","AEH_SF","AEH_MAC","AEH_LY"),
              level_row = rownames(dat),
              values = c(0,0.15,0.3,0.5,1))