首页 > 其他分享 >ChIP-seq 分析:基因集富集(11)

ChIP-seq 分析:基因集富集(11)

时间:2023-03-07 22:22:06浏览次数:51  
标签:11 seq msigdbr ChIP 基因 clusterProfiler 使用 GO 我们

动动发财的小手,点个赞吧!

1. 基因集检测

转录因子或表观遗传标记可能作用于按共同生物学特征(共享生物学功能、RNAseq 实验中的共同调控等)分组的特定基因组。

ChIPseq 分析中的一个常见步骤是测试常见基因集是否富含转录因子结合或表观遗传标记。

精心策划的基因集的来源包括 GO 联盟(基因的功能、生物过程和细胞定位)、REACTOME(生物途径)和 MsigDB(计算和实验衍生)。

2. ChIPseq 的基因集测试

可以对具有与其相关联的峰的基因集执行基因集富集测试。在这个例子中,我们将考虑峰值在基因 TSS 1000bp 以内的基因。

我们不会在测试中直接访问这些数据库库,但会使用广泛使用它们的其他 R/Bioconductor 库。

3. GO 和基因集测试

要在这里执行基因集测试,我们将使用 clusterProfiler 包。clusterProfiler 提供多种富集函数,允许将您的基因列表与已知(例如 GO、KEGG)或自定义基因集进行比较。

在这个例子中,我们使用我们发现与 Myc 峰重叠的所有 TSS 站点。落在 TSS 区域的峰将在我们带注释的 GRanges 对象的注释列中标记为“启动子”。

annotatedPeaksGR[1, ]

annotatedPeaksGR

我们可以通过对带注释的 GRanges 进行子集化并从 geneId 列中检索基因名称来提取 TSS 中具有峰的基因的唯一名称。

annotatedPeaksGR_TSS <- annotatedPeaksGR[annotatedPeaksGR$annotation == "Promoter",
    ]
genesWithPeakInTSS <- unique(annotatedPeaksGR_TSS$geneId)
genesWithPeakInTSS[1:2]

genesWithPeakInTSS

接下来,我们可以提取包含在 TxDb 对象中的所有基因,以用作我们用于通路富集的基因域。

allGeneGR <- genes(TxDb.Mmusculus.UCSC.mm10.knownGene)
allGeneGR[1:2, ]

allGeneGR

allGeneIDs <- allGeneGR$gene_id

一旦我们有了相同格式的基因列表和基因域,我们就可以在 enrichGO 函数中使用它们来执行基因本体分析

对于 ont 参数,我们可以在“BP”、“MF”和“CC”子本体之间进行选择,或者为所有三个选择“ALL”。

library(clusterProfiler)
library(org.Mm.eg.db)
GO_result <- enrichGO(gene = genesWithPeakInTSS, universe = allGeneIDs, OrgDb = org.Mm.eg.db,
    ont = "BP")

我们现在有一个 enrichResult 实例。从这个对象中,我们可以提取最丰富的基因本体类别的数据框。

GO_result_df <- data.frame(GO_result)
GO_result_df[1:5, ]

GO_result_df

可以使用 enrichplot 包从任何 enrichResult 对象生成网络图。

我们测量各种重要基因集之间的相似性并相应地对它们进行分组。 showCategory 参数指定要显示的顶级基因本体命中数。

library(enrichplot)
GO_result_plot <- pairwise_termsim(GO_result)
emapplot(GO_result_plot, showCategory = 20)

emapplot

除了基因本体之外,我们还可以使用 clusterProfiler enricher 函数针对我们作为 gmt 文件导入的自定义基因集测试我们的基因列表。类似于 enrichGO 函数,这将生成一个可用于可视化的 enrichResult 对象。

在这里,我们将使用 msigdbr 包从 MSigDB 获取基因集。我们还可以运行 msigdbr_collections 函数来查看将用于访问基因集的类别和子类别代码。

library(msigdbr)
msigdbr_collections()

msigdbr_collections

从上一张幻灯片的数据框中,我们可以识别我们想要的类别/子类别,并在 msigdbr 函数中使用它们。这里我们将使用“H”来访问 Hallmark 基因集,最后我们需要得到一个数据框,其中第一列包含基因集的名称,第二列包含基因 ID。

library(msigdbr)
msig_t2g <- msigdbr(species = "Mus musculus", category = "H", subcategory = NULL)
msig_t2g <- msig_t2g[, colnames(msig_t2g) %in% c("gs_name", "entrez_gene")]
msig_t2g[1:3, ]

msig_t2g

然后我们运行基因集富集,使用我们创建的术语到基因映射作为 enricher 函数中的 TERM2GENE 参数。

