search

HectorRDB / condiments

Trajectory inference across multiple conditions with condiments: differential topology, progression, differentiation, and expression
https://hectorrdb.github.io/condiments/
Other
24 stars 3 forks source link

Error with conditionTest #24

Closed jordiabante closed 10 months ago

jordiabante commented 1 year ago

Hi Hector and Nitesh,

Thanks for this package, it is really useful.

I am interested in running the different differential analyses available in condiments between two conditions. So far, I have been able to run the following ones:

I have also been able to fit a GAM model, regardless of the condition using (had to reduce the number of cells and genes though):

# select subset of cells for testing purposes (for testing purposes)
  sel_cells <- split(colnames(so@assays$RNA@data), so$cell_annot)
  sel_cells <- unlist(lapply(sel_cells, function(x) {
    return(sample(x, 50))
  }))

  # select subset of 500 most variable genes (for testing purposes)
  gv <- as.data.frame(na.omit(scran::modelGeneVar(so@assays$RNA@data[, sel_cells])))
  gv <- gv[order(gv$bio, decreasing = T), ]
  sel_genes <- sort(rownames(gv)[1:500])

  # design matrix
  des.mat <- model.matrix(
    ~ stage + genotype + stage:genotype,
    data=so[sel_genes,sel_cells]@meta.data[,c("stage","genotype")]
  )

  # fit GAM (for nknots see https://github.com/statOmics/tradeSeq/issues/121)
  sceGam <- fitGAM(
    counts = drop0(so@assays$RNA@data[sel_genes, sel_cells]),
    pseudotime = pseudotime[sel_cells,],
    cellWeights = cellWeights[sel_cells, ],
    U = des.mat, nknots = 6, verbose = T, parallel = T, sce = TRUE,
    BPPARAM = BiocParallel::MulticoreParam())

As well as a GAM model taking into account the condition (genotype):

# select subset of cells for testing purposes (for testing purposes)
  sel_cells <- split(colnames(so@assays$RNA@data), so$cell_annot)
  sel_cells <- unlist(lapply(sel_cells, function(x) {
    return(sample(x, 50))
  }))

  # select subset of N most variable genes (for testing purposes)
  gv <- as.data.frame(na.omit(scran::modelGeneVar(so@assays$RNA@data[, sel_cells])))
  gv <- gv[order(gv$bio, decreasing = T), ]
  sel_genes <- sort(rownames(gv)[1:500])

  # design matrix
  des.mat <- model.matrix(
    ~ stage + genotype + stage:genotype,
    data=so[sel_genes,sel_cells]@meta.data[,c("stage","genotype")]
  )

  # fit GAM (for nknots see https://github.com/statOmics/tradeSeq/issues/121)
  sceGamGt <- fitGAM(
    counts = drop0(so@assays$RNA@data[sel_genes, sel_cells]),
    pseudotime = pseudotime[sel_cells,],
    cellWeights = cellWeights[sel_cells,],
    conditions = factor(so@meta.data[sel_cells,]$genotype),
    nknots = 6, U=des.mat, verbose = T, parallel = T, sce=TRUE,
    BPPARAM = BiocParallel::MulticoreParam()
  )

However, when I try to extract the results of the latter (per lineage) I get the following error:

> condRes <- conditionTest(sceGamGt, global = FALSE, lineages=TRUE, pairwise=FALSE)
Error in if (r == 1) return(c(NA, NA)) : 
  missing value where TRUE/FALSE needed

If I simply run condRes <- conditionTest(sceGamGt) then it seems to work just fine, whereas any combination with lineages=TRUE seems to fail. However, we are very much interested in the individual results per lineage.

Any help would be greatly appreciate it. I would be happy to share the sceGamGtobject, which I have stored in an .rds file, if that can be of any help.

Thanks! jordi

> sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.5 (Green Obsidian)

Matrix products: default
BLAS:   /gpfs/apps/AMD/R/4.2.2/GCC/lib64/R/lib/libRblas.so
LAPACK: /gpfs/apps/AMD/R/4.2.2/GCC/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] condiments_1.5.0            viridis_0.6.2               viridisLite_0.4.1           rgl_1.1.3                   RANN_2.6.1                 
 [6] tradeSeq_1.12.0             slingshot_2.7.0             TrajectoryUtils_1.6.0       SingleCellExperiment_1.20.1 princurve_2.1.6            
[11] htmlwidgets_1.6.1           plotly_4.10.1               scales_1.2.1                cowplot_1.1.1               RColorBrewer_1.1-3         
[16] qusage_2.32.0               limma_3.54.2                SoupX_1.6.2                 fgsea_1.24.0                DESeq2_1.38.3              
[21] SummarizedExperiment_1.28.0 Biobase_2.58.0              MatrixGenerics_1.10.0       matrixStats_0.63.0          GenomicRanges_1.50.2       
[26] GenomeInfoDb_1.34.9         IRanges_2.32.0              S4Vectors_0.36.2            BiocGenerics_0.44.0         DoubletFinder_2.0.3        
[31] scCustomize_1.1.1           data.table_1.14.8           plyr_1.8.8                  stringr_1.5.0               cacoa_0.4.0                
[36] Matrix_1.5-1                pROC_1.18.0                 mixtools_2.0.0              patchwork_1.1.2             ggplot2_3.4.2              
[41] SeuratObject_4.1.3          Seurat_4.3.0                dplyr_1.1.1                

