2022-05-19 Peak注释学习总结
三个讲的很好的参考教程链接
学习一遍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)
