专栏首页生信修炼手册使用ChIPseeker进行peak注释

使用ChIPseeker进行peak注释

欢迎关注”生信修炼手册”!

ChIPseeker是使用的最广泛的peak注释软件之一,提供了以下多种功能

  1. peak在染色体和TSS位点附近分布情况可视化
  2. peak关联基因注释以及在基因组各种元件上的分布
  3. 获取GEO数据库中peak的bed文件
  4. 多个peak文件的比较和overlap分析

首先我们需要输入peak文件,支持两种格式,第一种是BED格式,最少只需要3列内容记录peak的染色体位置就可以了,示意如下

当然也可以有多余的列,只需要符合BED格式的标准即可;另外一种和MACS的peak calling输出结果类似,第一行为表头,示意如下

通过函数readPeaks读取peak文件,用法如下

peak <- readPeakFile("peak.bed")

函数根据文件名称的后缀来判断是否为bed格式,建议BED格式的输入文件后缀统一成.bed, 当然压缩文件也是支持的,比如.bed.gz;如果不是BED格式的输入,文件名称则不能使用BED格式对应的后缀。

下面来详细看下几个主要功能的代码和结果展示

1. peak 在染色体上的分布

用法如下

covplot(peak, chr = c("chr1", "chr2"))

输出结果示意如下

2. peak 在TSS位点附件的分布

用法如下

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
# 定义TSS上下游的距离
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

热图每一行代表一个基因,展示的是所有基因TSS两侧的分布,除了热图外,还可以对所有基因取均值,用折线图来展示TSS两侧分布情况,用法如下

plotAvgProf(
  tagMatrix,
  xlim=c(-3000, 3000),
  xlab="Genomic Region (5'->3')",
  ylab = "Read Count Frequency")

输出结果示意如下

3. pea关联基因注释

用法如下

peakAnno <- annotatePeak(
    peak,
    tssRegion = c(-3000, 3000),
    TxDb = txdb,
    annoDb = "org.Hs.eg.db")

write.table(
    as.data.frame(peakAnno),
    "peak.annotation.tsv",
    sep="\t",
    row.names = F,
    quote = F)

注释文件内容如下

给出了关联的基因以及对应的基因组区域的类别,根据这个结果,可以提取关联基因进行下游的功能富集分析,比如提取geneid这一列,用clusterProfiler进行GO/KEGG等功能富集分析。 注释的结果还提供了多种可视化方式,其中饼图最为常见,用法如下

plotAnnoPie(peakAnno)

输出结果示意如下

4. 下载GEO中的peak文件

以hg19为例,首先查询对应的GEO编号信息,用法如下

> hg19 <- getGEOInfo(genome="hg19", simplify=TRUE)
> head(hg19)
    series_id        gsm     organism
111  GSE16256  GSM521889 Homo sapiens
112  GSE16256  GSM521887 Homo sapiens
113  GSE16256  GSM521883 Homo sapiens
114  GSE16256 GSM1010966 Homo sapiens
115  GSE16256  GSM896166 Homo sapiens
116  GSE16256  GSM910577 Homo sapiens

由于列数太多,上述结果只展示了部分信息,对于每个bed文件,会列出对应的描述信息,方便筛选感兴趣的peak进行下载,可以根据GSM编号进行下载,用法如下

downloadGSMbedFiles("GSM521889", destDir="hg19")

也可以根据下载一个基因组对应的所有peak文件,用法如下

downloadGEObedFiles(genome="hg19", destDir="hg19")

5. peak间的overlap分析

peak的overlap分析不仅可以探究生物学重复样本间的一致性,还可以进一步识别多种蛋白或者转录因子在调控网络中的作用,如果两个蛋白的chip结果overlap显著,很可能这两个蛋白构成了复合体,或者两种蛋白具有相互作用,这对于探究其调控机制有相当大的帮助。用法如下

enrichPeakOverlap(
    queryPeak     = peak_setA,
    targetPeak    = c(peak_setB, peak_setC),
    TxDb          = txdb,
    pAdjustMethod = "BH",
    nShuffle      = 1000,
    chainFile     = NULL,
    verbose       = FALSE)

依次将query的peak与target中的每一个peak文件进行overlap分析,计算出一个p值代表两个peak之间overlap的程度,p值越小,overlap的程度越高。

ChIPseeker除了peak基因注释的基本功能外,整合了GEO的下载功能与peak的overlap分析,可以方便的将自己的chip_seq数据与GEO的公共数据集进行比较分析。

·end·

—如果喜欢,快分享给你的朋友们吧—

扫描关注微信号,更多精彩内容等着你!

本文分享自微信公众号 - 生信修炼手册(gh_0146e37a8a70),作者:lzyg

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-04-02

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • annoPeakR:一个peak注释的在线工具

    annoPeakR是一个peak注释工具,基于R语言中的shiny包开发出的web应用,网址如下

    生信修炼手册
  • peak差异分析的工具那么多,如何选择?

    对于ATAC_seq, chip_seq等抗体富集型文库而言,peak calling是分析的第一步。通过peak calling,可以得到抗体富集的区域,这些...

    生信修炼手册
  • 使用ChIPpeakAnno进行peak注释

    ChIPpeakAnno是一个bioconductor上的R包,针对peak calling之后的下游分析,提供了以下多种功能

    生信修炼手册
  • annoPeakR:一个peak注释的在线工具

    annoPeakR是一个peak注释工具,基于R语言中的shiny包开发出的web应用,网址如下

    生信修炼手册
  • peak注释信息揭秘

    在chip_seq数据分析中,peak calling是核心,得到peak区间之后,我们首先需要对peak进行注释。所谓的注释其实是一个比较宽泛的概念,其中包含...

    生信修炼手册
  • GTRD:最全面的人和小鼠转录因子chip_seq数据库

    GTRD从SRA, GEO, ENCODE等公共数据库中收集转录因子相关的chip_seq数据,采用标准流程进行peak calling分析,并基于已有的转录因...

    生信修炼手册
  • MySQL性能测试 : 新的InnoDB Double Write Buffer

    新的MySQL8.0.20版本重新设计了InnoDB Double Write(DBLWR),确实是一个大的历史烦人的事情。为什么在过去这么痛苦,让我们付出了这...

    田帅萌
  • 服务器配置优化

    query_cache_size:查询缓存 OLAP 类型数据库,需要重点加大此内存缓存. 但是一般不会超过 GB. 对于经常被修改的数据,缓存会立马失效。 我...

    Lemon黄
  • 树莓派Linux基础(四):修改文件权限与从属关系

    小雨编程
  • ReMap:人类Chip-seq数据大全

    ReMap收集来自GEO和Encode项目中人的chip_seq数据,对来自不同细胞系,不同类别转录因子的数据进行归类整理,网址如下

    生信修炼手册

扫码关注云+社区

领取腾讯云代金券