loaded via a namespace (and not attached):
  [1] estimability_1.4.1        ggprism_1.0.4             scattermore_0.8           ModelMetrics_1.2.2.2      spatstat.model_3.2-3     
  [6] Ecume_0.9.1               coda_0.19-4               tidyr_1.3.0               bit64_4.0.5               knitr_1.42               
 [11] irlba_2.3.5.1             multcomp_1.4-22           DelayedArray_0.24.0       rpart_4.1.19              doParallel_1.0.17        
 [16] KEGGREST_1.38.0           hardhat_1.2.0             RCurl_1.98-1.10           generics_0.1.3            ScaledMatrix_1.6.0       
 [21] TH.data_1.1-1             RSQLite_2.3.0             proxy_0.4-27              future_1.31.0             bit_4.0.5                
 [26] spatstat.data_3.0-1       lubridate_1.9.2           httpuv_1.6.9              gower_1.0.1               xfun_0.37                
 [31] promises_1.2.0.1          fansi_1.0.4               igraph_1.4.0              DBI_1.1.3                 geneplotter_1.76.0       
 [36] spatstat.geom_3.2-1       paletteer_1.5.0           purrr_1.0.1               ellipsis_0.3.2            annotate_1.76.0          
 [41] sparseMatrixStats_1.10.0  deldir_1.0-6              vctrs_0.6.1               ROCR_1.0-11               abind_1.4-5              
 [46] caret_6.0-94              cachem_1.0.6              withr_2.5.0               progressr_0.13.0          emmeans_1.8.5            
 [51] sctransform_0.3.5         scran_1.26.2              goftest_1.2-3             cluster_2.1.4             ape_5.7                  
 [56] segmented_1.6-2           lazyeval_0.2.2            crayon_1.5.2              spatstat.explore_3.1-0    labeling_0.4.2           
 [61] edgeR_3.40.2              recipes_1.0.5             pkgconfig_2.0.3           nlme_3.1-160              vipor_0.4.5              
 [66] transport_0.13-0          nnet_7.3-18               rlang_1.1.0               globals_0.16.2            lifecycle_1.0.3          
 [71] miniUI_0.1.1.1            sandwich_3.0-2            rsvd_1.0.5                ggrastr_1.0.1             polyclip_1.10-4          
 [76] lmtest_0.9-40             rngtools_1.5.2            zoo_1.8-11                base64enc_0.1-3           beeswarm_0.4.0           
 [81] ggridges_0.5.4            GlobalOptions_0.1.2       png_0.1-8                 bitops_1.0-7              KernSmooth_2.23-20       
 [86] Biostrings_2.66.0         DelayedMatrixStats_1.20.0 blob_1.2.3                doRNG_1.8.6               shape_1.4.6              
 [91] parallelly_1.34.0         spatstat.random_3.1-4     beachmat_2.14.0           sccore_1.0.3              memoise_2.0.1            
 [96] magrittr_2.0.3            ica_1.0-3                 distinct_1.10.2           zlibbioc_1.44.0           compiler_4.2.2           
[101] dqrng_0.3.0               fitdistrplus_1.1-8        snakecase_0.11.0          cli_3.6.0                 XVector_0.38.0           
[106] listenv_0.9.0             pbapply_1.7-0             MASS_7.3-58.1             mgcv_1.8-41               tidyselect_1.2.0         
[111] fftw_1.0-7                stringi_1.7.12            forcats_1.0.0             BiocSingular_1.14.0       locfit_1.5-9.7           
[116] ggrepel_0.9.3             grid_4.2.2                spatstat.linnet_3.1-0     fastmatch_1.1-3           tools_4.2.2              
[121] timechange_0.2.0          future.apply_1.10.0       parallel_4.2.2            circlize_0.4.15           rstudioapi_0.14          
[126] bluster_1.8.0             foreach_1.5.2             metapod_1.6.0             janitor_2.2.0             gridExtra_2.3            
[131] prodlim_2019.11.13        farver_2.1.1              Rtsne_0.16                digest_0.6.31             shiny_1.7.4              
[136] lava_1.7.2.1              Rcpp_1.0.10               scuttle_1.8.4             later_1.3.0               RcppAnnoy_0.0.20         
[141] httr_1.4.4                AnnotationDbi_1.60.0      kernlab_0.9-32            colorspace_2.1-0          XML_3.99-0.13            
[146] tensor_1.5                reticulate_1.28           splines_4.2.2             statmod_1.5.0             uwot_0.1.14              
[151] rematch2_2.1.2            spatstat.utils_3.0-3      scater_1.26.1             sp_1.6-0                  xtable_1.8-4             
[156] jsonlite_1.8.4            spatstat_3.0-5            timeDate_4022.108         ipred_0.9-14              R6_2.5.1                 
[161] pillar_1.8.1              htmltools_0.5.4           mime_0.12                 glue_1.6.2                fastmap_1.1.0            
[166] BiocParallel_1.32.5       BiocNeighbors_1.16.0      class_7.3-20              codetools_0.2-18          mvtnorm_1.1-3            
[171] utf8_1.2.3                lattice_0.20-45           spatstat.sparse_3.0-0     tibble_3.2.1              ggbeeswarm_0.7.1         
[176] leiden_0.4.3              survival_3.4-0            munsell_0.5.0             e1071_1.7-13              GenomeInfoDbData_1.2.9   
[181] iterators_1.0.14          reshape2_1.4.4            gtable_0.3.1
HectorRDB commented 1 year ago

Hi @jordiabante Sorry for the late response. Yes please do share the rds object at hector dot rouxdebezieux at berkeley.edu thanks

HectorRDB commented 10 months ago

Closed because of inactivity