R绘图实战|GSEA富集分析图

最新推荐文章于 2024-08-16 20:34:54 发布
最新推荐文章于 2024-08-16 20:34:54 发布
阅读量2.2w 收藏 231
点赞数 31
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_45822007/article/details/114989817
版权

在这里插入图片描述

写在前面

后台难得有读者私信，请教了下图中文章的GSEA图能不能用R来画，今天就来简单写个教学。

title figure_gsea

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

GSEA 和GO、KEGG pathway不同的地方在于，后两者会提前设定一个阈值，只关注差异变化大的基因（相当于重点班）。这样子容易遗漏部分差异表达不显著却有重要生物学意义的基因（成绩一般，但是很有天赋）。所以GSEA分析比较适用于，传统分析方法筛选后样本过少的数据集。

数据准备

直接用之前的转录组差异分析后的数据来演示，数据格式如下：

gene_diff

至于基因差异分析怎么l做，有好多种办法，常用的有DESeq2、EdgeR、limma，有问题的话可以私信~因为GSEA只需要SYMBOL（基因名）和foldchange （或logFC)两列，所以可以把不需要的删掉。（以上操作在EXCEL，或者用R的tidyverse（数据框处理可以参考文章： ）进行操作都可以，怎么方便怎么来）

gene_arrage

开始分析

安装并导入要用到的R包

BiocManager::install("clusterProfiler") #感谢Y叔的clusterprofiler包
BiocManager::install("enrichplot")  #画图需要
BiocManager::install("org.Hs.eg.db") #基因注释需要
library(org.Hs.eg.db)
library(clusterProfiler)
library(enrichplot)

导入数据

setwd("D:/Note/MZBJ/Q_A") #设置文件所在位置
df = read.table("gene_diff.txt",header = T) #读入txt
# df = read.csv("gene_diff.csv",header = T) #读入csv
head(df)#查看前面几行
dim(df)#数据总共几行几列

> head(df)#查看前面几行
        SYMBOL    logFC
1         CD74 41.99218
2      MAB21L3 35.00852
3     KCNQ1OT1 22.78417
4 RP3-323A16.1 22.25173
5    LINC00504 16.82801
6       MALAT1 16.64222
> dim(df)#数据总共几行几列
[1] 5057    2

转换基因ID

如基因名是symbol，需要将基因ID转换为Entrez ID格式。Entrez ID实际上是指的Entrez gene ID，是对应于染色体上一个gene location的。每一个发现的基因都会被编制一个统一的编号，而Entrez ID是指的来自于NCBI旗下的Entrez gene数据库所使用的编号。因为Entrez ID具有特异性，所以后续分析更适合用Entrez ID。

df_id<-bitr(df$SYMBOL, #转换的列是df数据框中的SYMBOL列
            fromType = "SYMBOL",#需要转换ID类型
            toType = "ENTREZID",#转换成的ID类型
            OrgDb = "org.Hs.eg.db")#对应的物种，小鼠的是org.Mm.eg.db
