Using addIterativeLSI with a PeakMatrix

[skwalker]

Hi!

Thanks for building such a great platform for analyzing scATAC-seq data ! I've been playing around with using different values for ceiling when running addPeakMatrix and seeing how they change the lower dimensionality space by running addIterativeLSI with useMatrix = "PeakMatrix". I've found that changing the value for ceiling, while it certainly changes the values in the Peak Matrix, does not seem to affect the lower dimensionality space; for example, my UMAP plots are exactly the same when I use ceiling=1 and ceiling=30, and in addition, the SVD matrix generated from the reduced dimensions is also unchanged. It seems a little bizarre to me that this wouldn't change anything since the Peak Matrix is quite different in these two situations. Do you have any idea of why updating the PeakMatrix like this wouldn't actually change the lower dimensionality space ? It is definitely using the PeakMatrix, right?

Thanks,
Sarah

[rcorces]

I dont see a clear reason for how this would happen. Could you show some of the code that you have used?

[skwalker]

I think perhaps the issue might have something to do with the size of the peakset (my original one had 52k peaks)?

I was able to reiterate my original issue on the tutorial dataset by using a subset of the peak set extracted from the tutorial dataset (with sample(peakset, 25000)). Then, even though the peak matrices are quite different when using a ceiling of 1 and 16, the lower dimensionality spaces are exactly the same.

[rcorces]

Can you provide the exact code to reproduce this from the tutorial data?

[skwalker]

Ah I'm sorry for the delay -- I thought I had already sent this message. Here is the code I used to show this:

library(ArchR)
set.seed(1)
addArchRThreads(threads = 1) 
inputFiles =   getTutorialData("Hematopoiesis")
addArchRGenome("hg19")

arrowfiles <- createArrowFiles(
  inputFiles = inputFiles,
  sampleNames = names(inputFiles),
  filterTSS = 4, #Dont set this too high because you can always increase later
  filterFrags = 1000, 
  addTileMat = TRUE,
  addGeneScoreMat = FALSE
)

proj <- ArchRProject(ArrowFiles = arrowfiles,
                     outputDirectory="archr_test_dir",
                     copyArrows=TRUE,)

# for adding clusters based on tiles
proj <- addIterativeLSI(
    ArchRProj = proj,
    useMatrix = "TileMatrix", 
    name = "IterativeLSI", 
    iterations = 2, 
    clusterParams = list( #See Seurat::FindClusters
        resolution = c(0.2), 
        sampleCells = 10000, 
        n.start = 10
    ), 
    varFeatures = 25000, 
    dimsToUse = 1:30
)

proj <- addClusters(
    input = proj,
    reducedDims = "IterativeLSI",
    method = "Seurat",
    name = "Clusters",
    resolution = 0.8
)

proj <- addGroupCoverages(ArchRProj = proj, groupBy = "Clusters")

# create peak set
pathToMacs2 <- findMacs2()
proj <- addReproduciblePeakSet(
    ArchRProj = proj, 
    groupBy = "Clusters",
    pathToMacs2 = pathToMacs2
)

proj = addPeakMatrix(proj, ceiling = 4)
max(assay(getMatrixFromProject(proj, useMatrix="PeakMatrix"))) # returns 4 

new_proj <- saveArchRProject(ArchRProj = proj, outputDirectory = "archr_test_dir_b", load = TRUE)
peakset <- getPeakSet(proj)

# create subsetted peak set
subsampled_peakset <- sample(peakset, 25000)
mcols(subsampled_peakset) <- NULL
subsampled_peakset #returns basic peak set with no meta data as expected

new_proj <- addPeakSet(ArchRProj = new_proj, peakSet = subsampled_peakset, force = TRUE)
getPeakSet(new_proj) #returns basic peak set with some meta data

# dimensionality reduction with ceiling of 1 for peak matrix
new_proj <- addPeakMatrix(new_proj, ceiling = 1)
max(assay(getMatrixFromProject(new_proj, useMatrix="PeakMatrix"))) ## returns 1, as expected

new_proj <- addIterativeLSI(
  ArchRProj = new_proj,
  useMatrix = "PeakMatrix", 
  name = "IterativeLSICeiling1", 
  iterations = 4, 
  clusterParams = list( 
    resolution = c(1.5), 
    sampleCells = 10000, 
    n.start = 10,
    maxClusters = 25
  ), 
  varFeatures = 52004, 
  dimsToUse = 1:30,
  force = TRUE,
  LSIMethod = 2,
)
ceiling1LSI <- new_proj@reducedDims$IterativeLSICeiling1$matSVD


# dimensionality reduction with ceiling of 13 for peak matrix
new_proj <- addPeakMatrix(new_proj, ceiling=13)
max(assay(getMatrixFromProject(new_proj, useMatrix="PeakMatrix"))) ## returns 13, as expected

new_proj <- addIterativeLSI(
  ArchRProj = new_proj,
  useMatrix = "PeakMatrix", 
  name = "IterativeLSICeiling13", 
  iterations = 4, 
  clusterParams = list( 
    resolution = c(1.5), 
    sampleCells = 10000, 
    n.start = 10,
    maxClusters = 25
  ), 
  varFeatures = 52004, 
  dimsToUse = 1:30,
  force = TRUE,
  LSIMethod = 2,
)
ceiling13LSI <- new_proj@reducedDims$IterativeLSICeiling13$matSVD

all.equal(ceiling1LSI, ceiling13LSI) # returns TRUE, not as expected

[rcorces]

I can confirm that I am able to recapitulate this behavior using the code you provided. I am also able to recapitulate this behavior when creating two new projects from a blank-slate tutorial project that lacks any PeakMatrix, meaning that I dont think this behavior has anything to do with project subsetting, weird issues with saveArchRProject, or issues with overwriting a peak matrix etc.

That being said, I havent yet figured out what the problem is because the actual PeakMatrix objects are indeed different when using different values for ceiling. I'll keep looking into this.

[rcorces]

I'm moving this to an Issue for easier tracking. Please follow #941 and post there for followup