使用ChIPpeakAnno进行peak注释
欢迎关注”生信修炼手册”!
ChIPpeakAnno属于生物信息学领域的R工具包，在生物可 accessibility平台（Bioconductor）上提供；它用于完成peak calling后的后续数据分析，并包括多种 downstream analysis功能。
搜索与peak区域直接相邻的基因，并且允许自定义查询功能，并且可以选择exons和miRNAs作为搜索对象
-
peak相邻基因的GO富集分析
-
提取peak及其周围区域的序列
在ChIPpeakAnno中, 包括_peak区间信息以及基因组注释信息均通过toGRanges方法被转化为了R语言中的GRanges对象. 以peaks为例, bed格式的具体内容如下
通过如下代码可以导入该信息
library(ChIPpeakAnno)
bed <- "peaks.bed"
gr <- toGRanges(bed, format="BED", header=FALSE)除了以BED格式存储的文件外（BED），该方法还同样支持通过指定GTF格式文件导入相关数据（GTF）。在完成peak信息与基因组注释数据的导入后（这些数据通常用于后续分析），就可以直接开始进行后续的数据处理工作了。
1. 进行peak之间的overlap分析
当导入了多个样本的peak信息时，可以进行venn分析，用法如下
# 导入A样本的peak
bedA <- "sampleA_peaks.bed"
sampleA <- toGRanges(bedA, format="BED", header=FALSE)
# 导入B样本的peak
bedB <- "sampleB_peaks.bed"
sampleB <- toGRanges(bedB, format="BED", header=FALSE)
# 求交集
ol <- findOverlapsOfPeaks(sampleA, sampleB)
# 绘制venn图
makeVennDiagram(ol)结果示意如下
当执行venn分析时，能够识别出所有venn图上区域的数量之和并不等于输入的所有peak区间数量，在默认情况下
2. 提取peak周围的序列
用法如下
library(BSgenome.Hsapiens.UCSC.hg19)
seq <- getAllPeakSequence(sampleA, upstream=20, downstream=20, genome=Hsapiens)
write2FASTA(seq, "sampleA.peaks.fa")3. 进行peak motif分析
提取到peak序列之后，可以进行motif分析，用法如下
# 用1号染色体的碱基分布当做背景
freqs <- oligoFrequency(Hsapiens$chr1, MarkovOrder=3)
# oligoLength规定了motif的长度
os <- oligoSummary(seq, oligoLength=6, MarkovOrder=3,
quickMotif=TRUE, freqs=freqs)
zscore <- sort(os$zscore)
# 绘制所有6个碱基组合的频率分布图
h <- hist(zscore, breaks=100, xlim=c(-50, 50), main="Histogram of Z-score")
# 频率最大的碱基组合即为motif的结果
text(zscore[length(zscore)], max(h$counts)/10,
labels=names(zscore[length(zscore)]), adj=1)结果示意如下
还可以通过motifStack这个R包绘制motif的sequence logo, 用法如下
library(motifStack)
pfms <- mapply(function(.ele, id)
new("pfm", mat=.ele, name=paste("SAMPLE motif", id)),
os$motifs, 1:length(os$motifs))
motifStack(pfms[[1]])输出结果示意如下
4. 进行peak注释
首先是peak在基因组各个特征区间的分布比例，用法如下
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
aCR<-assignChromosomeRegion(sampleA, nucleotideLevel=FALSE,
precedence=c("Promoters", "immediateDownstream",
"fiveUTRs", "threeUTRs",
"Exons", "Introns"),
TxDb=TxDb.Hsapiens.UCSC.hg19.knownGene)
barplot(aCR$percentage, las=3)输出结果如下所示
然后进行peak关联基因的注释，用法如下
# 准备基因组注释信息
library(EnsDb.Hsapiens.v75)
annoData <- toGRanges(EnsDb.Hsapiens.v75, feature="gene")
# 进行
overlaps.anno <- annotatePeakInBatch(sampleA,
AnnotationData=annoData,
output="nearestLocation"
)
library(org.Hs.eg.db)
overlaps.anno <- addGeneIDs(overlaps.anno,
"org.Hs.eg.db",
IDs2Add = "entrez_id")
pie1(table(overlaps.anno$insideFeature))输出结果示意如下
当调用annotatePeakInBatch进行标注操作时，默认会检索与peak位置最接近的基因。然而，在某些情况下还可以通过调整参数来实现特定功能。其中参数overlapping用于指定那些与peak区域存在重叠关系的具体条件。当设置相应的参数后，则会将所有与当前peak区间有重叠关系的基因标记为相关联。此外该方法还支持多种不同的取值设置，在不同的应用场景下都可以灵活应用。
5. 进行peak关联基因的富集分析
完成基因注释后就能获得与peak相关的基因信息，并可开展后续的功能富集分析，请参见下文以了解具体操作方法。
over <- getEnrichedGO(overlaps.anno, orgAnn="org.Hs.eg.db",
maxP=.05, minGOterm=10,
multiAdjMethod="BH", condense=TRUE)ChIPpeakAnno提供了这一项全面的peak下游功能包, 包括基因注释、富集统计以及motif检测等多方面内容, 这是一个高度集成的功能集合. 作为一项downstream analysis tool, 其基础操作相对简单明了, 但深入使用时还有很多细节值得探索, 可参考官方文档获取详细指导.
·end·
—如果喜欢，快分享给你的朋友们吧—
扫描关注微信号，更多精彩内容等着你！
