Error message with findNhoodGroupMarkers

[gianfilippo]

Error message with findNhoodGroupMarkers

[MikeDMorgan]

Hi @gianfilippo as stated on the github repo landing page please report the output of your sessionInfo() when reporting issues.

[gianfilippo]

sorry about it. Here it is
R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.9 (Maipo)

Matrix products: default
BLAS/LAPACK: /gpfs/ycga/apps/hpc/software/OpenBLAS/0.3.12-GCC-10.2.0/lib/libopenblas_haswellp-r0.3.12.so

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

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

other attached packages:
[1] miloR_1.5.0 edgeR_3.40.1
[3] limma_3.54.0 HGNChelper_0.8.1
[5] ProjecTILs_3.0.1 DoubletFinder_2.0.3
[7] scDblFinder_1.13.7 djvdj_0.0.0.9000
[9] scater_1.24.0 scuttle_1.6.3
[11] SingleCellExperiment_1.20.0 SummarizedExperiment_1.26.1
[13] Biobase_2.56.0 GenomicRanges_1.48.0
[15] GenomeInfoDb_1.32.4 IRanges_2.30.1
[17] S4Vectors_0.34.0 BiocGenerics_0.42.0
[19] MatrixGenerics_1.10.0 matrixStats_0.62.0
[21] muscat_1.10.1 CIPR_0.1.0
[23] ComplexHeatmap_2.12.1 gprofiler2_0.2.1
[25] ggplot2_3.4.0 cowplot_1.1.1
[27] future_1.28.0 purrr_1.0.1
[29] openxlsx_4.2.5 patchwork_1.1.1
[31] miQC_1.4.0 SeuratWrappers_0.3.1
[33] SeuratDisk_0.0.0.9020 glmGamPoi_1.8.0
[35] SeuratObject_4.1.3 Seurat_4.3.0
[37] dplyr_1.1.0

