- library("clusterProfiler")
- library("org.Hs.eg.db")
GO分析与KEGG分析
GO网站建设定制开发网站建设定制开发分析需要一个基因 symbol列表，网站建设定制开发列表中为差异表达基因。
一、读入数据
- result<- read.csv(file = "Results/gleason high vs low_DESeq2差异分析/gleason high vs low_result.csv", header=T, row.names=1,check.names=FALSE)
-
- t_index=result$Change %in% c('up','down')
- DEG_symbol=rownames(result)[t_index]
这里以gleason high vs low_result.csv为例，得到的DEG_symbol网站建设定制开发即为所需要的 symbol列表
View(result)
二、symbol转换为entrezid
1 转换
- DEG_entrezid = mapIds(x = org.Hs.eg.db,
- keys = DEG_symbol,
- keytype = "SYMBOL",
- column = "ENTREZID") #存在NA
转换后的DEG_entrezid是character vector，其中有NA值
2 去除DEG_entrezid中的NA值
DEG_entrezid=na.omit(DEG_entrezid) na.omit()网站建设定制开发函数能去除所有含有NA的行
三、GO分析与KEGG分析
- GO_BP = enrichGO(gene = DEG_entrezid,
- OrgDb = org.Hs.eg.db,
- keyType = "ENTREZID",
- ont = "BP", #'BP','CC','MF'
- pvalueCutoff = 0.5,
- qvalueCutoff = 0.5)
-
- KEGG <- enrichKEGG(DEG_entrezid,
- organism = 'hsa', ## hsa网站建设定制开发为人的简写
- keyType = 'kegg',
- pvalueCutoff = 0.05,
- pAdjustMethod = 'BH',
- minGSSize = 3,
- maxGSSize = 500,
- qvalueCutoff = 0.5,
- use_internal_data = FALSE)
四、可视化
- dotplot(GO_BP,
- showCategory=5, #展示5个通路
- title="EnrichmentGO_BP")
- barplot(GO_BP,
- showCategory=5,
- title="EnrichmentGO_BP")
GSEA分析
GSEA分析需要一个 gene_list vector，为foldchange（或logFC）的降序排列，同时标注ENTREZID
在Rstudio中显示如下：
- library(clusterProfiler)
- library(org.Hs.eg.db)
- library(enrichplot)
一、读入数据
- result<- read.csv(file = "Results/gleason high vs low_DESeq2差异分析/gleason high vs low_result.csv", header=T, row.names=1,check.names=FALSE)
-
- DEG_fc=data.frame('Symbol'=rownames(result),'LogFC'=result[,2],stringsAsFactors=F) #合并rownames和log2Foldchange列
View(DEG_fc)
二、转换ENTREZID，获取gene_list vector
- DEG_fc[,'ENTREZID'] = mapIds(x = org.Hs.eg.db,
- keys = DEG_fc$Symbol,
- keytype = "SYMBOL",
- column = "ENTREZID") #存在NA
- DEG_fc=na.omit(DEG_fc) #去除NA值
View(DEG_fc)
- gene_list=DEG_fc$LogFC #提取logFC列
- names(gene_list)=DEG_fc$ENTREZID #加上ENTREZID
- gene_list = sort(gene_list, decreasing = T) #降序排列
三、GSEA分析
- gsea_KEGG <- gseKEGG(gene_list,
- organism = "hsa",
- keyType = "kegg")
若将gsea_KEGG转化为dataframe格式，如下图：
四、可视化
1 path 为需要展示的pathway id，这里展示的是enrichment score最高的4条通路
- t_index=order(gsea_KEGG_d$enrichmentScore,decreasing = T)
- path=rownames(gsea_KEGG[t_index,])[1:4] #选择展示的pathway
2 作图
- gseaplot2(gsea_KEGG,
- path,
- subplots = 1:2, #展示前2个图（共3个）
- pvalue_table = T, #显示p值
- title = "Olfactory transduction", #设置title
- base_size = 10, #字体大小
- #color="red") #线条颜色可选
GO分析
1 GO分析分3个层面
- Cellular component，CC 细胞成分
- Biological process， BP 生物学过程
- Molecular function，MF 分子功能
- Cellular component解释的是基因存在在哪里，在细胞质还是在细胞核？如果存在细胞质那在哪个细胞器上？如果是在线粒体中那是存在线粒体膜上还是在线粒体的基质当中？这些信息都叫Cellular component。
- Biological process是在说明该基因参与了哪些生物学过程，比如，它参与了rRNA的加工或参与了DNA的复制，这些信息都叫Biological process
- Molecular function在讲该基因在分子层面的功能是什么？它是催化什么反应的？
So, we will have a gene annotation infarmation.
立足于这三个方面，我们将得到基因的注释信息。
2 结果解读
气泡图（dotplot）
- 横坐标是
GeneRatio，意思是说输入进去的基因，它每个term（纵坐标）占整体基因的百分之多少； - 圆圈的大小代表
gene count的多少，图中给出了最大的圆圈代表11个基因； - 圆圈的颜色代表
P-value；
总体来说，P-value越小gene count圈越大，富集就越可信
柱形图（barplot）
- 柱子长度代表
gene count的多少，图中给出了最长柱子代表11个基因； - 柱子的颜色代表
P-value；
总体来说，P-value越小gene count越大，富集就越可信
KEGG分析
1 KEGG富集分析，即Pathway富集分析
除了对基因本身功能的注释，我们也知道基因会参与人体的各个通路，基于人体通路而形成的数据库就是通路相关的数据库。而KEGG就是通路相关的数据库的一种。
2 结果解读
和GO分析的dotplot、barplot相同
GSEA分析
传统的基因功能富集分析的时候肯定遇到这样的情况，一条富集到的通路中，既有上调的差异表达基因，也有下调的差异表达基因，那么这条通路总体是被抑制还是被激活呢？那么这条通路中的基因表达水平在实验组相比于对照组究竟是上升了呢，还是下降了呢？
GSEA富集分析解决的就是 pathways 被上调或者下调的问题。
KEGG数据库
1 简介
KEGG是一个综合数据库，它们大致分为系统信息、基因组信息、化学信息和健康信息四大类。进一步可细分为15个主要的数据库。
2 查询 hsa pathways id 和 description
方法1
path:hsa00010 Glycolysis / Gluconeogenesis - Homo sapiens (human)path:hsa00020 Citrate cycle (TCA cycle) - Homo sapiens (human)path:hsa00030 Pentose phosphate pathway - Homo sapiens (human)path:hsa00040 Pentose and glucuronate interconversions - Homo sapiens (human)path:hsa00051 Fructose and mannose metabolism - Homo sapiens (human)path:hsa00052 Galactose metabolism - Homo sapiens (human)... ...
方法2
物种选择hsa，keywords输入description