Normalised Gene Activity matrix still affected by peak region fragments

[wudustan]

Following the vignette I generate my Gene Activity matrix from my scATAC dataset. I then NormalizeData() using the standard approach.

However, when I plot genes, I can clearly see that parts of the library with higher peak_region_fragments values also have higher gene activity scores after normalisation.

Is it possible to correct this?

Example:

Peak region fragments plot

Gene signature from normalized Gene Activity

[timoast]

It's difficult to fully remove these differences due to different sequencing depth using a simple method like that used by NormalizeData(). However, for most downstream applications, these differences are not usually an issue since the biological variation is typically a much stronger signal and are consistent across datasets.

[wudustan]

Do you have a recommendation for a package / method to use the peak_region_fragments as a latent variable to compare the gene expression .

[timoast]

Yes, Seurat. You can include a latent variable in a few of the DE methods in FindMarkers(), eg FindMarkers(..., test.use="LR", latent.vars="peak_region_fragments")