近来处理了一个22w个细胞的单细胞数据，数据过多，做GSVA直接卡住好几天，本来尝试每200个细胞这样计算然后拼接起来，函数也写好了，跑完后突然发现对于GSVA不能分割矩阵来计算，因为GSVA有一个估计核密度函数的过程，需要某个基因的所有样本，因此不能分割样本。具体原理参考前一篇文章：https://www.jianshu.com/p/47f9f0107ffe

library(reshape2)
library(ggpubr)
library(ggplot2)
library(GSVA)

ssgsea_by_steps <- function(object,by=200,genelist,group="",plot=T,test_method="wilcox.test"){
  results <- list()
  gsva_temp_list <- list()
  N = ncol(object)
  n=floor(N/by)
  if(n>0){
    pb_ssGSEA <- progress_bar$new(
      format = '完成 [:bar] :percent eta: :eta',
      total = n, clear = FALSE, width = 80
    )#进度条
    
    for(i in 1:n){
      exp_temp <- as.matrix(as.data.frame(object[,(by*(i-1)+1):(by*i)][["RNA"]]@counts))
      gsva_temp_list[[i]] <- gsva(exp_temp, genelist, method="ssgsea")
      
      pb_ssGSEA$tick()
      Sys.sleep(0.05)#进度条
    }
    if(floor(N/by)!=N/by){
      gsva_temp_list[[n+1]] <- gsva(as.matrix(as.data.frame(object[,(by*n+1):ncol(object)][["RNA"]]@counts)),
                                    genelist, method="ssgsea")}
    
    #合并所有的的ssgsea分数
    gsva_res <- t(gsva_temp_list[[1]])
    pb_ssGSEA <- progress_bar$new(
      format = '完成 [:bar] :percent eta: :eta',
      total = length(gsva_temp_list), clear = FALSE, width = 80
    )#
    
    for(i in 2:length(gsva_temp_list)){
      gsva_res <- rbind(gsva_res,t(gsva_temp_list[[i]]))
      
      pb_ssGSEA$tick()
      Sys.sleep(0.05)#进度条
    }
    gsva_res <- as.data.frame(gsva_res)
    results[["ssGSEA_Score"]] <- gsva_res}
  
  #n=0
  else{gsva_res <- gsva(as.matrix(as.data.frame(object[["RNA"]]@counts)),
                        genelist, method="ssgsea")
  gsva_res <- as.data.frame(t(gsva_res))
  results[["ssGSEA_Score"]] <- gsva_res
  }
  if(!group==""){
    group_data <- as.data.frame(object@meta.data[,group],row.names=rownames(object@meta.data))
    colnames(group_data) <- group
    group_data$barcode <- rownames(group_data)
    
    gsva_res$barcode <- rownames(gsva_res)
    gsva_res <- merge(gsva_res,group_data,by="barcode",all=F)
    gsva_res <- as.data.frame(gsva_res[,-which(colnames(gsva_res)=="barcode")],row.names = gsva_res[,"barcode"])
    results[["ssGSEA_Score_Group"]] <- gsva_res
    
    #ggplot2柱状图
    if(plot){
      gsva_res_melt <- melt(gsva_res,id.vars=group)
      colnames(gsva_res_melt)[which(colnames(gsva_res_melt)==group)]="group"
      results[["Boxplot"]] <- ggplot(data=gsva_res_melt,mapping=aes(x=variable,y=value,fill=group))+geom_boxplot()+
        stat_compare_means(data=gsva_res_melt,method = test_method, label = "p.signif",cex=3)+
        theme_bw()+#背景白色
        theme(axis.text.x=element_text(vjust=1,size=8,angle = 45,hjust = 1),)+#调整横轴基因名大小
        theme(panel.grid =element_blank())+#去除网格线
        xlab("EPIGENETIC PROCESS")+ylab("GSVA Score")}
    
  }
  
  return(results)
}
#object: seurat对象
#by: 每次计算多少个细胞,默认200个
#genelist: gsva函数支持的pathway文件，list格式
#group: 细胞的分组
#plot: 是否使用ggplot和ggpubr画柱状图
#test_method: 画柱状图时的检验方法
#最后结果返回一个list，包含ssgsea score、包含分组信息的ssgsea score以及plot.



#进度条
pb2 <- progress_bar$new(
  format = '完成 [:bar] :percent eta: :eta',
  total = 22, clear = FALSE, width = 80
)

ssgsea <- list()
#按照cluster分别计算
for(i in c(1:22)){
  cluster_temp <- subset(pbmc,seurat_clusters==i)
  ssgsea[[i]] <- ssgsea_by_steps(cluster_temp,genelist = genelist,group = "disease.severity",plot = T)
  pb2$tick()
  Sys.sleep(0.05)
}

#由于seurat包命名cluster是从0开始，但是list用[]取子集的时候不能取[[0]],所以单独算一次。
ssgsea[["0"]] <- ssgsea_by_steps(subset(pbmc,seurat_clusters==0),genelist = genelist,group = "disease.severity",plot = T)