loaded via a namespace (and not attached):
[1] rsvd_1.0.5 ica_1.0-2
[3] Rsamtools_2.12.0 foreach_1.5.2
[5] lmtest_0.9-40 crayon_1.5.1
[7] rbibutils_2.2.8 MASS_7.3-57
[9] nlme_3.1-157 backports_1.4.1
[11] rlang_1.0.6 XVector_0.36.0
[13] ROCR_1.0-11 irlba_2.3.5.1
[15] nloptr_2.0.0 xgboost_1.6.0.1
[17] BiocParallel_1.32.4 rjson_0.2.21
[19] bit64_4.0.5 glue_1.6.2
[21] sctransform_0.3.5 pbkrtest_0.5.1
[23] parallel_4.2.0 vipor_0.4.5
[25] spatstat.sparse_3.0-0 AnnotationDbi_1.60.0
[27] spatstat.geom_3.0-3 tidyselect_1.2.0
[29] fitdistrplus_1.1-8 variancePartition_1.26.0
[31] XML_3.99-0.9 tidyr_1.2.0
[33] zoo_1.8-10 GenomicAlignments_1.32.1
[35] xtable_1.8-4 magrittr_2.0.3
[37] Rdpack_2.3 cli_3.6.0
[39] zlibbioc_1.42.0 dbscan_1.1-11
[41] miniUI_0.1.1.1 sp_1.5-0
[43] aod_1.3.2 shiny_1.7.1
[45] BiocSingular_1.12.0 clue_0.3-60
[47] cluster_2.1.3 caTools_1.18.2
[49] tidygraph_1.2.1 KEGGREST_1.36.3
[51] tibble_3.1.8 ggrepel_0.9.1
[53] listenv_0.8.0 Biostrings_2.64.1
[55] png_0.1-7 withr_2.5.0
[57] bitops_1.0-7 ggforce_0.4.1
[59] plyr_1.8.7 pracma_2.4.2
[61] dqrng_0.3.0 coda_0.19-4
[63] pillar_1.8.1 gplots_3.1.3
[65] GlobalOptions_0.1.2 cachem_1.0.6
[67] multcomp_1.4-18 flexmix_2.3-17
[69] hdf5r_1.3.5 GetoptLong_1.0.5
[71] DelayedMatrixStats_1.18.0 vctrs_0.5.2
[73] ellipsis_0.3.2 generics_0.1.2
[75] tools_4.2.0 beeswarm_0.4.0
[77] munsell_0.5.0 tweenr_2.0.2
[79] emmeans_1.8.1-1 DelayedArray_0.22.0
[81] fastmap_1.1.0 compiler_4.2.0
[83] abind_1.4-5 httpuv_1.6.5
[85] rtracklayer_1.56.1 plotly_4.10.0
[87] GenomeInfoDbData_1.2.8 gridExtra_2.3
[89] glmmTMB_1.1.5 lattice_0.20-45
[91] deldir_1.0-6 utf8_1.2.2
[93] later_1.3.0 jsonlite_1.8.0
[95] scales_1.2.0 ScaledMatrix_1.4.0
[97] pbapply_1.5-0 sparseMatrixStats_1.8.0
[99] estimability_1.4.1 genefilter_1.78.0
[101] lazyeval_0.2.2 promises_1.2.0.1
[103] doParallel_1.0.17 R.utils_2.11.0
[105] goftest_1.2-3 spatstat.utils_3.0-1
[107] reticulate_1.26 sandwich_3.0-1
[109] blme_1.0-5 statmod_1.4.36
[111] Rtsne_0.16 uwot_0.1.14
[113] igraph_1.3.1 survival_3.3-1
[115] numDeriv_2016.8-1.1 yaml_2.3.5
[117] htmltools_0.5.2 memoise_2.0.1
[119] modeltools_0.2-23 BiocIO_1.6.0
[121] locfit_1.5-9.5 graphlayouts_0.8.0
[123] viridisLite_0.4.0 digest_0.6.29
[125] RhpcBLASctl_0.21-247.1 mime_0.12
[127] RSQLite_2.2.12 future.apply_1.9.0
[129] remotes_2.4.2 data.table_1.14.2
[131] blob_1.2.3 R.oo_1.24.0
[133] splines_4.2.0 RCurl_1.98-1.6
[135] broom_1.0.1 hms_1.1.1
[137] colorspace_2.0-3 BiocManager_1.30.19
[139] ggbeeswarm_0.6.0 shape_1.4.6
[141] nnet_7.3-17 Rcpp_1.0.8.3
[143] RANN_2.6.1 mvtnorm_1.1-3
[145] circlize_0.4.14 fansi_1.0.3
[147] tzdb_0.3.0 parallelly_1.32.1
[149] R6_2.5.1 ggridges_0.5.3
[151] lifecycle_1.0.3 zip_2.2.0
[153] bluster_1.6.0 abdiv_0.2.0
[155] minqa_1.2.4 leiden_0.4.3
[157] Matrix_1.5-0 RcppAnnoy_0.0.19
[159] TH.data_1.1-0 RColorBrewer_1.1-3
[161] iterators_1.0.14 spatstat.explore_3.0-5
[163] TMB_1.9.1 stringr_1.4.0
[165] htmlwidgets_1.5.4 beachmat_2.12.0
[167] polyclip_1.10-4 iNEXT_3.0.0
[169] globals_0.16.1 spatstat.random_3.0-1
[171] progressr_0.11.0 codetools_0.2-18
[173] metapod_1.4.0 gtools_3.9.2
[175] prettyunits_1.1.1 R.methodsS3_1.8.1
[177] gtable_0.3.0 DBI_1.1.2
[179] tensor_1.5 httr_1.4.2
[181] KernSmooth_2.23-20 stringi_1.7.6
[183] progress_1.2.2 farver_2.1.0
[185] reshape2_1.4.4 annotate_1.74.0
[187] viridis_0.6.2 boot_1.3-28
[189] BiocNeighbors_1.14.0 lme4_1.1-29
[191] restfulr_0.0.15 readr_2.1.2
[193] geneplotter_1.74.0 scattermore_0.8
[195] DESeq2_1.36.0 scran_1.24.1
[197] bit_4.0.4 spatstat.data_3.0-0
[199] ggraph_2.1.0 pkgconfig_2.0.3
[201] lmerTest_3.1-3

