最近学徒被委托模仿一个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图和热图的拼图展示。