label hoods in NhoodGraphDA

[Luuk024026]

Dear Dr. Morgan,

Thanks for this usefull package and clear documentation!

I am working with CD8 T cells from three groups (healthy control, kidney transplant patients and kidney transplant patients with a immune supressive drug). Per condition I have 7, 8 and 8 samples respectively. I have two main questions regarding DA testing>

Q1: To test DA between two groups I first subset my dataset to only include those two groups and then follow the rest of the vignette for DA testing. Is this the correct way to do this or is there a better way?

Q2: I found several hoods that are significant under FDR <1. I wanted to extract the cells that are in these hoods to check for markers and annotate them with a celltype. I was able to do this but I wanted to know which hoods came from which position on the Nhoodgraph, so my thought was to label the significant hoods. Then it occured to me that the results table from testNhoods names the hoods as indeces, but these indeces do not match with the hood names as displayed in the plotNhoodGraphDA (.ggraph.index). Is there a way to label the significant hoods in the NhoodGraphDA that correspond to the hoods from the testNhoods table?

Kind regard,
Luuk

The code I use:
milo.object <- DietSeurat(Seurat.object, assays = "RNA", dimreducs = c("umap", "pca"))

milo.object <- as.SingleCellExperiment(milo.object)
milo.meta <- milo.object@colData

milo.obj <- Milo(milo.object)

milo.obj <- milo.obj[ ,milo.obj@colData$group %in% c("KTR", "KTR_MMF")]
milo.meta <- milo.obj@colData

milo.obj <- buildGraph(milo.obj, k=20, d=30)
milo.obj <- makeNhoods(milo.obj, k=20, d=30, refined=TRUE, prop=0.1)
plotNhoodSizeHist(milo.obj)

milo.obj <- calcNhoodDistance(milo.obj, d=30)

milo.obj <- countCells(milo.obj, sample="orig.ident", meta.data=milo.meta)

milo.design <- as.data.frame(xtabs(~ group + orig.ident, data=milo.meta))
milo.design <- milo.design[milo.design$Freq > 0, ] # This gets rid of the group that is not there (HC in this case)
rownames(milo.design) <- milo.design$orig.ident # Make the row names the sample ids
milo.design <- milo.design[colnames(nhoodCounts(milo.obj)),] # This is to align and order your design matrix with your hoods stored in the milo object
milo.design$group <- factor(milo.design$group) # This is added to make sure that milo sees the groups as factors

milo.res <- testNhoods(milo.obj, design = ~group, design.df=milo.design)

significant_nhoods <- milo.res[milo.res$SpatialFDR < 0.1, ]

milo.obj <- buildNhoodGraph(milo.obj)
plotNhoodGraphDA(milo.obj, milo.res, alpha = 0.1)

nhoods <- milo.obj@nhoods
cell_names <- rownames(nhoods)
significant_hood <- cell_names[nhoods[, 136] == 1]

Result:

Session info:
R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 24.04 LTS

Matrix products: default
BLAS: /usr/local/R-4.3.3/lib/libRblas.so
LAPACK: /usr/local/R-4.3.3/lib/libRlapack.so; LAPACK version 3.11.0

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

time zone: Etc/UTC
tzcode source: system (glibc)

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

