modelGeneVarByPoisson error selecting bandwidth automatically
1
0
Entering edit mode
@lluis-revilla-sancho
Last seen 11 days ago
European Union

I was following the OSCA book for finding variable genes and while trying to calculate the variable genes with modelGeneVarByPoisson I ran into an error.

I have found only another mention to this function in the support forum but it was not related to this error AFAIK.

Code and error message + traceback and sessionInfo:

> var_poiss <- scran::modelGeneVarByPoisson(gex_filt, 
                                            size.factors = gex_filt$sum,
                                            block = gex_filt$type, subset_row = top_g, 
                                            BPPARAM = multi_core)
|==============================================================| 100%

|==============================================================| 100%

Error in density.default(x, adjust = adjust, from = min(x), to = max(x)) : 
  need at least 2 points to select a bandwidth automatically
In addition: Warning message:
  In seq.default(from = log(min(means)), to = log(max(means)), length = npts) :
  partial argument match of 'length' to 'length.out'
> traceback()
15: stop("need at least 2 points to select a bandwidth automatically")
14: density.default(x, adjust = adjust, from = min(x), to = max(x))
13: density(x, adjust = adjust, from = min(x), to = max(x)) at utils_variance.R#183
12: .inverse_density_weights(m, adjust = 1) at fitTrendVar.R#100
11: fitTrendVar(fm, fv, ...) at utils_variance.R#80
10: .decompose_log_exprs(x.stats$means, x.stats$vars, sim.out$means, 
                         sim.out$vars, x.stats$ncells, ...) at modelGeneVarByPoisson.R#123
9: .model_gene_var_by_poisson(assay(x, i = assay.type), ...) at modelGeneVarByPoisson.R#147
8: .local(x, ...)
7: .nextMethod(x = x, size.factors = size.factors, ...)
6: eval(call, callEnv)
5: eval(call, callEnv)
4: callNextMethod(x = x, size.factors = size.factors, ...) at modelGeneVarByPoisson.R#156
3: .local(x, ...)
2: scran::modelGeneVarByPoisson(gex_filt, size.factors = gex_filt$sum, 
                                block = gex_filt$type, subset_row = top_g, BPPARAM = multi_core) at modelGeneVarByPoisson.R#136
1: scran::modelGeneVarByPoisson(gex_filt, size.factors = gex_filt$sum, 
                                block = gex_filt$type, subset_row = top_g, BPPARAM = multi_core)
