R语言ssGSEA代码

以下是基于R语言实现的ssGSEA富集分析代码：

library(GSEABase)
library(ggplot2)
library(dplyr)

# 定义ssGSEA函数
ssgsea <- function(gene_exp, gene_sets, nperm=1000, weighted_score_type=1, permutation=TRUE,
                   min_size=1, max_size=5000, verbose=FALSE, seed=NULL) {

  # 设置随机数种子
  if (!is.null(seed)) set.seed(seed)

  # 准备基因表达矩阵
  gene_exp <- as.matrix(gene_exp)

  # 准备基因集列表
  gene_sets <- gene_sets[gene_sets$Size >= min_size & gene_sets$Size <= max_size, ]

  # 计算基因集得分
  pathways <- list()
  for (i in 1:nrow(gene_sets)) {
    # 获取基因集中的基因
    genes <- gene_sets$genes[[i]]

    # 计算每个样本的基因集得分
    gss <- apply(gene_exp, 1, function(x) {
      # 获取基因集中基因的表达值
      pathway_exp <- x[which(names(x) %in% genes)]

      # 计算基因集得分
      if (weighted_score_type == 0) sum(pathway_exp)
      else if (weighted_score_type == 1) mean(pathway_exp)
      else if (weighted_score_type == -1) mean(abs(pathway_exp))
    })

    # 计算富集得分和p值
    if (permutation) {
      null_gss <- replicate(nperm, {
        # 打乱基因表达值
        shuffle_exp <- apply(gene_exp, 2, sample)

        # 计算每个样本的基因集得分
        apply(shuffle_exp, 1, function(x) {
          # 获取基因集中基因的表达值
          pathway_exp <- x[which(names(x) %in% genes)]

          # 计算基因集得分
          if (weighted_score_type == 0) sum(pathway_exp)
          else if (weighted_score_type == 1) mean(pathway_exp)
          else if (weighted_score_type == -1) mean(abs(pathway_exp))
        })
      })

      null_gss <- apply(null_gss, 2, function(x) sort(x))
      null_es <- apply(null_gss, 2, function(x) (mean(x > gss) - mean(x < gss)) / sd(x))

      es <- (mean(gss) - mean(null_es)) / sd(null_es)
      pval <- mean(null_es < mean(gss))

    } else {
      es <- (mean(gss) - mean(gss)) / sd(gss)
      pval <- 1 - pnorm(es)
    }

    pathway <- gene_sets$Pathway[i]
    pathways[[pathway]] <- list(es=es, pval=pval)

    if (verbose) {
      cat(sprintf("%s: ES = %.3f, p-value = %.3f\n", pathway, es, pval))
    }
  }

  return(pathways)
}

# 使用示例
# 准备基因表达矩阵和基因集列表
data(geneExp)
data(geneSets)

# 进行ssGSEA富集分析
pathways <- ssgsea(gene_exp=geneExp, gene_sets=geneSets)

# 将结果转换为数据框并排序
pathways_df <- lapply(pathways, function(x) {
  data.frame(Pathway=names(x), ES=x$es, pval=x$pval)
}) %>% do.call(rbind, .)
pathways_df <- pathways_df[order(pathways_df$pval), ]

# 绘制p值分布图
ggplot(pathways_df, aes(x=pval)) +
  geom_density() +
  xlab("p-value") +
  ggtitle("Distribution of ssGSEA p-values")

该代码使用了GSEABase、ggplot2和dplyr等库进行计算和绘图。在使用时，需要将基因表达矩阵和基因集列表传递给`ssgsea`函数。此外，还可以设置其他参数，例如是否进行置换和置换次数等。函数返回一个包含富集分析结果的列表。可以将结果转换为数据框并排序，然后绘制p值分布图。