富集分析实战

最新推荐文章于 2024-05-30 09:12:31 发布

我是一只grizzly

最新推荐文章于 2024-05-30 09:12:31 发布

阅读量1.2k

点赞数 3

文章标签： r语言

本文链接：https://blog.csdn.net/weixin_49814193/article/details/131714811

版权

富集分析实战

富集是指将基因按照先验知识，也就是基因组注释信息，对基因进行分类的过程。基因经过分类后，能够帮助认知寻找到的基因是否具有某方面的共性(如功能、组成等等).

大家最常用或最常见的富集分析就是ORA（GO,KEGG）和GSEA富集分析，这节内容主要讲述如何对自己得到的基因集进行富集分析并可视化，以及如何选择自己想要的通路进行个性化的富集分析可视化。

ORA (Over Representation Analysis)

过表达分析

对差异的结果设置一个标准(通常也就是我们差异表达的阈值，比如（LogFC >1|P-Value <0.05))。如果达到这个标准的基因作为候选分析的基因。将基因名都输入到R语言中或者网站中，结合数据库就得到分析的结果了。
所以这个算法的主要输入条件其实就是基因名。通常来讲，ORA使用的背景数据集就是GO和KEGG。

做任何分析的第一步就是要整理好我们的数据，让数据成为程序可阅读的样子

首先我们拿到我们上一步差异分析得到的结果，以GSE7793为例，如图：

View(DEG_p_7793_ip_400)

可以看到这是我们差异分析之后通过阈值筛选得到的结果

首先是加载R包

library(tidyverse)
library(org.Hs.eg.db) #人源
library(org.Mm.eg.db) #鼠源
library(clusterProfiler) #对于基因名的转换，基因的富集分析等很多功能我们都会用到该包

对于GO和KEGG来说，我们只需要取基因名即可,但大部分人可能会分别做上调基因和下调基因的GO分析

DEG_df <- rownames_to_column(DEG_p_7793_ip_400,var = "symbol") #行名转列
#首先做GO分析
positive_rows <- DEG_df[DEG_df$logFC > 0,] #上调基因
gene.df_up <- bitr(positive_rows$symbol,
                fromType = "SYMBOL",
                toType = c("ENSEMBL", "ENTREZID"),
                OrgDb = org.Mm.eg.db) #注意这是小鼠的基因

gene_GO_up <- enrichGO(gene = gene.df_up$ENTREZID,
                   OrgDb = org.Mm.eg.db,ont = "ALL", 
                   pvalueCutoff = 0.05, 
                   qvalueCutoff = 0.05,
                   readable = TRUE)
#ont这里可以选择BP,CC,MF或者全选
p1 <- barplot(gene_GO_up, showCategory = 20,color = "pvalue")#展示前20个go结果
p11 <- dotplot(gene_GO_up, showCategory = 20,color = "pvalue")#展示前20个go结果
gene_GO

negative_rows <- DEG_df[DEG_df$logFC < 0,] #下调基因
gene.df_down <- bitr(negative_rows$symbol,
                fromType = "SYMBOL",
                toType = c("ENSEMBL", "ENTREZID"),
                OrgDb = org.Mm.eg.db) #注意这是小鼠的基因

gene_GO_down <- enrichGO(gene = gene.df_up$ENTREZID,
                   OrgDb = org.Mm.eg.db,ont = "ALL", 
                   pvalueCutoff = 0.05, 
                   qvalueCutoff = 0.05,
                   readable = TRUE)
#ont这里可以选择BP,CC,MF或者全选
p2 <- barplot(gene_GO_down, showCategory = 20,color = "pvalue")#展示前20个go结果
p22 <- dotplot(gene_GO_down, showCategory = 20,color = "pvalue")#展示前20个go结果

p3 <- p1|p2
p3

p33 <- p11|p22
p33

#KEGG
gene.df <- bitr(DEG_df$symbol,
                fromType = "SYMBOL",
                toType = c("ENSEMBL", "ENTREZID"),
                OrgDb = org.Mm.eg.db) #注意这是小鼠的基因

gene_kk <- gene.df$ENTREZID

kk <- enrichKEGG(gene = gene_kk, 
                 organism = "mmu", #注意物种是小鼠
                 pvalueCutoff = 1, 
                 qvalueCutoff =1)
p4 <- barplot(kk, showCategory = 20,color = "pvalue")#展示前20个go结果
p44 <- dotplot(kk, showCategory = 20,color = "pvalue")#展示前20个go结果
p4
p44

结果如图：

上调基因与下调基因的GO：

KEGG富集结果：

当然还有些朋友得到了GO，KEGG的结果之后，想要自己展示一些个性化的通路或者疾病富集分析，我们拿KEGG举例，当然GO也可以使用相同的操作

