关于程序员:ChIPseq-分析基因集富集11

38次阅读

共计 2773 个字符,预计需要花费 7 分钟才能阅读完成。

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

1. 基因集检测

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

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

精心策划的基因集的起源包含 GO 联盟(基因的性能、生物过程和细胞定位)、REACTOME(生物路径)和 MsigDB(计算和试验衍生)。

2. ChIPseq 的基因集测试

能够对具备与其相关联的峰的基因集执行基因集富集测试。在这个例子中,咱们将思考峰值在基因 TSS 1000bp 以内的基因。

咱们不会在测试中间接拜访这些数据库库,但会应用宽泛应用它们的其余 R/Bioconductor 库。

3. GO 和基因集测试

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

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

annotatedPeaksGR[1,]

咱们能够通过对带正文的 GRanges 进行子集化并从 geneId 列中检索基因名称来提取 TSS 中具备峰的基因的惟一名称。

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

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

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

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,]

能够应用 enrichplot 包从任何 enrichResult 对象生成网络图。

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

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

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

在这里,咱们将应用 msigdbr 包从 MSigDB 获取基因集。咱们还能够运行 msigdbr_collections 函数来查看将用于拜访基因集的类别和子类别代码。

library(msigdbr)
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,]

而后咱们运行基因集富集,应用咱们创立的术语到基因映射作为 enricher 函数中的 TERM2GENE 参数。

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

咱们在 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]

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

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

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

咱们能够应用与用于 clusterProfiler 的基因映射雷同的术语(只管它必须从 tibble 转换为 goseq 的数据框)来运行基因集富集测试。

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

本文由 mdnice 多平台公布

正文完
 0