数据分析:基于DESeq2的转录组功能富集分析

介绍

DESeq2常用于识别差异基因,它主要使用了标准化因子标准化数据,再根据广义线性模型判别组间差异(组间残差是否显著判断)。在获取差异基因结果后,我们可以进行下一步的富集分析,常用方法有基于在线网站DAVID以及脚本处理的两类,本文介绍基于fgsea的方法计算富集分析得分。

DESeq2差异分析

了解DESeq2如何标准化数据和识别差异基因。下面给出简要代码

library(DESeq2)
library(airway)
data("airway")
ddsSE <- DESeqDataSet(airway, design = ~ cell + dex)
ddsSE <- DESeq(ddsSE)
res <- results(ddsSE, tidy = TRUE) %>% na.omit() %>% as_tibble()head(res)
# A tibble: 6 x 7row             baseMean log2FoldChange  lfcSE   stat     pvalue      padj<chr>              <dbl>          <dbl>  <dbl>  <dbl>      <dbl>     <dbl>
1 ENSG00000000003    709.          0.381  0.101   3.79  0.000152   0.00128  
2 ENSG00000000419    520.         -0.207  0.112  -1.84  0.0653     0.197    
3 ENSG00000000457    237.         -0.0379 0.143  -0.264 0.792      0.911    
4 ENSG00000000460     57.9         0.0882 0.287   0.307 0.759      0.895    
5 ENSG00000000971   5817.         -0.426  0.0883 -4.83  0.00000138 0.0000182
6 ENSG00000001036   1282.          0.241  0.0887  2.72  0.00658    0.0328 

转换geneID

我们使用的MSigDB数据库的pathway 基因ID只有entrez和HGNC symbol两类,如果是ensemble id,需要转换。

library(org.Hs.eg.db)
library(tidyverse)
ens2symbol <- AnnotationDbi::select(org.Hs.eg.db,key=res$row, columns="SYMBOL",keytype="ENSEMBL")
ens2symbol <- as_tibble(ens2symbol)
head(ens2symbol)
# A tibble: 6 x 2ENSEMBL         SYMBOL  <chr>           <chr>   
1 ENSG00000000003 TSPAN6  
2 ENSG00000000419 DPM1    
3 ENSG00000000457 SCYL3   
4 ENSG00000000460 C1orf112
5 ENSG00000000971 CFH     
6 ENSG00000001036 FUCA2 
  • 合并数据;过滤NA值;去重;重复基因求stat(stat数据作为排序指标用于后续富集分析)
res2 <- inner_join(res, ens2symbol, by=c("row"="ENSEMBL")) %>% dplyr::select(SYMBOL, stat) %>% na.omit() %>% distinct() %>% group_by(SYMBOL) %>% summarize(stat=mean(stat))
head(res2 )
# A tibble: 6 x 2SYMBOL       stat<chr>       <dbl>
1 A1BG      0.680  
2 A1BG-AS1 -1.79   
3 A2M      -1.26   
4 A2M-AS1   0.875  
5 A4GALT   -4.14   
6 A4GNT     0.00777

构建fgsea输入数据

  • 基因排序值转换
library(fgsea)ranks <- deframe(res2)
head(ranks, 20)
        A1BG     A1BG-AS1          A2M      A2M-AS1       A4GALT        A4GNT         AAAS         AACS 0.679946437 -1.793291412 -1.259539478  0.875346116 -4.144839902  0.007772497  0.163986128  1.416071728 AADACL4        AADAT        AAGAB         AAK1        AAMDC         AAMP         AAR2        AARS1 
-1.876311694  3.079128034  1.554279946  1.141522348 -2.147527241 -3.170612332 -2.364380163  4.495474603 AARS2       AARSD1        AASDH     AASDHPPT 5.057470292  0.654208006  0.665531695 -0.353496148 
  • pathways的基因集合,上MSigDB下载基因集。演示使用KEGG基因集
