生物信息数据分析教程视频——16-单样本基因集富集分析(ssGSEA)用于肿瘤相关免疫细胞浸润水平评估

# setwd("H:/MedBioInfoCloud/analysis/TCGA/new/conventionalAnalysis")
options(stringsAsFactors = F)
library(TCGAbiolinks)
library(GSVA)
library(GSEABase)
FilePath <- dir("H:/MedBioInfoCloud/analysis/TCGA/new/processedTCGAdata/TCGA-STAR_Exp/",
                "STARdata.Rdata$",full.names = T)

opt <- "output/009-CellFractions/ssGSEA/"
ifelse(dir.exists(opt),FALSE,dir.create(opt,recursive = T))

source("H:/MedBioInfoCloud/analysis/TCGA/new/00-fun/filterGeneTypeExpr.R")
source("H:/MedBioInfoCloud/analysis/TCGA/new/00-fun/del_dup_sample.R")


###TCGA数据库中33中癌症类型
project <- getGDCprojects()$project_id
project <- project[grep("TCGA-",project)]
# proj = "TCGA-LUAD"

immune_cell_geneSet = getGmt("H:/MedBioInfoCloud/analysis/base_files/immue_BG_GeneSet/ssGSEA_immue_BG_GeneSet.gmt",
                             geneIdType=SymbolIdentifier())

# devtools::install_github("BioInfoCloud/csGeneset")
# library(csGeneset)
# gsImmCell <- gsGSVA[["ImmCell"]][["geneSet"]]
# proj = "TCGA-LUAD"

for(proj in project){
  message("===============================")
  message(proj)
  load(FilePath[grep(proj,FilePath)])#STARdata
  tpm <- STARdata[["tpm"]]
  tpm <- filterGeneTypeExpr(expr = tpm,
                            fil_col = "gene_type",
                            filter = FALSE)
  ##过滤不表达的基因
  tpm <- tpm[apply(tpm,1,var) !=0,]
  ##正常组织样本ID
  SamN <- TCGAquery_SampleTypes(barcode = colnames(tpm),
                                typesample = c("NT","NB","NBC","NEBV","NBM"))
  
  ##肿瘤组织样本ID
  SamT <- setdiff(colnames(tpm),SamN)
  
  ###去除重复样本
  tur_exp <- del_dup_sample(tpm[,SamT],col_rename = T)
  # ###long2转换
  # tur_exp <- log2(tur_exp + 1)
  ssGSEA_Score = gsva(as.matrix(tur_exp), 
                      immune_cell_geneSet,
                      method='ssgsea',
                      kcdf='Gaussian',
                      abs.ranking=TRUE)#ssGSEA计算
  normalize=function(x){
    return((x-min(x))/(max(x)-min(x)))}#定义ssGSEA_Score矫正函数
  
  norm_ssGSEA_Score <- normalize(ssGSEA_Score)#对ssGSEA_Score进行矫正
 
  file <- paste0(opt,proj,"-normalize_turmor_ssGSEA_Score.txt")
  write.table(norm_ssGSEA_Score,file = file,sep="\t",quote=F,col.names=F)#此处的输出文件即为ssGSEA富集得分文件
}