GSVA,GSEA,KEGG,GO学习

目录

GSVA

1:获取注释基因集

2:运行

GSEA

1,示例数据集

2,运行

GSEA_KEGG富集分析

GSEA_GO富集分析

DO数据库GSEA

MSigDB数据库选取GSEA

KEGG

1:运行

2:绘图

bar图

气泡图

绘图美化

GO


GSVA

1:获取注释基因集

2:运行

GSEA

1,示例数据集

2,运行

GSEA_KEGG富集分析

GSEA_GO富集分析

DO数据库GSEA

MSigDB数据库选取GSEA

KEGG

GO


GSVA

【精选】RNA 18. SCI 文章中基因集变异分析 GSVA_gsva分析-CSDN博客

RNA-seq入门实战(八):GSVA——基因集变异分析 - 知乎 (zhihu.com)

表达矩阵反映了样本和基因的关系,则GSVA将一个“样本×基因”的矩阵转化为“样本×通路”的矩阵,直接反映了样本和读者感兴趣的通路之间的联系。因此,如果用limma包做差异表达分析可以寻找样本间差异表达的基因,同样地,使用limma包对GSVA的结果(依然是一个矩阵)做同样的分析,则可以寻找样本间有显著差异的通路。这些“差异表达”的通路,相对于基因而言,更加具有生物学意义,更具有可解释性,是统计学与生物学成功结合后,对GSEA结果的一次升华,可以进一步用于肿瘤subtype的分型等等与生物学意义结合密切的探究。

1:获取注释基因集

可以用msigdbr包下载读取或者GSEA | MSigDB | Human MSigDB Collections (gsea-msigdb.org)下载整理。 按需下载整理

##进行数据读取
geneSets <- getGmt('c2.cp.kegg_medicus.v2023.2.Hs.symbols.gmt')    ###下载的基因集

##下载symbols的 注释文件,使用表达矩阵是省略了转换的麻烦
2:运行
rm(list = ls())  ## 魔幻操作,一键清空~
options(stringsAsFactors = F)

library(GSVA)
library(GSEABase)
library(clusterProfiler)

expr <- read.csv("easy_input_expr.csv", row.names = 1)#表达矩阵
geneSets <- getGmt('h.all.v2023.2.Hs.symbols.gmt')##symbols 较为方便

##运行
GSVA_hall <- gsva(expr=as.matrix(expr),#需要为matrix格式
                  gset.idx.list=geneSets, 
                  #   method="gsva", #c("gsva", "ssgsea", "zscore", "plage")
                  mx.diff=T, # 数据为正态分布则T,双峰则F
                  kcdf="Gaussian", #CPM, RPKM, TPM数据就用默认值"Gaussian", read count数据则为"Poisson",
                  parallel.sz=4,# 并行线程数目
                  min.sz=2) 

GSEA

快速拿捏KEGG/GO/Reactome/Do/MSigDB的GSEA富集分析! (qq.com)

【精选】RNA 11. SCI 文章中基因表达富集之 GSEA_gsea数据库-CSDN博客

GSEA(Gene Set EnrichmentAnalysis),即基因集富集分析,它的基本思想是使用预定义的基因,将基因按照在两类样本中的差异表达程度排序,然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。

1,示例数据集
rm(list = ls()) 
library(DOSE)##包里有测试数据
library(clusterProfiler)
require(enrichplot)
options(stringsAsFactors = F)

data(geneList, package = "DOSE")
head(geneList)

##DOSE提供的一个geneList,name是每一个entrez gene id, value是log2FoldChange值。
df <- as.data.frame(geneList)##查看:还原gene symbol
df$ENTREZID <- rownames(df)
#?bitr 函数查看
df1<-bitr(df$ENTREZID, #转换的列是df数据框中的SYMBOL列
          fromType = "ENTREZID",#需要转换ID类型
          toType = "SYMBOL",#转换成的ID类型
          OrgDb = "org.Hs.eg.db")#物种选择(小鼠的是org.Mm.eg.db)
