peak的可视化

peak在染色体上的分布

covplot可视化BED数据

ChIPseeker实现了covplot来可视化BED数据，covplot支持直接读文件出图：

library(ChIPseeker)
library(ggplot2)
files <- getSampleFiles()
covplot(files[[4]])

> files
[1]GSM1174480_ARmo_0M_peaks.bed.gz
[1]GSM1174481_ARmo_1nM_peaks.bed.gz
[1]GSM1174482_ARmo_100nM_peaks.bed.gz
[1]GSM1295076_CBX6_BF_ChipSeq_mergedReps_peaks.bed.gz
[1]GSM1295077_CBX7_BF_ChipSeq_mergedReps_peaks.bed.gz

GSM1295076_CBX6_BF_ChipSeq_mergedReps_peaks.bed.gz

支持GRanges对象，同时可以多个文件或者GRangesList

peak=GenomicRanges::GRangesList(CBX6=readPeakFile(files[[4]]),CBX7=readPeakFile(files[[5]]))
png(file="a.png")
covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)
dev.off()

支持可视化某个窗口

png(file="a.png")
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7)) + facet_grid(chr ~ .id)
dev.off()

拿到数据后，我们首先会可视化看一下数据，接下来就会想知道这些peak都和什么样的基因有关。

原文: https://mp.weixin.qq.com/s/3CMj0xejiV-FSMC-Vxd_-w

固定的窗口里peak

我们可以限定在一个固定的窗口里，然后把这些peak全部align起来，看看在某个窗口上的结合谱图。比如说启动子区域，使用转录起始位点，然后指定上下游，使用getPromoters函数，就可以为我们准备好窗口。再使用getTagMatrix函数，把peak比对到这个窗口，并生成矩阵，供我们分析、可视化。

library(ChIPseeker)
files <- getSampleFiles()
peak <- readPeakFile(files[[4]])
promoter <- getPromoters(TxDb=txdb, 
    upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, 
    windows=promoter)

这时候，我们可以用tagHeatmap来画热图，颜色随意指定：

png(file="a.png",width=300)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), 
    color="red")
dev.off()

另外你可以直接从文件出发，ChIPseeker提供了peakHeatmap函数，你给文件名，就可以直接出图，两个函数都支持多种数据同时展示。

require(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
files <- getSampleFiles()
png(file="a.png")
peakHeatmap(files, TxDb=txdb, 
            upstream=3000, downstream=3000, 
            color=rainbow(length(files)))
dev.off()

看看结合的强度

png(file="a.png")
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
            xlab="Genomic Region (5'->3')", 
            ylab = "Read Count Frequency")
dev.off()

上面这个图可以非常清晰地看到，这个结合位点是集中在转录起始位点的。另外ChIPseeker也提供了plotAvgProf2函数，支持以文件为输入一步做图。下面这个代码，会产生上面一样的图。

plotAvgProf2(files[[4]], TxDb=txdb, 
             upstream=3000, downstream=3000,
             xlab="Genomic Region (5'->3')", 
             ylab = "Read Count Frequency")

另外这个函数还支持使用bootstrap估计置信区间：

png(file="a.png")
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), 
            conf = 0.95, resample = 1000)
dev.off()

同样这个函数也支持多个数据的比较:

png(file="a.png")
tagMatrixList <- lapply(files, getTagMatrix, 
                        windows=promoter)
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))
dev.off()

plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), 
    conf=0.95,resample=500, facet="row")

这里我展示的都是启动子窗口，然而ChIPseeker支持的并不止是启动子，不单单是看转录调控，比如说你是研究可变剪切的，那么蛋白结合到外显子/内含子的起始位置，显然就更感兴趣了。所以除了getPromoters函数之外，ChIPseeker还提供了getBioRegion函数，你可以指定’exon’或’intron’，它会像getPromoters一样为你准备好数据，即使这不是你感兴趣的，你也依然可以用它来准备窗口，然后可视化看一下，或许你能意外地发现你研究的蛋白竟然在外显子/内含子起始位置上有很强的结合谱，做为前期data exploration