pathways.hallmark <- gmtPathways("../../Result/GeneID/msigdb.v7.1.symbols_KEGG.gmt")
pathways.hallmark %>% head() %>% lapply(head)
$KEGG_GLYCOLYSIS_GLUCONEOGENESIS
[1] "ACSS2" "GCK"   "PGK2"  "PGK1"  "PDHB"  "PDHA1"$KEGG_CITRATE_CYCLE_TCA_CYCLE
[1] "IDH3B" "DLST"  "PCK2"  "CS"    "PDHB"  "PCK1" $KEGG_PENTOSE_PHOSPHATE_PATHWAY
[1] "RPE"   "RPIA"  "PGM2"  "PGLS"  "PRPS2" "FBP2" $KEGG_PENTOSE_AND_GLUCURONATE_INTERCONVERSIONS
[1] "UGT1A10" "UGT1A8"  "RPE"     "UGT1A7"  "UGT1A6"  "UGT2B28"$KEGG_FRUCTOSE_AND_MANNOSE_METABOLISM
[1] "MPI"  "PMM2" "PMM1" "FBP2" "PFKM" "GMDS"$KEGG_GALACTOSE_METABOLISM
[1] "GCK"     "GALK1"   "GLB1"    "GALE"    "B4GALT1" "PGM2"
  • 运行
fgseaRes <- fgsea(pathways=pathways.hallmark, stats=ranks, nperm=1000)
head(fgseaRes[order(pval), ])
  • 从查看KEGG_REGULATION_OF_ACTIN_CYTOSKELETON富集分数分布
plotEnrichment(pathways.hallmark[["KEGG_REGULATION_OF_ACTIN_CYTOSKELETON"]],ranks) + labs(title="KEGG_REGULATION_OF_ACTIN_CYTOSKELETON")

  • 查看上下调通路结果
topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n=10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
plotGseaTable(pathways.hallmark[topPathways], ranks, fgseaRes, gseaParam=0.5)

  • 其他展示方式
fgseaResTidy <- fgseaRes %>%as_tibble() %>%arrange(desc(NES))# Show in a nice table:
fgseaResTidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% DT::datatable()ggplot(fgseaResTidy, aes(reorder(pathway, NES), NES)) +geom_col(aes(fill = padj<0.0001)) +coord_flip() +labs(x="Pathway", y="Normalized Enrichment Score",title="Hallmark pathways NES from GSEA") + theme_minimal()

查看通路的基因

res_temp <- inner_join(res, ens2symbol, by=c("row"="ENSEMBL"))
pathways.hallmark %>% enframe("pathway", "SYMBOL") %>% unnest(cols = c(SYMBOL)) %>% inner_join(res_temp , by="SYMBOL") %>%head()
# A tibble: 6 x 9pathway                         SYMBOL row             baseMean log2FoldChange lfcSE   stat pvalue   padj<chr>                           <chr>  <chr>              <dbl>          <dbl> <dbl>  <dbl>  <dbl>  <dbl>
1 KEGG_GLYCOLYSIS_GLUCONEOGENESIS ACSS2  ENSG00000131069    669.         -0.269  0.114 -2.35  0.0188 0.0756
2 KEGG_GLYCOLYSIS_GLUCONEOGENESIS GCK    ENSG00000106633     28.8         0.305  0.374  0.815 0.415  0.662 
3 KEGG_GLYCOLYSIS_GLUCONEOGENESIS PGK1   ENSG00000102144   7879.         -0.300  0.353 -0.850 0.395  0.642 
4 KEGG_GLYCOLYSIS_GLUCONEOGENESIS PDHB   ENSG00000168291    648.         -0.257  0.102 -2.52  0.0117 0.0521
5 KEGG_GLYCOLYSIS_GLUCONEOGENESIS PDHA1  ENSG00000131828    651.         -0.0744 0.104 -0.715 0.475  0.710 
6 KEGG_GLYCOLYSIS_GLUCONEOGENESIS PGM2   ENSG00000169299    302.         -0.315  0.136 -2.33  0.0201 0.0797

其他用法

  • miR targets
fgsea(pathways=gmtPathways("msigdb/c3.mir.v6.2.symbols.gmt"), ranks, nperm=1000) %>% as_tibble() %>% arrange(padj)
  • GO annotations
fgsea(pathways=gmtPathways("msigdb/c5.all.v6.2.symbols.gmt"), ranks, nperm=1000) %>% as_tibble() %>% arrange(padj)
  • 非人物种
