Hi there.
I am using two tools to annotate peaks: annotatepeak from the ChIPseeker library and bedmap.
The codes I use are these. For annotatepeak:
#### I load the .bed file that has 3 columns: Chr, start, end for the start and end positions. My bed file is called bed
peak = GRanges(seqnames = bed$Chr, ranges = IRanges(start = bed$start, end = bed$end)) # convert the bed database in granges object
library(GenomicFeatures)
txdb = makeTxDbFromGFF(file = "my_gtf.gtf", format = "gtf", dataSource = "Annotation file",taxonomyId = 123456, organism = "R Phil")
peakAnno <- annotatePeak(peak = peak, tssRegion=c(-1000, 1000), TxDb=txdb)
For bedmap:
bedmap --echo --echo-map-id --range 1000 Mypeaks.bed my_gtf.gtf.sorted.bed > Mypeaks.wGenes.bed
Judging by the two outputs, it seems that bedmap gives an output that makes more sense because if I look at the coordinates within the .gtf some genes that are correctly named with bedmap (because they fall within the gene coordinates of that specific gene), with annotatepeak are labeled distalintergenic and are given the ID of a gene which is much more distant... So I am not sure how annotatepeak is choosing how to annotate genes?
Any help would be much appreciated
Thanks
Luca
Hi there.
I am using two tools to annotate peaks:
annotatepeakfrom the ChIPseeker library andbedmap.The codes I use are these. For annotatepeak:
#### I load the .bed file that has 3 columns: Chr, start, end for the start and end positions. My bed file is called bedpeak = GRanges(seqnames = bed$Chr, ranges = IRanges(start = bed$start, end = bed$end)) # convert the bed database in granges objectlibrary(GenomicFeatures)txdb = makeTxDbFromGFF(file = "my_gtf.gtf", format = "gtf", dataSource = "Annotation file",taxonomyId = 123456, organism = "R Phil")peakAnno <- annotatePeak(peak = peak, tssRegion=c(-1000, 1000), TxDb=txdb)For bedmap:
bedmap --echo --echo-map-id --range 1000 Mypeaks.bed my_gtf.gtf.sorted.bed > Mypeaks.wGenes.bedJudging by the two outputs, it seems that
bedmapgives an output that makes more sense because if I look at the coordinates within the .gtf some genes that are correctly named withbedmap(because they fall within the gene coordinates of that specific gene), withannotatepeakare labeled distalintergenic and are given the ID of a gene which is much more distant... So I am not sure how annotatepeak is choosing how to annotate genes?Any help would be much appreciated
Thanks
Luca