三个讲的很好的参考教程链接

学习一遍ChIPseeker的使用 https://www.jianshu.com/p/c76e83e6fa57

2021-05-23 ChIP-seq数据从头到尾比对及分析汇总（个人分析记录贴）https://www.jianshu.com/p/96688fecd864

使用R包ChIPseeker实现对ChIP-seq peaks的注释（注释ChIP-seq的峰位置（多个bed文件的批量处理））https://www.jianshu.com/p/0c71f1495dee

软件背景介绍

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

输入文件之BED文件

全称是：Browser Extensible Data，为基因组浏览器而生

包括3个必须字段和9个可选字段：

3个必须

• 1 chrom - 染色体名字

• 2 chromStart - 染色体起始位点（起始于0，而不是1）许多软件忽略了这一点，存在一个碱基的位移（如peakAnalyzer, ChIPpeakAnno存在这个问题），Homer、ChIPseeker没有

• 3 chromEnd - 染色体终止位点

9个可选

• 4 name - 名字

• 5 score - 分值(0-1000), 用于genome browser展示时上色。

• 6 strand - 正负链，对于ChIPseq数据来说，一般没有正负链信息。

• 7 thickStart - 画矩形的起点

• 8 thickEnd - 画矩形的终点

• 9 itemRgb - RGB值

• 10 blockCount - 子元件（比如外显子）的数目

• 11 blockSizes - 子元件的大小

• 12 blockStarts - 子元件的起始位点

一般只用前5个足矣（MACS的输出结果也是前5个字段）

分析代码实例

###加载R包

library(ChIPseeker)

library(ggplot2)

###读取数据

peak <- readPeakFile("ChIP-seq1.bed", as="GRanges")

#使用的是3.5中通过macs3得到的NAME_peaks.narrowPeak文件，as是指读取后转化成的对象 （Granges是genomic ranges的简称,它是bioconductor中的R包存储genomic intervals的数据结构）

###推荐手动构建TxDb注释文件

library(GenomicFeatures)

txdb <- makeTxDbFromGFF('XXX.gff3')

###peaks在基因组上的分布

pdf("Figure1.sample_distribution_of_peaks_in_chrs.pdf")

covplot(peak, weightCol="V5")

#covplot(peak, weightCol="V5", chrs=c("chr1", "chr2"), xlim=c(4.5e7, 5e7)) #chr1 chr2上的peaks

dev.off()

###多个文件同时画图

pdf("多个文件同时画染色体分布图")

peak=GenomicRanges::GRangesList(CBX6=readPeakFile("~/3_macs/sample1_peaks.narrowPeak", as="GRanges"),CBX7=readPeakFile("~/3_macs/sample2_peaks.narrowPeak", as="GRanges"))

covplot(peak, weightCol="V5") + facet_grid(chr ~ .id)

dev.off()

###TSS区域，上下游2.5kb，可自行设置

pdf("Figure2.sample_TSS_peaks_heatmap and plot.pdf")

##热图

tss <- getPromoters(TxDb=txdb, upstream=2500, downstream=2500)

tagMatrix <- getTagMatrix(peak, windows=tss)

#tagHeatmap(tagMatrix, xlim=c(-2500, 2500), color="red")

peakHeatmap(peak, TxDb=txdb, upstream=2500, downstream=2500, color="red") #等效语句

##峰图

plotAvgProf(tagMatrix, xlim=c(-2500, 2500),

            xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")

#含置信区间峰图

plotAvgProf(tagMatrix, xlim=c(-2500, 2500), conf = 0.95, resample = 1000)

dev.off()

###peaks注释,tss选取上下游2.5Kb，可自行调整，这里默认promoter是上游2.5Kb

peakAnno <- annotatePeak(peak, tssRegion=c(-2500, 2500), TxDb=txdb)

peakAnno

#在注释时，有的peak可能同时落在两个或者更多的gene feature里（例如是一个基因的外显子而同时又是另一个基因的内含子），但只能注释其中一个。

#默认按照Promoter、5’ UTR、3’ UTR、Exon、Intron、Downstream、Intergenic顺序先后注释。

#

### 保存注释结果

peakAnno.df <- as.data.frame(peakAnno) #保存为数据框

peakAnno.gr <- as.GRanges(peakAnno) #保存为GenomicRanges

head(peakAnno.gr, 3)

write.table(peakAnno.df,file="sample.annotation.xls",quote=F,sep="\t")

###peaks注释可视化(多种图形及最近基因)

pdf("Figure_3.sample_peaks_annotation_pie_bar_venn_upset.pdf")

plotAnnoPie(peakAnno)

plotAnnoBar(peakAnno)

vennpie(peakAnno)

upsetplot(peakAnno, vennpie=TRUE)

dev.off()

------------------------------------------------------------------------------------------------------------------------------------------

------------------------------------------------------------------------------------------------------------------------------------------

注释ChIP-seq的峰位置（多个bed文件的批量处理）

#读取chipseq峰的bed文件

peak1 <- readPeakFile('ChIP-seq1.bed')

peak2 <- readPeakFile('ChIP-seq2.bed')

peaks <- list(peak1 = peak1, peak2 = peak2)

#注释

peakAnnoList <- lapply(peaks, annotatePeak, TxDb = txdb, tssRegion = c(-3000, 3000), addFlankGeneInfo = TRUE, flankDistance = 5000)

#分别输出结果

write.table(peakAnnoList[1], file = 'peak1.txt',sep = '\t', quote = FALSE, row.names = FALSE)

write.table(peakAnnoList[2], file = 'peak2.txt',sep = '\t', quote = FALSE, row.names = FALSE)

#可视化，两个可以放在一起比较

plotAnnoBar(peakAnnoList)

vennpie(peakAnnoList[[1]])

plotAnnoPie(peakAnnoList[[1]])

plotDistToTSS(peakAnnoList)