GSEAvis富集分析可视化

最新推荐文章于 2024-06-26 09:53:25 发布

医学和生信笔记

最新推荐文章于 2024-06-26 09:53:25 发布

阅读量1k

点赞数

分类专栏：生信数据挖掘文章标签： r语言数据挖掘生信 TCGA

本文链接：https://blog.csdn.net/Ayue0616/article/details/132692045

版权

生信数据挖掘专栏收录该内容

54 篇文章 199 订阅

订阅专栏

本文介绍了如何使用R包GseaVis进行GSEA富集分析可视化，包括彩色RES图、generank图、多图展示和火山图等，并展示了如何准备数据、进行差异分析和对接GSEA桌面软件结果。

摘要由CSDN通过智能技术生成

今天继续介绍富集分析可视化哦~往期推文链接：

今天要介绍的老俊俊的R包：GseaVis，专门用于GSEA富集分析可视化，相比于enrichplot，增加了很多好用的功能，很多功能让我直呼泰裤辣！！

本期目录：

文章目录

准备数据

用gse87466这个GEO的数据做演示，下载整理的过程这次就不演示了。数据可在粉丝QQ群免费下载。

library(easyTCGA)
load(file = "G:/easyTCGA_test/gse87466.Rdata")

这是一个炎症性肠病的数据集，一共108个样本，21个normal，87个uc（ulcerative colitis）。

探针注释我已经提前做好了，但是有一些探针对应多个symbol，为了方便我这里直接删掉了：

exprSet[1:4,1:4]
##                                           GSM2332098 GSM2332099 GSM2332100
## IGK@ /// IGKC                               13.86197   13.76880   13.95740
##                                             13.95740   13.92619   13.79664
## IGL@                                        13.73797   13.61266   13.86197
## IGH@ /// IGHA1 /// IGHA2 /// LOC100126583   13.79664   13.16844   13.76880
##                                           GSM2332101
## IGK@ /// IGKC                               13.95740
##                                             13.86197
## IGL@                                        13.76880
## IGH@ /// IGHA1 /// IGHA2 /// LOC100126583   13.73797
group <- factor(group_list,levels = c("normal","UC"))
table(group)
## group
## normal     UC 
##     21     87

首先对这个数据做下差异分析，也是用easyTCGA包，1行代码即可，基因芯片数据也是支持的，并且它会自动检测需不需要进行log2转换，如果是count矩阵，会自动使用DESeq2、limma、edgeR进行差异分析，如果不是，会自动进行wilcoxon和limma的差异分析：

library(easyTCGA)

diff_res <- diff_analysis(exprset = exprSet
                          , group = group
                          , is_count = F
                          )
## log2 transform not needed
## => Running limma
## => Running wilcoxon test
## => Analysis done.

# limma的结果
diff_limma <- diff_res$deg_limma

# 多个gene symbol的直接删除，方便演示
diff_limma <- diff_limma[!grepl("/",diff_limma$genesymbol),]
head(diff_limma)
##               logFC   AveExpr         t      P.Value    adj.P.Val        B
## SLC6A14    5.024103  9.413107  21.56440 4.104849e-41 8.514279e-37 82.58182
## LOC389023 -3.550396  5.541681 -21.01057 4.054400e-40 4.204818e-36 80.36199
## SLC23A1   -2.473180  5.649224 -17.88487 3.378001e-34 2.335550e-30 67.08748
## DUOX2      4.911030  9.916299  17.37129 3.569259e-33 1.850839e-29 64.78265
## DPP10     -1.910958  3.991413 -16.98863 2.113068e-32 7.304876e-29 63.04259
## TIMP1      2.125930 11.402645  16.88534 3.425860e-32 1.015131e-28 62.56956
##           genesymbol
## SLC6A14      SLC6A14
## LOC389023  LOC389023
## SLC23A1      SLC23A1
## DUOX2          DUOX2
## DPP10          DPP10
## TIMP1          TIMP1

GSEA富集分析

富集分析首选clusterProfiler，没有之一！简单，好用！

clusterProfiler为我们提供了非常好用的ID转换函数，这里的ID转换和上面说的探针注释并不是一回事：

suppressMessages(library(clusterProfiler))

gene_entrezid <- bitr(geneID = diff_limma$genesymbol
                         , fromType = "SYMBOL" # 从symbol
                         , toType = "ENTREZID" # 转成ENTREZID
                         , OrgDb = "org.Hs.eg.db"
                         )
## 
## 'select()' returned 1:many mapping between keys and columns
head(gene_entrezid)
##    SYMBOL ENTREZID
## 1 SLC6A14    11254
## 3 SLC23A1     9963
## 4   DUOX2    50506
## 5   DPP10    57628
## 6   TIMP1     7076
## 7    LCN2     3934

富集分析最好用ENTREZID进行，关于多种不同的ID，在曾老师的书中都有详细介绍，强烈推荐初学者一定要看：生信初学者基础知识资源推荐

做GSEA分析对数据格式有要求，之前也说过，需要是一个有序的数值型向量，其名字是基因的ID

gene_entrezid <- merge(gene_entrezid,diff_limma,by.x = "SYMBOL", by.y = "genesymbol")
genelist <- gene_entrezid$logFC
names(genelist) <- gene_entrezid$ENTREZID
genelist <- sort(genelist,decreasing = T)

head(genelist)
##     4314    11254    50506     1673     1116     6279 
## 5.123666 5.024103 4.911030 4.608619 4.552790 4.256463

我们使用msigdbr包从msigdb数据库下载人类的C5注释集，大家常用的GO、KEGG的数据其实都是包括在msigdb数据库中的。

library(msigdbr)

