Identifying clusters on scATAC-seq data integrated with Harmony

[zrcjessica]

Hello, I have integrated scATAC-seq data from several individuals using Harmony based on a unified set of all peaks observed across all my samples. I want to identify clusters in my data next so that I can perform label transfer from my scRNA-seq data and measure differential peak accessibility in each cell type. I noticed that in the Signac vignettes for merging and integrating datasets, as well as the Seurat vignette for integration, only RunUMAP() is applied to the merged/integrated dataset. My questions are as follows:

• If I want to identify clusters for cell type-identification, do I also need to run FindNeighbors() and FindClusters()? If so, should I do these steps before/after RunUMAP() on the merged/integrated object, or does the order not matter?
• If I am running FindNeighbors() on an object that was integrated with Harmony, should I set reduction = 'harmony' or reduction = 'lsi'?

I've tried running FindNeighbors() on my Harmony integrated object with both of the options described in the previous bullet point and I've noticed more clusters which are more clearly distinct from one another when I set reduction = 'harmony'. There are fewer clusters when I set reduction = 'lsi', but it looks like the clusters are messier and cells from different clusters colocalize.

Thank you for your help!

[timoast]

If you want to cluster using the harmony-corrected LSI, you can run FindNeighbors() with reduction="harmony", followed by FindClusters(). You can do this before or after computing the UMAP, they're independent.