Difficult error message - bug or documentation issue? or problem on my object?

[skilpinen]

Hi,

I wonder is there bug in the Signac code or in data objects, but I get following error:

peaks.to.genes <- LinkPeaks(seurat.object, peak.assay = "peaks", expression.assay = "RNA", n_sample=10, genes.use = G.a[1:10])
Testing 3 genes and 196807 peaks
Error in if (sum(x = links.keep) == 0) { :
missing value where TRUE/FALSE needed
In addition: Warning message:
In cor(x = as.matrix(x = bg.access), y = as.matrix(x = gene.expression), :
the standard deviation is zero

This is really hard to traceback why this error happens? Where exactly it finds Na where logical is needed? That seurat.object is otherwise working fine. And error is not dependent on that gene amount, I just dropped it to very small to quickly confirm repeat of this problem.

[timoast]

Does this still happen if you increase n_sample? It looks like the standard deviation for the background peak set is zero. To compute the z-score we divide by the sd, which would result in a division by zero in this case and give Inf.

I wouldn't recommend running this function with a very low n_sample value, I'd suggest 50 at a minimum.

[skilpinen]

Yes, it does happen. I just dropped n_sample to run it quick. However, I have here a seurat object with which this this function works and another with which it doesn't. Difference is that the one that does not work is subset of the one that works.

Works:
peaks.to.genes <- LinkPeaks(seurat.object, peak.assay = "peaks", expression.assay = "RNA", n_sample=10, genes.use = G.a[1:10])

Doesn't work:
P3GABA<-subset(seurat.object, idents=c("P3 GABAergic neurons","P3 GABAergic precursors","GABA/Ach", "FB GABA"))
peaks.to.genes <- LinkPeaks(P3GABA, peak.assay = "peaks", expression.assay = "RNA", n_sample=10, genes.use = G.a[1:10])

[timoast]

How many cells are in the subset data?

[skilpinen]

869 in subset, 3577 in original

[timoast]

If you're able to share the subset object that recreates this error I can look into it further. You can email a download link to the object to tstuart@nygenome.org

[skilpinen]

Yes, done.

[skilpinen]

You can also directly email to Sami.Kilpinen@helsinki.fi

[timoast]

This error was due to incorrect feature metadata being stored, which could happen if you run FindTopFeatures and then subset the object. This would cause the total counts for each peak to change, since cells are removed. Since LinkPeaks decides which peaks to use based on total accessibility, this cause problems where it was choosing peaks that actually had zero total counts in the dataset. I've pushed a fix to the develop branch that solves this issue by re-computing the feature metadata when an object is subset. You should be able to run:

devtools::install_github("timoast/signac", ref = "develop")
library(Signac)

obj <- FindTopFeatures(obj, min.cutoff = NA) # re-compute feature metadata
obj <- LinkPeaks(obj, ...)

[skilpinen]

Yes, I confirm that this solved the issue. Thanks.