library(biomaRt)
mart <- useDataset("mmusculus_gene_ensembl", mart=useMart("ensembl"))
bm <- getBM(attributes=c("ensembl_gene_id", "hsapiens_homolog_associated_gene_name"), mart=mart) %>%distinct() %>%as_tibble() %>%na_if("") %>% na.omit()
bm

参考

  1. Fast Gene Set Enrichment Analysis

  2. DESeq results to pathways in 60 Seconds with the fgsea package

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

如若内容造成侵权/违法违规/事实不符,请联系一条长河网进行投诉反馈,一经查实,立即删除!

相关文章

## CSDN创作活动:缓解工作压力:程序员的健康之道

缓解工作压力&#xff1a;程序员的健康之道 在当今快节奏的社会中&#xff0c;程序员作为一个高度专业化和技术密集的群体&#xff0c;往往需要面对持续的工作压力和创新挑战。在如此高强度的工作环境下&#xff0c;如何有效缓解工作压力&#xff0c;保持工作效率和个人健康成…

Java线程池的七大参数说明

线程池中的七大参数如下&#xff1a; &#xff08;1&#xff09;corePoolSize&#xff1a;线程池中的常驻核心线程数。 &#xff08;2&#xff09;maximumPoolSize&#xff1a;线程池能够容纳同时执行的最大线程数&#xff0c;此值大于等于1。 &#xff08;3&#xff09;keepAl…

Recruit App

招聘类APP小程序

【Vulhub靶场】Nginx 漏洞复现

Nginx 漏洞复现 一、Nginx 文件名逻辑漏洞&#xff08;CVE-2013-4547&#xff09;1、影响版本2、漏洞原理3、漏洞复现 二、Nginx 解析漏洞1、版本信息&#xff1a;2、漏洞详情3、漏洞复现 一、Nginx 文件名逻辑漏洞&#xff08;CVE-2013-4547&#xff09; 1、影响版本 Nginx …

HSDB使用教程

HSDB&#xff1a;Hostspot Debugger&#xff0c;JVM内置的工具&#xff0c;用于深入分析JVM运行时的内部状态 启动HSDB java -cp D:/tools/jdk-1.8/lib/sa-jdi.jar sun.jvm.hotspot.HSDB 获取进程id jps 连接到指定进程 查找类 通过查询查找对象 输入查询语句 select d from …

SQL 基础 | AS 的用法介绍

SQL&#xff08;Structured Query Language&#xff09;是一种用于管理和操作数据库的标准编程语言。 在SQL中&#xff0c;AS关键字有几种不同的用法&#xff0c;主要用于重命名表、列或者查询结果。 以下是AS的一些常见用法&#xff1a; 重命名列&#xff1a;在SELECT语句中&a…

Python中的Higher库详解

Python的Higher库是一个用于对深度学习模型进行微分和元学习的强大工具。这个库允许用户在不修改底层代码的情况下对任何可微模型进行高阶梯度计算。在本文中&#xff0c;我们将详细介绍Higher库的基本用法&#xff0c;展示一些实际的代码示例&#xff0c;并探讨其在模型训练中…

【RBAC鉴权】node-casbin基础教程

一、RBAC概述 RBAC鉴权&#xff0c;完整的英文描述是&#xff1a;Role-Based Access Control&#xff0c;中文意思是&#xff1a;基于角色&#xff08;Role&#xff09;的访问控制。这是一种广泛应用于计算机系统和网络安全领域的访问控制模型。 简单来说&#xff0c;就是通过…

为人处事电影解说,全新升级瀚海跑道一分钟一条视频,全平台可推广,轻轻松松日入1000

自古以来&#xff0c;我国流行的一种现象是&#xff0c;大多数人都会与领导或上司打交道。由于某些话题不宜公开讨论&#xff0c;因此出现了许多含蓄的表达方式。随着年龄的增长&#xff0c;人们的态度也发生了变化&#xff0c;从最初的轻视到现在的重视。 下 载 地 址&#…

38-4 Web应用防火墙 - WAF的使用及规则

准备:38-3 Web应用防火墙 - 安装配置WAF-CSDN博客 WAF的使用 启动 Nginx /usr/local/nginx/sbin/nginx 为了测试未启动 ModSecurity 时的访问效果,我们可以模拟攻击。要查看当前虚拟机的 IP 地址,可以使用命令 ifconfig 浏览器中访问ip,如果要在真实机中访问就需要关闭…

