Signac/DeepTools data exchange

Question

[andrea.rubiop88]

Hi all,

I'm trying to analyze some DOGMA single cell data (GEX+ ATAC). I have been following the tutorials from Seurat + Signac packages to analyze the ATAC part of the data (this one and this one), and so far everything is working fine.
My problem is that I want to plot the differentially accessible peaks in my conditions in a heatmap like the ones output by deepTools. In Signac, they offer a similar thing by computing a RegionMatrix, which counts fragments within a set of regions for different conditions. It then outputs one matrix per condition with the counts per region. This seems to be the same thing done by deepTools function computeMatrix. The thing is I can't find an R approach to plot that data in a deepTools "copycat" way. The function offered by Signac to do that (RegionHeatmap) doesn't offer as much versatility and is definitely not as beautiful/informative. The only solution I can think of is trying to tweak the matrix output by Signac so that it can be an input for deepTools, but first I want to make sure I understood the outputs correctly. In the deepTools documentation for their function they speak of "scores" while in Signac they say they use fragment counts. Are those the same thing? And if not, is it a way to transform one to the other?
I'm also concerned about the conversion itself, do you think it is feasible to add the necessary headers etc to a fragment count matrix file so that deepTools can plot it?
As you may have guessed, I am a complete novice in analyzing ATAC data, which is part of the reason I used SIgnac for the analysis, as it was similar to a standard Seurat analysis, with which I am very familiar.
Any info or insight will be appreciated. Thanks!

[DARLOR]

Hi, did you find a more effiecient solution than RegionMatrix used by Signac? However, I could suggest to you to convert the matrix produced by RegionMatrix() into a matrix for EnrichedHeatmap package.

library(EnrichedHeatmap)
try = RegionMatrix(seu_obj, regions = peaks_gr, assay="ATAC", group.by = "condition_of_interest",
                     idents = c("idents_of_interest1", "idents_of_interest_2"), key = "peaks", upstream = 1000, downstream = 1000, verbose = TRUE)


b=try@assays$ATAC@positionEnrichment$peaks$`idents_of_interest1`
mat= as.normalizedMatrix(b, k_upstream = 1001, k_downstream = 1000, k_target = 0, extend = 1000, smooth = TRUE, keep = c(0, 0.95))
EnrichedHeatmap(mat, name = "title", column_title = "Enrichment of ....")

that works fine producing nice heatmaps. The problem is RegionMatrix doesn't work if you put the whole gene body. I don't really understand the mechanism that they use to create the regionamtrix.