跟着顶刊学习单细胞CD8T细功能基因集评分！

# 加载必要的库
library(Seurat)
library(tidyverse)
library(pheatmap)
library(scales) # 用于rescale函数
##这里加载自己的CD8T细胞处理好的Seurat对象即可
load("CD8_Obj.Rdata")

# 设置输出路径
figurePath <- "./figures" # 输出图表的目录
dir.create(figurePath, showWarnings = FALSE, recursive = TRUE)
##将seurat_clusters设置为因子类型
CD8_Obj$seurat_clusters=as.factor(CD8_Obj$seurat_clusters)
Idents(CD8_Obj) <- CD8_Obj$seurat_clusters

# 读取标记基因 - 正确处理宽格式CSV
CD8_markers <- read_csv("CD8-markers.csv", skip = 1)  # 如果第一行是标题行，跳过它

# 将宽格式转换为标记基因列表
marker.list <- list()
names(marker.list[1])
for(col_name in colnames(CD8_markers)) {
  # 提取非NA的基因
  genes <- CD8_markers[[col_name]]
  genes <- genes[!is.na(genes) & genes != ""]
  
  # 替换列名中的特殊字符
  clean_name <- gsub("/", ":", col_name)
  
  # 添加到列表
  if(length(genes) > 0) {
    marker.list[[clean_name]] <- genes
  }
}

CD8_Obj <-  AddModuleScore(CD8_Obj,
                           features = marker.list,
                           ctrl = 5,
                           name = "FunctionScore")
for(i in 1:length(marker.list)){
  colnames(CD8_Obj@meta.data)[colnames(CD8_Obj@meta.data) == paste0("FunctionScore", i)] <- names(marker.list)[i]
}
Idents(CD8_Obj) <- CD8_Obj$seurat_clusters
Differentiation <- c("Naive", "Activation:Effector function", "Exhaustion")
Function <- c("TCR Signaling", "Cytotoxicity", "Cytokine:Cytokine receptor",
              "Chemokine:Chemokine receptor", "Senescence", "Anergy",
              "NFKB Signaling", "Stress response", "MAPK Signaling", "Adhesion",
              "IFN Response")
Metabolism <- c("Oxidative phosphorylation", "Glycolysis", "Fatty acid metabolism")
Apoptosis <- c("Pro-apoptosis", "Anti-apoptosis")
MarkerNameVector <- c(Differentiation, Function, Metabolism, Apoptosis)
FunctionScoreMatrix <- matrix(0,
                              ncol = length(unique(CD8_Obj$seurat_clusters)),
                              nrow = length(marker.list))
colnames(FunctionScoreMatrix) <- paste0("CD8_", 0:10)
rownames(FunctionScoreMatrix) <- MarkerNameVector
for(ci in 1:ncol(FunctionScoreMatrix)){
  for(ri in 1:nrow(FunctionScoreMatrix)){
    FunctionVec <- as_tibble(CD8_Obj@meta.data) %>% pull(MarkerNameVector[ri])
    fv <- mean(FunctionVec[CD8_Obj$seurat_clusters == levels(CD8_Obj$seurat_clusters)[ci]])
    FunctionScoreMatrix[ri, ci] <- fv
  }
}
FunctionScoreMatrix <- t(apply(FunctionScoreMatrix, 1, rescale, to=c(-1, 1)))
orderC =c(0,1,2,3,4,5,6,7,8,9,10)  # 指定因子水平的顺序 
orderC_with_prefix = paste0("CD8_", orderC)  # 将 "HCC_" 前缀添加到每个元素
FunctionScoreMatrix <- FunctionScoreMatrix[,orderC_with_prefix]

my.breaks <- c(seq(-1, 0, by=0.1), seq(0.1, 1, by=0.1))
my.colors <- c(
  colorRampPalette(colors = c("#6DCCFD", "white"))(length(my.breaks)/2),
  colorRampPalette(colors = c("white", "#FD9AA0"))(length(my.breaks)/2))
## cellType_col <- data.frame(cell.type = CD8_Obj_CellType)
## rownames(cellType_col) <- colnames(FunctionScoreMatrix)
signatureType_row <- data.frame(Signature.type = c(
  rep("Differentiation", length(Differentiation)),
  rep("Function", length(Function)),
  rep("Metabolism", length(Metabolism)),
  rep("Apoptosis", length(Apoptosis))))
rownames(signatureType_row) <- MarkerNameVector
pheatmap(FunctionScoreMatrix,
         show_colnames = T,
         show_rownames = T,
         ## annotation_col = cellType_col,
         annotation_row = signatureType_row,
         gaps_row = c(3, 14, 17),
         cluster_rows = F,
         cluster_cols = F,
         breaks = my.breaks,
         color = my.colors,
         border_color = "NA",
         fontsize = 8,
         width = 5,
         height = 3.8,
         filename = file.path(figurePath, paste0("CD8_FunctionScore_heatmap.pdf")))