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
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包,然后看中文的原版说明书还是比较轻松的~~
网友评论