GEO生信数据挖掘(八)富集分析(GO 、KEGG、 GSEA 打包带走)

第六节,我们使用结核病基因数据,做了一个数据预处理的实操案例。例子中结核类型,包括结核,潜隐进展,对照和潜隐,四个类别。第七节延续上个数据,进行了差异分析。 本节对差异基因进行富集分析。

目录

数据展示

GO富集分析 -对基因名称映射基因ID

GO富集分析 -从org.Hs.eg.db库中去匹配基因

KEGG富集分析 (不详细讲了看注释)

GSEA 富集分析

更多复杂的图(关联网络图、八卦图 、弦图)


数据展示

差异基因计算完毕的指标如下图所示

差异基因筛选后表达矩阵

GO富集分析 -对基因名称映射基因ID

加载数据


#&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#+&&&&&&&&&&&&&&&&&&加载数据&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
load( "DEG_TB_LTBI_step13.Rdata")
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#+&&&&&&&&&&&&&&&&&&加载数据&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&library(clusterProfiler)
library(org.Hs.eg.db)#增加基因名
all_diff$SYMBOL=rownames(all_diff)#基因名称转换注释
gene_ids_DEG_TB_LTBI = bitr(geneID = rownames(dataset_TB_LTBI_DEG),fromType="SYMBOL",toType = c("ENTREZID","ENSEMBL","SYMBOL"),OrgDb = 'org.Hs.eg.db',drop =TRUE)#合并 增加logFC 为后续GSEA富集分析所需数据准备
gene_ids_DEG_TB_LTBI <- merge(gene_ids_DEG_TB_LTBI,all_diff,by="SYMBOL")
#观察
dim(gene_ids_DEG_TB_LTBI)
head(gene_ids_DEG_TB_LTBI)#获取基因ID ENSEMBL
gene_ENSEMBL <- gene_ids_DEG_TB_LTBI$ENSEMBL
gene_ENTREZID <- gene_ids_DEG_TB_LTBI$ENTREZID
gene_SYMBOL<- gene_ids_DEG_TB_LTBI$SYMBOL

经过映射,2813个差异基因得到2551个基因ID,下图为三种不同形式的基因名称,富集分析时,按需进行转换。

GO富集分析 -从org.Hs.eg.db库中去匹配基因

#Go富集分析,从库中去匹配
go <- enrichGO(gene_SYMBOL,OrgDb = org.Hs.eg.db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'SYMBOL')#进行GO富集,确定P值与Q值得卡值并使用BH方法对值进行调整。
#查看富集结果
dim(go)
#导出GO富集的结果
write.csv(go,file="go1.csv")

绘制气泡图

#绘制气泡图
pdf(file="15aGO富集分析step15.pdf", width = 9, height = 6)
dotplot(go,showCategory=20,label_format = 80)#气泡图dev.off()

