RNA-seq(9):功能富集分析

这部分开始进行基本的富集分析,两类

  • A:差异基因富集分析(不需要表达值,只需要gene name)

  • B: 基因集(gene set)富集分析(不管有无差异,需要全部genes表达值)

############################################################

A:差异基因富集分析(不需要表达值,只需要gene name)

############################################################

-----------先说富集什么-----------

  • 最常用的基因注释工具是GOKEGG注释,这基本上是差异基因分析一定做的两件事。GO可以在GO:BP(生物过程),GO:MF(分子功能),GO:CC(细胞组分)三个方面分别进行注释,用的比较多的是GO:BP,但其他两方面也很重要。
  • 外还有一个软件不得不提,那就是IPA(Ingenuity pathway analysis),这是一个收费软件,有基本版和高级版。我个人觉得它的upstream regulator analysis还是很不错的。分子激活功能等也可以用用。另外一个就是它内置的热图功能。高级版我没用过,但是知道可以导出一些数据等。

-----------什么是富集(原理)-----------

富集的统计学基础是超几何分布,简单来说根据下面的Fisher精确检验(Fisher exact test)公式,对每个GO或KEGG term计算一个p值
p=(M/K)[(N-M)/(n-k)]/(N/n),其中
N:所有gene总数
n:N中差异表达gene的总数
M:N中属于某个GO term的gene个数
k: n中属于某个GO term的gene个数
p:表示差异表达gene富集到这个GO term上的可信程度

  • 当p<0.05或0.01,则认为差异表达gene显著到这个GO term上(自己定义p值)
  • 意义:提供的信息更集中,更有意义

---------------拿什么来富集---------------

得到的差异表达基因列表就可以,也就是说不需要其他的值

---------------用什么工具富集--------------

只能说实在是太多太多了。。。。但是用的时候要小心,因为你多用几个工具,即使设定同样的p值也会发现结果有出入,有时还差异挺大。

1 按使用方式来说(简单度)有3种

  • (1)在线版:最主流的就是DAVID,各种级别杂志总见其身影,使用非常简单,不再赘述。另外还有GatherGOrilla,revigo,还有很多很多我就不在贴了。网页版有网页版的好处,可以先大概看下自己筛选的genes。另外很多工具有很好的可视化功能,自己一一去探索吧。
  • (2)客户端版:IPA(IPA不是用的GO和KEGG数据库)和FUNRICH,后者更新速度很慢,但里面有好玩又实用的功能,并且可以加载自己的数据。
  • (3)R包:介绍一个就行了,那就是Y叔的clusterProfiler,我论文中的富集功能很多都是用这个包做的(还有的用了IPA)。
    ##########################################################

B: 基因集(gene set)富集分析(不管有无差异,需要全部genes表达值)

##########################################################

  • 好处:可以发现被差异基因舍弃的genes可能参与了某重要生理过程或信号通路(稍后我会把以前手写的GSEA原理和结果意义解读发上来)
  • 工具:GSEA
  • 使用方法:R(还是clusterProfiler)或客户端

-------------------具体操作---------------------

#enrichment analysis using clusterprofiler package created by yuguangchuang
library(clusterProfiler)
library(DOSE)
library(org.Mm.eg.db)
#get the ENTREZID for the next analysis
setwd("F:/rna_seq/data/matrix")
sig.gene<-read.csv(file="DEG_treat_vs_control.csv")
head(sig.gene)
gene<-sig.gene[,1]
head(gene)
gene.df<-bitr(gene, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Mm.eg.db)

head(gene.df)

GO富集分析:

#Go classification
#Go enrichment
ego_cc<-enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Mm.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
ego_bp<-enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Mm.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")+ 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(ego_bp) str_wrap(ego_bp,width = 25))
gobp.jpeg
#KEGG enrichment
install.packages("stringr")
library(stringr)
kk<-enrichKEGG(gene      =gene.df$ENTREZID,
               organism = 'mmu',
               pvalueCutoff = 0.05)
kk[1:30]
barplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
dotplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
kegg.jpeg

keggdot.jpeg
# Gene Set Enrichment Analysis(GSEA)
# 获取按照log2FC大小来排序的基因列表
genelist <- sig.gene$log2FoldChange
names(genelist) <- sig.gene[,1]
genelist <- sort(genelist, decreasing = TRUE)
# GSEA分析
gsemf <- gseGO(genelist,
               OrgDb = org.Mm.eg.db,
               keyType = "ENSEMBL",
               ont="BP"
)
# 查看信息
head(gsemf)
# 画出GSEA图
gseaplot(gsemf, geneSetID="GO:0001819")
gsea.jpeg

后记:做完这部分富集分析,接着按我的流程进入下一部分分析RNA-seq(10):KEGG通路可视化,因为直接用到这部分数据,

参考Y叔的包说明,里面写的特别详细
还有lxmic的简书

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

推荐阅读更多精彩内容