> sessionInfo()
R version 4.4.0 (2024-04-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_ES.UTF-8        LC_COLLATE=en_GB.UTF-8    
 [5] LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=es_ES.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Madrid
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices datasets  utils    
[7] methods   base     

other attached packages:
 [1] readxl_1.4.3                rutils_0.0.1.9004          
 [3] PCAtools_2.16.0             ggrepel_0.9.5              
 [5] ensembldb_2.28.0            AnnotationFilter_1.28.0    
 [7] GenomicFeatures_1.56.0      forcats_1.0.0              
 [9] dplyr_1.1.4                 here_1.0.1                 
[11] ggplot2_3.5.1               DelayedMatrixStats_1.26.0  
[13] DelayedArray_0.30.1         SparseArray_1.4.8          
[15] S4Arrays_1.4.1              abind_1.4-5                
[17] Matrix_1.7-0                robustbase_0.99-2          
[19] GO.db_3.19.1                AnnotationDbi_1.66.0       
[21] BiocSingular_1.20.0         patchwork_1.2.0            
[23] scDblFinder_1.18.0          scds_1.20.0                
[25] BiocParallel_1.38.0         scuttle_1.14.0             
[27] biomaRt_2.60.0              AnnotationHub_3.12.0       
[29] BiocFileCache_2.12.0        dbplyr_2.5.0               
[31] DropletUtils_1.24.0         SingleCellExperiment_1.26.0
[33] SummarizedExperiment_1.34.0 Biobase_2.64.0             
[35] GenomicRanges_1.56.0        GenomeInfoDb_1.40.1        
[37] IRanges_2.38.0              S4Vectors_0.42.0           
[39] BiocGenerics_0.50.0         MatrixGenerics_1.16.0      
[41] matrixStats_1.3.0           targets_1.7.0              

loaded via a namespace (and not attached):
  [1] BiocIO_1.14.0            bitops_1.0-7            
  [3] filelock_1.0.3           tibble_3.2.1            
  [5] R.oo_1.26.0              cellranger_1.1.0        
  [7] pROC_1.18.5              XML_3.99-0.16.1         
  [9] lifecycle_1.0.4          httr2_1.0.1             
 [11] edgeR_4.2.0              rprojroot_2.0.4         
 [13] processx_3.8.4           lattice_0.22-6          
 [15] MASS_7.3-60.2            backports_1.5.0         
 [17] magrittr_2.0.3           limma_3.60.2            
 [19] yaml_2.3.8               metapod_1.12.0          
 [21] cowplot_1.1.3            DBI_1.2.3               
 [23] zlibbioc_1.50.0          purrr_1.0.2             
 [25] R.utils_2.12.3           RCurl_1.98-1.14         
 [27] rappdirs_0.3.3           GenomeInfoDbData_1.2.12 
 [29] irlba_2.3.5.1            dqrng_0.4.1             
 [31] codetools_0.2-20         xml2_1.3.6              
 [33] tidyselect_1.2.1         farver_2.1.2            
 [35] UCSC.utils_1.0.0         ScaledMatrix_1.12.0     
 [37] viridis_0.6.5            GenomicAlignments_1.40.0
 [39] jsonlite_1.8.8           BiocNeighbors_1.22.0    
 [41] scater_1.32.0            tools_4.4.0             
 [43] progress_1.2.3           Rcpp_1.0.12             
 [45] glue_1.7.0               gridExtra_2.3           
 [47] xfun_0.44                HDF5Array_1.32.0        
 [49] withr_3.0.0              BiocManager_1.30.23     
 [51] fastmap_1.2.0            rhdf5filters_1.16.0     
 [53] bluster_1.14.0           fansi_1.0.6             
 [55] callr_3.7.6              digest_0.6.35           
 [57] rsvd_1.0.5               secretbase_0.5.0        
 [59] R6_2.5.1                 colorspace_2.1-0        
 [61] scattermore_1.2          RSQLite_2.3.7           
 [63] R.methodsS3_1.8.2        utf8_1.2.4              
 [65] generics_0.1.3           renv_1.0.7.9000         
 [67] data.table_1.15.4        rtracklayer_1.64.0      
 [69] prettyunits_1.2.0        httr_1.4.7              
 [71] pkgconfig_2.0.3          gtable_0.3.5            
 [73] blob_1.2.4               XVector_0.44.0          
 [75] base64url_1.4            ProtGenerics_1.36.0     
 [77] scales_1.3.0             png_0.1-8               
 [79] scran_1.32.0             knitr_1.47              
 [81] rstudioapi_0.16.0        reshape2_1.4.4          
 [83] rjson_0.2.21             curl_5.2.1              
 [85] cachem_1.1.0             rhdf5_2.48.0            
 [87] stringr_1.5.1            BiocVersion_3.19.1      
 [89] parallel_4.4.0           vipor_0.4.7             
 [91] restfulr_0.0.15          pillar_1.9.0            
 [93] grid_4.4.0               vctrs_0.6.5             
 [95] beachmat_2.20.0          cluster_2.1.6           
 [97] beeswarm_0.4.0           cli_3.6.2               
 [99] locfit_1.5-9.9           compiler_4.4.0          
 [ reached getOption("max.print") -- omitted 24 entries ]

Sorry that I don't provide a reproducible example, I'm not sure if it is worth to create it as I'm also checking the variable genes via modelGeneVar.

scran • 562 views
ADD COMMENT
0
Entering edit mode
Peter Hickey ▴ 740
@petehaitch
Last seen 4 weeks ago
WEHI, Melbourne, Australia

Without a reproducible example or access to the data, it's a bit hard to say. I assume the examples work for you (example("modelGeneVarByPoisson", "scran"))?

ADD COMMENT
0
Entering edit mode

Yes, the examples run fine. I only get some warnings (I activated that check in my environment):

Warning messages:
1: In seq.default(from = log(min(means)), to = log(max(means)), length = npts) :
  partial argument match of 'length' to 'length.out'
2: In seq.default(from = log(min(means)), to = log(max(means)), length = npts) :
  partial argument match of 'length' to 'length.out'

Sorry the data is not public and I am not sure if it is worth for me to chase this error down. I hoped that with the error and the traceback it could help to find the part of the code related to the error. I'll check if there is some tweaking in fitTrendVar that could help work around this error.

ADD REPLY
0
Entering edit mode

Afraid I can't tell much from the traceback. Providing dim(gex_filt), summary(gex_filt$sum), table(gex_filt$type), str(top_g) and multi_core might shed a little more light, but may still be in the too-hard basket for me.

ADD REPLY
0
Entering edit mode

Thank you very much for following this. Here is the information about the input values:

dim(gex_filt)
[1]  32298 191382
> summary(gex_filt$sum)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    500    3315    4429    5262    5997   70637 
> table(gex_filt$type)

   CSF   PBMC 
 78012 113370 
> str(top_g)
 chr [1:1876] "ENSG00000255823" "ENSG00000019582" "ENSG00000204287" ...
> multi_core
class: MulticoreParam
  bpisup: FALSE; bpnworkers: 5; bptasks: 2147483647; bpjobname: BPJOB
  bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
  bpRNGseed: ; bptimeout: NA; bpprogressbar: TRUE
  bpexportglobals: TRUE; bpexportvariables: FALSE; bpforceGC: FALSE
  bpfallback: TRUE
  bplogdir: NA
  bpresultdir: NA
  cluster type: FORK
ADD REPLY
0
Entering edit mode

Hmm, I'm afraid nothing jumps out at me.

ADD REPLY

Login before adding your answer.

Traffic: 456 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6