How does SnapAtac2 snap.tl.merge_peaks work?

[1010shane]

Hello,

I have been playing around a bit with snapatac2 and have so far found it to be a very good tool for scATACseq analysis. However, I am struggling to understand how peaks called at the cluster level are merged to generate a consensus peak set.

Most notably, it seems that snap.tl.macs3 returns variable peak lengths when peak calling, but snap.tl.merge_peaks returns fixed width peaks the same size used in genome tiling (500 bp in my case).

I have created a small reprex demonstrating this which can be found here.

Can someone please shed some light on how snap.tl.merge_peaks works?

[kaizhang]

This function initially expands the summits of identified peaks by half_width on both sides. Following this expansion, it addresses the issue of overlapping peaks through an iterative process. The procedure begins by prioritizing the most significant peak, determined by the smallest p-value. This peak is retained, and any peak that overlaps with it is excluded. Subsequently, the same method is applied to the next most significant peak. This iteration continues until all peaks have been evaluated, resulting in a final list of non-overlapping peaks, each with a fixed width determined by the initial extension.

[1010shane]

Thanks! Very insightful. Quick follow up question–what is the rationale behind converting the variable-width peaks identified by macs3 into fixed width peaks?

By taking the approach you've just explained and artificially expanding the width of a considerable number of peaks called by macs3 (see density plot in the .ipynb file I've linked), I would think you would increase your rate of overlap and thus increase the rate of potentially valuable data being excluded from the consensus peak set. What is the benefit or justification of this?

[kaizhang]

Fixed-width non-overlapping peak list is usually desirable for downstream analyses, e.g., motif enrichment analysis. I agree with your concern that merging may cause some genuine peaks been removed. But you can mitigate this by setting a smaller peak width. There are also other strategies for peak merging. I'm open to implement different methods in SnapATAC2, as there isn't a universal method that works in all cases.

[TingTingShao]

I still don't understand why there is merge step, as at the beginning, the genome has already been segmented into tile matrixes, which for my understanding has already achieved fixed-width of non-overlapping peaks? and downstream analysis such as feature selection is based on the tile matrix.

Could someone explain this to me, thanks!

tingting

[adamklie]

Hi TingTing,

Here we are referring to the process of determining a uniform peak set after calling peaks at the pseudobulk level. Once we have determined how we want to group cells (via clustering of the tile matrix or other methods), we can export the fragments for each grouping of cells and call peaks on those. This will give us a set of genomic regions that are highly accessible in each group of cells we are interested in.

But now we run into the problem of having different sets of regions for each grouping of cells, and often times many of these peaks overlap between groups. So we need a procedure to create a uniform peak set that applies to all the cells. Simply merging all peaks leads to the problem of "daisy-chaining" peaks together that end up as super wide peaks in the final set (which is often undesirable). So we do the iterative procedure Kai described instead!

I like to think of this process as another version of the feature selection step that is less biased towards things that are accessible in the majority cell type(s) (i.e. what you are more prone to getting if you just look at the most accessible features in pp.select_features). Depending on your data, this may lead to better resolving rarer cell identities/states in clustering.

Hope this helps!

[TingTingShao]

Hi,

Thanks so much for this detailed explanation. I think now I have a much clearer view on this!

Regards,
tingting