[MikeDMorgan]

Ok - first things first - make sure you are using the most up to date version of Milo (1.6) from Bioconductor (3.16).

[gianfilippo]

Hi,

I updated miloR to 1.6.0, but I see the same issue.
I have the following

da_results_Ann is the output from the following steps: testNhoods, groupNhoods and annotateNhoods
tmiloOBJ is the original after logNormCounts and gene filtering

isig = da_results_Ann$SpatialFDR < 0.05
subnhoods = da_results_Ann[isig,"Nhood"]
subgroups = da_results_Ann[isig,"NhoodGroup"]

I tried the following combinations

1. Specifying both subset.groups and subset.nhoods
   da_nhood_markers <- findNhoodGroupMarkers(tmiloOBJ, da_results_Ann, "counts", aggregate.samples=T, sample_col="SampleName", gene.offset=T, subset.groups=subgroups, subset.nhoods=subnhoods)
   Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) :
   contrasts can be applied only to factors with 2 or more levels

2. Specifying only subset.nhoods
   da_nhood_markers <- findNhoodGroupMarkers(tmiloOBJ, da_results_Ann, "counts", aggregate.samples=T, sample_col="SampleName", gene.offset=T, subset.nhoods=subnhoods)
   Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) :
   contrasts can be applied only to factors with 2 or more levels

3. Specifying only subset.groups
   da_nhood_markers <- findNhoodGroupMarkers(traj_milo_da, da_results_Ann, "counts", aggregate.samples=T, sample_col="SampleName", gene.offset=T, subset.groups=subgroups)
   Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) :
   contrasts can be applied only to factors with 2 or more levels

4. Specifying neither subset.groups or subset.nhoods, but this is taking a very long time.

Thanks for your help!

sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.9 (Maipo)

Matrix products: default
BLAS/LAPACK: /gpfs/ycga/apps/hpc/software/OpenBLAS/0.3.12-GCC-10.2.0/lib/libopenblas_haswellp-r0.3.12.so

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

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

other attached packages:
[1] HGNChelper_0.8.1 ProjecTILs_3.0.1
[3] DoubletFinder_2.0.3 scDblFinder_1.13.7
[5] djvdj_0.0.0.9000 scater_1.24.0
[7] scuttle_1.6.3 SingleCellExperiment_1.20.0
[9] SummarizedExperiment_1.26.1 Biobase_2.56.0
[11] GenomicRanges_1.48.0 GenomeInfoDb_1.32.4
[13] IRanges_2.30.1 S4Vectors_0.34.0
[15] BiocGenerics_0.42.0 MatrixGenerics_1.10.0
[17] matrixStats_0.62.0 muscat_1.10.1
[19] CIPR_0.1.0 ComplexHeatmap_2.12.1
[21] gprofiler2_0.2.1 ggplot2_3.4.0
[23] cowplot_1.1.1 future_1.28.0
[25] purrr_1.0.1 openxlsx_4.2.5
[27] patchwork_1.1.1 miQC_1.4.0
[29] SeuratWrappers_0.3.1 SeuratDisk_0.0.0.9020
[31] glmGamPoi_1.8.0 SeuratObject_4.1.3
[33] Seurat_4.3.0 dplyr_1.1.0
[35] miloR_1.6.0 edgeR_3.40.1
[37] limma_3.54.0