df2<-merge(df,df1,by="ENTREZID",all=F)##进行合并

测试数据查看

2,运行
GSEA_KEGG富集分析
rm(list = ls()) 
library(DOSE)##包里有测试数据
library(clusterProfiler)
require(enrichplot)
options(stringsAsFactors = F)
##input需要的是entrez ID+log2fc文件(包里的文件即为entrez ID+log2fc)
##如果不是entrez ID+log2fc,需要进行转换

##数据集查看和整理
data(geneList, package = "DOSE")
head(geneList)
#4312     8318    10874    55143    55388      991 
#4.572613 4.514594 4.418218 4.144075 3.876258 3.677857 
df <- as.data.frame(geneList)##查看:还原gene symbol
df$ENTREZID <- rownames(df)
#?bitr 函数查看
df1<-bitr(df$ENTREZID, #转换的列是df数据框中的SYMBOL列
          fromType = "ENTREZID",#需要转换ID类型
          toType = "SYMBOL",#转换成的ID类型
          OrgDb = "org.Hs.eg.db")#物种选择(小鼠的是org.Mm.eg.db)
df2<-merge(df,df1,by="ENTREZID",all=F)##进行合并
df3 <-df2[,-1] 
head(df3)#自己分析常见的格式
#geneList SYMBOL
#1 -0.28492113   NAT2

##以df3 为input进行转换:添加entrez ID列,symbol转entrez ID
##注意还原结果少几十个基因:因为一个entrez ID可能对应多个symbol(一基因多symbol)
dat<-bitr(df3$SYMBOL, 
          fromType = "SYMBOL", #现有的ID类型
          toType = "ENTREZID",#需转换的ID类型
          OrgDb = "org.Hs.eg.db")#物种
head(dat)##转换时部分SYMBOL会转换失败

dat1<-merge(df3,dat,by="SYMBOL",all=F)##进行合并
#按照foldchange排序
sortdf<-dat1[order(dat1$geneList, decreasing = T),]#这里geneList其实是logFC值
head(sortdf)

geneList1 <- sortdf$geneList##先把foldchange按照从大到小提取出来
names(geneList1) <- sortdf$ENTREZID###给上面提取的foldchange加上对应上ENTREZID
head(geneList1 )
#4312     8318    10874    55143    55388      991 
#4.572613 4.514594 4.418218 4.144075 3.876258 3.677857  

#GSEA_KEGG富集分析:
KEGG_ges <- gseKEGG(
  geneList = geneList1,#input
  organism = "hsa")#物种

#按照enrichment score从高到低排序,便于查看富集通路
sortKEGG_ges<-KEGG_ges[order(KEGG_ges$enrichmentScore, decreasing = T),]
#sortKEGG_ges <- KEGG_ges@result

说明:在GSEA中,基因集的富集分数可以为正或负,表示该基因集在对应生物条件下的富集程度。富集分数的绝对值越大,表示富集程度越高。 对于AB两个亚型的差异基因,在GSEA富集分析中,结果按照富集分数排序。通常情况下,富集分数为正最大的前几个表示在A亚型中富集的集中的通路,而富集分数负值最大的几个表示在B亚型中富集的集中的通路。 需要注意的是,富集分数的正负并不代表富集的方向,而是表示富集程度的大小。因此,正值最大的前几个通路表示在A亚型中富集程度最高的通路,负值最大的几个通路表示在B亚型中富集程度最高的通路。 这样的排序方式可以帮助我们理解不同亚型或条件下基因集的富集模式,以及与这些通路相关的生物学过程和功能。

#进行绘图
gseaplot2(KEGG_ges, row.names(sortKEGG_ges)[1:5])##可以自行选择通路

