数据分析:基于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.mzph.cn/news/831946.shtml

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

相关文章

银行智能化数据安全分类分级实践分享

文章目录 前言一、数据安全智能分类分级平台建设背景二、数据安全分类分级建设思路和实践1、做标签– 数据安全标签体系2、打标签– 鹰眼智能打标平台 3.03、用标签– 全行统一“数据安全打标签结果”服务提供前言 随着国家对数据安全的高度重视,以及相关法律法规的出台,数据…

python数据分析中数据可视化简单入门

1.折线图表 首先引入相关包pyecharts&#xff0c;如果没下载可以先下载 pip install pyecharts from pyecharts.charts import Lineline Line() # 添加x轴 line.add_xaxis([呱了个呱,羊村,牟多,蜂地,喵帕斯]) # 添加y轴 line.add_yaxis("GDP",[50,30,40,34,63,22])…

01 JVM --

JVM (Java Virtual Machine) 是一个虚拟机HotSpot 是 JVM 概念的一个实现。HotSpot 虚拟机通过即时编译 (JIT) 技术将 Java 字节码转换为本地机器码&#xff0c;以提高程序的执行效率。OpenJDK 是一个项目名&#xff0c;它在 HotSpot 的基础上开发了 HotSpot 的开源实现方法区是…

Redhat Linux忘记密码解决方案

1、重启系统 2、将光标移动到要启动的内核 3、按e编辑当前条目 4、将光标移动以Linux开头的行&#xff0c;此为内核命令行 5、在末尾添加人的rd.break 6、按ctrlx继续启动&#xff0c;如果发现输入的rd.break不能进入到伪系统&#xff0c;那么改为 rd.break consoletty0 (按ctr…

设计模式-03 设计模式-工厂模式factory-内部工厂

设计模式-03 设计模式-工厂模式factory-内部工厂 目录 设计模式-03 设计模式-工厂模式factory-内部工厂 1.定义 2.内涵 3.案例对比 4.特点 4.总结 1.定义 内部工厂模式是一种创建类对象的方式&#xff0c;其中工厂方法被封装在类内部&#xff0c;客户端只能通过类的公共…

epoll的LT和ET模式介绍

