Reproducing CoveragePlot results during differential peak analysis

[si3]

Hello

Thanks for developing this toolkit! I am trying to identify differentially accessible chromatin regions between cell subpopulation of an integrated Seurat object. The subpopulations are made up of cells from different scATAC-seq experiments and have different depths. Per the documentation, CoveragePlot displays aggregated pseudobulk values of Tn5 insertion sites on the X axis which are normalized using a "per-group scaling factor computed as the number of cells in the group multiplied by the mean sequencing depth for that group of cells".

You can see an example here, where clusters 0 and 8 show clear differences in peak height:

SOX2_Cov_plot.pdf

While I see clear differences in peak height in this plot, the same region does not appear as a differential peak using the standard FindMarkers workflow. In fact, I obtain almost no differentially accessible peaks in the subpopulation of interest (cluster 8) across the entire genome. I've read here on github and elsewhere that sequencing depth is a large factor influencing the number of peaks recovered, but that doesn't seem to be an issue when using the CoveragePlot function.

My question is the following:
How can I pseudobulk values and normalize values in such a way as to yield results that are consistent with what is obtained using CoveragePlot, even for cell clusters that have low coverage? I've tried using AggregateExpression and AverageExpression without success. There doesn't seem to be a vignette that specifically addresses how to do this type of pseudobulk analysis using the same normalization approach as the CoveragePlot function.