m_t2g <- msigdbr(species = "Homo sapiens", category = "C5") %>% 
  dplyr::select(gs_name, entrez_gene)
head(m_t2g)
## # A tibble: 6 × 2
##   gs_name                                          entrez_gene
##   <chr>                                                  <int>
## 1 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS       60496
## 2 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS       10840
## 3 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS      160428
## 4 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS        4522
## 5 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS       25902
## 6 GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS      441024

然后是进行GSEA分析：

gsea_res <- GSEA(genelist, 
                 TERM2GENE = m_t2g,
                 minGSSize = 10,
                 maxGSSize = 500,
                 pvalueCutoff = 0.05,
                 pAdjustMethod = "BH",
                 seed = 456
                 )
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...

gsea_res_symbol <- setReadable(gsea_res,"org.Hs.eg.db",keyType = "ENTREZID")

富集分析得到的结果是一个对象，关于这个对象包括哪些东西，如何对它进行各种操作，我们在之前的推文都介绍过了，这里就不多说了~

GseaVis包的功能很强大。

library(GseaVis)

彩色RES图

Running Enrichment Score的线条变成了彩色，注意geneSetID只能接受字符串，不可以直接使用数字。

gseaNb(gsea_res, geneSetID = "GOBP_HUMORAL_IMMUNE_RESPONSE")

plot of chunk unnamed-chunk-9

通路名截断

还会很贴心的帮你把下划线去掉：

gseaNb(gsea_res, 
       geneSetID = "GOBP_ANTIMICROBIAL_HUMORAL_IMMUNE_RESPONSE_MEDIATED_BY_ANTIMICROBIAL_PEPTIDE",
       termWidth = 30
       )

plot of chunk unnamed-chunk-10

新的gene rank图

彩色的RES以及gene rank图形：

gseaNb(gsea_res,
       geneSetID = 'GOBP_HUMORAL_IMMUNE_RESPONSE',
       newGsea = T,
       addPoint = T,
       newHtCol = c("#336699", "white", "#993399")
       )

plot of chunk unnamed-chunk-11

添加NES信息

enrichplot需要自己DIY才能添加NES信息，这里提供了另一种思路：

gseaNb(gsea_res,
       geneSetID = 'GOBP_HUMORAL_IMMUNE_RESPONSE',
       addPval = T,
       pCol = "steelblue",
       pvalX = 0.65, # 位置
       pvalY = 0.7,
       pHjust = 0, # 对齐方式
       nesDigit = 4, # 小数点位数
       pDigit = 4
       )

在这里插入图片描述

多个gseaplot拼图

enrichplot的gseaplot2函数的结果是gglist，所以不能直接拼图，但是GseaVis可以：

terms <- gsea_res@result$ID[1:4]

gseaplot_list <- lapply(terms, function(x){
  gseaNb(object = gsea_res,
         geneSetID = x,
         termWidth = 30,
         addPval = T,
         pvalX = 0.75,
         pvalY = 0.6
         )
})

# 可以直接拼
cowplot::plot_grid(plotlist=gseaplot_list, ncol = 2)

plot of chunk unnamed-chunk-13

多条通路另类展示

不得不说颜值上升了好几个高度！

terms <- gsea_res@result$ID[1:3]

gseaNb(object = gsea_res,
       geneSetID = terms,
       subPlot = 2,
       termWidth = 35,
       legend.position = c(0.8,0.8),
       addPval = T,
       pvalX = 0.05,pvalY = 0.05)

plot of chunk unnamed-chunk-14

dotplot

可以分开展示上调和下调的通路。

dotplotGsea(data = gsea_res,
            topn= 10,
            str.width = 20 # 折叠通路名，我这里不能用...
            )

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

还可以画成棒棒糖图：

dotplotGsea(data = gsea_res,
            topn= 10,
            order.by = "NES",
            add.seg = T,
            line.col = 'orange',
            line.type = 1
            )

plot of chunk unnamed-chunk-16

火山图

火山图展示GSEA的结果，泰裤辣！再也不用自己提取结果了。

volcanoGsea(data = gsea_res,
            nudge.y = c(-0.8,0.8)
            )

plot of chunk unnamed-chunk-17

添加热图

提供表达矩阵即可，但是第一列需要是基因名字！

library(tibble)
exprset <- rownames_to_column(exprSet,var = "gene")
exprset[1:4,1:4]
##                                        gene GSM2332098 GSM2332099 GSM2332100
## 1                             IGK@ /// IGKC   13.86197   13.76880   13.95740
## 2                                             13.95740   13.92619   13.79664
## 3                                      IGL@   13.73797   13.61266   13.86197
## 4 IGH@ /// IGHA1 /// IGHA2 /// LOC100126583   13.79664   13.16844   13.76880

无需手动匹配基因，会自动帮你进行整理：

gseaNb(object = gsea_res_symbol,
       geneSetID = 'GOBP_HUMORAL_IMMUNE_RESPONSE',
       add.geneExpHt = T,
       exp = exprset)

plot of chunk unnamed-chunk-19

对接GSEA桌面软件的结果

gseaNb 可以直接传入输出文件路径 (filePath) 自动进行分析, 你也可以先使用 readGseaFile 函数读取,再传给 filePath 参数,这样画图就不用每次重新分析整合数据了。

我并不用那个GSEA的软件，所以这里就不展示了~大家可以参考俊俊的推文：GseaVis 一键对接 GSEA 软件结果并可视化

参考资料

大家可以去老俊俊的github查看更加详细的内容：https://github.com/junjunlab/GseaVis

给他点个小星星~

医学和生信笔记

关注

0
点赞
踩
6

收藏

觉得还不错? 一键收藏
0
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫

专栏目录