Problem generating observational weights (zinbwave_1.8.0)

[bgrone]

I have a problem in which I am trying to use the zinbwave package for differential expression analysis from snRNAseq data, but I cannot get the necessary observational weights.

I am using: R version 3.6.1 zinbwave_1.8.0

I created a SummarizedExperiment object by using data and metadata from a Seurat object:

se <- SummarizedExperiment(seurat@assays$RNA@counts, colData = seurat@meta.data)

After filtering the se_object by gene counts, I use the 'zinbwave' command to generate a SingleCellObject

singlecell <- zinbwave(se, K=2, X="~genotype", epsilon=1000)

The SingleCellExperiment object ('singlecell') is generated, but after creating the object there is no "weights" data in the object.

Trying to access weights data:

weights <- assay(singlecell, "weights")

results in this error message:

Error in assay(singlecell, "weights") : 'assay(, i="character", ...)' invalid subscript 'i' 'weights' not in names(assays())

The output of 'assayNames(singlecell)' is: NULL

The output of 'singlecell@assays' is: An object of class "SimpleAssays" Slot "data": List of length 1

In the documentation, it says that "Since version 1.1.5, zinbwave computes the observational weights by default." However, I also tried specifying "observationalWeights = TRUE" and this did not appear to lead to any weights being generated, either.

Why are the weights not generated by the zinbwave function? Could it be because of a problem with how I set up the SummarizedExperiment? Or is it more likely related to the way I have used the zinbwave command?

[drisso]

Hi @bgrone

I'm not sure why I say otherwise, but the default value is actually observationalWeights = FALSE as you can see in the man page of the zinbwave function. Where in the documentation did you see that its default value is TRUE? I should correct it.

Hence, the default behavior that you observe is expected. What is not expected of course is that the weights are not computed when you specify observationalWeights = TRUE. When I try this myself everything was as expected, see below for a reproducible example.

library(zinbwave)
library(scRNAseq)

se = ReprocessedFluidigmData(assays = "tophat_counts")
se = as(se, "SummarizedExperiment")

filter <- rowSums(assay(se[,1:30]))>100
se = se[filter,1:30]

sce = zinbwave(se[1:100,], K=2, epsilon=1000)
sce

In this case no weights assay:

class: SingleCellExperiment 
dim: 100 30 
metadata(3): sample_info clusters which_qc
assays(1): tophat_counts
rownames: NULL
rowData names(0):
colnames(30): SRR1275356 SRR1274090 ... SRR1275363 SRR1275280
colData names(28): NREADS NALIGNED ... Cluster1 Cluster2
reducedDimNames(1): zinbwave
spikeNames(0):
altExpNames(0):

But when I specify that I want the weights:

sce = zinbwave(se[1:100,], K=2, epsilon=1000, observationalWeights = TRUE)
sce

class: SingleCellExperiment 
dim: 100 30 
metadata(3): sample_info clusters which_qc
assays(2): tophat_counts weights
rownames: NULL
rowData names(0):
colnames(30): SRR1275356 SRR1274090 ... SRR1275363 SRR1275280
colData names(28): NREADS NALIGNED ... Cluster1 Cluster2
reducedDimNames(1): zinbwave
spikeNames(0):
altExpNames(0):

Would you be able to create a reproducible example that reproduces the behavior that you observed?

[bgrone]

Hi Davide,

Thank you very much for your reply. I found that statement about the observational weight default: "Since version 1.1.5, zinbwave computes the observational weights by default. See the man page of zinbwave."

at the zinbwave vignette, in section 8: https://bioconductor.org/packages/release/bioc/vignettes/zinbwave/inst/doc/intro.html#differential-expression

Using your reproducible example, I was able to get the same result that you showed.

I then tried to make a reproducible example using a Seurat object as the starting point. The "pbmc3k_final.rds" Seurat object is available from this vignette: https://satijalab.org/seurat/v3.0/interaction_vignette.html

library(zinbwave)
library(scRNAseq)

pbmc <- readRDS(file = "pbmc3k_final.rds")

#convert Seurat object to SummarizedExperiment
pbmc_se <- SummarizedExperiment(pbmc@assays$RNA@counts,
                           colData = pbmc@meta.data)

#removing those genes that do not have more than 5 reads in more than 5 samples
filter <- rowSums(assay(pbmc_se)>5)>5
table(filter)

pbmc_se  <- pbmc_se [filter,]

#zinbwave on subset of cells
pbmc_se_singlecell <- zinbwave(pbmc_se[1:100,], K=2, epsilon=1000, observationalWeights = TRUE)

pbmc_se_singlecell

These commands led to errors:

pbmc_se_singlecell <- zinbwave(pbmc_se[1:100,], K=2, epsilon=1000, observationalWeights = TRUE) Error in t.default(x) : argument is not a matrix In addition: Warning message: In .local(Y, ...) : No assay named counts, using first assay.Use assay to specify a different assay.

pbmc_se_singlecell Error: object 'pbmc_se_singlecell' not found

So, because your example works but mine does not, it seems likely that there is some problem with how the SummarizedExperiment object was generated from the Seurat object, rather than with the way I am using the zinbwave command. Is that correct?

Is there a proper way to extract a SummarizedExperiment object from Seurat so that it will work with zinbwave?

[bgrone]

I also tried to generate the SummarizedExperiment object a different way, but still had problems.

#convert Seurat object to SummarizedExperiment
pbmc_se <- as.SingleCellExperiment(pbmc)
pbmc_se <- as(pbmc_se, "SummarizedExperiment")

which led to errors:

> pbmc_se_singlecell <- zinbwave(pbmc_se[1:100,], K=2, epsilon=1000, observationalWeights = TRUE)
Error in t.default(x) : argument is not a matrix
> pbmc_se_singlecell
Error: object 'pbmc_se_singlecell' not found

There may be a problem because the Seurat object has a formal class 'dgCMatrix', instead of a normal matrix?

[vincr04]

Maybe you could try

pbmc_se <- SummarizedExperiment(as.matrix(pbmc@assays$RNA@counts),
                           colData = pbmc@meta.data)

I had the same issue and it has solved it for me...

[bgrone]

Yes, vincr04, this appears to be a solution! I tried it with my dataset and zinbwave was able to run and generate the observational weights. Thank you.