Harmony batch correction for integration of snATAC-seq and scRNA-seq data

[jacklyn-ashmus]

Hi,
I've been having some issues trying to get Harmony to work correctly using Seurat. I can get it to batch correct our RNA data, but when I use Harmony through Seurat, I somehow can't get it to correct the batch effects at all. The resulting atac-seq UMAP plots look kind of like spaghetti and meatballs. This problem persists even after I transfer anchors and merge the two data sets, and then call Harmony. I have tried both
hm_atac <- HarmonyMatrix(Embeddings(seurat_ATAC, reduction='lsi'), seurat_ATAC$library, do_pca=FALSE)
and
atac_adata <- RunHarmony(atac_adata, group.by.vars='library', reduction.save='harmony')
after calling normalization and PCA, and after making sure that the correct variables are in the 'library' metadata column and that I have the correct assay active.

I should probably also mention that I've successfully run Harmony on this ATAC data set outside of the context of Seurat before, where the only difference in commands was to include the parameter tau=2. I tried that here and it didn't make a difference.

I've included my code and some figures in this post, although I'm fairly new to github, so I apologize if anything is wrong with this.

Edit:
OK, here are the fragments of code that I have been using that I think are involved in the issue:

reading in the h5ad files (done within a loop, for every sample):

atac_seurat_list <- list()

loop:

    h5FileName <- sprintf('/dev/ATAC_MTX/%s.raw_peak_bc_matrix.h5', sample)
    h5 <- Read10X_h5(h5FileName, use.names = TRUE, unique.features = TRUE)
    pbmc_adata_atac_single <- CreateSeuratObject(counts = h5, assay = 'ATAC', min.cells = 3, min.features = 200)
    pbmc_adata_atac_single[['percent.mt']] <- PercentageFeatureSet(pbmc_adata_atac_single, pattern = '^MT-')
    DefaultAssay(pbmc_adata_atac_single) <- 'ATAC'
    atac_seurat_list[[sample]]<- pbmc_adata_atac_single

merging the obects together:

mergeLeftList <- c(atac_seurat_list$pbmc2, atac_seurat_list$pbmc3, atac_seurat_list$pbmc4, 
                  atac_seurat_list$pbmc5, atac_seurat_list$pbmc6, atac_seurat_list$pbmc8, 
                   atac_seurat_list$pbmc9, atac_seurat_list$pbmc10, atac_seurat_list$pbmc12)
atac_adata <- merge(
  x = atac_seurat_list$pbmc1,
  y = mergeLeftList
) 

normalization and running Harmony, etc:

atac_adata <- RunTFIDF(atac_adata)
atac_adata <- FindTopFeatures(atac_adata, min.cutoff = "q0")
atac_adata <- RunSVD(atac_adata)
atac_adata <- RunHarmony(atac_adata, group.by.vars='library', reduction.save = "harmony", assay.use="ATAC")
atac_adata <- RunUMAP(atac_adata, reduction = "lsi", dims = 2:30, reduction.name = "umap.atac", reduction.key='atacUMAP_')
atac_adata <- FindNeighbors(object = atac_adata, reduction = 'harmony', dims = 2:30)
 atac_adata <- FindClusters(object = atac_adata, algorithm=4, verbose=FALSE, method="igraph", resolution=0.5)

[timoast]

Hi, could you provide the minimal code required to reproduce the issue you're having, and put it in the post itself rather than a pdf?

[jacklyn-ashmus]

Hi, I am really sorry about that. I've fixed it to only include the relevant bits of code now.

[timoast]

I would recommend quantifying the same peak set across the different samples you have before merging them. See the examples given here: https://satijalab.org/signac/articles/merging.html