hallmark <- enricher(gene = genesWithPeakInTSS, universe = allGeneIDs, TERM2GENE = msig_t2g)
hallmark_df <- data.frame(hallmark)
hallmark_df[1:3, ]

hallmark_df

我们在RNAseq课程中了解到了goseq包,这是另一个类似于clusterProfiler的功能标注包,在这里,我们对 MSigDB Hallmark 基因集执行相同的富集测试。

对于 goseq,我们需要所有基因(宇宙)的命名向量,其中 1 或 0 代表基因是否在 TSS 中达到峰值。我们可以简单地使用 as.integer 函数将逻辑向量转换为 1 表示 TRUE 和 0 表示 FALSE。

allGenesForGOseq <- as.integer(allGeneIDs %in% genesWithPeakInTSS)
names(allGenesForGOseq) <- allGeneIDs
allGenesForGOseq[1:3]

allGenesForGOseq

首先,我们必须使用 nullp 函数构建一个 nullp data.frame 以便在 goseq 中使用,并提供我们的命名向量、要使用的基因组和使用的基因标识符。

nullp 函数试图纠正我们在基因集测试中可能看到的基因长度偏差。也就是说,较长的基因可能有更多机会在其中出现峰值。

library(goseq)
pwf = nullp(allGenesForGOseq, "mm10", "knownGene", plot.fit = FALSE)

pwf

我们可以使用与用于 clusterProfiler 的基因映射相同的术语(尽管它必须从 tibble 转换为 goseq 的数据框)来运行基因集富集测试。

Myc_hallMarks <- goseq(pwf, "mm10", "knownGene", gene2cat = data.frame(msig_t2g))
Myc_hallMarks[1:3, ]

Myc_hallMarks

本文由mdnice多平台发布

标签:11,seq,msigdbr,ChIP,基因,clusterProfiler,使用,GO,我们
From: https://www.cnblogs.com/swindler/p/17189954.html

相关文章

  • PAT 甲级 1011 World Cup Betting(20)
    Withthe2010FIFAWorldCuprunning,footballfanstheworldoverwerebecomingincreasinglyexcitedasthebestplayersfromthebestteamsdoingbattlesfor......
  • 「解题报告」CF1178B WOW Factor
    ¿题目链接这是一道非常富有启发性的题目,值得一做,闪耀着人类和机器人的智慧光辉的绝佳题目.首先注意到(vv)o(vv)的结构,可以考虑枚举中间的o,这样只需要算两边的选法......
  • 没等到Windows 11 Moment 2更新?可直接下载更新包升级
    没等到Windows11Moment2更新?可直接下载更新包升级2023-03-0102:50·cnBetaWindows11KB5022913正在作为运行22H2的PC上的可选更新推出。它是备受期待的"Moment2"版......
  • WIN11无法访问win7的共享打印机,(操作不能完成(错误0x00000709)),WIN10能正常访问该共享打
    1、问题描述:WIN11无法访问win7的共享打印机(操作不能完成(错误0x00000709)),win10可以访问。三台电脑都在同一个局域网内,分别为win7,win10,win11。WIN7系统为旗舰版,作为共......
  • 116、tail+grep命令——2023年3月7日10:01:06
    2023年2月20日14:50:371、tail基本命令tail命令.因为查看日志通常从后面最新的日志去看,tail命令就是从后往前找.比如下述命令会显示access.log的最后10行的内......
  • 删除win11右键一级菜单的AMD驱动栏
    1、按Win+R快捷键,输入"regedit”打开注册表。2、地址栏输入HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\Windows\CurrentVersion\ShellExtensions回车。3、点击选择Sh......
  • L11_买东西询问是否有某个东西
    概述在日语中,买东西是想问店里有没有某个东西,可以采用:物品名称はありますか的句式。おまもりはありますか有护身符吗?动画会话A:このTシャツ見て......
  • 每日小结(11)
    今天写了一个小测试,关于数组的子数组之和最大的。今天我学习了一种解决数组问题的算法,即统计数组中子数组之和最大的问题。该算法可以在O(n)的时间复杂度内解决这个问题,因......
  • 11_Redis
    Redis【简介&安装篇】-知乎(zhihu.com)【1】什么是redis,谈谈你对redis的理解redis就是一个数据库,不过与传统数据库不同的是,redis的数据是key-value存储在内存中的,......
  • 3.36每日总结11
    今天利用不到两个小时的时间进行了第一次作业界面的设计以及后台部分代码的设计,在这期间遇到了布局中按钮位置不能改变的问题,然后经过百度查找到了改变线性布局的位置......