[R]bioconductor之ChIPseeker学习
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
-
plotDistToTSScan 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包，然后看中文的原版说明书还是比较轻松的~~
