先加載R包
> library(GSEABase)
> library(limma)
> library(clusterProfiler)
> library(enrichplot)
第一步：差異分析
首先是用limma包做差異分析的過程
#差異分析：
dat <- exprSet
group_list <- c(rep("Control",3),
rep("Combine",3))
group_list <- factor(group_list,ordered = F)
design <- model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg)
這里由于做的時候是FPKM數(shù)據(jù)，所以參考了生信技能樹的教程：
https://mp.weixin.qq.com/s/BEvFqu-BNA9lWFNaE9Jkdw
第二步：獲取基因集
MSigDB可以方便地獲取GMT文件，讀進來
> geneset <- read.gmt("./h.all.v7.5.1.symbols.gmt")
第三步：進行GSEA
> geneList <- deg$logFC
> names(geneList) <- toupper(rownames(deg))
> geneList <- sort(geneList,decreasing = T)
GSEA輸入的geneList是排序后的gene名提茁，因此這里需要對geneList進行一個排序饺窿。實際上得到差異分析的結(jié)果后鸿捧，只需要logFC和gene名就可以構(gòu)建geneList進行GSEA富集分析了嗤军。
> gsea_results <- GSEA(
+ geneList = geneList,
+ TERM2GENE = geneset,
+ verbose = F,
+ eps=1e-10
+ )
Warning messages:
1: In serialize(data, node$con) :
'package:stats' may not be available when loading
2: In serialize(data, node$con) :
'package:stats' may not be available when loading
3: In fgseaMultilevel(...) :
For some pathways, in reality P-values are less than 1e-10. You can set the `eps` argument to zero for better estimation.
P值太小時可能會報warning靶瘸，不過不影響寝优√跆颍可以把eps這個參數(shù)設置成0。得到的結(jié)果是一個gseaResult對象
> class(gsea_results)
[1] "gseaResult"
attr(,"package")
[1] "DOSE"
可視化
接下來可以對這個gseaResult對象進行可視化乏矾。
> for (i in 1:length(list)) {
+ p <- gseaplot2(x=gsea_results,geneSetID=paste("HALLMARK_",list[i],sep=""))
+ d <- paste("./",list[i],".pdf",sep="")
+ pdf(file=d,
+ family = "Times",width=10,height = 6)
+ print(p)
+ dev.off()
+ }
這里循環(huán)讀取感興趣的每一個基因集進行作圖孟抗，并輸出為PDF。
GSEA
當然也可以多個圖放在一起
> p <- gseaplot2(gsea_results,
+ gsea_results@result[["ID"]][1:3],
+ pvalue_table = TRUE)
比如這個樣子
gseaplot2
