chipseeker annotation issue

Question

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
library(clusterProfiler)
library(org.Hs.eg.db)

peaks <- readPeakFile("MACS_folder/pol2_chip/Pol_2pausing/pol2_diff_peak.bed")

print(peaks)


#peakAnno <- annotatePeak("merge/blast.bed", tssRegion=c(-3000, 3000),
#                        TxDb=txdb, annoDb="org.Hs.eg.db")

peakAnno <- annotatePeak(peaks, tssRegion=c(-1000, 1000),
                         TxDb=txdb, annoDb="org.Hs.eg.db")

So this is basically what im doing to annotate the peaks.

Now the issue comes with annotation for example I see TLR4 annotated to this

118082896   118083202   307 *   shp1_R1_peak_679    Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9   117708197   117708979   783 1   7099    ENST00000490685.1   374699  ENSG00000136869 TLR4

Now The peak is 679 as you can see from the IGV plot as well as annotated file. But when i see in IGV browser TLR4 is far away from the peak where i get the annotation.

So now is it an issue with the code im using or due to the inherent nature of the data.

Any help or suggestion would be really appreciated.

Answer 1

I had a similar case just recently.

The tssRegion argument ChIPseeker is used to define what is considered a promoter, which is frequently the question at hand. However, ChIPseeker will not perform any kind of filtering on the annotated peak file.

In the results there is a column annotation that contains information what kind of genomic region the peak is in. If it is within your tssRegion window it will be labeled as a promoter.
In your case the annotation is 1st Exon = Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9. This tells you that the peak is not part of the promoter region you defined and located within the first exon of the indicated transcript. The peak is really far away from TLR4 (374699 bp) but I think all other sequences nearby are pseudogenes and this is why TLR4 is selected as the closest gene.

You can get an overview over the annotations using plotAnnoPie(peakAnno)

tl;dr: you need to filter the results of the ChIPseeker annotation to get the genomic regions you want.