dev.off()
#个性化展示(选取结果中的ID)
p1 <- gseaplot2(KEGG_ges,
                geneSetID = c("hsa03030","hsa03050","hsa04710","hsa00350"),#通路
                color = c("#003399", "#FFFF00", "#FF6600","black"),#颜色
                pvalue_table = TRUE,#显示P值
                ES_geom = "line")#"dot"将线转换为点
p1

GSEA_GO富集分析

主要是函数的不同 gseGO 函数

#GSEA
rm(list = ls()) 
library(DOSE)##包里有测试数据
library(clusterProfiler)
require(enrichplot)
options(stringsAsFactors = F)

##数据集查看和整理
data(geneList, package = "DOSE")
head(geneList)

##GSEA_GO富集分析:
GO_ges <- gseGO(geneList = geneList,
                OrgDb = "org.Hs.eg.db",
                ont = "CC", #one of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
                minGSSize = 10,
                maxGSSize = 500,
                pvalueCutoff = 0.05,
                eps = 0,
                verbose = FALSE)
#res <- GO_ges@result
res1<-GO_ges[order(GO_ges$enrichmentScore, decreasing = T),]
DO数据库GSEA

需要library(DOSE)   DO(Disease Ontology)数据库GSEA

#GSEA
rm(list = ls()) 
library(DOSE)##包里有测试数据
library(clusterProfiler)
require(enrichplot)
options(stringsAsFactors = F)

##数据集查看和整理
data(geneList, package = "DOSE")
head(geneList)

##GSEA_DO(Disease Ontology)富集分析:
DO_ges <- gseDO(geneList,
                minGSSize = 10,
                maxGSSize = 500,
                pvalueCutoff = 0.05,
                pAdjustMethod = "BH",
                verbose = FALSE,
                eps = 0)

#res <- DO_ges@result
res1<-DO_ges[order(DO_ges$enrichmentScore, decreasing = T),]
MSigDB数据库选取GSEA

msigdf + clusterProfiler全方位支持MSigDb (guangchuangyu.github.io)

Q&A | 如何使用clusterProfiler对MSigDB数据库进行富集分析 - 知乎 (zhihu.com)

#GSEA
rm(list = ls()) 
library(DOSE)##包里有测试数据
library(clusterProfiler)
require(enrichplot)
library(msigdbr)
options(stringsAsFactors = F)

##数据集查看和整理
data(geneList, package = "DOSE")
head(geneList)

##msigdbr 提取注释自己所需的注释基因集
H <- msigdbr(species = "Homo sapiens", 
              category = "H") %>% 
  dplyr::select(gs_name, entrez_gene)

#C2all <- msigdbr(species = "Homo sapiens", 
#              category = "H")#完整的注释

##富集分析
H_ges <- GSEA(geneList,
               TERM2GENE = H,##注释基因集
               minGSSize = 10,
               maxGSSize = 500,
               pvalueCutoff = 0.05,
               pAdjustMethod = "BH",
               verbose = FALSE,
               eps = 0)

#res <- H_ges@result
res1<-H_ges[order(H_ges$enrichmentScore, decreasing = T),]

KEGG

KEGG富集分析及可视化,一把子拿捏! (qq.com)

RNA 10. SCI 文章中基因表达富集之 KEGG 注释_kegg中qvalue_桓峰基因的博客-CSDN博客

  在线KEGGDAVID Functional Annotation Bioinformatics Microarray Analysis (ncifcrf.gov)

DAVID 在线数据库进行 GO/ KEGG 富集分析_david数据库go富集分析-CSDN博客

1:运行
##差异基因KEGG富集分析
rm(list = ls())  ## 魔幻操作,一键清空~
options(stringsAsFactors = F)
library(dplyr)#数据清洗
library(org.Hs.eg.db)#ID转换
library(clusterProfiler)#富集分析
library(ggplot2)#绘图
library(RColorBrewer)#配色调整
library(DOSE)##包里有测试数据

