Strange behavior with scATAC integration using Harmony

[jeremymsimon]

Hi @timoast
I'm currently using Harmony to integrate multiple scATACseq samples together, but I'm observing there are several instances where clusters are miles apart in UMAP space (despite all the caveats of that observation) but don't seem to have any peaks that distinguish them; nCount_peaks and nFeature_peaks however are quite different between them. Look at e.g. clusters 6 and 21 below.

I started to wonder whether I should be modifying the dims.use parameter in RunHarmony to omit the first PC, like we do in other steps as described in the scATAC integration vignette. The vignette describes Seurat-based integration, but the section on integrating LSI embeddings suggests dims = 1:30 rather than 2:30. Is this correct, or a typo? If correct, is there a rationale for using all 30 PCs for integration but omitting the first for subsequent steps?

Or is there something else here I've missed that could be causing this behavior?

My process is as follows:

[...]
chromap.combined.singlets <- RunTFIDF(chromap.combined.singlets, assay = "peaks")
chromap.combined.singlets <- FindTopFeatures(chromap.combined.singlets, min.cutoff = "q0", assay = "peaks")
chromap.combined.singlets <- RunSVD(chromap.combined.singlets, assay = "peaks", reduction.name = "lsi_singlets")
chromap.combined.singlets <- RunUMAP(chromap.combined.singlets, assay = "peaks", reduction = "lsi_singlets", dims = 2:30, reduction.name = "umap_singlets")

integrated <- RunHarmony(
  object = chromap.combined.singlets,
  group.by.vars = 'Sample',
  reduction.use = 'lsi_singlets',
  assay.use = 'peaks',
  project.dim = FALSE
)
integrated <- RunUMAP(integrated, dims = 2:30, reduction = 'harmony', reduction.name = 'harmony.umap')
integrated <- FindNeighbors(integrated, reduction = 'harmony', dims = 2:30)
integrated <- FindClusters(integrated, verbose = FALSE, algorithm = 3, resolution = 0.8)
DimPlot(integrated, reduction = "harmony.umap", label = TRUE)

I generated the following dotplot by running FindAllMarkers as follows:

markers <- FindAllMarkers(integrated,
                      assay="peaks",
                      min.pct = 0.2,
                      test.use = 'LR',
                      latent.vars = 'nCount_peaks',
                      only.pos = T, 
                      logfc.threshold = 0.25)

then plotting the top ~10 per cluster