other attached packages:
[1] loupeR_1.1.4 scales_1.3.0 circlize_0.4.16 Cairo_1.6-2 immunarch_0.9.1
[6] patchwork_1.3.0 data.table_1.17.0 dtplyr_1.3.1 remotes_2.5.0 hdf5r_1.3.12
[11] scater_1.30.1 scuttle_1.12.0 Azimuth_0.5.0 shinyBS_0.61.1 pak_0.8.0.1
[16] scRepertoire_2.3.2 ZINBMM_0.1.0 devtools_2.4.5 usethis_3.0.0 SingleCellExperiment_1.24.0
[21] miloR_1.10.0 edgeR_4.0.16 limma_3.58.1 openxlsx_4.2.8 org.Mm.eg.db_3.18.0
[26] AnnotationDbi_1.64.1 clusterProfiler_4.10.1 msigdbr_7.5.1 fgsea_1.28.0 scCustomize_3.0.1
[31] celldex_1.12.0 SingleR_2.4.1 SummarizedExperiment_1.32.0 Biobase_2.62.0 GenomicRanges_1.54.1
[36] GenomeInfoDb_1.38.8 IRanges_2.36.0 S4Vectors_0.40.2 BiocGenerics_0.48.1 MatrixGenerics_1.14.0
[41] matrixStats_1.5.0 dplyr_1.1.4 clustree_0.5.1 ggraph_2.2.1 ggplot2_3.5.1
[46] Seurat_5.2.1 SeuratObject_5.0.2 sp_2.2-0 BiocManager_1.30.23