data(geneList, package = "DOSE")
head(geneList)
df <- as.data.frame(geneList)##查看:还原gene symbol
df$ENTREZID <- rownames(df)
#?bitr 函数查看
df1<-bitr(df$ENTREZID, #转换的列是df数据框中的SYMBOL列
          fromType = "ENTREZID",#需要转换ID类型
          toType = "SYMBOL",#转换成的ID类型
          OrgDb = "org.Hs.eg.db")#物种选择(小鼠的是org.Mm.eg.db)
df2<-merge(df,df1,by="ENTREZID",all=F)##进行合并

##选择logFC>1.5的基因:我们以这个筛选的差异基因集为input测试
df3 <- df2[abs(df2$geneList)>1.5,]##abs() 表示绝对值
##自己使用时进行自定义

KEGG_diff <- enrichKEGG(gene = df3$ENTREZID,
                        organism = "hsa",#物种,Homo sapiens (human)
                        pvalueCutoff = 0.05,
                        qvalueCutoff = 0.05)

KEGG_result <- KEGG_diff@result

#保存富集结果:
save(KEGG_diff,KEGG_result,file = c("KEGG_diff.Rdata"))
#write.csv(KEGG_result,file = "KEGG_result.csv")

ID:pathway的ID名;GeneRatio:差异基因中富集到该pathway的基因数目/富集到所有pathway的总差异基因数目;BgRatio:所有背景基因中富集到该pathway的基因数目/总背景基因数目;

Count:富集到该pathway的基因数目

2:绘图
bar图
#绘图
library(enrichplot)
library(stringr)
library(cowplot)
library(ggplot2)
barplot(KEGG_diff,
        x = "Count", #or "GeneRatio"
        color = "pvalue", #or "p.adjust" and "qvalue"
        showCategory = 20,#显示前top20
        font.size = 12,
        title = "KEGG enrichment barplot",
        label_format = 30 #超过30个字符串换行
)

气泡图
dotplot(
  KEGG_diff,
  x = "GeneRatio",
  color = "p.adjust",
  title = "Top 20 of Pathway Enrichment",
  showCategory = 20,
  label_format = 30
)

绘图美化
#将pathway按照p值排列
KEGG_top20 <- KEGG_result[1:20,]
KEGG_top20$pathway <- factor(KEGG_top20$Description,levels = rev(KEGG_top20$Description))
p2 <- ggplot(data = KEGG_top20,
             aes(x = Count,y = pathway))+
  geom_point(aes(size = Count,
                 color = -log10(pvalue)))+ # 气泡大小及颜色设置
  theme_bw()+
  scale_color_distiller(palette = "Spectral",direction = 1) +
  labs(x = "Gene Number",
       y = "",
       title = "Dotplot of Enriched KEGG Pathways",
       size = "Count")
p2

GO

GO富集分析及可视化,一把子拿捏! (qq.com)

【精选】RNA 9. SCI 文章中基因表达之 GO 注释_consider increasing max.overlaps_桓峰基因的博客-CSDN博客

就是函数的差别

GO_MF_diff <- enrichGO(gene = diff_entrez$ENTREZID, #用来富集的差异基因
OrgDb = org.Hs.eg.db, #指定包含该物种注释信息的org包
ont = "MF", #可以三选一分别富集,或者"ALL"合并
pAdjustMethod = "BH", #多重假设检验矫正方法
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE) #是否将gene ID映射到gene name
#提取结果表格:
GO_MF_result <- GO_MF_diff@result
View(GO_MF_result)

感谢上面的许多教程:更详细的大家可以去学习!

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

springBoot 配置druid多数据源 MySQL+SQLSERVER

1:pom 文件引入数据 <dependency> <groupId>com.alibaba</groupId> <artifactId>druid-spring-boot-starter</artifactId> <version>1.1.0</version> </dependency>…

mysql服务器数据同步

在Linux和Windows之间实现MySQL服务器数据的同步。下面是一些常见的方法和工具&#xff1a; 复制&#xff08;Replication&#xff09;&#xff1a;MySQL复制是一种常见的数据同步技术&#xff0c;可用于将一个MySQL服务器的数据复制到其他服务器。您可以设置主服务器&#xff…

