33 ATAC Batch Correction (Harmony)
ATAC batch effects are assessed and corrected in this chapter. The ATAC modality cannot use the same CCA approach as RNA because peak matrices are sparser and lack the gene-level biological annotations needed to define integration anchors. Instead, Harmony is applied to the LSI embedding to align donors and samples in the latent space while preserving biological variation.
33.1 TF-IDF, LSI and ATAC UMAP
Run TF-IDF normalisation, select top features, compute SVD to generate the LSI embedding, and produce an uncorrected ATAC UMAP. The first LSI dimension is excluded from all downstream steps because it correlates with sequencing depth rather than biological variation.
combined <- LoadSeuratRds("data/seurat_object_donor_assigned.rds")
DefaultAssay(combined) <- "ATAC"
combined <- RunTFIDF(combined)
combined <- FindTopFeatures(combined, min.cutoff = 20)
combined <- RunSVD(combined)
combined <- RunUMAP(combined, dims = 2:50,
reduction = 'lsi',
reduction.name = "umap.atac",
reduction.key = "atacUMAP_")33.2 Uncorrected ATAC UMAP
Assess whether sample or donor effects are visible before applying batch correction.
DimPlot(combined, reduction = "umap.atac", group.by = "sample") +
ggtitle("ATAC UMAP (uncorrected) — by sample")
DimPlot(combined, reduction = "umap.atac", group.by = "sample_donor") +
ggtitle("ATAC UMAP (uncorrected) — by sample-donor")
DimPlot(combined, reduction = "umap.atac", group.by = "individual_condition") +
ggtitle("ATAC UMAP (uncorrected) — by condition")
33.3 Harmony batch correction on LSI
Harmony is applied to the LSI embedding (dims 2:50) to correct for
sample and donor batch effects. The corrected embedding is stored as
harmony_atac and a new UMAP (umap.atac.harmony) is computed from it.
library(harmony)
combined <- RunHarmony(combined,
group.by.vars = "sample_donor",
reduction = "lsi",
dims.use = 2:50,
reduction.save = "harmony_atac",
project.dim = FALSE)
n_harmony_atac <- ncol(Embeddings(combined, "harmony_atac"))
combined <- RunUMAP(combined, dims = 2:n_harmony_atac,
reduction = "harmony_atac",
reduction.name = "umap.atac.harmony",
reduction.key = "atacHarmonyUMAP_")
combined <- FindNeighbors(combined, reduction = "harmony_atac", dims = 2:n_harmony_atac,
graph.name = c("atac.harmony_nn", "atac.harmony_snn"))
combined <- FindClusters(combined, graph.name = "atac.harmony_snn",
resolution = seq(0.1, resolution, 0.1))33.4 Harmony-corrected ATAC UMAP
DimPlot(combined, reduction = "umap.atac.harmony", group.by = "sample") +
ggtitle("ATAC UMAP (Harmony) — by sample")
DimPlot(combined, reduction = "umap.atac.harmony", group.by = "sample_donor") +
ggtitle("ATAC UMAP (Harmony) — by sample-donor")
DimPlot(combined, reduction = "umap.atac.harmony", group.by = "individual_condition") +
ggtitle("ATAC UMAP (Harmony) — by condition")
33.5 Conclusion
Harmony correction on the ATAC LSI embedding, grouped by sample_donor, was evaluated but not adopted for the final WNN integration. Two considerations motivated this decision. First, the uncorrected LSI embedding (dimensions 2–50, excluding the depth-correlated first component) already showed substantial mixing of cells from Sample B, C, and D, indicating that sample-of-origin was not a dominant source of technical variance in this dataset. Second, and more critically, the sample_donor grouping variable is confounded with individual_condition, since each sample-donor combination corresponds to a single individual with a fixed PASC or Healthy status. Harmony correction on this variable therefore risks removing genuine disease-associated chromatin accessibility variation rather than purely technical batch effects. Consistent with this concern, the Harmony-corrected embedding produced a fragmented UMAP structure with numerous small, disconnected clusters not present in the uncorrected embedding, suggestive of overcorrection. The uncorrected lsi reduction (dims 2:50) was therefore retained for weighted nearest-neighbor integration with the integrated RNA PCA space.