loaded via a namespace (and not attached):
[1] igraph_2.1.4 Formula_1.2-5 ica_1.0-3 plotly_4.10.4
[5] rematch2_2.1.2 fpc_2.2-13 zlibbioc_1.48.2 doParallel_1.0.17
[9] tidyselect_1.2.1 bit_4.6.0 lattice_0.22-6 rjson_0.2.23
[13] evmix_2.12 factoextra_1.0.7 blob_1.2.4 stringr_1.5.1
[17] urlchecker_1.0.1 S4Arrays_1.2.1 parallel_4.3.3 seqLogo_1.68.0
[21] png_0.1-8 cli_3.6.4 ggplotify_0.1.2 ProtGenerics_1.34.0
[25] goftest_1.2-3 gargle_1.5.2 BiocIO_1.12.0 kernlab_0.9-33
[29] purrr_1.0.4 BiocNeighbors_1.20.2 Signac_1.14.0 stringdist_0.9.15
[33] uwot_0.2.3 shadowtext_0.1.4 curl_6.2.1 mime_0.12
[37] evaluate_1.0.3 tidytree_0.4.6 stringi_1.8.4 backports_1.5.0
[41] gsl_2.1-8 XML_3.99-0.18 lubridate_1.9.4 httpuv_1.6.15
[45] flexmix_2.3-20 paletteer_1.6.0 magrittr_2.0.3 rappdirs_0.3.3
[49] splines_4.3.3 mclust_6.1.1 RcppRoll_0.3.1 DT_0.33
[53] sctransform_0.4.1 ggbeeswarm_0.7.2 sessioninfo_1.2.2 DBI_1.2.3
[57] withr_3.0.2 class_7.3-22 enrichplot_1.22.0 lmtest_0.9-40
[61] tidygraph_1.3.1 rtracklayer_1.62.0 htmlwidgets_1.6.4 fs_1.6.5
[65] biomaRt_2.58.2 ggrepel_0.9.6 SparseArray_1.2.4 DEoptimR_1.1-3-1
[69] cellranger_1.1.0 diptest_0.77-1 annotate_1.78.0 reticulate_1.41.0.1
[73] zoo_1.8-13 JASPAR2020_0.99.10 XVector_0.42.0 knitr_1.49
[77] TFBSTools_1.40.0 TFMPvalue_0.0.9 timechange_0.3.0 foreach_1.5.2
[81] caTools_1.18.3 grid_4.3.3 rhdf5_2.46.1 ggtree_3.10.1
[85] R.oo_1.27.0 quantreg_6.1 poweRlaw_1.0.0 RSpectra_0.16-2
[89] irlba_2.3.5.1 ggrastr_1.0.2 fastDummies_1.7.5 gridGraphics_0.5-1
[93] ellipsis_0.3.2 lazyeval_0.2.2 yaml_2.3.10 survival_3.7-0
[97] scattermore_1.2 BiocVersion_3.18.1 crayon_1.5.3 RcppAnnoy_0.0.22
[101] RColorBrewer_1.1-3 tidyr_1.3.1 progressr_0.15.1 tweenr_2.0.3
[105] later_1.4.1 ggridges_0.5.6 codetools_0.2-20 GlobalOptions_0.1.2
[109] profvis_0.3.8 rlist_0.4.6.2 KEGGREST_1.42.0 Rtsne_0.17
[113] shape_1.4.6.1 Rsamtools_2.18.0 filelock_1.0.3 pkgconfig_2.0.3
[117] ggpubr_0.6.0 xml2_1.3.7 spatstat.univar_3.1-2 GenomicAlignments_1.38.2
[121] aplot_0.2.5 iNEXT_3.0.1 evd_2.3-7.1 spatstat.sparse_3.1-0
[125] BSgenome_1.70.2 ape_5.8-1 viridisLite_0.4.2 xtable_1.8-4
[129] car_3.1-3 plyr_1.8.9 httr_1.4.7 prabclus_2.3-4
[133] tools_4.3.3 globals_0.16.3 pkgbuild_1.4.4 broom_1.0.7
[137] beeswarm_0.4.0 nlme_3.1-165 HDO.db_0.99.1 dbplyr_2.5.0
[141] ExperimentHub_2.10.0 shinyjs_2.1.0 MatrixModels_0.5-3 assertthat_0.2.1
[145] digest_0.6.37 Matrix_1.6-5 farver_2.1.2 tzdb_0.4.0
[149] AnnotationFilter_1.26.0 reshape2_1.4.4 yulab.utils_0.2.0 viridis_0.6.5
[153] DirichletMultinomial_1.44.0 glue_1.8.0 cachem_1.1.0 BiocFileCache_2.10.2
[157] polyclip_1.10-7 UpSetR_1.4.0 generics_0.1.3 Biostrings_2.70.3
[161] ggalluvial_0.12.5 googledrive_2.1.1 presto_1.0.0 parallelly_1.42.0
[165] pkgload_1.4.0 statmod_1.5.0 RcppHNSW_0.6.0 ScaledMatrix_1.10.0
[169] carData_3.0-5 pbapply_1.7-2 spam_2.11-1 gson_0.1.0
[173] graphlayouts_1.2.2 gtools_3.9.5 readxl_1.4.5 ggsignif_0.6.4
[177] gridExtra_2.3 shiny_1.10.0 GenomeInfoDbData_1.2.11 rhdf5filters_1.14.1
[181] R.utils_2.13.0 RCurl_1.98-1.16 memoise_2.0.1 rmarkdown_2.29
[185] pheatmap_1.0.12 R.methodsS3_1.8.2 googlesheets4_1.1.1 future_1.34.0
[189] RANN_2.6.2 spatstat.data_3.1-4 rstudioapi_0.16.0 cluster_2.1.6
[193] janitor_2.2.1 spatstat.utils_3.1-2 hms_1.1.3 fitdistrplus_1.2-2
[197] munsell_0.5.1 cowplot_1.1.3 colorspace_2.1-1 quadprog_1.5-8
[201] rlang_1.1.5 ggdendro_0.2.0 shinythemes_1.2.0 DelayedMatrixStats_1.24.0
[205] sparseMatrixStats_1.14.0 shinydashboard_0.7.2 truncdist_1.0-2 dotCall64_1.2
[209] ggforce_0.4.2 xfun_0.51 ggseqlogo_0.2 CNEr_1.38.0
[213] modeltools_0.2-23 iterators_1.0.14 abind_1.4-8 GOSemSim_2.28.1
[217] tibble_3.2.1 interactiveDisplayBase_1.40.0 EnsDb.Hsapiens.v86_2.99.0 treeio_1.26.0
[221] Rhdf5lib_1.24.2 readr_2.1.5 bitops_1.0-9 promises_1.3.2
[225] scatterpie_0.2.4 RSQLite_2.3.9 qvalue_2.34.0 DelayedArray_0.28.0
[229] GO.db_3.18.0 compiler_4.3.3 forcats_1.0.0 prettyunits_1.2.0
[233] beachmat_2.18.1 SparseM_1.84-2 listenv_0.9.1 BSgenome.Hsapiens.UCSC.hg38_1.4.5
[237] Rcpp_1.0.14 AnnotationHub_3.10.1 BiocSingular_1.18.0 tensor_1.5
[241] MASS_7.3-60.0.1 progress_1.2.3 uuid_1.2-1 BiocParallel_1.36.0
[245] cubature_2.1.1 babelgene_22.9 spatstat.random_3.3-2 R6_2.6.1
[249] rstatix_0.7.2 fastmap_1.2.0 fastmatch_1.1-6 vipor_0.4.7
[253] ensembldb_2.26.0 ROCR_1.0-11 SeuratDisk_0.0.0.9021 nnet_7.3-19
[257] rsvd_1.0.5 gtable_0.3.6 phangorn_2.12.1 KernSmooth_2.23-24
[261] miniUI_0.1.1.1 deldir_2.0-4 htmltools_0.5.8.1 bit64_4.6.0-1
[265] spatstat.explore_3.3-4 lifecycle_1.0.4 ggprism_1.0.5 Brobdingnag_1.2-9
[269] zip_2.3.1 restfulr_0.0.15 vctrs_0.6.5 robustbase_0.99-4-1
[273] VGAM_1.1-13 spatstat.geom_3.3-5 snakecase_0.11.1 DOSE_3.28.2
[277] ggfun_0.1.8 SeuratData_0.2.2.9001 future.apply_1.11.3 pillar_1.10.1
[281] GenomicFeatures_1.54.4 locfit_1.5-9.10 jsonlite_1.9.1

