Error in FindClustersRecall: missing value required for `stats::runif < prop.zero`

[wong-nw]

When working with a previous object containing RNA, SCT, and integrated assays, the object is set up to analyze clusters on the integrated assay. Running FindClustersRecall produces the following output:

so_recall = FindClustersRecall(so)
Pulling data from Seurat object
Computing artificial features
Error in if (stats::runif(1) < prop.zero) { : 
  missing value where TRUE/FALSE needed

Do you know if this is an internal bug or an issue with the data set itself?

[agdenadel]

@wong-nw are you able to share your data with me (even privately)? I am a bit busy because I am defending my PhD on Thursday, but after that I would be happy to look into this. I worked directly with RNA assays during development so I am not too surprised that something like this happened. It's most likely a bug and not an issue with the data.

Also, recall doesn't yet have any built-in support for integrated data so I'm not sure how the integrated result will turn out. Which integration method did you use?

[agdenadel]

@wong-nw circling back to this, I am happy to look into this for you. Can you help me reproduce this error?

[wong-nw]

Hi @agdenadel, sorry about the delay. This particular dataset used the "integrate" function from Seurat, which, if I recall correctly, is CCA. I don't have permission to share this data, as it is still unpublished, but I can try to use some public datasets to reproduce the error.

From what I can interpret, it looks like the error starts with the insertion of the pseudodata (my terminology is probably wrong), which would make sense, since the integrated assay only contains the variable genes (~3000).

I'll reach out with a test dataset/code soon.

Good luck on your dissertation defense! Been there, done that, and eventually got the degree in the end

[agdenadel]

@wong-nw if you can share the code that you are using to integrate the data (even with a column of artificial labels on the pbmc3k data) I will happily look into this. If you use the SeuratData package it will be as easy as:

InstallData("pbmc3k")
data("pbmc3k")
pbmc3k

I am reasonably confident that any integrated dataset will reproduce this error since the method was only tested using RNA assays when I was working on the publication. I can also try out the seurat integration myself in a few days.

[wong-nw]

@agdenadel I'm currently testing with the latest version of Seurat, which does not create an integrated assay. It actually raises a new issue, in that the preserved data is restricted to the reductions slot of the Seurat object. Can we add a parameter to FindClustersRecall that allows a different reduction to be used? This would likely solve this issue.

As an additional question, I was testing with the SCTransform assay, and it consistently ends after one iteration. Do you have any insight?

Code is attached

#Data and processing copied from https://satijalab.org/seurat/articles/integration_introduction

library(Seurat)
library(SeuratData)
library(recall)
library(patchwork)
library(presto)

InstallData("ifnb")

# load dataset
ifnb_orig <- LoadData("ifnb")
# split the RNA measurements into two layers one for control cells, one for stimulated cells

ifnb_orig[["RNA"]] <- split(ifnb_orig[["RNA"]], f = ifnb_orig$stim)
ifnb_orig

# Initial analysis workflow with log normalization
# run standard analysis
ifnb <- NormalizeData(ifnb_orig)
ifnb <- FindVariableFeatures(ifnb)
ifnb <- ScaleData(ifnb)
ifnb <- RunPCA(ifnb)
ifnb <- FindNeighbors(ifnb, dims = 1:30, reduction = "pca")
ifnb <- RunUMAP(ifnb, dims = 1:30, reduction = "pca")

#Integration on log-normalized data
ifnb_int <- IntegrateLayers(object = ifnb, method = CCAIntegration, orig.reduction = "pca", new.reduction = "integrated.cca",
                        verbose = TRUE)

# re-join layers after integration
ifnb_int[["RNA"]] <- JoinLayers(ifnb_int[["RNA"]])
ifnb_int <- FindNeighbors(ifnb_int, reduction = "integrated.cca", dims = 1:30)
ifnb_int <- RunUMAP(ifnb_int,reduction="integrated.cca",dims=1:30)

#Analysis workflow with SCTransform normalization
ifnb_sct <- SCTransform(ifnb_orig)
ifnb_sct <- RunPCA(ifnb_sct)
ifnb_sct <- FindNeighbors(ifnb_sct, dims = 1:30, reduction = "pca")
ifnb_sct <- RunUMAP(ifnb_sct,dims=1:30)

#integration with SCTransform normalized data
ifnb_int_sct <- IntegrateLayers(object = ifnb_sct, method = CCAIntegration, normalization.method = "SCT", verbose = T,new.reduction="integrated.cca")
ifnb_int_sct <- FindNeighbors(ifnb_int_sct, reduction = "integrated.cca", dims = 1:30)
ifnb_int_sct = RunUMAP(ifnb_int_sct,dims=1:30,reduction="integrated.cca")


#Run "recall" on uncorrected lognorm data
ifnb[["RNA"]] = JoinLayers(ifnb[["RNA"]]) #JoinLayers reserved for this step to allow other processing steps
ifnb = FindClustersRecall(ifnb,algorithm = "louvain",reduction_percentage = 0.1,dims=1:30,assay="RNA")
ifnb$recall_clusters = Idents(ifnb)
DimPlot(ifnb,group.by=c("stim","recall_clusters"))
        
#Run "recall" on batch corrected lognorm data
ifnb_int = FindClustersRecall(ifnb_int,algorithm="louvain",reduction_percentage=0.1, dims=1:30,assay="RNA")
ifnb_int$recall_clusters = Idents(ifnb_int)
DimPlot(ifnb_int,group.by=c("stim","recall_clusters"))

#Run "recall" on uncorrected SCTransform data
ifnb_sct = FindClustersRecall(ifnb_sct,algorithm = "louvain",reduction_percentage = 0.1,dims=1:30,assay="SCT")
ifnb_sct$recall_clusters = Idents(ifnb_sct)
DimPlot(ifnb,group.by=c("stim","recall_clusters"))

#Run "recall" on batch corrected SCTransform data
ifnb_int_sct = FindClustersRecall(ifnb_int_sct,algorithm="louvain",reduction_percentage=0.1, dims=1:30,assay="SCT")
ifnb_int_sct$recall_clusters = Idents(ifnb_int_sct)
DimPlot(ifnb_int_sct,group.by=c("stim","recall_clusters"))

[agdenadel]

@wong-nw I don't have a fix yet but I want to let you know that I am looking into this today.

[agdenadel]

@wong-nw I was able to get multiple iterations to occur on the SCT data by using a higher resolution, but I am hesitant to encourage the use of recall with SCTransform normalized (and also integrated data) since the method has only been validated using raw counts.

The code you attached mostly ran, but the initial integration step failed (not sure if that is a Seurat versioning issue or something). I will take a look at the initial issue that you reported as well.