RNA富集分析

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

  • A：差异基因富集分析（不需要表达值，只需要gene name）
  • B: 基因集(gene set)富集分析（不管有无差异，需要全部genes表达值）

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

A：差异基因富集分析（不需要表达值，只需要gene name）

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

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

  • 最常用的基因注释工具是GO和KEGG注释，这基本上是差异基因分析一定做的两件事。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,各种级别杂志总见其身影，使用非常简单，不再赘述。另外还有Gather，GOrilla,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的简书

赞 (0)

相关推荐