loaded via a namespace (and not attached):
[1] rsvd_1.0.5 ica_1.0-2
[3] Rsamtools_2.12.0 foreach_1.5.2
[5] lmtest_0.9-40 crayon_1.5.1
[7] rbibutils_2.2.8 MASS_7.3-57
[9] nlme_3.1-157 backports_1.4.1
[11] rlang_1.0.6 XVector_0.36.0
[13] ROCR_1.0-11 irlba_2.3.5.1
[15] nloptr_2.0.0 xgboost_1.6.0.1
[17] BiocParallel_1.32.4 rjson_0.2.21
[19] bit64_4.0.5 glue_1.6.2
[21] sctransform_0.3.5 pbkrtest_0.5.1
[23] parallel_4.2.0 vipor_0.4.5
[25] spatstat.sparse_3.0-0 AnnotationDbi_1.60.0
[27] spatstat.geom_3.0-3 tidyselect_1.2.0
[29] fitdistrplus_1.1-8 variancePartition_1.26.0
[31] XML_3.99-0.9 tidyr_1.2.0
[33] zoo_1.8-10 GenomicAlignments_1.32.1
[35] xtable_1.8-4 magrittr_2.0.3
[37] Rdpack_2.3 cli_3.6.0
[39] zlibbioc_1.42.0 dbscan_1.1-11
[41] miniUI_0.1.1.1 sp_1.5-0
[43] aod_1.3.2 shiny_1.7.1
[45] BiocSingular_1.12.0 clue_0.3-60
[47] cluster_2.1.3 caTools_1.18.2
[49] tidygraph_1.2.1 KEGGREST_1.36.3
[51] tibble_3.1.8 ggrepel_0.9.1
[53] listenv_0.8.0 Biostrings_2.64.1
[55] png_0.1-7 withr_2.5.0
[57] bitops_1.0-7 ggforce_0.4.1
[59] plyr_1.8.7 pracma_2.4.2
[61] dqrng_0.3.0 coda_0.19-4
[63] pillar_1.8.1 gplots_3.1.3
[65] GlobalOptions_0.1.2 cachem_1.0.6
[67] multcomp_1.4-18 flexmix_2.3-17
[69] hdf5r_1.3.5 GetoptLong_1.0.5
[71] DelayedMatrixStats_1.18.0 vctrs_0.5.2
[73] ellipsis_0.3.2 generics_0.1.2
[75] tools_4.2.0 beeswarm_0.4.0
[77] munsell_0.5.0 tweenr_2.0.2
[79] emmeans_1.8.1-1 DelayedArray_0.22.0
[81] fastmap_1.1.0 compiler_4.2.0
[83] abind_1.4-5 httpuv_1.6.5
[85] rtracklayer_1.56.1 plotly_4.10.0
[87] GenomeInfoDbData_1.2.8 gridExtra_2.3
[89] glmmTMB_1.1.5 lattice_0.20-45
[91] deldir_1.0-6 utf8_1.2.2
[93] later_1.3.0 jsonlite_1.8.0
[95] scales_1.2.0 ScaledMatrix_1.4.0
[97] pbapply_1.5-0 sparseMatrixStats_1.8.0
[99] estimability_1.4.1 genefilter_1.78.0
[101] lazyeval_0.2.2 promises_1.2.0.1
[103] doParallel_1.0.17 R.utils_2.11.0
[105] goftest_1.2-3 spatstat.utils_3.0-1
[107] reticulate_1.26 sandwich_3.0-1
[109] blme_1.0-5 statmod_1.4.36
[111] Rtsne_0.16 uwot_0.1.14
[113] igraph_1.3.1 survival_3.3-1
[115] numDeriv_2016.8-1.1 yaml_2.3.5
[117] htmltools_0.5.2 memoise_2.0.1
[119] modeltools_0.2-23 BiocIO_1.6.0
[121] locfit_1.5-9.5 graphlayouts_0.8.0
[123] viridisLite_0.4.0 digest_0.6.29
[125] RhpcBLASctl_0.21-247.1 mime_0.12
[127] RSQLite_2.2.12 future.apply_1.9.0
[129] remotes_2.4.2 data.table_1.14.2
[131] blob_1.2.3 R.oo_1.24.0
[133] labeling_0.4.2 splines_4.2.0
[135] RCurl_1.98-1.6 broom_1.0.1
[137] hms_1.1.1 colorspace_2.0-3
[139] BiocManager_1.30.19 ggbeeswarm_0.6.0
[141] shape_1.4.6 nnet_7.3-17
[143] Rcpp_1.0.8.3 RANN_2.6.1
[145] mvtnorm_1.1-3 circlize_0.4.14
[147] fansi_1.0.3 tzdb_0.3.0
[149] parallelly_1.32.1 R6_2.5.1
[151] ggridges_0.5.3 lifecycle_1.0.3
[153] zip_2.2.0 bluster_1.6.0
[155] abdiv_0.2.0 minqa_1.2.4
[157] leiden_0.4.3 Matrix_1.5-0
[159] RcppAnnoy_0.0.19 TH.data_1.1-0
[161] RColorBrewer_1.1-3 iterators_1.0.14
[163] spatstat.explore_3.0-5 TMB_1.9.1
[165] stringr_1.4.0 htmlwidgets_1.5.4
[167] beachmat_2.12.0 polyclip_1.10-4
[169] iNEXT_3.0.0 globals_0.16.1
[171] spatstat.random_3.0-1 progressr_0.11.0
[173] codetools_0.2-18 metapod_1.4.0
[175] gtools_3.9.2 prettyunits_1.1.1
[177] R.methodsS3_1.8.1 gtable_0.3.0
[179] DBI_1.1.2 tensor_1.5
[181] httr_1.4.2 KernSmooth_2.23-20
[183] stringi_1.7.6 progress_1.2.2
[185] reshape2_1.4.4 farver_2.1.0
[187] annotate_1.74.0 viridis_0.6.2
[189] boot_1.3-28 BiocNeighbors_1.14.0
[191] lme4_1.1-29 restfulr_0.0.15
[193] readr_2.1.2 geneplotter_1.74.0
[195] scattermore_0.8 DESeq2_1.36.0
[197] scran_1.24.1 bit_4.0.4
[199] spatstat.data_3.0-0 ggraph_2.1.0
[201] pkgconfig_2.0.3 lmerTest_3.1-3

