ssGSEA

一腔孤勇，携我至此；一把年纪，还是幼稚；

Goals:

文献为《Local mutational diversity drives intratumoral
immune heterogeneity in non-small cell lung cancer》
要复现的图片为

Clues:

a.主要是两个操作：ssGSEA和ggplot（点图+cor.test结果呈现）;

b.大家看问题角度不同，其实会比较有趣，比如，用R包提取pdf的做法（因为我当时比较有时间，所以就不断地折腾，其实时间不够的话，我会复制粘贴到excel中再读进去）；比如，我的代码被不断说‘太丑’（还好有把参照的其他代码放进去，那你说丑就丑吧；我学习一下漂亮的；）；

c.对我来说，需要考虑的地方有两点，（1）如果你读了文章的话，会看到，绘制热图有提到nomalized to 0-1，后面搜索之后，才有了normalize那个function；（2）duplicated gene_name的过滤，要依据一定的规则会比较好，网上没有查到相关的有让我信服的推荐，后来，采用了jimmy的median排序进行操作（算是没有标准的标准）；这里需要了解，GSVA函数里的两个输入，data和gene_set，既然，gene_set是gene-name，data指定需要是gene-name为主体；

Steps:

1.从supplementarytable中将数据集读出来，这个在文章中有提到

虽然我不咋想承认，但故事总得圆圆满满，请点击学习这部分漂亮的代码；
附件

rm(list=ls())
library(pdftools)
options(stringsAsFactors = F)
b <- pdf_text('SupplementaryTables.pdf')
geneset_substract<- function(tmp){split_to_line<- gsub('\r','',strsplit(tmp,split = '\n')[[1]])
              gene_name<- apply(data.frame(split_to_line),1,function(x){ line<- strsplit(x,split=' ')[[1]]
                                                   pos<- grep('[A-Za-z\\d+]|\\d+',line)
                                                   res <- line[pos[1]]})
              cell_type<- apply(data.frame(split_to_line),1,function(x){ line<- strsplit(x,split=' ')[[1]]
                                                           pos<- grep('[A-Za-z\\d+]|\\d+',line)
                                                           res <- line[pos]
                                                           res <- res[c(-1,-length(res))]
                                                           s <- ''
                                                           for (i in 1:length(res)){
                                                             s<- paste(s,res[i])}
                                                           return(s)})
              result<- data.frame(gene_name,cell_type)
              return(result)
              }
gene_set <- data.frame()
for(i in 20:36){
  gene_set<- rbind(gene_set,geneset_substract(b[i]))
}
gene_set <- gene_set[c(-1,-2),]
list <- list()
for(i in 1:length(unique(gene_set$cell_type))){
  list[[i]] <- gene_set$gene_name[gene_set$cell_type== (unique(gene_set$cell_type)[i])]
}
names(list)<- unique(gene_set$cell_type)
save(list,file='gene_set.Rdata')

漂亮的代码：

trimws去掉字符串前后的空格（\n\r\t）

由于pdf_text读到的内容直接是一整个字符串，所以先按\n拆分（因为每行是一个单元，用换行符就比较容易理解了），再用trimws对拆出来以行为单位去掉前后的空格；

paste collapse参数可以将一个向量里的所用元素用collapse设置的内容连接成单个字符；

do.call 对以list形式反馈的参数，执行相应的功能；

rm(list=ls())
library(pdftools)
options(stringsAsFactors = F)

b <- pdf_text('SupplementaryTables.pdf')

tmp = unlist(lapply(20:36, function(i){
  trimws(strsplit(b[[i]],split = '\n')[[1]])
}))
tmp=tmp[-c(1,2)]
library(stringr)
tmp=do.call(rbind,lapply(str_split(tmp,' '), function(x){
  x=x[-length(x)]
  gene_name<- x[1]
  cell_type<- trimws(paste(x[-1],collapse = ' '))
  return(c(gene_name,cell_type))
}))
immune_list <- split(tmp[,1],tmp[,2])

2.提取矩阵和表型信息，需要手动从GEO下载，试了就知道为啥（如果试了不知道为啥，就留言）；进行ssGSEA分析，只是用到了处理后的矩阵和基因集两个内容；对score结果归一化后进行热图绘制；

rm(list=ls())
###矩阵信息提取
a <- read.table('GSE112996_merged_fpkm_table.txt.gz',
                header = T,
                row.names=1)
raw_data<- a[,-1]
###表型信息提取
pheno <- read.csv(file = 'GSE112996_series_matrix.txt')
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]))
####数据过滤
data<- a[!apply(raw_data,1,sum)==0,]
####去除重复基因名的行，归一化
data$median=apply(data[,-1],1,median)
data=data[order(data$GeneName,data$median,decreasing = T),]
data=data[!duplicated(data$GeneName),]
rownames(data)=data$GeneName
uni_matrix <- data[,grep('\\d+',colnames(data))]
uni_matrix <- log2(uni_matrix+1)
colnames(uni_matrix)<- gsub('X','',gsub('\\.','\\-',colnames(uni_matrix)))
uni_matrix<- uni_matrix[,order(colnames(uni_matrix))]
library(genefilter)
library(GSVA)
library(Biobase)
load('gene_set.Rdata')
gsva_matrix<- gsva(as.matrix(uni_matrix), list,method='ssgsea',kcdf='Gaussian',abs.ranking=TRUE)
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)))}
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.png')
save(gsva_matrix,gsva_matrix1,pheno,file = 'score.Rdata')

ssgsea.png

3.计算score加和后，ggplot2进行绘图；

rm(list=ls())
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')
load('score.Rdata')
anti<- as.data.frame(gsva_matrix1[gsub('^ ','',rownames(gsva_matrix1))%in%anti_tumor,])
pro<- as.data.frame(gsva_matrix1[gsub('^ ','',rownames(gsva_matrix1))%in%pro_tumor,])
anti_n<- apply(anti,2,sum)
pro_n<- apply(pro,2,sum)
patient <- pheno$num2[match(colnames(gsva_matrix1),pheno$num1)]
library(ggplot2)
data <- data.frame(anti=anti_n,pro=pro_n,patient=patient)
anti_pro<- cor.test(anti_n,pro_n,method='pearson')
gg<- ggplot(data,aes(x = anti, y = pro),color=patient) + 
  xlim(-20,15)+ylim(-15,10)+
  labs(x="Anti-tumor immunity", y="Pro-tumor suppression") +
  geom_point(aes(color=patient),size=3)+geom_smooth(method='lm')+
  annotate("text", x = -5, y =7.5,label=paste0('R=',round(anti_pro$estimate,4),'\n','p<0.001'))
ggsave(gg,filename = 'cor.png')

参考内容：
1.GSVA: The Gene Set Variation Analysis package for microarray and RNA-seq data

课程分享
生信技能树全球公益巡讲
（https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g）
B站公益74小时生信工程师教学视频合辑
（https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw）
招学徒：
（https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw）