在之前通过macs2得到了差异peaks的数据，现在利用ChIPseeker包进行注释。
参考：CS2: BED文件；CS3: peak注释；CS番外3：upsetplot + vennpie
- 加载R包
library(ChIPseeker)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(clusterProfiler)
#library(annotables)
library(org.Mm.eg.db)
library(ggplot2)
1. peak在整个染色体上的分布
#读取narrowPeak文件，默认GRanges格式：readPeakFile(peakfile, as = "GRanges", ...)
peak1 <- readPeakFile("D:\\1A\\macs2\\A_macs2_q0.1.peaks_peaks.narrowPeak")
peak2 <- readPeakFile("D:\\1A\\macs2\\B_macs2_q0.1.peaks_peaks.narrowPeak")
#covplot(peak1)仅显示分布，跟二维码一样，weightCol="V5"显示峰高
covplot(peak1, weightCol="V5")
#chr = c()选择展示哪些染色体
covplot(peak1, chr = c("chr1","chr2","chr3","chr4","chr5","chr6","chr7","chr8")
#多个文件绘制成一张图
peak=GenomicRanges::GRangesList(B=peak2,A=peak1)
covplot(peak, chr = c("chr3","chr4","chr5","chr6","chr7","chr8")) + facet_grid(chr ~ .id)
#, xlim=c(4e7, 5e7))展示哪块窗口
2. peak注释
require(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb = TxDb.Mmusculus.UCSC.mm10.knownGene
peakAnno = annotatePeak(peak1, tssRegion=c(-1000, 1000), TxDb=txdb,
#peak上下游某个范围内（比如说-5k到5k的距离）都有什么基因
addFlankGeneInfo=TRUE, flankDistance=5000,
#ID转换成SYMBOL
annoDb="org.Mm.eg.db")
peakAnno
as.GRanges(peakAnno) %>% head(3)
write.table(as.data.frame(peakAnno), file="D:\\E_macs_annotation.txt", sep="\t", quote=F, row.names=F)
3. 可视化
peakAnno <- annotatePeak(peak2, tssRegion=c(-3000, 3000), TxDb=txdb)
upsetplot(peakAnno, vennpie=T)
我们知道上游很重要啊，因为可能会调控转录，但注释的时候，没有上游这个东西，为什么呢？因为转录起始位点TSS的上下游被定义为promoter，所以啊上游被包括在promoter中，也就没有上游这个category了。
那么转录终止位点TTS的上下游呢？上游还在基因主体里，它可以是外显子、内含子、3'UTR这些，优先拿这些来注释，而下游呢？基因间区！基因间区就是各种不编码蛋白的区域，当然也可能编码一些非编码RNA之类的，这一块从基因的角度来看，是比较‘没用’的。但对于TTS后面紧接着的基因间区，它可能对基因的转录还是有些影响的，所以单独拿出来，就是这里要说的downstream了。
所以一个基因主体的immediate upstream，包含在promter里，而immediate downstream，我们也单独拿出来注释为downstream，这两块其实都在基因间区，但被我们拿出来了，因为和基因直接连接，很近的区域，可以说这是近端基因间区。而其它的基因间区，我们称之为远端基因间区，distal intergenic.