[MikeDMorgan]

This will happen where the grouping contains insufficient nhoods - how many nhood groups do you have, and how many of them have ≤ 2 nhoods as members?

[gianfilippo]

Hi,

because of the issue I reported related to groupNhoods, I run the groupNhoods function with a very stringent max.lfc.delta=0.1, which probably ends up with various nhood groups with a single nhood.

What I see within the function is when building the fake.meta data.frame, at the line below, as I reported previously.

fake.meta[nhood.gr.cells, "Nhood.Group"] <- ifelse(is.na(fake.meta[nhood.gr.cells, "Nhood.Group"]), nhood.gr[i], NA)

How does this changes with groups having >= 2 nhoods ?

Thanks

[MikeDMorgan]

The error tells you the problem - you can't compare 2 groups if there are no observations in one of those groups. This will happen if you have nhood groups with very similar nhoods, i.e. many overlapping cells. The assignment of single-cells to nhood groups inside findNhoodMarkers is greedy so, because a cell cannot be assigned to multiple nhood groups for the aggregation step prior to marker identification. If you have fewer larger nhood groups then this should not be a problem.

[gianfilippo]

I see. I suspected that the reason the line of code was replacing the nhoodgroup, causing the problem, was because of overlapping cells.
Now this is more clear.
I will follow your suggestions for the other issue and, if that solves the grouping issue, than this should be solve as well.

Thanks!!