三种不同类别的合并的气泡图(#CC细胞组件,MF分子功能,BP生物学过程)

pdf(file="15bGO富集分析三组step15.pdf", width = 9, height = 6)#CC细胞组件,MF分子功能,BP生物学过程
goCC <- enrichGO(gene = gene_ENTREZID,  #基因列表(转换的ID)keyType = "ENTREZID",  #指定的基因ID类型,默认为ENTREZIDOrgDb=org.Hs.eg.db,  #物种对应的org包ont = "CC",   #CC细胞组件,MF分子功能,BP生物学过程pvalueCutoff = 0.05,  #p值阈值pAdjustMethod = "fdr",  #多重假设检验校正方式minGSSize = 1,   #注释的最小基因集,默认为10maxGSSize = 500,  #注释的最大基因集,默认为500qvalueCutoff = 0.2,  #q值阈值readable = TRUE)  #基因ID转换为基因名goBP <- enrichGO(gene_ENTREZID,OrgDb = org.Hs.eg.db, ont='BP',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')goMF <- enrichGO(gene_ENTREZID,OrgDb = org.Hs.eg.db, ont='MF',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')
#通过ggplot2将BP、MF、CC途径的富集结果挑选前8条绘制在一张图上
barplot(go,label_format=100, split="ONTOLOGY")+ facet_grid(ONTOLOGY~.,scale="free")dev.off()

KEGG富集分析 (不详细讲了看注释)

#+=================================================================
#============================================================
#+========KEGG富集分析 气泡图step16===================
#+==========================================
#+================================#KEGG富集分析
pdf(file="16KEGG富集分析step16.pdf", width = 9, height = 6)
kegg<- enrichKEGG(gene = gene_ENTREZID,   #基因列表(ENTREZID ID: 54490,51144,31,3906) organism = "hsa",  #物种keyType = "kegg",  #指定的基因ID类型,默认为keggminGSSize = 3, maxGSSize = 500,pvalueCutoff = 0.05,  pAdjustMethod = "fdr", # pAdjustMethod = 'BH'qvalueCutoff = 0.02)
#观察
dim(kegg)
#绘制气泡图
dotplot(kegg)
dev.off()#kegg 增加可读性,对基因ID 转基因名
kegg_enrich_results <- DOSE::setReadable(kegg, OrgDb="org.Hs.eg.db", keyType='ENTREZID') #ENTREZID to gene Symbol#保存kegg结果
write.csv(kegg_enrich_results@result,'KEGG_gene_up_enrichresults.csv') 
#save(kegg_enrich_results, file ='KEGG_gene_up_enrichresults.Rdata')##查看与选择所需通路
kegg_enrich_results@result$Description[1:10] #查看前10通路###选择所需通路的ID号
i=1
select_pathway <- kegg_enrich_results@result$ID[i] #选择所需通路的ID号#&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
save(gene_ids_DEG_TB_LTBI,go,keggfile ="15_gene_ids_DEG_TB_LTBI.Rdata")
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

GSEA 富集分析


#+=================================================================
#============================================================
#+========GSEA 富集分析 气泡图step17===================
#+==========================================
#+================================# GSEA 分析
#需要把多个方法取并集
#该方法的输入需要基因和 logFC 排序后的结果
#不同方法 相同基因的的logFC值不一样,直接保留第一个重复基因library(stringr)## 去重  #去除NA值
dim(gene_ids_DEG_TB_LTBI)
colnames(gene_ids_DEG_TB_LTBI)gene_list_df = gene_ids_DEG_TB_LTBI[,c('ENTREZID','logFC')]
gene_list_df_na <- na.omit(gene_list_df)gene_ids_TB_LTBI_distinct <- dplyr::distinct(gene_list_df_na,ENTREZID,.keep_all=TRUE) dim(gene_ids_TB_LTBI_distinct)
gene_list=gene_ids_TB_LTBI_distinct$logFC   #提取logFC列
names(gene_list)=gene_ids_TB_LTBI_distinct$ENTREZID             #加上ENTREZID
gene_list_gsea = sort(gene_list, decreasing = T)  #降序排列gsea_KEGG <- gseKEGG(gene_list_gsea,organism = "hsa",keyType = "kegg")gsea_KEGG_d <- as.data.frame(gsea_KEGG)gsea_KEGG_d
#path 为需要展示的pathway id,这里展示的是enrichment score最高的4条通路t_index=order(gsea_KEGG_d$enrichmentScore,decreasing = T)
path=rownames(gsea_KEGG[t_index,]) #选择展示的 pathwayrownames(gsea_KEGG[t_index,]) [1:4]#作图 
pdf(file="17GSEA富集分析step17.pdf", width = 9, height = 6)
gseaplot2(gsea_KEGG,path, subplots = 1:2, #展示前2个图pvalue_table = T, #显示p值title = "Olfactory transduction",  #设置titlebase_size = 10, #字体大小color="red") #线条颜色可选
dev.off()#&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
save(gene_ids_DEG_TB_LTBI,go,kegg,gene_list_gsea,gsea_KEGG,file ="17_gene_ids_DEG_TB_LTBI.Rdata")
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

更多复杂的图(关联网络图、八卦图 、弦图)

参考 最全的GO, KEGG, GSEA分析教程(R),你要的高端可视化都在这啦![包含富集圈图] - 糖糖家的老张的文章 - 知乎 https://zhuanlan.zhihu.com/p/377356510

原文链接:https://blog.csdn.net/qq_50898257/article/details/120588222

#+=================================================================
#============================================================
#+========富集分析 更多的图step18===================
#+==========================================
#+================================
library(clusterProfiler)
library(enrichplot)
#+富集基因与所在功能集/通路集的关联网络图:
enrichplot::cnetplot(go,circular=FALSE,colorEdge = TRUE)#基因-通路关联网络图
enrichplot::cnetplot(kegg,circular=FALSE,colorEdge = TRUE)#circluar为指定是否环化,基因过多时建议设置为FALSEGO2 <- pairwise_termsim(go)
KEGG2 <- pairwise_termsim(kegg)
enrichplot::emapplot(GO2,showCategory = 50, color = "p.adjust", layout = "kk")#通路间关联网络图
enrichplot::emapplot(KEGG2,showCategory =50, color = "p.adjust", layout = "kk")write.table(kegg$ID, file = "KEGG_IDs.txt", #将所有KEGG富集到的通路写入本地文件查看append = FALSE, quote = TRUE, sep = " ",eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE, qmethod = c("escape", "double"),fileEncoding = "")
#打印几条通路名称看看
kegg$ID[1:3]
#打开浏览器观察通路
browseKEGG(kegg,"hsa04660")#选择其中的hsa05166通路进行展示#富集弦图
genedata<-data.frame(ID=gene_ids_DEG_TB_LTBI$SYMBOL  ,logFC=gene_ids_DEG_TB_LTBI$logFC)
write.table(go$ONTOLOGY, file = "GO_ONTOLOGYs.txt", #将所有GO富集到的基因集所对应的类型写入本地文件从而得到BP/CC/MF各自的起始位置如我的数据里是1,2103,2410append = FALSE, quote = TRUE, sep = " ",eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE, qmethod = c("escape", "double"),fileEncoding = "")'''
根据计算出的go 文件数量,调整
'''
GOplotIn_BP<-go[1:178,c(2,3,7,9)] #提取GO富集BP的前10行,提取ID,Description,p.adjust,GeneID四列
GOplotIn_CC<-go[179:194,c(2,3,7,9)]#提取GO富集CC的前10行,提取ID,Description,p.adjust,GeneID四列
GOplotIn_MF<-go[195:209,c(2,3,7,9)]#提取GO富集MF的前10行,提取ID,Description,p.adjust,GeneID四列library(stringr)
GOplotIn_BP$geneID <-str_replace_all(GOplotIn_BP$geneID,'/',',') #把GeneID列中的’/’替换成‘,’
GOplotIn_CC$geneID <-str_replace_all(GOplotIn_CC$geneID,'/',',')
GOplotIn_MF$geneID <-str_replace_all(GOplotIn_MF$geneID,'/',',')
names(GOplotIn_BP)<-c('ID','Term','adj_pval','Genes')#修改列名,后面弦图绘制的时候需要这样的格式
names(GOplotIn_CC)<-c('ID','Term','adj_pval','Genes')
names(GOplotIn_MF)<-c('ID','Term','adj_pval','Genes')
GOplotIn_BP$Category = "BP"#分类信息
GOplotIn_CC$Category = "CC"
GOplotIn_MF$Category = "MF"BiocManager::install('GOplot')
library(GOplot)
circ_BP<-GOplot::circle_dat(GOplotIn_BP,genedata) #GOplot导入数据格式整理
circ_CC<-GOplot::circle_dat(GOplotIn_CC,genedata) 
circ_MF<-GOplot::circle_dat(GOplotIn_MF,genedata) chord_BP<-chord_dat(data = circ_BP,genes = genedata) #生成含有选定基因的数据框
chord_CC<-chord_dat(data = circ_CC,genes = genedata) 
chord_MF<-chord_dat(data = circ_MF,genes = genedata) 
'''
> chord_CC<-chord_dat(data = circ_CC,genes = genedata) 
Error in `[<-`(`*tmp*`, g, p, value = ifelse(M[g] %in% sub2$genes, 1,  : subscript out of bounds我去检查了go和genelist的数据结构发现,genelist里的gene用的是gene名,而go里的基因用的是基因ID,不一样了,所以跑不出结果,所以我把genelist的gene换成了基因ID,就能跑出来了。作者:ff的小世界勿扰
链接:https://www.jianshu.com/p/ee4012fd253f
来源:简书
著作权归作者所有。商业转载请联系作者获得授权,非商业转载请注明出处。
'''#可以画 数量太多了
GOChord(data = chord_BP,#弦图title = 'GO-Biological Process',space = 0.01,#GO Term间距limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,lfc.col = c('red','white','blue'), #上下调基因颜色process.label = 10) #GO Term字体大小
GOChord(data = chord_CC,title = 'GO-Cellular Component',space = 0.01,limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,lfc.col = c('red','white','blue'), process.label = 10) 
GOChord(data = chord_MF,title = 'GO-Mollecular Function',space = 0.01,limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,lfc.col = c('red','white','blue'), process.label = 10)
'''
Warning messages:
1: Using size for a discrete variable is not advised. 
2: Removed 15 rows containing missing values (`geom_point()`). 
'''

富集分析完毕!

回顾我们用到方法,差异分析后进行富集分析,理论基础实际上就是简单的找不同,分析。

实际应用种,由于基因之间存在关联,另一套分析理论考虑的是基因之间的相互作用,下节,我们来看非常火的WGCNA 共表达加权网络进行基因分析。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/111940.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

大咖云集,智慧碰撞|第 18 届 CLK 大会完整议程揭晓(内附报名通道)

自 2006 年以来&#xff0c;在国内 Linux 技术爱好者和行业公司的鼎力支持下&#xff0c;中国 Linux 内核开发者大会已走过 17个年头&#xff0c;是中国 Linux 内核领域最具影响力的峰会之一。今年的中国内核开发者大会依然秉承历届理念&#xff0c;以“自由、协作、创新”为理…

Kotlin协程简介

文章目录 官网地址1 来源2 出现原因3 特点1&#xff09; 轻量2&#xff09;内存泄漏更少3&#xff09;内置取消支持4&#xff09;Jetpack 集成 4 依赖项信息5 在后台线程中执行6 使用协程确保主线程安全7 处理异常 官网地址 https://developer.android.google.cn/kotlin/corou…

挑战分布式架构,解密Java在业务场景下的高效应用面试题目介绍常用的通信方式有哪些请分别说明各自的特点和适用场景

本人详解 作者:王文峰,参加过 CSDN 2020年度博客之星,《Java王大师王天师》作者 公众号:山峯草堂,非技术多篇文章,专注于天道酬勤的 Java 开发问题、中国国学、传统文化和代码爱好者的程序人生,期待你的关注和支持!本人外号:神秘小峯 转载说明:务必注明来源(注明:…

Python-Python高阶技巧:闭包、装饰器、设计模式、多线程、网络编程、正则表达式、递归

版本说明 当前版本号[20231018]。 版本修改说明20231018初版 目录 文章目录 版本说明目录Python高阶技巧闭包简单闭包修改外部函数变量的值实现以下atm取钱的闭包实现了闭包注意事项 装饰器装饰器的一般写法&#xff08;闭包写法&#xff09;装饰器的语法糖写法 设计模式单例…

x86 架构的机载计算机,它来了!

Allspark 2-x86采用Intel酷睿11代或12代CPU&#xff0c;x86架构&#xff0c;适用于无人机等机器人运行SLAM、VIO等复杂逻辑和高精度的机器视觉任务。预装 Ubuntu 22.04或Windows 11&#xff0c;满足多种使用场景。 市面上现有的一些NUC产品&#xff0c;不仅没有针对移动机器人使…

【数据结构】线性表(三)循环链表的各种操作(创建、插入、查找、删除、修改、遍历打印、释放内存空间)

目录 线性表的定义及其基本操作&#xff08;顺序表插入、删除、查找、修改&#xff09; 四、线性表的链接存储结构 1. 单链表 2. 循环链表 a. 循环链表节点结构体 b. 创建新节点 c. 在循环链表末尾插入节点 d. 删除循环链表中指定值的节点 e. 在循环链表中查找指定值的…

Java 中用的是值传递还是引用传递?

值传递&#xff08;Pass by Value&#xff09;和引用传递&#xff08;Pass by Reference&#xff09;是两种参数传递的方式。 值传递是指在调用函数或方法时&#xff0c;将参数的值复制给一个临时变量然后传递给函数或方法。在函数或方法内部&#xff0c;对参数进行修改不会影…

HZOJ-271: 滑动窗口

题目描述 ​ 给出一个长度为 N&#xfffd; 的数组&#xff0c;一个长为 K&#xfffd; 的滑动窗口从最左移动到最右&#xff0c;每次窗口移动&#xff0c;如下图&#xff1a; 找出窗口在各个位置时的极大值和极小值。 输入 ​ 第一行两个数 N,K&#xfffd;,&#xfffd;。 …

linux部署gitlab

1. 配置yum源&#xff1a; vim /etc/yum.repos.d/gitlab-ce.repo [gitlab-ce] nameGitlab CE Repository baseurlhttps://mirrors.tuna.tsinghua.edu.cn/gitlab-ce/yum/el$releasever/ gpgcheck0 enabled1 2. 更新本地缓存 sudo yum install -y gitlab-ce 3. 安装相关依赖 yum …

白水三佳电脑ERP部署

安装宝塔面板&#xff0c;有这个方便很多&#xff0c;可以省下3天的环境部署时间。 移动端&#xff0c; 先取移动版的压缩包&#xff0c;上传至服务器/www/wwwroot/目录下面&#xff0c;直接解压到当前目录后会生成/www/wwwroot/m/的目录&#xff0c;移动版就在这里面了。以下…

Python3----------抽象(多态、封装、继承等)

术语多态&#xff08;polymorphism&#xff09;源自希腊语&#xff0c;意思是“有多种形态”。这大致意味着即便你不知道变量指向的是哪种对象&#xff0c;也能够对其执行操作&#xff0c;且操作的行为将随对象所属的类型&#xff08;类&#xff09;而异。 封装&#xff08;en…

基于深度优先搜索的图遍历

这里写目录标题 基于深度优先搜索的无向图遍历算法流程图Python实现Java实现 基于深度优先搜索的有向图遍历Python实现 基于深度优先搜索的无向图遍历 使用深度优先搜索遍历无向图&#xff0c;将无向图用邻接表存储&#xff1a; 算法流程图 初始化起点 source&#xff0c;当…

RAID2.0优势

一、定义 RAID2.0技术将硬盘域中的硬盘空间切分成固定大小的物理空间-CK&#xff08;Chunk64M&#xff09;&#xff0c;实现底层虚拟化&#xff0c;不同硬盘的多个CK组成存储池&#xff0c;相同类型的CK按照RAID策略组成&#xff08;CKG&#xff09;&#xff0c;CKG将再次切分成…

机器人控制算法——TEB算法—Obstacle Avoidance and Robot Footprint Model(避障与机器人足迹模型)

1.How Obstacle Avoidance works 1.1处罚条款 避障是作为整体轨迹优化的一部分来实现的。显然&#xff0c;优化涉及到找到指定成本函数&#xff08;目标函数&#xff09;的最小成本解&#xff08;轨迹&#xff09;。简单地说&#xff1a;如果一个计划的&#xff08;未来&…

Jmeter测试关联接口

Jmeter用于接口测试时&#xff0c;后一个接口经常需要用到前一次接口返回的结果&#xff0c;本文主要介绍jmeter通过正则表达式提取器来实现接口关联的方式&#xff0c;可供参考。 一、实例场景&#xff1a; 有如下两个接口&#xff0c;通过正则表达式提取器&#xff0c;将第一…

【ROS 2 基础-常用工具】-6 Rviz基础使用

所有内容请查看&#xff1a;博客学习目录_Howe_xixi的博客-CSDN博客

chromium 52 chrome 各个版本发布功能列表(58-84)

chromium Features 58-84 From https://chromestatus.com/features chromium58 Features:41 ‘allow-top-navigation-by-user-activation’ <iframe sandbox> keyword Adds a new keyword named “allow-top-navigation-by-user-activation” for iframe sandbox, wh…

整理uvc驱动相关函数的调用流程

目录 1、uvc_video.c初始化函数的调用关系 2、uvc_queue.c3、uvc_v4l2.c4、v4l2-core5、数据传输1、分配一个gadget请求2、请求一个queue 1、uvc_video.c // uvc_video.c uvc_video_encode_header uvc_video_encode_data uvc_video_encode_bulk uvc_video_encode_isoc uvcg_vi…

【算法|前缀和系列No.4】leetcode238. 除自身以外数组的乘积

个人主页&#xff1a;兜里有颗棉花糖 欢迎 点赞&#x1f44d; 收藏✨ 留言✉ 加关注&#x1f493;本文由 兜里有颗棉花糖 原创 收录于专栏【手撕算法系列专栏】【leetcode】 &#x1f354;本专栏旨在提高自己算法能力的同时&#xff0c;记录一下自己的学习过程&#xff0c;希望…

Linux UWB Stack实现——MCPS帧处理

MCPS帧处理 用于处理IEEE 802.15.4中的相关帧&#xff0c;Frame Processing&#xff0c;简写为&#xff1a;fproc。 在实现中&#xff0c;维护了关于帧处理的有限状态机(FSM)。本文从帧处理的数据结构和部分典型处理实现上进行简要的介绍。 1. 数据结构定义 关于帧处理状态…