简单介绍一下KEGG

KEGG包含了六类生物学板块，有时候可能你只想给老板汇报一下哪些通路被富集到了。

write.csv(kk,file = "C:\\Users\\78467\\Desktop\\kk.csv") #先保存下来看看

得到了很多生物学板块的富集结果

我们选取其中五个来作图吧，分别是MAPK，TNF，cAMP ，Wnt和Rap1信号通路

kk <- read.csv(file = "C:\\Users\\78467\\Desktop\\kk.csv") #加载进来
selected_paths <- kk[c(1,2,7:9),]#这是这五个通路的位置
# 条形图
p5 <- ggplot(selected_paths, aes(x =-log10(pvalue) , y = reorder(Description, pvalue,decreasing = TRUE), fill = pvalue)) +
  geom_bar(stat = "identity") +
  scale_fill_gradient(low = "red", high = "blue") +
  labs(x = "Description", y = "-log10(p-value)", fill = "p-value") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))
#注意P值是由小到大排序的
# 点图
p55 <- ggplot(selected_paths, aes(x = reorder(Description, pvalue, decreasing = TRUE), y = -log10(pvalue))) +
  geom_point(aes(fill = pvalue, size = Count), shape = 21, color = "NA", stroke = 0.5) +
  scale_fill_gradient(low = "red", high = "blue") +
  scale_size(range = c(5, 10)) +
  labs(x = "Description", y = "-log10(p-value)", fill = "p-value", size = "Count") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  coord_flip()
p6 <- p5|p55

个性化通路结果：

GSEA (Gene Set Enrichment Analysis)

基因集富集分析

ORA富集分析时，我们只需要一个差异基因的列表，根本不关心这个差异基因究竟是上调还是下调。这是因为，传统的富集分析根本不需要考虑基因表达量的变化趋势，其算法的核心只关注这些差异基因的分布是否和随机抽样得到的分布一致，即使后期在可视化时，我们在通路图上用不同颜色标记了上下调的基因，但是由于没有采用有效的统计学手段去分析这条通路下所有差异基因的总体变化趋势，这使得传统的富集分析结果无法回答上述的问题。

相较于ORA一刀切的方式来选择输入基因，GSEA的算法认为，虽然有一些基因不满足人为设定的的筛选标准，但是也可能在生物学过程中发挥重要作用。

GSEA通过整个基因组表达的情况来评估主要是哪些通路有意义。所以我们输入的数据类型要包括测序或芯片数据的所有基因以及基因在两个group之间的差异(一般我们使用logFC)。

注意如果你的测序或芯片数据是小鼠的基因表达矩阵，你需要将基因名映射为人类的，因为目前GSEA这个函数只支持人类的GSEA富集分析！

View(DEG_7793_ip_400) #所有基因的一个表达结果

数据处理与富集分析：

expr <- DEG_7793_ip_400[,c(1,2)] #我们只需要基因名和logFC
library(homologene) #这个包用于基因名同源转换
homogenes <- homologene(expr$Gene, inTax = 10090, outTax = 9606) #9606是人类的ENTREZID
colnames(homogenes)[1] = "Gene"#名字改为一样的才可以取交集
genes_co <- inner_join(expr,homogenes,by = "Gene")
GSEA <- genes_co[,c(5,2)] #我们只需要基因名和logFC
GSEA<-na.omit(GSEA) #去除缺失值，一般来说没有
GSEA$logFC<-sort(GSEA$logFC,decreasing = T)#一定要排序
geneList = GSEA[,2]
names(geneList) = as.character(GSEA[,1])
str(geneList)
head(geneList) #因为GSEA分析的数据是具有命名的数值向量，所以要这一步是必要的
# 富集分析
Go_gseresult <- gseGO(geneList, 'org.Hs.eg.db', keyType = "ENTREZID", ont="all", nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=1) #这里一定要注意种属

KEGG_gseresult <- gseKEGG(geneList, nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=1)

Go_Reactomeresult <- gsePathway(geneList, nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=1)
dim(KEGG_gseresult)
[1] 338  11

结果可视化：

head(KEGG_gseresult,n = 20)
library(enrichplot)
gseaplot2(KEGG_gseresult, "hsa04370", color = "firebrick", rel_heights=c(1, .2, .6)) #改变更多参数，为了美观

ult)
[1] 338 11


**结果可视化：**

```R
head(KEGG_gseresult,n = 20)
library(enrichplot)
gseaplot2(KEGG_gseresult, "hsa04370", color = "firebrick", rel_heights=c(1, .2, .6)) #改变更多参数，为了美观