Error in .local(x, ...) : size factors should be positive #97
Description
Hi,
Thanks for a cool package!
I am running into a problem with size factors when I try to run velociraptor::scvelo(). Let me describe my situation: I have three datasets (seurat objects) that have been constructed from inex (intron and exon counts) and processed through the same (fairly standard seurat) pipeline. The inex counts come from the .rds output from the zUMIs pipeline (smart-seq3 data).
I then read in the exon and intron matrices from the corresponding .rds files and append to the seurat object.
I then convert the seurat object to SingleCellExperiment like this:
#as.SingleCellExperiment() is not treating the assays in a satisfactory way so I do
#it manually like this instead:
#Extract exon and intron counts directly from the Seurat assays
exo_mat <- GetAssayData(seurat_obj[["exon"]], layer = "counts")
int_mat <- GetAssayData(seurat_obj[["intron"]], layer = "counts")
#sct_mat <- GetAssayData(seurat_obj[["SCT"]], layer = "data")
cnt_mat <- GetAssayData(seurat_obj[["RNA"]], layer = "counts")
#Restrict all to shared genes
shared_genes <- Reduce(intersect, list(
rownames(exo_mat),
rownames(int_mat),
#rownames(sct_mat),
rownames(cnt_mat)
))
exo_mat <- exo_mat[shared_genes, ]
int_mat <- int_mat[shared_genes, ]
#sct_mat <- sct_mat[shared_genes, ]
cnt_mat <- cnt_mat[shared_genes, ]
#Optional: Add metadata and dimensionality reductions from Seurat
cell_metadata <- seurat_obj@meta.data
umap_coords <- Embeddings(seurat_obj, "umap")
pca_coords <- Embeddings(seurat_obj, "pca")
#Create SCE with spliced/unspliced and metadata
SCE_obj <- SingleCellExperiment(
assays = list(spliced = exo_mat,
unspliced = int_mat,
#SCT = sct_mat
counts = cnt_mat
),
colData = cell_metadata,
reducedDims = SimpleList(X_umap = umap_coords,
X_pca = pca_coords)
)
I then run scvelo like this:
#Extract the precomputed variable genes
top.hvgs <- VariableFeatures(seurat_obj)
top.hvgs <- top.hvgs[top.hvgs %in% rownames(SCE_obj)]
# Log normalize
SCE_obj <- scuttle::logNormCounts(SCE_obj, assay.type = "counts")
# Run scVelo
velo.out <- velociraptor::scvelo(
SCE_obj,
subset.row = top.hvgs,
sf.X = librarySizeFactors(SCE_obj, assay.type = "logcounts"),
assay.X = "logcounts",
assay.spliced = "spliced",
assay.unspliced = "unspliced",
mode = "dynamical",
scvelo.params = list(# set n cores (mode = "dynamical can't parallelize)
velocity_graph = list(n_jobs = 16L)))
For two out of three of the datasets this seems to work. However for one of them I get this error message:
> velo.out <- velociraptor::scvelo(
+ SCE_obj,
+ subset.row = top.hvgs,
+ sf.X = librarySizeFactors(SCE_obj, assay.type = "logcounts"),
+ assay.X = "logcounts",
+ assay.spliced = "spliced",
+ assay.unspliced = "unspliced",
+ mode = "dynamical",
+ scvelo.params = list(# set n cores (mode = "dynamical can't parallelize)
+ velocity_graph = list(n_jobs = 16L)))
Error in .local(x, ...) : size factors should be positive
Can you spot what I am doing wrong?
Here is some info about the SCE_obj:
> librarySizeFactors(SCE_obj) %>% summary()
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.001632 0.793102 0.990877 1.000000 1.200908 3.042436
> SCE_obj
class: SingleCellExperiment
dim: 13039 1261
metadata(0):
assays(4): spliced unspliced counts logcounts
rownames(13039): A2ML1 AAAS ... ZZEF1 ZZZ3
rowData names(0):
colnames(1261): P30915_105_AACCATCGGCCAAGTGTAGC P30915_105_AACCATCGGCCCAATTAGCA ... P32812_127_TTGGCCACGATTGCCTTATC
P32812_127_TTGGCCACGATTGCGACCAC
colData names(30): orig.ident nCount_RNA ... nFeature_intron sizeFactor
reducedDimNames(2): X_umap X_pca
mainExpName: NULL
altExpNames(0):
I am verry grateful for any help I can get!
Best,
Robin