projectBulkATAC not working with Harmony?

[cschmidl]

Describe the bug
Dear Ryan and Jeff, I was running into a problem with the projectBulkATAC function (which I find really useful!!).

I use the recent version from githuub and commented out
.loadUWOT (line 168) as it gave an error Error ('unused argument (embedding$params$nc)')
Also, we added the force option.

The function runs perfect with the IterativeLSI dimredduction in the testdataset. However, when I run the function with the testdataset on the Object after Harmony and the corresponding UMAP I get the error:

Error in normalizeDoubleBracketSubscript(i, x, exact = exact, allow.NA = TRUE, :
attempt to extract less than one element
9.
stop("attempt to extract less than one element")
8.
normalizeDoubleBracketSubscript(i, x, exact = exact, allow.NA = TRUE,
allow.nomatch = TRUE)
7.
getListElement(x@listData, i, exact = exact)
6.
getListElement(x@listData, i, exact = exact)
5.
getListElement(x, i, ...)
4.
getListElement(x, i, ...)
3.
rD[[grep("Features", names(rD))]]
2.
rD[[grep("Features", names(rD))]]
1.
projectBulkATAC(ArchRProj = proj2, seATAC = rse_list$GSE179593,
reducedDims = "Harmony", embedding = "HarmonyUMAP", n = 100,
verbose = TRUE, threads = 1, logFile = createLogFile("projectBulkATAC"))

I attached the logfile!

Thank you very much!
ArchR-projectBulkATAC-a4b71339b1a8-Date-2022-04-07_Time-16-06-36.log

[rcorces]

Hi @cschmidl! Thanks for using ArchR! Please make sure that your post belongs in the Issues section. Only bugs and error reports belong in the Issues section. Usage questions and feature requests should be posted in the Discussions section, not in Issues.

Before we help you, you must respond to the following questions unless your original post already contained this information:
1. If you've encountered an error, have you already searched previous Issues to make sure that this hasn't already been solved?
2. Can you recapitulate your error using the tutorial code and dataset? If so, provide a reproducible example.
3. Did you post your log file? If not, add it now.
4. Remove any screenshots that contain text and instead copy and paste the text using markdown's codeblock syntax (three consecutive backticks). You can do this by editing your original post.

[rcorces]

theres some active development on .loadUWOT that I'm not up to date on (6e7abec). But I imagine this is related to your problem. Can you try updating to the latest release_1.0.2 and trying again.

For starters, it isnt surprising that commenting out loadUWOT results in downstream problems. So the error you are reporting is not what needs to be addressed per se. The function should work without commenting out code.

[cschmidl]

Thank you for the fast response, using simply the 1.0.2 version yields this error now even with the LSI dimreduction from the testdataset:
ArchR logging to : ArchRLogs/ArchR-projectBulkATAC-c01820d594a8-Date-2022-04-07_Time-18-06-22.log
If there is an issue, please report to github with logFile!
2022-04-07 18:06:24 : Overlap Ratio of Reduced Dims Features = 0.71432
Error in rep(xi, length.out = nvar) :
attempt to replicate an object of type 'S4'
ArchR-projectBulkATAC-c01820d594a8-Date-2022-04-07_Time-18-06-22.log

[rcorces]

can you provide the output of traceback()

[cschmidl]

4: rbind(deparse.level, ...)
3: rbind(mat, matNotIn)
2: .safeSubset(mat = .getAssay(subATAC, "counts"), subsetRows = paste0("f",
seq_along(rDGR)))
1: projectBulkATAC(ArchRProj = proj, seATAC = rse_list$GSE179593,
reducedDims = "IterativeLSI", embedding = "UMAP", n = 100,
verbose = TRUE, threads = 1, logFile = createLogFile("projectBulkATAC"))

[rcorces]

@cschmidl could you try this new branch dev_projectBulkATAC which includes this change 3a07120

[cschmidl]

2022-04-08 07:51:20 : Overlap Ratio of Reduced Dims Features = 0.71432
Error in dimnames<-.data.frame(*tmp*, value = list(n)) :
invalid 'dimnames' given for data frame

traceback()
11: stop("invalid 'dimnames' given for data frame")
10: dimnames<-.data.frame(*tmp*, value = list(n))
9: dimnames<-(*tmp*, value = list(n))
8: as.array.default(.getAssay(subATAC, "counts"))
7: as.array(.getAssay(subATAC, "counts"))
6: rownames(mat)
5: match(x, table, nomatch = 0)
4: subsetRows %ni% rownames(mat)
3: which(subsetRows %ni% rownames(mat))
2: .safeSubset(mat = as.array(.getAssay(subATAC, "counts")), subsetRows = paste0("f",
seq_along(rDGR)))
1: projectBulkATAC(ArchRProj = proj, seATAC = rse_list$GSE179593,
reducedDims = "IterativeLSI", embedding = "UMAP", n = 100,
verbose = TRUE, threads = 1, force = TRUE, logFile = createLogFile("projectBulkATAC"))
ArchR-projectBulkATAC-108802d6bed81-Date-2022-04-08_Time-07-56-07.log

Sorry I was also wring about the commenting out loadUWOT, it was just the second argument that gave an error
'unused argument (embedding$params$nc)', so this soved it:
umapModel <- .loadUWOT(embedding$params$uwotModel) #, embedding$params$nc)

I attach the function that is working on my computer (a colleague fixed it and I think it was mainly removing embedding$params$nc and enabling force=TRUE, but the could not fix the issue with projectBulkATAC & Harmony, and also not the other issue I opened with my real data that reports inconsistencies with the reduces dimensions also in the LSI domreduction. Thanks again for your effort
projectBulkATAC.txt
!

[rcorces]

I'm having some trouble following the issue.

A reproducible example would be useful. You mentioned that this is failing with the tutorial data set so perhaps all you would need to provide would be the object being passed to seATAC.

[cschmidl]

I'm sorry for the confusion, here is my se ("rse_list"):
https://myfiles.uni-regensburg.de/filr/public-link/file-download/044787a07ff55ccb01800a7226bf2c1e/80511/-5564395770685507458/giles_2022_GSE179593_GSE179606_ranged_se_list2.RDS

rse_list$GSE179593 is he one to run:

res=projectBulkATAC(
ArchRProj = proj,
seATAC = rse_list$GSE179593,
reducedDims = "IterativeLSI",
embedding = "UMAP",
n = 100,
verbose = TRUE,
threads = 1,
force = TRUE,
logFile = createLogFile("projectBulkATAC")
)

Thank you!

[rcorces]

Long story short, you cannot currently project into a Harmony reduced dimension space.

I'm not sure the details because I'm not that familiar with the projection code but I made a change to the .safeSubset() function on a new branch called dev_safeSubset that should solve your problem
05c1378

I think there is some incompatibility between rbinding a matrix and a sparse matrix and if I had to guess projectBulkATAC() seems to want a sparse matrix as the counts matrix in seATAC. Either way, the dev_safeSubsetbranch should work for you withIterativeLSI`

[cschmidl]

Thank you for your answer and updated code!

It runs now through the samples of the se, but in the end yields an error with my data, not with the test data (logfile attached):
...
2022-04-13 11:59:06 : Projecting Sample (217 of 217), 2.068 mins elapsed.
2022-04-13 11:59:06 : Error incosistency found with matching LSI dimensions to those used in addEmbedding
Returning with simulated reduced dimension coordinates...

In addition, will you implement project into a Hermony reduced dimenstion space? I think that would be super useful not only for me!

Thanks a lot!

ArchR-projectBulkATAC-1521abc9a3e-Date-2022-04-13_Time-11-57-02.log

[rcorces]

I think there is something wrong with the logic in this section of the code but I'm really not sure.

In this section, there are some nested if statements, two of which check for the same exact thing which doesnt make sense to me (L139 and L155).

ArchR/R/BulkProjection.R

Lines 139 to 166 in 3af7c30

	if(embedding$params$nc != ncol(simRD)){
	if(is.null(dimsToUse)){
	dimsToUse <- seq_len(ncol(rD[[1]]))
	}
	if(!is.null(corCutOff)){
	if(scaleDims){
	corToDepth <- rD$corToDepth$scaled
	dimsToUse <- dimsToUse[corToDepth < corCutOff]
	}else{
	corToDepth <- rD$corToDepth$none
	dimsToUse <- dimsToUse[corToDepth < corCutOff]
	}
	}
	if(embedding$params$nc != ncol(simRD)){
	.logMessage("Warning! Inconsistency found with matching LSI dimensions to those used in addUMAP or addTSNE",
	"\nReturning with simulated reduced dimension coordinates...", verbose = TRUE, logFile = logFile)
	out <- SimpleList(
	simulatedReducedDims = simRD
	)
	return(out)
	}
	simRD <- simRD[, dimsToUse, drop = FALSE]
	}

It seems like this if statement is checking whether the number of dimensions used in the original addUMAP() call match the simulated data but this would only be true if you did not_ use the dimsToUse param when originally adding your embedding. If you did happen to subset the dimensions to use, then it seems like this code would always fail. Maybe I'm misinterpreting this though.

Could you confirm the values of corCutOff, dimsToUse, and scaleDims when you run:

  embedding <- getEmbedding(ArchRProj = ArchRProj, embedding = embedding, returnDF = FALSE)
  corCutOff <- embedding$params$corCutOff
  dimsToUse <- embedding$params$dimsToUse
  scaleDims <- embedding$params$scaleDims