Harmony correction results in no UMAP separation

[RegnerM2015]

Hi @rcorces and @jgranja24,

Please note that I am using ArchR version 1.0.1

I have 17 samples profiled with 10x Genomics scATAC-seq. To mitigate batch/patient of origin effects, I decided to run Harmony after adding the LSI embeddings. There are 122,493 cells after quality control (minTSS = 6, minFrags = 1000, maxFrags =1e+05, doublets removed using filterDoublets() with default parameters). These seem like reasonable quality control cutoffs, but when I visualize the Harmony-based UMAP, I observe a large blob of cells with little separation. Please see the UMAP plots attached, highlighted by sample, by TSSEnrichment, and by log10(nFrags).

Which seems to be the most reasonable path forward?

1. Change quality control cutoffs to filter out even more "low-quality" barcodes (more stringent QC).
2. Explore iterative LSI and UMAP parameters until I find an "optimal" parameter set.
3. Explore UMAP visualizations without Harmony correction.

Below are the PDF files and code used to generate these plots:

Plot-UMAP-TSS.pdf
Plot-UMAP-Sample.pdf
Plot-UMAP-logFrags.pdf

library(scater)
library(dplyr)
library(Seurat)
library(patchwork)
library(SingleCellExperiment)
library(scater)
library(ComplexHeatmap)
library(ConsensusClusterPlus)
library(msigdbr)
library(fgsea)
library(dplyr)
library(tibble)
library(Signac)
library(ggplot2)
library(stringr)
library(EnsDb.Hsapiens.v86)
library(Seurat)
library(Signac)
library(ggplot2)
library(ensembldb)
library(EnsDb.Hsapiens.v86)
library(ArchR)
library(SingleR)
library(viridis)

set.seed(1234)
addArchRThreads(threads =64) 
addArchRGenome("hg38")

# Set up directories and file variables:
##################################################################################################
SAMPLE.ID <- "All"

output.dir <- "."
####################
setwd(output.dir)
####################

inputFiles <- list.files(pattern = "\\.gz$")

names <- str_remove(inputFiles,pattern = "fragments_")
names <- str_remove(names,pattern = ".tsv.gz")

names[c(1:11,16:17)] <- paste0("HT",names[c(1:11,16:17)])
names[12:15] <- paste0("HN",names[12:15])
sampleNames <- names

# Make patient sample metadata and color assignments 

sampleColors.1 <- RColorBrewer::brewer.pal(9,"Set3")
sampleColors.2 <- RColorBrewer::brewer.pal(8,"Set1")
sampleColors <- c(sampleColors.1,sampleColors.2)

# Create Arrow and ArchR project
##########################################################################
ArrowFiles <- createArrowFiles(
  inputFiles =inputFiles,
  sampleNames = sampleNames,
  outputNames = sampleNames,
  validBarcodes = NULL,
  geneAnnotation = getGeneAnnotation(),
  genomeAnnotation = getGenomeAnnotation(),
  minTSS = 6,
  minFrags = 1000,
  maxFrags =1e+05,
  QCDir = "QualityControl",
  nucLength = 147,
  promoterRegion = c(2000, 100),
  TSSParams = list(),
  excludeChr = c("chrM", "chrY"),
  nChunk = 5,
  bcTag = "qname",
  gsubExpression = NULL,
  bamFlag = NULL,
  offsetPlus = 4,
  offsetMinus = -5,
  addTileMat = TRUE,
  TileMatParams = list(),
  addGeneScoreMat = F,
  GeneScoreMatParams = list(),
  force = FALSE,
  threads = getArchRThreads(),
  parallelParam = NULL,
  subThreading = TRUE,
  verbose = TRUE,
  cleanTmp = TRUE,
  logFile = createLogFile("createArrows"),
  filterFrags = NULL,
  filterTSS = NULL
)


doubScores <- addDoubletScores(
  input = ArrowFiles,
  useMatrix = "TileMatrix",
  k = 10,
  nTrials = 5,
  dimsToUse = 1:30,
  LSIMethod = 1,
  scaleDims = FALSE,
  corCutOff = 0.75,
  knnMethod = "UMAP",
  UMAPParams = list(n_neighbors =30, min_dist = 0.3, metric = "cosine", verbose =FALSE),
  LSIParams = list(outlierQuantiles = NULL, filterBias = FALSE),
  threads = getArchRThreads(),
  force = F,
  parallelParam = NULL,
  verbose = TRUE,
  logFile = createLogFile("addDoubletScores")
)

proj <- ArchRProject(
  ArrowFiles = ArrowFiles,
  outputDirectory = "All",
  copyArrows = T #This is recommened so that you maintain an unaltered copy for later usage.
)
saveRDS(proj@cellColData,"proj_barcodes-meta.rds")

proj <- filterDoublets(proj,filterRatio = 1,cutEnrich = 1,cutScore = -Inf)
saveRDS(proj@cellColData,"proj_postfilter_doublets-meta.rds")

saveArchRProject(
  ArchRProj = proj,
  outputDirectory = getOutputDirectory(proj),
  overwrite = TRUE,
  load = TRUE,
  dropCells = T,
  logFile = createLogFile("saveArchRProject"),
  threads = getArchRThreads()
)


