fgsea单细胞转录组GSEA富集分析

Gene Set Enrichment Analysis(GSEA),Subramanian et al. 20051年发表在Proc Natl Acad Sci USA.提出来的一个基因功能富集方法。

fgsea能够快速对预选基因集进行富集分析,预选基因集可以是自己设定,一般使用MSigdb(Molecular Signatures Database)数据库,同样由提出GSEA方法团队提供。该数据库包含了以下9种不同基因的基因,可供下载以及R软件包载入。

  • H:癌症特征基因集合,共50 gene sets,最常用。
  • C1:染色体位置基因集合,共299 gene sets,用的很少。
  • C2:包含了已知数据库,文献和专家支持的基因集信息,包含5529 gene sets。
  • C3:调控靶基因集合,包括miRNA-target以及转录因子-target调控模式,3735 gene sets。
  • C4:计算机软件预测出来的基因集合,主要是和癌症相关的基因,858 gene sets。
  • C5:Gene Ontology对应的基因集合,10192 gene sets。
  • C6:致癌基因集合,大部分来源于NCBI GEO发表芯片数据,189 gene sets。
  • C7:免疫相关基因集,4872 gene sets。
  • C8:single cell identitiy gene sets, 302 gene sets

载入相关R包;
载入数据库基因集;
载入单细胞转录组Seurat对象,获取marker基因,以0群为例,相当于0群与其余各群的差异分析,非单细胞测序数据直接做差异分析。

library(msigdbr) #提供MSigdb数据库基因集
library(fgsea)
library(Seurat)
class(PBMC)
[1] "Seurat"

## 预定义基因集
mdb_c2 <- msigdbr(species = "Homo sapiens", category = "C2")
mdb_kegg = mdb_c2 [grep("^KEGG",mdb_c2 $gs_name),]
fgsea_sets<- mdb_kegg %>% split(x = .$gene_symbol, f = .$gs_name)

##  Seurat获取0群marker基因
cmarkers <- FindMarkers(PBMC,ident.1 = 0, only.pos = TRUE, 
                        min.pct = 0.1, logfc.threshold = 0)

head(cmarkers)
                 p_val avg_logFC pct.1 pct.2     p_val_adj
NR4A1     0.000000e+00 1.1942463 0.977 0.599  0.000000e+00
SERPINF1  0.000000e+00 1.1067154 0.939 0.385  0.000000e+00
FOSB     5.600685e-304 0.8970092 0.999 0.893 1.734196e-299
GEM      2.809004e-301 1.1730836 0.739 0.230 8.697800e-297
ANTXR1   2.237076e-300 1.0842559 0.825 0.303 6.926881e-296
C1R      2.444552e-282 0.9066228 0.988 0.496 7.569310e-278

GSEA富集分析(fgsea),使用预定义的基因集(通常来自功能注释或先前实验的结果),将基因按照在两类样本中的差异表达程度排序,然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。

#基因按logFC排序
cmarkers$genes = rownames(cmarkers)
cluster0.genes<- cmarkers %>% arrange(desc(avg_logFC)) %>% dplyr::select(genes,avg_logFC)
ranks<- deframe(cluster0.genes)

#fgsea
fgseaRes<- fgsea(fgsea_sets, stats = ranks, nperm = 1000)

画感兴趣的通路

plotEnrichment(fgsea_sets[["KEGG_P53_SIGNALING_PATHWAY"]],
               ranks) + labs(title="KEGG_P53_SIGNALING_PATHWAY")

如下图1,展示了其中KEGG_P53_SIGNALING_PATHWAY通路结果,横坐标从大到小排序好的所有基因,靠近0的基因代表表达量与表型呈现正相关,远离0的为负相关基因。纵坐标纵坐标是富集打分ES,这个富集打分是一个动态的值,通常用偏离0最远的值作为富集打分。 中间的竖线为感兴趣基因集所处的位置。


图1 KEGG_P53_SIGNALING_PATHWAY

以P值小于0.05筛选有显著统计学意义的通路,以富集分数NES排序筛选前20个通路,绘制条形图。

ggplot(fgseaRes %>% as_tibble() %>% arrange(desc(NES)) %>% filter(pval < 0.05) %>% 
         head(n= 20), aes(reorder(pathway, NES), NES)) +
  geom_col(aes(fill= NES)) +
  coord_flip() +
  labs(x="KEGG", y="Normalized Enrichment Score",title="KEGG gene sets NES from GSEA")

如图2所示,NES为每个基因子集ES根据基因集大小标准化的富集分数,ES反应基因集在排序列表两端的富集程度。


图2 GSEA通路富集分析

关于GSEA结果解读见另一篇笔记https://www.jianshu.com/p/f7f89d7f1cb7

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 201,312评论 5 473
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,578评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 148,337评论 0 333
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,134评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,161评论 5 363
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,303评论 1 280
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,761评论 3 393
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,421评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,609评论 1 295
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,450评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,504评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,194评论 3 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,760评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,836评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,066评论 1 257
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,612评论 2 348
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,178评论 2 341