【原】GSEA确实搭配热图后更直观易懂

最近学徒被委托模仿一个2022的发表在Biomaterials杂志的文章《Engineering dual catalytic nanomedicine for autophagy-augmented and ferroptosis-involved cancer nanotherapy》的转录组两分组标准分析以及出图，我们默认流程是差异分析和生物学功能数据库注释的富集分析。

其中生物学功能数据库注释目前最稳妥的是GSEA方法，但是文章在标准的gsea图下面加上了一个热图，蛮有意思的：

gsea图下面加上了一个热图

因为纯粹的GSEA方法出图后很多人都没办法理解。

假设芯片或者其它测量方法测到了2万个基因，那么这两万个基因在case和control组的差异度量 (其实有六种差异度量，默认是signal 2 noise，GSEA官网有提供公式，也可以选择大家熟悉的foldchange) 肯定不一样， 那么根据它们的差异度量，就可以对它们进行排序，并且Z-score标准化，在GSEA图的最底端展示的就是排序后的基因列表 ：

纯粹的GSEA方法出图

那么图中间，就是我们每个gene set里面的基因（通常是几十个基因甚至上百个基因）在所有的2万个排序好基因的位置，如果gene set里面的基因集中在2万个基因的前面部分，就是在case里面富集，如果集中在后面部分，就是在control里面富集着。

而最上面的那个ES score的算法，大概如下：

ES score的算法

仔细看，其实还是能看明白的，排序好的2万多个基因里面的每个基因在每次gsea分析的每个gene set里面的ES score取决于这个基因是否属于该gene set，还有就是它的差异度量，上图的差异度量就是FC（foldchange）, 对每个gene set来说，所有的基因的ES score都要一个个加起来，叫做running  ES score，在加的过程中，什么时候ES score达到了最大值，就是这个gene set最终的ES score ！

所以实际上大家只需要看每个gene set最终的ES score即可，比如：

library(clusterProfiler)
library(org.Hs.eg.db)
library(DOSE)
data(geneList)
head(geneList) #排序好的基因序列，而且是entrezeID的形式
R.utils::setOption( "clusterProfiler.download.method",'auto' )
kkgsea <- gseKEGG(geneList     = geneList,
               organism     = 'hsa', 
               minGSSize    = 10,
               maxGSSize    = 500,
               pvalueCutoff = 1,
               pAdjustMethod = "none" ) #进行gseKEGG富集分析 
kkgsea=setReadable(kkgsea,keyType = 'ENTREZID',OrgDb = 'org.Hs.eg.db')
tmp = kkgsea@result[,c(1:5,11)]
save(kkgsea,file = 'gsea_kk.Rdata')

如下所示可以看到不同通路确实里面的基因是不一样的，而且每个gene set最终的ES score 有正负值以及数值的大小差异 ：

image-20221030153240324

我们常规的可视化代码就是enrichplot包的gseaplot2函数，代码如下所示：

library(enrichplot)
gseaplot2(kkgsea,geneSetID = head(kkgsea@result$ID,12))

显示12个通路，这样可以看到正负值的区别， 可视化效果；

看到正负值的区别

上面的gsea和分析和可视化是直接使用了DOSE里面的data(geneList) ，这个排序好的基因序列，而且是entrezeID的形式。

实际上我们自己的 基因排序，往往是来自于自己的差异分析后的：

rm(list = ls())
options(stringsAsFactors = F)
load(file = "./data/Step03-DESeq2_nrDEG.Rdata")
load(file = "./data/Step01-airwayData.Rdata")
head(deg_anno)

library(clusterProfiler)
library(org.Hs.eg.db)
# 转成ENTREZID
df <- bitr(unique(deg_anno$SYMBOL), fromType = "SYMBOL",
           toType = c( "ENTREZID"),
           OrgDb = org.Hs.eg.db) 
DEG=merge(deg_anno,df,by='SYMBOL')
colnames(DEG)
data_all_sort <- DEG %>%  #排序
  arrange(desc(log2FoldChange))
geneList = data_all_sort$log2FoldChange #把foldchange按照从大到小提取出来
names(geneList) <- data_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(geneList) #排序好的基因序列，而且是entrezeID的形式

因为我们的基因排序，来自于差异分析的结果，就有表达量矩阵和分组信息，而我们的gsea分析里面的也是有对应的通路的里面的基因名字。就可以把每个通路的里面的基因拿去热图可视化，参考代码是：

# 其中 exprSet 是前面的转录组测序后的counts矩阵
# group_list 是矩阵里面的每个样品的分组信息
# up_kegg 是自己挑选好的通路 

pro = 'up'
print(pro)
dir.create('output/plot/up_kk_gse_heatmap')
savDir <- 'output/plot/up_kk_gse_heatmap/'

lapply(1:nrow(up_kegg), function(i){ 
  # i=12
  cg = strsplit(up_kegg$core_enrichment[i],'/')[[1]]
  cg
  length(cg)
  dat=log2(edgeR::cpm(exprSet)+1)
  #dat[1:4,1:4] 
  table(group_list)
  pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出，组合起来为一个新的表达矩阵
  n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
  n[n>2]=2 
  n[n< -2]= -2
  #  n[1:4,1:4]
  pheatmap(n,show_colnames =F,show_rownames = F)
  ac=data.frame(group=group_list)
  rownames(ac)=colnames(n) 
  pheatmap(n,show_colnames =F,show_rownames = T,
           annotation_col=ac)
  
  heatmap.title = paste0('Genes of ', up_kegg$Description[i])
  filename =  paste0(savDir, pro,'_',
                     gsub('/','-',up_kegg$Description[i]),
                     '.pdf')
  if(length(cg) > 15){plot.height <- length(cg)/5}
  if(length(cg) %in% c(3:15)){plot.height <- length(cg)*0.5}
  if(length(cg) %in% c(1:2)){plot.height <- length(cg)*1.5}
  pheatmap(n,main = heatmap.title,
           show_colnames =T,show_rownames = T,
           annotation_col=ac,
           height = plot.height,
           filename = filename) 
}) 

这样的热图并不是最开始文章里面的那样，因为它展现了每个通路里面的全部的基因，而我们的默认gsea分析后返回的是core_enrichment里面的基因，不一样。

学徒作业

拿最经典的airway这个转录组测序数据集里面的表达量矩阵和分组信息，走标准的差异分析后，对基因进行logFC的排序，然后走kegg数据库的gsea注释，选取上下调通路各一条进行 gsea图和热图的拼图展示。