【C++】vector类的增删改查模拟实现(图例超详细解析!!!)

目录 一、前言 二、源码引入 三、vector的模拟实现 ✨实现框架 ✨前情提要 ✨Member functions —— 成员函数 ⚡构造函数 ⭐无参构造 ⭐迭代器区间构造 ⭐n个值构造 ⚡拷贝构造 ⚡运算符赋值重载 ⚡析构函数 ✨Element access —— 元素访问 ⚡operator[ ] …

C语言 | Leetcode C语言题解之第63题不同路径II

题目&#xff1a; 题解&#xff1a; int uniquePathsWithObstacles(int** obstacleGrid, int obstacleGridSize,int* obstacleGridColSize) {int n obstacleGridSize, m obstacleGridColSize[0];int f[m];memset(f, 0, sizeof(f));f[0] (obstacleGrid[0][0] 0);for (int i…

融创共赢,算网领航 | 移动云智能芯片开放实验室正式发布

4月29日上午&#xff0c;中国移动算力网络大会“融创共赢&#xff0c;算网领航-算网原生技术”分论坛在苏州金鸡湖国际会议中心顺利召开&#xff0c;中国移动云能力中心副总经理吴世俊出席论坛并发表致辞。大会举行了智能芯片开放实验室发布仪式&#xff0c;同时发布移动云最新…

服务器IP选择

可以去https://ip.ping0.cc/查看IP的具体情况 1.IP位置--如果是国内用&#xff0c;国外服务器的话建议选择日本&#xff0c;香港这些比较好&#xff0c;因为它们离这里近&#xff0c;一般延时低&#xff08;在没有绕一圈的情况下&#xff09;。 不过GPT的话屏蔽了香港IP 2. 企…

2024年税务师报名照片处理工具使用详细教程✅

第1️⃣步&#xff1a;登陆税务师考试报名入口 第2️⃣步&#xff1a;登陆税务师考试报名入口 第3️⃣步&#xff1a;下载税务师报名照片处理工具 第4️⃣步&#xff1a;双击打开照片检测工具 第5️⃣步&#xff1a;点击“打开照片文件”上传合规照片 第6️⃣步&#xff1a;在对…

JAVA面试题---WEB部分

网络通讯 TCP与UDP TCP(Transmission Control Protocol 传输控制协议)是一种面向连接(连接导向)的、 可靠的、 基于 IP 的传输层协议。 UDP 是 User Datagram Protocol 的简称&#xff0c;中文名是用户数据报协议&#xff0c;是 OSI 参考模 型中的传输层协议&#xff0c;它是…

用WebUI生成一个毛茸茸的图标教程

&#x1f680;内容概要 好久没有更新webui的教程了&#xff0c;五一期间刷新liblib&#xff0c;发现更新了几个有趣的lora&#xff0c;正好这段时间看到有朋友发了毛绒图标&#xff0c;所以这里做了一个简单的webui 教程&#xff0c;教你如何一步生成毛绒图标&#xff0c;就像…

Unity SteamVR入门

概述 VR项目现在在当前已经是非常热门的技术&#xff0c;可以给玩家身临其境的感觉&#xff0c;接下来让我们学习这部分的内容吧&#xff01; SteamVR Input SteamVR绑定流程&#xff0c;在Windows窗口的点击SteamVR-input&#xff0c;图1&#xff0c;在这里可以选择你需要绑定…

『跨端框架』Flutter环境搭建

『跨端框架』Flutter环境搭建 资源网站简介跨平台高性能发展历程跨平台框架的比较成功案例 环境搭建&#xff08;windows&#xff09;基础环境搭建Windows下的安卓环境搭建Mac下的安卓环境配置资源镜像JDKAndroid StudioFlutter SDK问题一问题二问题三修改项目中的Flutter版本 …

HANA小知识点记录:SQL JOIN 条件中的条件判断(CASE WHEN )

今天写sql的时候要用到这个&#xff0c;查到其他数据库是这么写的&#xff1a; 在hana里试了下上面那样写不行&#xff0c;试了一下可以这么写&#xff0c;满足需求 LEFT JOIN "SAPHANADB"."/BI0/PCUSTOMER" AS F-- 通过附加客户关联客户主数据 ON CASE W…