
rm(list=ls())
library(GSVA)
library(GSEABase)
library(limma)
library(tidyverse)

setwd("01 raw data/")
## tpm
names = dir("E:\\Nutsore\\001-待处理\\001 文章项目处理\\胶质瘤原始数据\\01 raw data",
            pattern = "tpm")

## for recle
res.out = data.frame()
for (i in 1:5) {
  rt = read.table(paste0("E:\\Nutsore\\001-待处理\\001 文章项目处理\\胶质瘤原始数据\\01 raw data\\", names[i]),
                   sep = "\t", header = T, check.names = F)
  row.names(rt) = rt[,1]
  rt = rt[,-1]
  rt = round(rt, 2)
  dimnames=list(rownames(rt),colnames(rt))
  mat=matrix(as.numeric(as.matrix(rt)),nrow=nrow(rt),dimnames=dimnames)
  mat=avereps(mat)
  mat=mat[rowMeans(mat)>0,]
  
  ## 读取特征基因集
  gmtFile="hallmark.gmt"  
  geneSet=getGmt(gmtFile, geneIdType=SymbolIdentifier())
  
  #ssgsea
  ssgseaScore=gsva(mat, geneSet, method='ssgsea', 
                   kcdf='Gaussian',abs.ranking=TRUE)
  
  #ssGSEA score
  normalize=function(x){
    return((x-min(x))/(max(x)-min(x)))}
  
  ssgseaOut=normalize(ssgseaScore)
  ssgsea = as.data.frame(t(ssgseaOut))
  res.out = rbind(res.out, ssgsea)
}

## result out
res.out = round(res.out, 3)
res.out = tibble::rownames_to_column(res.out, var = "ID")
write.csv(res.out, file="01 hallmarks.ssGSEA.csv", row.names=F) 