CMSIS-RTOS在stm32使用

目录&#xff1a; 一、安装和配置CMSIS_RTOS.1.打开KEIL工程&#xff0c;点击MANAGE RUN-TIME Environment图标。2.勾选CMSIS CORE和RTX.3.配置RTOS 时钟频率、任务栈大小和数量&#xff0c; 软件定时器. 二、CMSIS_RTOS内核启动和创建线程。1.包含头文件。2.内核初始化和启动。…

C#,数值计算——插值和外推,曲线插值(Curve_interp)的计算方法与源程序

1 文本格式 using System; namespace Legalsoft.Truffer { /// <summary> /// Object for interpolating a curve specified by n points in dim dimensions. /// </summary> public class Curve_interp { private int dim { get; s…

openGauss通过VIP实现的故障转移

&#x1f4e2;&#x1f4e2;&#x1f4e2;&#x1f4e3;&#x1f4e3;&#x1f4e3; 哈喽&#xff01;大家好&#xff0c;我是【IT邦德】&#xff0c;江湖人称jeames007&#xff0c;10余年DBA及大数据工作经验 一位上进心十足的【大数据领域博主】&#xff01;&#x1f61c;&am…

VisualGDB 6.0 R2 Crack

轻松跨平台"VisualGDB 使 Visual Studio 的跨平台开发变得简单、舒适。它支持&#xff1a; 准系统嵌入式系统和物联网模块&#xff08;查看完整列表&#xff09; C/C Linux 应用程序 本机 Android 应用程序和库 Raspberry Pi 和其他Linux 板 Linux 内核模块&#xff08;单…

【PTA题目】6-13 求叠数(递归版) 分数 10

6-13 求叠数(递归版) 分数 10 全屏浏览题目 切换布局 作者 李祥 单位 湖北经济学院 请编写递归函数&#xff0c;生成叠数。 例如&#xff1a;Redup(5,8)88888 函数原型 long long Redup(int n, int d); 说明&#xff1a;参数 n 为重复次数(非负整数)&#xff0c;d 为数字…

未来科技中的云计算之路

随着科技的不断发展&#xff0c;云计算已经不再是一个陌生的词汇&#xff0c;而是我们日常生活中不可或缺的一部分。从智能家居到无人驾驶&#xff0c;再到虚拟现实和人工智能&#xff0c;云计算在这些领域都扮演着至关重要的角色。在这篇博客中&#xff0c;我们将一同探索云计…

【如何学习Python自动化测试】—— 页面元素定位

接上篇自动化测试环境搭建&#xff0c;现在我们介绍 webdriver 对浏览器操作的 API。 2、 页面元素定位 通过自动化操作 web 页面&#xff0c;首先要解决的问题就是定位到要操作的对象&#xff0c;比如要模拟用户在页面上的输入框中输入一段字符串&#xff0c;那就必须得定位到…

UiPath Studio 2023.10 Crack

UiPath Studio是一款功能强大且用户友好的集成开发环境 (IDE)&#xff0c;专为机器人流程自动化 (RPA) 设计。它由自动化技术领域的领先公司UiPath开发。 以下是 UiPath Studio 的一些主要功能和组件&#xff1a; 图形用户界面 (GUI)&#xff1a;UiPath Studio 具有直观且用户友…

RT-Thread STM32F407 BMI088--SPI

BMI088是一款高性能6轴惯性传感器&#xff0c;由16位数字三轴24g加速度计和16位数字三轴2000/ s陀螺仪组成。 这里用SPI来驱动BMI088进行数据解读 第一步&#xff0c;首先在 RT-Thread Settings中进行配置 第二步&#xff0c;退出RT-Thread Settings&#xff0c;进入board.h…

数模建模竞赛——写作手三天速成(文末领取)

目录 第一天&#xff1a;准备论文模板&#xff0c;学习各类基础画图技巧 1、论文模板 2、基础画图能力 第二天&#xff1a;看按模型算法分类的优秀论文&#xff0c;学习其模型的写作方式 第三天&#xff1a;配合团队完成真题练习 第一天&#xff1a;准备论文模板&#xff…

【网络通信】探索UDP与TCP协议、IP地址和端口号的奥妙

&#x1f33a;个人主页&#xff1a;Dawn黎明开始 &#x1f380;系列专栏&#xff1a;网络奇幻之旅 ⭐每日一句&#xff1a;往前走&#xff0c;朝着光 &#x1f4e2;欢迎大家&#xff1a;关注&#x1f50d;点赞&#x1f44d;评论&#x1f4dd;收藏⭐️ 文章目录 &#x1f4cb;前…

嵌入式 Linux 移植与系统启动方法

1、Linux系统启动与U-Boot 所谓移植就是把程序代码从一种运行环境转移到另一种运行环境。对于内核移植来说&#xff0c;主要是从一种硬件平台转移到另一种硬件平台上运行。 体系结构级别的移植是指在不同体系结构平台上Linux内核的移植&#xff0c;例如&#xff0c;在ARM、MI…

【2023春李宏毅机器学习】生成式学习的两种策略

文章目录 1 各个击破2 一步到位3 两种策略的对比 生成式学习的两种策略&#xff1a;各个击破、一步到位 对于文本生成&#xff1a;把每一个生成的元素称为token&#xff0c;中文当中token指的是字&#xff0c;英文中的token指的是word piece。比如对于unbreakable&#xff0c;他…

【docker】iptables实现NAT

iptables是一个Linux内核中的防火墙工具&#xff0c;可以被用来执行各种网络相关的任务&#xff0c;如过滤、NAT和端口转发等&#xff0c;可以监控、过滤和重定向网络流量。 iptables可以用于以下应用场景&#xff1a; 网络安全&#xff1a;iptables可以过滤网络流量&#xf…

潇洒郎: 小白一次性成功——小米红米手机解BL锁+ ROOT-刷面具

一、账号与设备绑定 手机登录账号,绑定账号,使用手机卡流量,等待7天后解BL锁。 二、解BL锁 下载工具 申请解锁小米手机 (miui.com) https://www.miui.com/unlock/index.html 1、登录账号-与绑定的账号一样 2、驱动检测安装 驱动安装进入Fastboot模式后,会自动识别已连接…

【数据结构】树与二叉树(二十):树获取大儿子、大兄弟结点的算法(GFC、GNB)

文章目录 5.1 树的基本概念5.1.1 树的定义5.1.2 森林的定义5.1.3 树的术语 5.2 二叉树5.3 树5.3.1 树的存储结构1. 理论基础2. 典型实例3. Father链接结构4. 儿子链表链接结构5. 左儿子右兄弟链接结构 5.3.2 获取结点的算法1. 获取大儿子结点的算法&#xff08;GFC&#xff09;…

Linux-top命令解释

Linux-top命令解释 常用参数查看所有逻辑核的运行情况&#xff1a;1查看指定进程的情况&#xff1a;-p pid显示进程的完整命令&#xff1a;-c 面板指标解释第一行top第二行tasks第三行%Cpu第四行Mem第五行Swap第六行各进程监控PID&#xff1a;进程IDUSER&#xff1a;进程所有者…

“流量为王”的时代一去不返!如何押注互联网下一个黄金十年

目录 1“流量为王”的时代一去不返&#xff01;如何押注互联网下一个黄金十年 2AI夺走的第一份工作竟是OpenAI CEO&#xff1f;阿尔特曼被“扫地出门”&#xff0c;网友热评&#xff1a;是被GPT-5取代了吗&#xff1f;马斯克更“毒”&#xff0c;挂出求职申请链接 3GPT-4V新玩…