这篇文章主要讲解了“如何使用ChIPpeakAnno进行peak注释”，文中的讲解内容简单清晰，易于学习与理解，下面请大家跟着小编的思路慢慢深入，一起来研究和学习“如何使用ChIPpeakAnno进行peak注释”吧！
ChIPpeakAnno是一个bioconductor上的R包，针对peak calling之后的下游分析，提供了以下多种功能
查找与peak区域最相邻的基因, 也支持自定义查找的特征，可以是exon,miRNA等
peak相邻基因的GO富集分析
提取peak及其周围区域的序列
在ChIPpeakAnno中，无论是peak区间信息还是基因组的注释信息，都通过toGRanges方法转化为R语言中的GRanges对象，以peak为例，bed格式的内容如下
通过如下代码可以导入该信息
library(ChIPpeakAnno)bed <- "peaks.bed"
gr <- toGRanges(bed, format="BED", header=FALSE)
除了BED格式外，该方法也支持导入GTF格式的信息，只需要修改format参数即可。导入peak信息和基因组注释信息后就可以进行后续分析了。
1. 进行peak之间的overlap分析当导入了多个样本的peak信息时，可以进行venn分析，用法如下
# 导入A样本的peakbedA <- "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]])
输出结果示意如下
首先是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))
输出结果示意如下
郑重声明：本文版权归原作者所有，转载文章仅为传播更多信息之目的，如作者信息标记有误，请第一时间联系我们修改或删除，多谢。
