How to integrate more than two scATAC-seq datasets?

[zrcjessica]

I was wondering how to best integrate more than two datasets. I've got multiple scATAC-seq datasets, each from a different individual classified as either one of two categories (i.e. different reactions when exposed to the same conditions). I want to make comparisons between the two categories.

The vignette for integration uses the findOverlaps() function to identify intersecting regions between the two datasets in the vignette and then creates a new assay for the sci-ATAC-seq dataset using only the peaks from that dataset that overlapped any of the intersecting regions (intersecting.regions). How would you advise efficiently finding overlaps between more than two datasets? Should I simply merge them and use the reduce() function as used in the vignette for merging datasets to generate the combined.peaks object?

Thank you for your input!

[timoast]

I would recommend first creating a unified set of peaks across all of the datasets, then quantifying the unified peak set in each dataset, and merging the resulting Seurat objects together. You can use the GenomicRanges::reduce() function for this, and there's an example in the merge vignette: https://satijalab.org/signac/articles/merging.html

If you find that you see a batch effect (cells separate by both cell state and dataset of origin) after merging the objects, then you could also apply data integration methods to remediate this. For >2 scATAC datasets, I'd recommend trying Harmony (example here). We're working on some updates to the Seurat integration to better support single-cell chromatin data, but these updates are not available yet.

[Dragonmasterx87]

So @timoast This is what I do for 4 datasets. So far it works alright, the issue is you need an HPC if you don't subset the Features or Peaks so I downsample using FindTopFeatures and then just use intersections amongst all VariableFeatures. Works fine, and I get good batch correction and identity preservation across clusters and origin. As always feel free to be brutal if I did something stupid.

# As outlined in #35 Signac
# convert to genomic ranges
peaks1 <- StringToGRanges(regions = rownames(cm1_atac))
peaks2 <- StringToGRanges(regions = rownames(cm2_atac))
peaks3 <- StringToGRanges(regions = rownames(cf1_atac))
peaks4 <- StringToGRanges(regions = rownames(aam1_atac))

# Create a unified set of peaks to quantify in each dataset
plan("multiprocess", workers = 12)
combined.peaks <- reduce(x = c(peaks1, peaks2, peaks3, peaks4))

# Newcounts based on unions
newcounts1 <- FeatureMatrix(fragments = Fragments(cm1_atac), features = combined.peaks, cells = colnames(cm1_atac))
newcounts2 <- FeatureMatrix(fragments = Fragments(cm2_atac), features = combined.peaks, cells = colnames(cm2_atac))
newcounts3 <- FeatureMatrix(fragments = Fragments(cf1_atac), features = combined.peaks, cells = colnames(cf1_atac))
newcounts4 <- FeatureMatrix(fragments = Fragments(aam1_atac), features = combined.peaks, cells = colnames(aam1_atac))

cm1_atac[['peakunion']] <- CreateAssayObject(counts = newcounts1)
cm2_atac[['peakunion']] <- CreateAssayObject(counts = newcounts2)
cf1_atac[['peakunion']] <- CreateAssayObject(counts = newcounts3)
aam1_atac[['peakunion']] <- CreateAssayObject(counts = newcounts4)

### ENV saved: atac_env4_R3 ####
DefaultAssay(cm1_atac) <- "peakunion"
DefaultAssay(cm2_atac) <- "peakunion"
DefaultAssay(cf1_atac) <- "peakunion"
DefaultAssay(aam1_atac) <- "peakunion"
cm1_atac <- RunTFIDF(cm1_atac)
cm1_atac <- FindTopFeatures(cm1_atac, min.cutoff = 20)
cm2_atac <- RunTFIDF(cm2_atac)
cm2_atac <- FindTopFeatures(cm2_atac, min.cutoff = 20)
cf1_atac <- RunTFIDF(cf1_atac)
cf1_atac <- FindTopFeatures(cf1_atac, min.cutoff = 20)
aam1_atac <- RunTFIDF(aam1_atac)
aam1_atac <- FindTopFeatures(aam1_atac, min.cutoff = 20)

unintegrated <- merge(
  x = cm1_atac,
  y = list(cm2_atac, cf1_atac, aam1_atac),
  add.cell.ids = c("cm1", "cm2", "cf1", "aam1")
)

DefaultAssay(unintegrated) <- "peakunion"
unintegrated <- RunTFIDF(unintegrated)
unintegrated <- FindTopFeatures(unintegrated, min.cutoff = 50)
unintegrated <- RunSVD(unintegrated)
unintegrated <- RunUMAP(unintegrated, reduction = 'lsi', dims = 2:30)
p1 <- DimPlot(unintegrated, group.by = 'sample', pt.size = 0.1) + ggplot2::ggtitle("Unintegrated")
p1

#Integration
peaks.use <- Reduce(intersect, list(VariableFeatures(cm1_atac), VariableFeatures(cm2_atac), VariableFeatures(cf1_atac), VariableFeatures(aam1_atac)))
peaks.use <- sample(peaks.use, 20000)

# find integration anchors between 10x and sci-ATAC
anchors <- FindIntegrationAnchors(
  object.list = list(cm1_atac, cm2_atac, cf1_atac,aam1_atac),
  anchor.features = peaks.use,
  assay = c('peakunion', 'peakunion', 'peakunion', 'peakunion'),
  k.filter = NA
)


# integrate data and create a new merged object
integrated <- IntegrateData(
  anchorset = anchors,
  weight.reduction = cm1_atac[['lsi']],
  dims = 2:30,
  preserve.order = TRUE
)

integrated <- IntegrateData(anchorset = anchors, dims = 2:30, preserve.order = TRUE)

Also sorry I forgot to email you back that this seems to work. I'm not a big fan of harmony, somehow the clusters 'Meld' into one another more than Signac. I'm looking into ArchR these days, it looks like it handles multi-dataset integration well. Really looking forward to your update mate, because I will need it to integrate ~9 datasets soon, and merging doesn't cut it.

Cheers for the help.
🐉

[timoast]

This looks reasonable, if you're having issues with RAM limits you could also remove any unneeded assays (eg, the original ATAC assay before you defined common peaks) by setting obj[[assay]] <- NULL. You can also try using Harmony rather than Seurat integration, we're working on more scalable integration methods for scATAC that will hopefully be available soon

[Dragonmasterx87]

I think co-embedding with scRNA is great tbh, but reviewers these days..... anyway thanks for the protip, also please incorporate control over RAM with Future or something, so that we can control analysis on non-HPCs. Will try ArchR and Harmony. Have a smooth weekend!

Cheers!
🐉

[zrcjessica]

Thanks to you both for your responses! @timoast, if I were to use Harmony, do I need to run RunTFIDF, FindTopFeatures, RunSVD, and RunUMAP on each individual dataset prior to merging and running Harmony on the merged dataset?

[timoast]

@zrcjessica You'd need to merge all the datasets, then process until you have the LSI, then you can run Harmony and run UMAP on the harmony-corrected LSI coordinates