fgsea做GSEA

2020-07-21 10:00:13 浏览数 (2)

代码语言:javascript复制
library(fgsea)

1.导入测试数据,fgesa的examplePathways,exampleRanks测试数据分别是通路的list和经过fold change排序的基因。

代码语言:javascript复制
data(examplePathways)
data(exampleRanks)

2.运行 fgsea

The resulting table contains enrichment scores and p-values:

代码语言:javascript复制
fgseaRes <- fgsea(pathways = examplePathways, 
                  stats = exampleRanks,
                  minSize=15,
                  maxSize=500,
                  nperm=10000)
head(fgseaRes[order(pval), ])

3.计算padj<0.01的个数

代码语言:javascript复制
sum(fgseaRes[, padj < 0.01])

4.画特定的通路

代码语言:javascript复制
plotEnrichment(examplePathways[["5991130_Programmed_Cell_Death"]],
               exampleRanks)   labs(title="Programmed Cell Death")

5.选择上调下调的top10的pathtway

代码语言:javascript复制
topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n=10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
plotGseaTable(examplePathways[topPathways], exampleRanks, fgseaRes, 
              gseaParam = 0.5)

From the plot above one can see that there are very similar pathways in the table (for example 5991502_Mitotic_Metaphase_and_Anaphase and 5991600_Mitotic_Anaphase). To select only independent pathways one can use collapsePathways function:

代码语言:javascript复制
collapsedPathways <- collapsePathways(fgseaRes[order(pval)][padj < 0.01], 
                                      examplePathways, exampleRanks)
mainPathways <- fgseaRes[pathway %in% collapsedPathways$mainPathways][
                         order(-NES), pathway]
plotGseaTable(examplePathways[mainPathways], exampleRanks, fgseaRes, 
              gseaParam = 0.5)

最后保存结果。

代码语言:javascript复制
library(data.table)
fwrite(fgseaRes, file="fgseaRes.txt", sep="t", sep2=c("", " ", ""))

0 人点赞