>'select()' returned 1:many mapping between keys and columns
Warning message:
In bitr(df$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db") :
  7.87% of input gene IDs are fail to map...  #7.87%没有比对到就是没有转换成功

把两个数据框df 和 df_id根据SYMBOL列合并。

df_all<-merge(df,df_id,by="SYMBOL",all=F)#使用merge合并
head(df_all) #再看看数据
dim(df_all) #因为有一部分没转换成功，所以数量就少了。

> head(df_all)
    SYMBOL        logFC ENTREZID
1      A2M -0.713519723        2
2     AAK1 -0.089497971    22848
3     AAMP -0.014536797       14
4    AARS2  0.077105219    57505
5 AASDHPPT -0.000560858    60496
6    ABCA1  0.436678052       19
> dim(df_all)
[1] 4660    3

GAEA

df_all_sort <- df_all[order(df_all$logFC, decreasing = T),]#先按照logFC降序排序
gene_fc = df_all_sort$logFC #把foldchange按照从大到小提取出来
head(gene_fc)
names(gene_fc) <- df_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(gene_fc)

> head(gene_fc)
[1] 41.99218 35.00852 22.78417 16.82801 16.64222 15.33221
> head(gene_fc)
     972   126868    10984   201853   378938     3514 
41.99218 35.00852 22.78417 16.82801 16.64222 15.33221

准备以上的东西，接下来一行代码解决。

#以KEGG Pathway示例
KEGG <- gseKEGG(gene_fc, organism = "hsa") #具体参数在下面

> KEGG <- gseKEGG(gene_fc, organism = "hsa")
Reading KEGG annotation online:

Reading KEGG annotation online:

preparing geneSet collections...
GSEA analysis...
leading edge analysis...
done...
Warning messages:
1: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
  There are ties in the preranked stats (0.13% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.
2: In serialize(data, node$con) : 载入时'package:stats'可能无用
3: In serialize(data, node$con) : 载入时'package:stats'可能无用
4: In serialize(data, node$con) : 载入时'package:stats'可能无用
5: In serialize(data, node$con) : 载入时'package:stats'可能无用
6: In serialize(data, node$con) : 载入时'package:stats'可能无用
7: In serialize(data, node$con) : 载入时'package:stats'可能无用
8: In fgseaMultilevel(...) :
  For some pathways, in reality P-values are less than 1e-10. You can set the `eps` argument to zero for better estimation.

如果要做GO富集呢？

#GO富集
GO <- gseGO(
  gene_fc, #gene_fc
  ont = "BP",# "BP"、"MF"和"CC"或"ALL"
  OrgDb = org.Hs.eg.db,#人类注释基因
  keyType = "ENTREZID",
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",#p值校正方法
)
#KEGG富集
gseKEGG(
  geneList,
  organism = "hsa",
  keyType = "kegg",
  exponent = 1,
  minGSSize = 10,
  maxGSSize = 500,
  eps = 1e-10,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  verbose = TRUE,
  use_internal_data = FALSE,
  seed = FALSE,
  by = "fgsea",
  ...
)

head(KEGG)#看一下这个文件
> head(KEGG)
               ID                    Description setSize enrichmentScore       NES
hsa03010 hsa03010                       Ribosome      99      -0.8707285 -2.370839
hsa05152 hsa05152                   Tuberculosis      87       0.8678558  1.786981
hsa05171 hsa05171 Coronavirus disease - COVID-19     142      -0.5976011 -1.704522
hsa04512 hsa04512       ECM-receptor interaction      19      -0.8866402 -1.913989
               pvalue     p.adjust      qvalues rank                  leading_edge
hsa03010 0.0000000001 0.0000000257 2.431579e-08  289 tags=65%, list=6%, signal=62%
hsa05152 0.0002124294 0.0272971804 2.582695e-02  279 tags=30%, list=6%, signal=29%
hsa05171 0.0004376904 0.0290749106 2.750893e-02  289 tags=46%, list=6%, signal=45%
hsa04512 0.0004525278 0.0290749106 2.750893e-02  250 tags=58%, list=5%, signal=55%
                                                                                                                                                                                                                                                                                                                                        core_enrichment
hsa03010           6231/6193/4736/6235/2197/6218/6166/6167/6157/3921/6129/140801/6152/6125/6169/6124/9349/6141/6138/6187/6228/6144/6135/6202/6155/6154/6132/6160/6159/6147/6156/6210/6230/6175/6122/6128/11224/23521/9045/25873/6161/6201/6208/6189/6181/6188/6133/6165/6194/6139/6168/6224/6143/6142/6222/6164/6176/6232/6206/6223/6171/6233/6134/6137
hsa05152

简单解释一下结果的意思：

  • ID 代表KEGG中的信号通路
  • Description 对信号通路的描述
  • setSize 该信号通路的基因个数
  • enrichmentScore 富集分数，也就是ES
  • NES 标准化以后的ES，全称normalized enrichment score、
  • qvalues ，或者说FDR q-val（false discovery rate）错误发现率
  • rank 排名
  • core_enrichment，富集该目的通路的基因列表。
sortKEGG<-KEGG[order(KEGG$enrichmentScore, decreasing = T),]#按照enrichment score从高到低排序
head(sortKEGG)
dim(sortKEGG)
write.table(sortKEGG,"gsea_sortKEGG.txt") #保存结果

结果可视化

#gseaplot2用法
gseaplot2(
  x, #gseaResult object，即GSEA结果
  geneSetID,#富集的ID编号
  title = "", #标题
  color = "green",#GSEA线条颜色
  base_size = 11,#基础字体大小
  rel_heights = c(1.5, 0.5, 1),#副图的相对高度
  subplots = 1:3, #要显示哪些副图 如subplots=c(1,3) #只要第一和第三个图，subplots=1#只要第一个图
  pvalue_table = FALSE, #是否添加 pvalue table
  ES_geom = "line" #running enrichment score用先还是用点ES_geom = "dot"
)

先来个常规的

gseaplot2(KEGG, "hsa05152", color = "firebrick", rel_heights=c(1, .2, .6))
hsa05152

再来个文章里的：

paths <- c("hsa03010", "hsa05152", "hsa05171", "hsa04512")#选取你需要展示的通路ID
gseaplot2(KEGG,paths, pvalue_table = TRUE)
pathsGSEA 
gseaplot2(KEGG,paths,color = colorspace::rainbow_hcl(4),subplots=c(1,2), pvalue_table = TRUE)
#换个颜色，只显示上面两个副图
pathsGSEA_rainbow

剩下的一点细节用AI处理一下就可以~

以上教程为本人的粗浅记录，如有错误还请各位读者指正如果觉得有用，希望大家多多点赞分享

代码领取：公众号后台私信“GSEA"即可领取全部代码。

往期内容：

《R数据科学》学习笔记|Note14:向量

《R数据科学》学习笔记|Note13:函数

《R数据科学》学习笔记|Note12:使用magrittr进行管道操作

木舟笔记
关注
r语言进行go富集分析_生信实操|如何利用R语言进行GSEA分析
weixin_32491317的博客
12-31 1万+
看图说话栏目曾介绍过GSEA的原理(看图说话|GSEA分析--教你解锁高级的富集分析)，今天我们来看一下如何利用R语言进行GSEA分析。如果你有RNA-seq的数据，就可以这样做，先把数据整成这样，一共两列，一列是SYMBOL，一列是foldChange。然后输入代码：>setwd("E:\\20201214GSEA分析怎么做")>#if(!requireNamesp...
GSEA富集分析可视化
医学和生信笔记的博客
09-05 6641
上一次推文已经给大家介绍了常见的富集分析类型以及如何使用全能的R包相必大家也知道的富集结果可以画出很多漂亮的图，其中enrichplot是专门用于对接的富集结果的可视化R包。但是大家有没有这样的困惑呢？ehrichplot今天就先给大家详细介绍enrichplot包。当然了，富集分析的可视化还有许多其他R包，我们会在后续的推文中继续介绍。
一文解决R语言GSEA分析及可视化
qq_42090739的博客
07-08 5403
富集分析是组学研究中绕不开的，常见的有GO、KEGG、GSEA等。GO、KEGG分析是基于差异基因的，也就是人为设定阈值，选定基因去做富集分析，其实这有很大的主观性，GO、KEGG分析及可视化我们之前讲过，可参考（转录组不求人系列(十三): GO、KEGG富集个性化作图）。GSEA则避免了这种主观性，是基于所有的基因进行分析的，将基因按照变化倍数从大到小排列，评估基因分布趋势，及其对表型的贡献。GSEA分析一般使用GSEA官网：GSEA提供的GSEA软件：GSEA | Login但是根据我自己的体验，GSE
单细胞irGSEA分析:整合多种富集分析方式的R包
最新发布
zfyyzhys的博客
08-16 769
irGSEA整合了多种基于单个细胞表达等级的富集分析方法(AUCell、UCell、singscore、ssGSEA、JASMINE和Viper)，并通过秩聚合算法(robust rank aggregation, RRA)对差异分析的结果进行评估，筛选出在这种几种方法中表现出相似的富集程度的差异基因集。通俗的来说，这个R包就是，获取在。
GSEA富集分析 - 界面操作
悟道西方
11-17 3万+
欢迎关注微信公众号生信宝典：http://mp.weixin.qq.com/s/3Nd3urhfRGkw-F0LGZrlZQGSEA定义Gene Set Enrichment Analysis (基因集富集分析)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的贡献。其输入数据包含两部分，一是已知功能的基因集 (可以是GO注释、MsigDB的注释或其它符合格
一文掌握GSEA富集分析-最详细教程
热门推荐
悟道西方
04-26 5万+
之前总结了一篇关于GSEA富集分析的推文——《GSEA富集分析 - 界面操作》，大略介绍了GSEA的定义、GSEA原理、GSEA分析、Leading-edge分析等，不太了解的朋友可以点击阅读先理解下概念。 最近用自己数据实战分析时用到了该方法，故将一些之前遗漏的点补充整理出来分享给大家。 从前文中我们了解到GSEA分析的目的是要判断S集基因（基于先验知识的基因注释信息）中的基因是随机分布还是聚集...
R语言实现GSEA分析（1）
weixin_49320263的博客
07-03 4908
R语言实现GSEA分析（1）
R语言：GSEA分析
weixin_69558614的博客
05-13 861
安装软件包> if (!#加载软件包#设置变量> head(rt)head(rt))\\-(.*?)\\-(.*?#排序#GESA分析dev.off()dev.off()一起学习交流。
R语言GSEA分析
weixin_47634487的博客
06-12 674
在R语言中进行基因集富集分析（Gene Set Enrichment Analysis, GSEA），你可以使用包，这是一个强大的工具，用于富集分析。以下是一个简单的例子，演示如何使用包进行GSEA分析。
GSEA绘图——R
ziixiaoshenwang的博客
09-24 1948
rm(list=ls()) options(stringsAsFactors = F) source("gseaplot_modified.r") path <- "E:/Users/Administrator/Desktop/GSEA/out/HT_LT_GO14.Gsea.1632210488415" replotGSEA(path,"CALMODULIN BINDING")
GSEA | 基因富集分析
生信练习生
04-24 1511
GSEA富集分析
gsea:用于基因组富集分析的R包
05-20
基因组富集分析（GSEA）是一种广泛应用的生物信息学方法，它可以帮助研究者理解基因表达数据中的生物学意义。在R编程环境中，`gsea`包提供了强大的工具来执行这种分析。本文将深入探讨`gsea`包的核心功能、安装与...
R语言绘制SCI科研多GSEA富集图源代码.zip
05-01
本资源提供了使用R语言绘制SCI科研多GSEA富集图的源代码，这对于进行基因功能注释和通路分析的研究人员来说是极有价值的工具。下面我们将详细讨论如何理解和应用这些源代码。 首先，了解R语言是至关重要的。R是一种...
多GSEA富集图_R语言绘制SCI图的输入文件及代码.rar
07-30
在本压缩包文件“多GSEA富集图_R语言绘制SCI图的输入文件及代码.rar”中，包含了使用R语言绘制多组GSEA富集结果的输入文件和源代码，这对于科研人员发表科学论文时展示数据分析结果非常有帮助。 首先，我们需要理解...
一文掌握单基因GSEA富集分析
BioinfoDu
02-03 1463
关于GSEA分析，我们在前期的教程单基因GSEA富集分析 | 20220404有出过类似的分享。今天，我们也结合相关的资源整理出一篇关于GSEA的教程及出图教程。每个方法的教程很多，我们大家结合自己的需求进行分析即可。以及，对于目前知识分享博主很多，只要你自己动手搜索，基本可以找到你的需求。更新！对于GSEA的教程原计划是在2月2日发表，但是由于有预约被占用了，因此这个教程也就是往后推迟。我在2月1日将我们的教程发在社群中，**也有同学提出疑问：**单基因也可以做GSEA分析，以及给出自己的疑问？
gsea结果分析图怎么看_利用GSEA对基因表达数据做富集分析
weixin_30905655的博客
12-23 2995
imageGene Set Enrichment Analysis (GSEA) is a computational method that determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biologic...
差异基因富集分析（R语言——GO&KEGG&GSEA）
m0_73145613的博客
07-03 2583
通路富集分析气泡图
GSEA软件基因富集分析
Eric_blog的博客
02-27 4751
快速目录链接GSEA分析简介分析步骤1、准备数据（1）表达数据文件（2）表型标签文件2、数据导入3、参数设置及运行 GSEA分析简介 基因集富集分析（GSEA）是一种计算方法，用于确定一组定义好的基因是否在两种生物状态（如表型）之间显示出统计上显著的一致性差异。 分析步骤 首先在官网下载软件：软件下载 1、准备数据 使用GSEA时，可以提供四个数据文件：表达数据集文件、表型标签文件、基因集文件和芯...
基因集富集分析(GSEA)
lazysnake666的博客
06-11 1万+
GSEA简介首先简单介绍一下GSEA，它是2005年在PNAS上发扬光大的方法，沿用至今，目的是看差异表达的基因在哪些基因集中富集。相比于Over-representation只关注显著差异表达的基因，GSEA分析纳入所有基因，将一些微弱但不显著的效应考虑在内。假设做了AHBA的10000个基因表达和大脑表征相关，其中300个基因经过检验是显著相关的，ORA富集分析这30...
r语言gsea生信分析代码
07-15
### 回答1： GSEA（基因集富集分析）是一种常用的生物信息学分析方法，用于研究基因集在基因表达谱中的富集情况。下面是使用R语言进行GSEA生信分析的代码示例： 1. 首先，需要安装和加载必要的R包，例如GSEA包和其他必要的依赖包。 ```R install.packages("GSEA") library(GSEA) ``` 2. 加载基因表达数据集，通常是一个包含基因表达矩阵的数据文件。假设文件名为"expression_data.txt"，其中包含基因表达矩阵和对应的样本信息。 ```R expression_matrix <- read.table("expression_data.txt", header = TRUE) ``` 3. 定义基因集，可以是预定义的基因集数据库（例如MSigDB）中的基因集，也可以是自定义的基因集。 ```R gene_sets <- c("GO_Biological_Process", "KEGG_Pathways", "Custom_Gene_Set") ``` 4. 进行GSEA分析，使用`gsea()`函数。其中，`gene_expr_matrix`参数为基因表达矩阵，`gene_sets`参数为基因集，`class_vector`参数为样本类别信息向量。 ```R gsea_results <- gsea(gene_expr_matrix = expression_matrix, gene_sets = gene_sets, class_vector = sample_classes) ``` 5. 分析结果包括富集分数（Enrichment Score）、正负富集基因集和富集图谱等。可以通过可视化方法进一步探索和解释这些结果。 ```R enrichment_score <- gsea_results$es positive_sets <- gsea_results$pos_sets negative_sets <- gsea_results$neg_sets gene_set_plot <- plot(gsea_results) ``` 以上是使用R语言进行GSEA生信分析的基本代码示例。根据具体的研究问题和分析目标，还可以进行更多的数据预处理和可视化分析。 ### 回答2： GSEA（Gene Set Enrichment Analysis）是一种生物信息学分析工具，可用于确定基因集在给定基因表达数据中的富集程度。下面是R语言中实现GSEA分析的示例代码。 首先，需要安装并加载GSEABase、clusterProfiler和enrichplot等相关的R包。 ```R install.packages("GSEABase") install.packages("clusterProfiler") install.packages("enrichplot") library(GSEABase) library(clusterProfiler) library(enrichplot) ``` 接下来，准备基因表达数据和基因集数据。假设基因表达数据保存在一个矩阵中，行表示基因，列表示样本；基因集数据保存在GMT格式文件中，每行包含一个基因集的名称、描述和基因列表。 ```R expression_data <- read.table("expression_data.txt", header = TRUE, row.names = 1) gmt_file <- system.file("extdata", "c2.cp.kegg.v7.4.symbols.gmt", package = "DOSE") gene_sets <- readGMT(gmt_file) ``` 然后，进行GSEA分析。可以选择使用差异表达基因列表作为输入，或者将基因表达数据与基因集数据一起传递。以下是基于基因表达数据进行GSEA分析的示例。 ```R gene_rank <- computeGeneRank(expression_data, method = "t.test") result <- enrichGSEA(gene_sets, gene_rank) ``` 最后，可以使用enrichplot包中的函数绘制GSEA结果的可视化，例如绘制富集图和基因集热图。 ```R dotplot(result, showCategory = 20) gene_heatmap(result, top = 10) ``` 通过这些代码，我们可以使用R语言实现GSEA生信分析，从而确定基因集在给定基因表达数据中的富集程度，并可视化展示分析结果。 ### 回答3： GSEA (基因集富集分析) 是一种用于分析生物学实验数据的生物信息学工具，它可以确定在给定条件下，特定基因集中的基因与实验结果相关性的显著性。下面是一个用R语言进行GSEA生信分析的代码示例： 1. 导入所需的R包。 ```R library(clusterProfiler) ``` 2. 导入基因表达数据。 ```R expression_data <- read.table("expression_data.txt", header = TRUE, sep = "\t") ``` 3. 根据实验分组信息创建一个分组向量。 ```R group <- c(rep("Group A", 3), rep("Group B", 3)) ``` 4. 根据基因的符号名称创建一个基因符号向量。 ```R gene_symbols <- c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5", "Gene6") ``` 5. 创建一个基因集对象。 ```R gene_set <- list( GroupA_genes = c("Gene1", "Gene2", "Gene3"), GroupB_genes = c("Gene4", "Gene5", "Gene6") ) ``` 6. 运行GSEA分析。 ```R gsea_result <- gseGO(expression_data, geneSet = gene_set, nPerm = 1000, minGSSize = 3, maxGSSize = 500, pvalueCutoff = 0.05) ``` 7. 查看GSEA结果。 ```R print(gsea_result) ``` 这段代码中，首先导入了clusterProfiler包，它包含了进行GSEA分析所需的函数。然后，基因表达数据被读入到一个名为expression_data的数据框中。接下来创建了一个分组向量，它指定了每个样品所属的实验组。然后，基因符号向量被创建，其中包含了基因的符号名称。根据实验组信息和基因符号，一个基因集对象被创建。最后，调用gseGO函数运行GSEA分析，其中包括参数，如基因集、置换次数、最小/最大基因集大小和显著性阈值。最后，打印GSEA分析的结果。
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值