ATAC seq peak annotation : Which peak to choose for downstream analysis

Question

After doing peak call I'm annotating the peaks using chipseeker tool, which I want to take further downstream analysis. Now unlike RNA seq where i have single gene and their respective expression. But here the peaks for a single gene multiple peaks. As an example have filtered promoter peaks.

seqnames  start    end width strand     V4 annotation geneChr geneStart geneEnd geneLength geneStrand geneId
1     chr1 826797 828101  1305      *  Peak1   Promoter       1    826832  852225      25394          1 643837
2     chr1 869647 870324   678      *  Peak4   Promoter       1    868240  870201       1962          2 284593
3     chr1 876462 877140   679      *  Peak5   Promoter       1    868071  876903       8833          2 284593
4     chr1 877219 878105   887      *  Peak6   Promoter       1    874529  877234       2706          2 284593
5     chr1 923459 924318   860      * Peak14   Promoter       1    923928  939291      15364          1 148398
6     chr1 924379 926062  1684      * Peak15   Promoter       1    925150  935793      10644          1 148398
       transcriptId distanceToTSS         ENSEMBL    SYMBOL                                    GENENAME
1 ENST00000623808.3             0 ENSG00000228794 LINC01128 long intergenic non-protein coding RNA 1128
2 ENST00000432963.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
3 ENST00000446136.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
4 ENST00000635557.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
5 ENST00000420190.6             0 ENSG00000187634    SAMD11    sterile alpha motif domain containing 11
6 ENST00000437963.5             0 ENSG00000187634    SAMD11    sterile alpha motif domain containing 11

Now if we see first 6 rows we have two genes with three different peaks. When i load and see in the igv browser I do see the there are subtle differences in the accessibility signals.

How to address this issue ,in other words which peak to consider or take all the peak and consider their accessibility as one.

Any suggestion or help would be really appreciated.

Answer 1

I'll consider all peaks around promoter of particular gene as one peak. Moreover I'll interpret this as an open chromatin region = active gene. Why? Why there's no one wide peak which represents open chromatin in promoter of particular gene? Why is divided into several narrower peaks? There might be some transcription factors or other proteins binded to DNA, hence you see peak (open, free DNA) - break (valley, TF binded) - peak(again free DNA where Tn5 cut DNA and incorporate adapters) in IGV. In my study I compared ATAC-seq and H3k27Ac. In such case usually you may notice that peaks from chip and atac never literally overlap. Because ATAC represents open chromatin (DNA without any binded proteins) and chip histone modification (DNA is actually wound on histone proteins), hence we have something like this figure C, first we see high ATAC peak - open chromatin, then H3k27ac peak - lysine acetylated = gene activated. In other words ATAC-seq won't show you peaks around histones, binded polymerase, or transcription factors.

Answer 2

You can see from your table that the different peaks correspond to different promoters for different transcripts, for example, FAM41C has ENST00000432963.1, ENST00000446136.1, ENST00000635557.1. And the fact that they are all within the TSS (see distanceToTSS) suggest these peaks are actually falling onto actually TSS.

The underlying issue is are these TSS or do will we mis classify the peaks that are inside genes or other non-coding elements as TSS peaks. Nothing to do with TF, histones etc.

I cannot see the number of reads in those peaks, but if one of them is really strong, and the rest almost non-existent it would make sense.

How can we improve this, I think it depends on the aim of your analysis. Since you are looking at TSS, I guess it makes sense to be able to link it to something that is actually expressed. You can do:

1. To use an updated annotation (i see that some of the transcripts come from old hg37 annotation)

2. you have prior knowledge which transcripts might be expressed from RNA seq data, I would use that to reduce the number of transcripts