美文网首页rna_seq
[R]bioconductor之ChIPseeker学习

[R]bioconductor之ChIPseeker学习

作者: 小贝学生信 | 来源:发表于2020-10-28 10:19 被阅读0次

    ChIPseeker包南方医科大学Y叔大牛写的许多有名的生信R包之一,其最初设计用于chip-seq的macs peak calling结果分析以及可视化,后来逐渐也适用于相关的peak分析。
    参考链接:https://www.bioconductor.org/packages/release/bioc/vignettes/ChIPseeker/inst/doc/ChIPseeker.html
    以及Y叔自己的微信公众号教程:https://mp.weixin.qq.com/mp/appmsgalbum?__biz=MzI5NjUyNzkxMg==&action=getalbum&album_id=1300625300497268737&scene=173&from_msgid=2247488238&from_itemidx=1&count=3#wechat_redirect

    1、关于ChIP-seq

    详见之前的笔记

    • 之前在学习macs文章时,有了解过;这里再简单学习一下。
    • 如下图,DNA上的蛋白结合位点往往是基因表达调控的关键位置,ChIP技术就是针对性的挑选出这些位置。
    • DNA和蛋白质交联(cross-linking),超声(sonication)将染色体随机切割,利用抗原抗体的特异性识别(IP),把目标蛋白相结合的DNA片段沉淀下来,反交联释放DNA片段,最后是测序(sequencing)。
    • MACS软件通过一定原理算法,对测序比对结果识别出有意义的peak。ChIPseeker包就是衔接这一步之后开始做的。


      ChIP

    2、准备工作 preparation

    • 因为macs结果为bed输出格式,所以需要了解bed,即bedtools软件
    • 了解GRangesTxDb这两种常见的生信基础数据对象
    library(TxDb.Hsapiens.UCSC.hg19.knownGene)
    txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
    
    • 安装R包,找到示例数据
    BiocManager::install("ChIPseeker")
    library(ChIPseeker)
    files <- getSampleFiles()
    print(files)
    #bed转为Granges对象
    peak <- readPeakFile(files[[4]])
    peak
    
    • 如下图,即为ChIPseeker包分析所需的peak GRange对象。
      分割线左边三列分别为所在染色体信息,起止位点,正负链情况;
      右边两列分别为peak name与score(我认为可以理解与reads数正相关)


      peak

    3、ChIPseeker基础peak可视化

    3.1、概况covplot()

    观察所有peak在染色体的分布、表达情况

    #依据第五列score,表明峰的高低情况
    covplot(peak, weightCol="V5")
    
    covplot all peak&all chromesome
    #筛选指定染色体的指定区域的分布情况
    covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4.5e7, 5e7))
    
    covplot some area

    3.2、针对某一feature的分布情况

    • heatmap
      常见的分析是观察不同peak分布在TSS的promoter区域情况
    #自己定义promoter区域,上下游3000bp
    promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
    #不理解这个函数也没关系,是为下一步做热图提供matrix
    tagMatrix <- getTagMatrix(peak, windows=promoter)
    tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
    

    如下图结果,每一行代表一个promoter区域,红线的即为peak分布


    tagheatmap
    #一键绘图,效果同上
    peakHeatmap(peak, TxDb=txdb, upstream=3000, downstream=3000, color="red")
    
    • 峰图
      上面的热图是描绘了所有的promoter情况,可以绘制一个峰图描述所有分布的平均情况。
    plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
                xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")
    
    plotAvgProf
    #加一个置信区间
    plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)
    
    plotAvgProf with conf

    we developed getBioRegion function to support centering all peaks to the start region of Exon/Intron. Users can also create heatmap or average profile of ChIP peaks binding to these regions.

    4、ChIPseeker peak annotation

    4.1 what's peak annotation

    • 简单理解peak 注释就是peak落在染色体的哪一个位置上。常见的基因结构组成如下图所示。


      basic structure of gene
    • 此外ChIPseeker的peak注释时还提供另外一种注释方法,具体在注释结果时再具体了解(nearest gene annotation)。

    4.2 annotatePeak()

    (1)just do it
    • ChIPseeker包主要用annotatePeak()注释peak。需要提供两个文件:一是peak文件,可以是bed或者Granges;另一个是对应物种的TxDb对象(提供原始注释信息)
    • 此外promoter的区间可以自己定义,默认设置为TSS上下游3k区域
    peak
    #共计1331个peak
    txdb
    peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000), TxDb=txdb)
    peakAnno
    

    如下图,如果在R里直接观察结果,它会告诉我们ChIPseq的位点落在基因组上什么样的区域,分布情况如何。(即第一种注释方法genomic annotation)


    genomic annotation

    在注释时,有的peak可能同时落在两个或者更多的gene feature里(例如是一个基因的外显子而同时又是另一个基因的内含子),但只能注释其中一个。默认按照Promoter、5’ UTR、3’ UTR、Exon、Intron、Downstream、Intergenic顺序先后注释。

    • 一般会将上述的结果输出为GRanges格式、或者data.frame格式;便于查看,同时也能了解到annotatePeak第二种nearest gene annotation结果。
    class(peakAnno)
    peakAnno.df <- as.data.frame(peakAnno)
    peakAnno.gr <- as.GRanges(peakAnno)
    head(peakAnno.gr, 3)
    

    如下图,右上角为genomic annotation结果、下面为nearest gene annotation结果。

    • nearest gene annotation最近基因注释:是peak相对于转录起始位点的距离,不管这个peak是落在内含子或者别的什么位置上,即使它落在基因间区上,我都能够找到一个离它最近的基因(即使它可能非常远)。
    • 如果peak和TSS有overlap,genomic annotation就是promoter,距离就是0,而最近基因也是同一个,所以在这种情况下,两种注释都指向同一个基因。
    • 最近基因的注释信息虽然是以基因为单位给出,但我们针对的是转录起始位点来计算距离,针对于不同的转录本,一个基因可能有多个转录起始位点,所以注释是在转录本的水平上进行的,我们可以看到输出有一列是transcriptId.


      head(peakAnno.gr, 3)

    另外一种思路:注意上述nearest gene annotation默认找的是最近的TSS,即first anno与second anno对应的可以不是同一个基因。如果我想说只要和基因有overlap就是最近基因,那么这两种注释的基因应该是一致的,只需把overlap="TSS"(default)设置为overlap="all"

    5、ChIPseeker基于注释的peak可视化

    (1)genomic annotation可视化
    • 饼图或柱状图可视化组成比例
    plotAnnoPie(peakAnno)
    plotAnnoBar(peakAnno)
    
    pie chart
    • 考虑到注释到多种feature的可能
    vennpie(peakAnno)
    
    venn + pie
    upsetplot(peakAnno, vennpie=TRUE)
    

    如下图可以清楚地看到绝大多数的peak同时落到了多种feature里


    upsetplot
    (2)nearest gene annotation结果可视化
    • 可视化the distance from the peak (binding site) to the TSS of the nearest gene
    • plotDistToTSS can calculate the percentage of binding sites upstream and downstream from the TSS of the nearest genes, and visualize the distribution.
    plotDistToTSS(peakAnno,
                  title="Distribution of transcription factor-binding loci\nrelative to TSS")
    
    plotDistToTSS

    ChIPseeker包暂时先学习到这里,还有很多深入的功能,比如富集分析等,之后有机会再学习。
    感觉到国人写的R包,然后看中文的原版说明书还是比较轻松的~~

    相关文章

      网友评论

        本文标题:[R]bioconductor之ChIPseeker学习

        本文链接:https://www.haomeiwen.com/subject/yjefvktx.html