目录 1.epoll的LT和ET模式介绍 2.epoll的ET模式如何处理 2.1 epoll的ET模式编程读取数据的处理方式 2.2 将描述符设置为非阻塞模式的方法 3.ET模式的总结 4.epoll的LT模式和ET模式总结 5.IO复用总结 1.epoll的LT和ET模式介绍 epoll对文件描述符有两种操作模式: LT(Leve…

Linux(openEuler、CentOS8)常用的IP修改方式(文本配置工具nmtui+配置文件+nmcli命令)

----本实验环境为openEuler系统<以server方式安装>&#xff08;CentOS类似&#xff0c;可参考本文&#xff09;---- 一、知识点 &#xff08;一&#xff09;文本配置工具nmtui(openEuler已预装) nmtui&#xff08;NetworkManager Text User Interface&#xff09;是一…

关于apache+php用户验证

一.直接在apache配置配置用户信息 1.apache配置可以参考外部文档 https://developer.aliyun.com/article/507049 2.上面配置好在php获取用户信息&#xff08;登录apache会拦截&#xff09; $userName $_SERVER[PHP_AUTH_USER];$password $_SERVER[PHP_AUTH_PW]; 二.上面直…

在线OJ——链表经典例题详解

引言&#xff1a;本篇博客详细讲解了关于链表的三个经典例题&#xff0c;分别是&#xff1a;环形链表&#xff08;简单&#xff09;&#xff0c;环形链表Ⅱ&#xff08;中等&#xff09;&#xff0c;随机链表的复制&#xff08;中等&#xff09;。当你能毫无压力地听懂和成功地…

再学Java基础——final关键字

在Java中&#xff0c;final关键字是一个修饰符&#xff0c;它可以应用于类、方法、变量和常量。以下是final关键字在Java中的不同用途及其作用&#xff1a; final类&#xff1a; 当一个类被声明为final时&#xff0c;它不能被继承。这意味着没有其他类可以扩展这个类。作用&…

【莫比乌斯变换-02】关于莫比乌斯变换属性梳理

文章目录 一、说明二、多视角看莫比乌斯变换2.1 从几何角度2.2 复分析中的莫比乌斯变换2.3 莫比乌斯变换运算组合2.4 莫比乌斯变换的不动点2.5 三个点决定一个莫比乌斯变换2.6 交叉比2.7 莫比乌斯变换的逆变换 三 莫比乌斯变换性质证明3.1 证明1&#xff1a;莫比乌斯变换将圆变…

【莫比乌斯变换-03】python实现圆对圆的变换

文章目录 一、说明二、python实现复平面的莫比乌斯变换三、线的变换四、画笑脸 一、说明 我们在前面的文章中&#xff0c;叙述了莫比乌斯变换的复数分析&#xff0c;以及种种几何属性&#xff0c;本篇中叙述如何程序地实现&#xff1a;复平面上的圆在莫比乌斯变换下的图像是另…

基于高德 API 的自动获取气候数据的 Python 脚本

文章目录 高德申请 Key脚本介绍运行结果示例 源代码&#xff1a; https://github.com/ma0513207162/PyPrecip。pyprecip\reading\read_api.py 路径下。 项目介绍&#xff1a;PyPrecip 是一个专注于气候数据处理的 Python 库&#xff0c;旨在为用户提供方便、高效的气候数据处理…

linux基本操作

vim的基本操作 正常模式&#xff1a;启动vim后默认处于正常模式。不论位于什么模式&#xff0c;按下Esc建都会进入正常模式。 插入模式&#xff1a;在正常模式中按下i&#xff0c;l&#xff0c;a&#xff0c;A等键&#xff0c;会进入插入模式。现在只用记住按i键会进行插入模…

文本内容如何扫码展示?文本转二维码常用技巧

随着网络的不断发展&#xff0c;现在很多人会通过二维码的方式来作为内容的载体&#xff0c;用来给其他人提供需要的信息或者其他内容。其中&#xff0c;将文本制作二维码就是很常用的一种方式&#xff0c;通过扫码展示文字内容&#xff0c;比如信件、物品信息、联系方式等类型…

python学习笔记----函数进阶(七)

一、函数多个返回值 def test_return():x 1y "hello"z Truereturn x, y, za, b, c test_return()print(a) print(b) print(c)二、函数多种传参方式 2.1 函数参数种类 根据使用方式上的不同&#xff0c;函数有4种常见参数的使用方式&#xff1a; 位置参数关键…

EPAI手绘建模APP资源管理和模型编辑器2

g) 矩形 图 26模型编辑器-矩形 i. 修改矩形的中心位置。 ii. 修改矩形的长度和宽度。 h) 正多边形 图 27模型编辑器-内接正多边形 图 28模型编辑器-外切正多边形 i. 修改正多边形的中心位置。 ii. 修改正多边形中心距离端点的长度。 iii. 修改正多边形的阶数。阶数为3&…

Linux Multipath极简安装配置

1.安装软件 #yum install -y device-mapper* 2.开启并设置开机启动 #systemctl start multipathd && systemctl enable multipathd 3.连接存储 在存储上对主机进行卷映射 4.修改配置文件 注&#xff1a;此配置文件可能存在&#xff0c;也可能不存在&#xff0c…

禅道项目管理系统身份认证绕过漏洞复现(QVD-2024-15263)

0x01 产品简介 禅道项目管理软件是国产的开源项目管理软件,专注研发项目管理,内置需求管理、任务管理、bug管理、缺陷管理、用例管理、计划发布等功能,完整覆盖了研发项目管理的核心流程。 0x02 漏洞概述 2024年4月,互联网披露禅道系统存在身份认证绕过的漏洞情报。经分…

【Unity 协程】

Unity中的协程&#xff08;Coroutine&#xff09;是一种编程结构&#xff0c;它允许你以一种看似同步的方式编写可能需要异步执行的代码。协程特别适用于需要在一定时间后执行操作&#xff0c;或者在循环执行某段代码直到某个条件满足时的场景。 协程使用IEnumerator委托来实现…