set.seed(1234)
addArchRThreads(threads =64) 
addArchRGenome("hg38")

# Set up directories and file variables:
##################################################################################################
SAMPLE.ID <- "All"

output.dir <- "."
####################
setwd(output.dir)
####################
addArchRThreads(threads =8)

proj <- loadArchRProject(path = "./All")
# # Perform LSI reduction and clustering with ATAC data only
# #######################################################################
# Add LSI dimreduc
proj <- addIterativeLSI(
  ArchRProj = proj,
  useMatrix = "TileMatrix",
  name = "IterativeLSI",
  iterations = 2,
  clusterParams = list(resolution = c(2), sampleCells = 10000, maxClusters = 6, n.start= 10),
  firstSelection = "top",
  depthCol = "nFrags",
  varFeatures = 25000,
  dimsToUse = 1:30,
  LSIMethod = 2,
  scaleDims = TRUE,
  corCutOff = 0.75,
  binarize = TRUE,
  outlierQuantiles = c(0.02, 0.98),
  filterBias = TRUE,
  sampleCellsPre = 10000,
  projectCellsPre = FALSE,
  sampleCellsFinal = NULL,
  selectionMethod = "var",
  scaleTo = 10000,
  totalFeatures = 5e+05,
  filterQuantile = 0.995,
  excludeChr = c(),
  saveIterations = TRUE,
  UMAPParams = list(n_neighbors =30, min_dist = 0.3, metric = "cosine", verbose =FALSE),
  nPlot = 10000,
  outDir = getOutputDirectory(proj),
  threads = getArchRThreads(),
  seed = 1,
  verbose = TRUE,
  force = T,
  logFile = createLogFile("addIterativeLSI")
)

proj <- addHarmony(
  ArchRProj = proj,
  reducedDims = "IterativeLSI",
  name = "Harmony",
  force = T,
  corCutOff = 0.75,
  groupBy = "Sample"
)

# Add UMAP based on LSI dims
proj <- addUMAP(proj,nNeighbors = 30,minDist = 0.3,
                dimsToUse = 1:30,metric = "cosine",
                force = T,reducedDims = "Harmony",
                corCutOff = 0.75)

p1 <- plotEmbedding(ArchRProj = proj, colorBy = "cellColData", name = "Sample", embedding = "UMAP")
plotPDF(p1, name = "Plot-UMAP-Sample.pdf", ArchRProj = proj, addDOC = FALSE, width = 5, height = 5)

p1 <- plotEmbedding(ArchRProj = proj, colorBy = "cellColData", name = "TSSEnrichment", embedding = "UMAP")
plotPDF(p1, name = "Plot-UMAP-TSS.pdf", ArchRProj = proj, addDOC = FALSE, width = 5, height = 5)

p1 <- plotEmbedding(ArchRProj = proj, colorBy = "cellColData", name = "log10(nFrags)", embedding = "UMAP")
plotPDF(p1, name = "Plot-UMAP-logFrags.pdf", ArchRProj = proj, addDOC = FALSE, width = 5, height = 5)
##################################################################################################

[rcorces]

There is no right answer here. I would try all of the things that you listed. In terms of QC, just make sure that you are inspecting the reads x TSS QC plots and selecting cutoffs that make sense. In my opinion, its ok to use a different cutoff for number of frags for different samples in order to exclude junk (but I would keep the TSS cutoff the same whenever possible).

[RegnerM2015]

Sounds like a plan! Thank you so much for your help.

In the past, I implemented a Gaussian Mixture Model from mclust to identify cellular and non-cellular barcodes in each of these quality control metrics. I think Cusanovich et al. performed a similar approach to identifying high quality barcodes. However, in my experience not every sample will have a biomodal distribution in TSS enrichment or log10(nFrags) so I decided to try hard thresholds this time. I keep brainstorming and provide updates for the ArchR community!

In response to, "make sure that you are inspecting the reads x TSS QC plots and selecting cutoffs that make sense":

Is it correct to assume these QC plots are generated after the creation of Arrow Files but before Doublet removal? These QC plots are very helpful, but they visualize the barcodes after the removal of outliers (according to thresholds). Would it be more reasonable to view these plots using all raw barcodes before QC and determine QC thresholds from there?

[rcorces]

Sometimes it is easy to draw QC cutoffs but often times the default thresholds are not sufficient. For example, in this plot, it is clear that the log10 frags threshold is way too low:

In this example, one might even think that the large orange blob are real cells but they are not:

[RegnerM2015]

Is it correct to assume these QC plots are generated after the creation of Arrow Files but before Doublet removal? These QC plots are very helpful, but they visualize the barcodes after the removal of outliers (according to thresholds). Would it be more reasonable to view these plots using all raw barcodes before QC and determine QC thresholds from there?

[rcorces]

Well, the plots I attached above are actually before removal of any cells. But by default ArchR makes these QC plots during Arrow File creation after removal of cells based on log10 frags. I believe this was done to eliminate the bottom left corner of the plot from dominating the density coloration (as it does in the above plots, though these arent the best examples). You can always create Arrow files with floored QC thresholds and ArchR will make you these plots without filtering any cells. Its a little time consuming but accomplishes what you asked for.

[RegnerM2015]

Sounds good! Thanks for your help and advice.