[MikeDMorgan]

Hi @Luuk024026
It is possible to perform a DA test comparing just two groups using contrasts: https://www.bioconductor.org/packages/release/bioc/vignettes/miloR/inst/doc/milo_contrasts.html

This would avoid doing any subsetting and re-running of whole Milo workflows.

Re: extracting cells - the nhood names are in the same order as the columns of the cell x nhood matrix from the nhoods() function.

[Luuk024026]

Thanks mr Morgan!

I have one additional question.

In the milo publication is written:

" To control the spatial FDR in the KNN graph, we apply a weighted version of the Benjamini-Hochberg (BH) method, where P values are weighted by the reciprocal of the neighbourhood connectivity"

So, if the neighborhood is "dense", then the Eucledian distance to the K-nearest neighbor is small

Is there any way you can explain this better. I understand BH, and how you calculate the connectivity,
But how do you then include those connectivities to obtain weights in BH?

Many thanks

[MikeDMorgan]

https://academic.oup.com/biostatistics/article/18/1/91/2555340

[Luuk024026]

Thank you for your reference to the weighted FDR paper. I understand that in dense areas of the "single cell space" you would allow more false discoveries than in less dense areas. Can you give me the rationale behind that?
Are those areas more interesting because they are more densely sampled and thus more reflective of the biology?

[MikeDMorgan]

The weight is inversely proportional to amount of overlap with adjacent neighbourhoods. So if there is a densely connected area in the graph, with lots of overlap, then these are down-weighted to reduce the number of false positives.

[Luuk024026]

Thanks for the answer, just to be sure I understand this correctly. In the MILO paper you write:

"To control the spatial FDR in the KNN graph, we apply a weighted version of the Benjamini–Hochberg method, where P values are weighted by the reciprocal of the neighborhood connectivity. As a measure of neighborhood connectivity, we use the Euclidean distance to the kth nearest neighbor of the index cell for each neighborhood."

You use the Eucledian distance to the kth nearest neighbor to measure connectivity, so if that distance is large do you expect little or much overlap with adjacent neighborhoods?

Many thanks Luuk