Entering edit mode
Hello,
I tried to use "TxDb.Hsapiens.UCSC.hg38.knownGene" to annotate my peaks by ChIPpeakAnno 3.16.1.
My question is how I can annotate each peaks whether it is exon, intron, 3'UTR, 5'UTR or in promoter region?
Thanks!
yue
The following code snippets might fit your needs.
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
TxDb <- TxDb.Hsapiens.UCSC.hg38.knownGene
TSS <- genes(TxDb)
exons <- exons(TxDb)
fiveUTRs <- unique(unlist(fiveUTRsByTranscript(TxDb)))
threeUTRs <- unique(unlist(threeUTRsByTranscript(TxDb)))
introns <- unique(unlist(intronsByTranscript(TxDb)))
assuming peaksGR is a GRanges containing a set of peaks.
ann.TSS.peaks <- annotatePeakInBatch(peaksGR, AnnotationData = TSS, output= "upstream")
ann.exons.peaks.overlap <- annotatePeakInBatch(peaksGR, AnnotationData = exons, output= "overlapping", maxgap = -1)
ann.introns.peaks.overlap <- annotatePeakInBatch(peaksGR, AnnotationData = introns, output= "overlapping", maxgap = -1)
ann.utr5.peaks.overlap <- annotatePeakInBatch(peaksGR, AnnotationData = fiveUTRs, output= "overlapping", maxgap = -1)
ann.utr3.peaks.overlap <- annotatePeakInBatch(peaksGR, AnnotationData = threeUTRs, output= "overlapping", maxgap = -1)
Please note that you can set the output and maxgap to meet your needs. For example, if you would like to get peaks strictly inside the exons, then you would specify output = "inside".
ann.exons.peaksInsideExons <- annotatePeakInBatch(peaksGR, AnnotationData = exons, output= "inside)
ann.exons.peaksInsideIntrons <- annotatePeakInBatch(peaksGR, AnnotationData = introns, output= "inside)
For a complete set of options and documentation, please type the following commands in a R console.
args(annotatePeakInBatch)
help(annotatePeakInBatch)
If you are interested in the distribution of your peaks among different features, you could use the function
assignChromosomeRegion.
If you would like to find more information about this function, please type help(assignChromosomeRegion).