R语言---生信分析---ssGSEA基因集富集分析、免疫浸润

    setwd('D:\ Desktop\ GEO分析\ GSE1991')
    rm(list=ls())

    a <- read.table('data_tpms.xls',
                header = T)
                # ,row.names=1)#读取已经下载好的补充的表达矩阵压缩文件
    a[1:4,1:4]
    
    #去掉基因名，得到纯粹的表达矩阵raw_data
    raw_data<- a[,-1]
    raw_data[1:4,1:4]
    
    ###读取表型信息
    pheno <- read.csv(file = 'GSE1991_series_matrix.txt')
    # 获取Sample_title Sample_characteristics, 标本名称、样本编号 进行处理，替换为P
    # pheno <- data.frame(num1 = strsplit(as.character(pheno[42,]),split='\t')[[1]][-1],
    #                     num2 = gsub('patient: No.','P',strsplit(as.character(pheno[51,]),split='\t')[[1]][-1]))

    ####数据过滤，把表达量和为0的基因去掉（去O)
    data<- a[!apply(raw_data,1,sum)==0,]
    ####去除重复基因名的行，归一化
    data$median=apply(data[,-1],1,median)#计算每行的中位数，添加到 data数据中
    data[1:4,1:4]
    
    data=data[order(data$gene_name, data$median,decreasing = T),]#排序
    data=data[!duplicated(data$gene_name),]#去除重复的基因名
    rownames(data)=data$gene_name#把基因名变成行名
    uni_matrix <- data[,grep('\ d+',colnames(data))]
    #grep('\ d+',colnames(data)）识别有数字的位置，然后挑出来
    uni_matrix <- log2(uni_matrix+1)#把表达矩阵log
    colnames(uni_matrix)<- gsub('DSP.','',gsub('\ .','\ -',colnames(uni_matrix)))#处理行名
    uni_matrix<- uni_matrix[,order(colnames(uni_matrix))]#根据行名排序表达矩阵
    ### 中间保存
    save(uni_matrix, file = 'uni_matrix.Rdata')##把处理好的数据存好
    write.table(uni_matrix, 'log2_tpm_uni_matrix.xls', sep='\t')
    uni_matrix[1:4,1:4]

    #======== 进行ssGSEA分析 =============
    ## 
    #if (!require("BiocManager", quietly = TRUE))
    #  install.packages("BiocManager")
    
    #BiocManager::install("GSVA")
    
    rm(list=ls())
    #加载包
    {
      library(genefilter)
      library(GSVA)
      library(Biobase)
      library(stringr)
    }
    ##载入数据
    load('uni_matrix.Rdata')

    # 读取基因列表，得到免疫细胞对应的特异的基因
    gene_set<- read.csv('../TILs_marker.csv',
                    header = T)##读取已经下载好的免疫细胞和对应基因列表，来源见文献附件
    gene_set<-gene_set[, 1:2]#选取特异基因和对应的免疫细胞两行
    head(gene_set)

    ### 读取基因列表，得到免疫细胞对应的特异的基因
    list<- split(as.matrix(gene_set)[,1], gene_set[,2])
    
    #评估Estimates GSVA enrichment scores.
    gsva_matrix<- gsva(as.matrix(uni_matrix), list,method='ssgsea',
                   kcdf='Gaussian',abs.ranking=TRUE)
    save(gsva_matrix, file = 'gsva_matrix.Rdata')##把处理好的数据存好
    write.table(gsva_matrix, 'gsva_matrix.xls', sep='\t')

    ##======   富集分数画热图   =========
    load('gsva_matrix.Rdata')
    # install.packages("pheatmap")
    library(pheatmap)
    gsva_matrix1<- t(scale(t(gsva_matrix)))#归一化
    gsva_matrix1[gsva_matrix1< -2] <- -2
    gsva_matrix1[gsva_matrix1>2] <- 2
    anti_tumor <- c('Activated CD4 T cell', 'Activated CD8 T cell', 'Central memory CD4 T cell', 'Central memory CD8 T cell', 'Effector memeory CD4 T cell', 'Effector memeory CD8 T cell', 'Type 1 T helper cell', 'Type 17 T helper cell', 'Activated dendritic cell', 'CD56bright natural killer cell', 'Natural killer cell', 'Natural killer T cell')
    pro_tumor <- c('Regulatory T cell', 'Type 2 T helper cell', 'CD56dim natural killer cell', 'Immature dendritic cell', 'Macrophage', 'MDSC', 'Neutrophil', 'Plasmacytoid dendritic cell')
    anti<- gsub('^ ','',rownames(gsva_matrix1))%in%anti_tumor
    pro<- gsub('^ ','',rownames(gsva_matrix1))%in%pro_tumor
    non <- !(anti|pro)##设定三种基因
    gsva_matrix1<- rbind(gsva_matrix1[anti,],gsva_matrix1[pro,],gsva_matrix1[non,])#再结合起来，使图分成三段
    normalization<-function(x){
      return((x-min(x))/(max(x)-min(x)))}#设定normalization函数
    nor_gsva_matrix1 <- normalization(gsva_matrix1)
    
    # annotation_col = data.frame(patient=pheno$num2)#加上病人编号
    # rownames(annotation_col)<-colnames(uni_matrix)#使编号能互相对应
    
    # bk = unique(c(seq(0,1, length=100)))#设定热图参数
    pheatmap(nor_gsva_matrix1,
         show_colnames = F,
         cluster_rows = F,cluster_cols = F,
         # annotation_col = annotation_col,
         # breaks=bk,
         # cellwidth=5,cellheight=5,
         # fontsize=5,
         gaps_row = c(12,20),)
         # filename = 'ssgsea.pdf',width = 8)#画热图
    # 保存数据
    # save(gsva_matrix,gsva_matrix1,pheno,file = 'score.Rdata')