這篇文章中有一張圖很有趣，如下：

作者使用Hallmarks通路進(jìn)行GSEA富集分析，共發(fā)現(xiàn)26條通路顯著與兩種表型相關(guān)，與stemness表型相關(guān)的有16條通路，與cancer表型相關(guān)的有10條通路。

本次演練中，我們選擇MSIDB數(shù)據(jù)庫(kù)中的50條Hallmarks通路進(jìn)行示例，通路信息下載鏈接：http://www.gsea-msigdb.org/gsea/downloads.jsp

下面我們用我們自己的數(shù)據(jù)來(lái)做一下這張圖：

rm(list = ls()) library(edgeR) library(DESeq2) library(fgsea) library(clusterProfiler) library(enrichplot) library(ggplot2) load('data.rda')?##?加載我們的數(shù)據(jù) 包括臨床數(shù)據(jù)clin和表達(dá)數(shù)據(jù)expr head(expr)

表達(dá)矩陣為 raw read count data

head(clin)

第一列為expr每一列對(duì)應(yīng)的ID，第二列為分組信息。

table(clin$Group)

0和1分別有223個(gè)和135個(gè)

##構(gòu)建分組信息 group?<-?factor(rep(c('1','0'),times=c(135,223))) colData <- data.frame(row.names=rownames(clin),group) ##保留在50%以上的樣本中count>=1的基因 keep?<-?rowSums(expr>=1)?>=?ncol(expr)*0.5 table(keep) cc?<-?expr[keep,] ##差異分析 dds <- DESeqDataSetFromMatrix(round(cc), colData, design= ~group) dds <- DESeq(dds) res<-?results(dds,contrast=c("group","1","0"),independentFiltering=FALSE) ##差異分析結(jié)果 alldiff <- as.data.frame(res)%>%na.omit() alldiff$type <- ifelse(alldiff$padj>0.05,'No-Sig',ifelse(alldiff$log2FoldChange>1,'Up',ifelse(alldiff$log2FoldChange< -1,'Down','No-Sig')))table(alldiff$type)

##?順手畫(huà)個(gè)火山圖 ggplot(alldiff,aes(log2FoldChange,-log10(padj),fill=type))+geom_point(shape=21,aes(size=-log10(padj),color=color))+scale_fill_manual(values=c('seagreen','gray','orange'))+scale_color_manual(values=c('gray60','black'))+geom_vline(xintercept=c(-1,1),lty=2,col="gray30",lwd=0.6) +geom_hline(yintercept = -log10(0.05),lty=2,col="gray30",lwd=0.6)+theme_bw(base_rect_size = 1)+theme(axis.title = element_text(size = 15),axis.text = element_text(size = 12),legend.title = element_blank(),legend.text = element_text(size = 12),panel.grid = element_blank(),plot.title = element_text(family = 'regular',hjust = 0.5),legend.position = c(0.5, 1),legend.justification = c(0.5, 1),legend.key.height = unit(0.5,'cm'),legend.background = element_rect(fill = NULL, colour = "black",size = 0.5))+xlim(-4,4)+guides(size=F,color=F)+ylab('-log10 (FDR)')+xlab('log2 (Fold Change)')

接下來(lái)開(kāi)始重頭戲，開(kāi)始畫(huà)GSEA table

## 根據(jù)logfc降序排列基因 alldiff?<-?alldiff[order(alldiff$log2FoldChange,decreasing?=?T),] ##?fgsea中輸入的關(guān)鍵基因信息 id <- alldiff$log2FoldChange names(id)?<-?rownames(alldiff) ##?fgsea中輸入的關(guān)鍵通路信息 gmtfile <- "./h.all.v7.4.symbols.gmt" hallmark <- read.gmt(gmtfile) hallmark$term <- gsub('HALLMARK_','',hallmark$term) hallmark.list <- hallmark %>% split(.$term) %>% lapply( "[[", 2) ##?Perform?the?fgsea analysis fgseaRes <- fgsea(pathways = hallmark.list, stats = id,minSize=1,maxSize=10000,nperm=10000) sig <- fgseaRes[fgseaRes$padj<0.05,] sig <- sig[order(sig$NES,decreasing = T),] ##?最后一步 開(kāi)始繪圖 plotGseaTable(hallmark.list[sig$pathway],id, fgseaRes,gseaParam = 0.5)

這樣子基本上畫(huà)完了，但是貌似不是很好看，可以保存為PPT格式再處理一下

library(export) graph2ppt(file = 'GSEA-table.pptx',height = 7,width = 8.5)

其中每個(gè)元素都可以調(diào)整哦

調(diào)整完之后如下所示：

加編者微信入群 "生信交流群-醫(yī)學(xué)僧"

加微信時(shí)請(qǐng)備注 "學(xué)校-專業(yè)-姓名"

往期精品(點(diǎn)擊圖片直達(dá)文字對(duì)應(yīng)教程)

機(jī)器學(xué)習(xí)

后臺(tái)回復(fù)“生信寶典福利第一波”或點(diǎn)擊閱讀原文獲取教程合集

總結(jié)

以上是生活随笔為你收集整理的这篇Cell里面的GSEA展示很不错！的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。