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
| chrom | chromStart | chromEnd | name | score |
|---|---|---|---|---|
| chr1 | 1243287 | 1244338 | MACS_peak_2 | 63.19 |
| chr1 | 2979976 | 2981228 | MACS_peak_3 | 100.16 |
| chr1 | 3566181 | 3567876 | MACS_peak_4 | 558.89 |
| chr1 | 3816545 | 3818111 | MACS_peak_5 | 57.57 |
| chr1 | 6304864 | 6305704 | MACS_peak_6 | 54.57 |
| chr1 | 10850776 | 10853929 | MACS_peak_7 | 160.10 |
| chr1 | 15899600 | 15900536 | MACS_peak_8 | 67.20 |
| chr1 | 21415285 | 21416273 | MACS_peak_9 | 50.85 |
| chr1 | 21671455 | 21673405 | MACS_peak_10 | 51.76 |
支持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