第4步用bedtools得到了hnRNPK与CBFB公有peak，结果为merge后的一致性的hnRNPK_CBFB_FLAG_antibody.bed文件，接下来就是对peaks的注释。
这篇主要用Y叔的R包ChIPseeker对peaks的位置（如peaks位置落在启动子、UTR、内含子等）以及peaks临近基因的注释。
数据准备
在用ChIPseeker包进行注释前，需要准备两个文件：
1. 注释peaks的文件
该文件需要满足BED格式。BED格式文件至少得有chrom（染色体名字），chromStart（染色体起始位点）和chromEnd（染色体终止位点），其它信息如name，score，strand等可有可无。一般情况而言，可直接用做peak calling的MACS输出文件（以_peaks.bed结尾文件）。这里用的是merge后的一致性的hnRNPK_CBFB_FLAG_antibody.bed文件
2.TxDb文件
TxDb.Hsapiens.UCSC.hg19.knownGene
> setwd("E:/ZYH/处理流程/RIP-Seq/练习1")
> library(TxDb.Hsapiens.UCSC.hg19.knownGene)
> txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
载入各种包
library("ChIPseeker")
library(clusterProfiler)
library("org.Hs.eg.db")
读入peaks文件
bedfile <- readPeakFile("hnRNPK_CBFB_FLAG_antibody.bed")
注释peaks
这里需要注意的是，启动子区域是没有明确的定义的，所以你可能需要指定tssRegion，把基因起始转录位点的上下游区域来做为启动子区域。
有了这两个输入（BED文件和TxDb对象），你就可以跑注释了，然后就可以出结果了。
peakAnnoList <- annotatePeak(bedfile, tssRegion=c(-5000, 5000), TxDb=txdb, annoDb="org.Hs.eg.db")
可视化,展示peak在基因组特征区域的分布以及转录因子在TSS区域的结合
library(ggplot2)
#饼图
plotAnnoPie(peakAnnoList)
#nearest gene annotation
plotDistToTSS(peakAnnoList,title="Distribution of TF-binding loci relative to TSS")
#covplot函数可以直接接受macs2-callpeak的输出bed文件进行画图。
#函数返回的图中，每个竖线表示一个peaks在染色体的位置，线的宽度表示peaks在染色体的分布范围。
covplot(bedfile)
#具体到指定染色体上的指定范围
covplot(bedfile, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4e7, 5e7))+ facet_grid(chr ~ .id)
#covplot选项：
#chrs 指定想要展示染色体，默认为全部
#xlim 指定x轴坐标范围，即染色体的范围
#如果需要ggplot刻面展示，不同样本是用id来表示
image.png
image.png
image.png
网友评论