MarioniLab / miloR

R package implementation of Milo for testing for differential abundance in KNN graphs
https://bioconductor.org/packages/release/bioc/html/miloR.html
GNU General Public License v3.0
316 stars 20 forks source link

plotNhoodGraphDA does not display logFC data #309

Open tnystul opened 4 months ago

tnystul commented 4 months ago

Hello, I am running this pipeline on four datasets (2 replicates x 2 genotypes, ~50k cells total) that I merged, SC transformed in Seurat and then converted to an SingleCellExperiment. Running through the code as it is in the tutorial (except with prop = 0.2 and k = 30), the first error I ran into was that the testNhoods reeturned NAs for the SpatialFDR column. The solution posted here fixed this issue for me, and I can confirm that I see values in both the SpatialFDF and logFC columns. However, when I got to the step of actually plotting the results, the logFC is 0 for all neighborhoods. It seems to be an issue with the graphing function because it will assign 0 for any column I chose (by setting res_column equal to another column). Can you please advise? Output

Combined.meta <- combined@meta.data
Combined.SCE <- as.SingleCellExperiment(combined)
colData(Combined.SCE) <- DataFrame(Combined.meta)

logcounts(Combined.SCE) <- log(counts(Combined.SCE) + 1)
Combined.SCE <- runPCA(Combined.SCE, ncomponents=30)
Combined.SCE <- runUMAP(Combined.SCE)
plotUMAP(Combined.SCE)

Comb_milo <- Milo(Combined.SCE)
reducedDim(Comb_milo, "UMAP") <- reducedDim(Combined.SCE, "UMAP")
Comb_milo

Comb_milo <- buildGraph(Comb_milo, k = 30, d = 30)
Comb_milo <- makeNhoods(Comb_milo, prop = 0.2, k = 30, d=40, refined = TRUE, refinement_scheme="graph")

plotNhoodSizeHist(Comb_milo)

Comb_milo <- countCells(Comb_milo, meta.data = data.frame(colData(Comb_milo)), sample="dataset")
head(nhoodCounts(Comb_milo))

Comb_design <- data.frame(colData(Comb_milo))[,c("dataset", "genotype")]
Comb_design <- distinct(Comb_design)
Comb_design

Comb_milo <- calcNhoodDistance(Comb_milo, d=30)

rownames(Comb_design) <- Comb_design$dataset
da_results <- testNhoods(Comb_milo, design = ~ genotype, design.df = Comb_design, fdr.weighting="graph-overlap")

da_results %>%
  arrange(- SpatialFDR) %>%
  head()

Comb_milo <- buildNhoodGraph(Comb_milo)

plotUMAP(Comb_milo) + plotNhoodGraphDA(Comb_milo, da_results, alpha=0.05, res_column = "logFC") +
  plot_layout(guides="collect")

Session info R version 4.2.2 (2022-10-31) Platform: aarch64-apple-darwin20 (64-bit) Running under: macOS Monterey 12.1

Matrix products: default LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages: [1] patchwork_1.1.2 dplyr_1.1.0 scater_1.26.1 ggplot2_3.4.0
[5] scuttle_1.8.4 SingleCellExperiment_1.20.1 SummarizedExperiment_1.28.0 Biobase_2.58.0
[9] GenomicRanges_1.50.2 GenomeInfoDb_1.34.7 IRanges_2.32.0 S4Vectors_0.36.1
[13] BiocGenerics_0.44.0 MatrixGenerics_1.10.0 matrixStats_0.63.0 miloR_1.99.9
[17] edgeR_3.40.2 limma_3.54.1 SeuratObject_4.1.3 Seurat_4.3.0

loaded via a namespace (and not attached): [1] plyr_1.8.8 igraph_1.3.5 lazyeval_0.2.2 sp_1.6-0
[5] splines_4.2.2 BiocParallel_1.32.5 listenv_0.9.0 scattermore_0.8
[9] digest_0.6.31 htmltools_0.5.4 viridis_0.6.2 fansi_1.0.4
[13] magrittr_2.0.3 ScaledMatrix_1.6.0 tensor_1.5 cluster_2.1.4
[17] ROCR_1.0-11 graphlayouts_0.8.4 globals_0.16.2 R.utils_2.12.2
[21] spatstat.sparse_3.0-0 colorspace_2.1-0 ggrepel_0.9.5 xfun_0.36
[25] RCurl_1.98-1.10 crayon_1.5.2 jsonlite_1.8.4 progressr_0.13.0
[29] spatstat.data_3.0-0 survival_3.5-0 zoo_1.8-11 glue_1.6.2
[33] polyclip_1.10-4 gtable_0.3.1 zlibbioc_1.44.0 XVector_0.38.0
[37] leiden_0.4.3 DelayedArray_0.24.0 BiocSingular_1.14.0 future.apply_1.10.0
[41] abind_1.4-5 scales_1.2.1 DBI_1.1.3 spatstat.random_3.1-3
[45] miniUI_0.1.1.1 Rcpp_1.0.10 viridisLite_0.4.1 xtable_1.8-4
[49] reticulate_1.31 rsvd_1.0.5 htmlwidgets_1.6.1 httr_1.4.4
[53] FNN_1.1.3.1 RColorBrewer_1.1-3 ellipsis_0.3.2 ica_1.0-3
[57] pkgconfig_2.0.3 R.methodsS3_1.8.2 farver_2.1.1 sass_0.4.5
[61] uwot_0.1.14 deldir_1.0-6 locfit_1.5-9.7 utf8_1.2.2
[65] tidyselect_1.2.0 labeling_0.4.2 rlang_1.0.6 reshape2_1.4.4
[69] later_1.3.0 munsell_0.5.0 tools_4.2.2 cachem_1.0.6
[73] cli_3.6.0 generics_0.1.3 ggridges_0.5.4 evaluate_0.20
[77] stringr_1.5.0 fastmap_1.1.0 yaml_2.3.7 goftest_1.2-3
[81] knitr_1.42 tidygraph_1.2.3 fitdistrplus_1.1-8 purrr_1.0.1
[85] RANN_2.6.1 ggraph_2.1.0 sparseMatrixStats_1.10.0 pbapply_1.7-0
[89] future_1.30.0 nlme_3.1-161 mime_0.12 ggrastr_1.0.2
[93] R.oo_1.25.0 compiler_4.2.2 rstudioapi_0.14 beeswarm_0.4.0
[97] plotly_4.10.1 png_0.1-8 spatstat.utils_3.0-1 statmod_1.5.0
[101] tweenr_2.0.2 tibble_3.1.8 bslib_0.4.2 stringi_1.7.12
[105] lattice_0.20-45 Matrix_1.5-3 vctrs_0.5.2 pillar_1.8.1
[109] lifecycle_1.0.3 spatstat.geom_3.0-6 lmtest_0.9-40 jquerylib_0.1.4
[113] BiocNeighbors_1.16.0 RcppAnnoy_0.0.20 bitops_1.0-7 data.table_1.14.6
[117] cowplot_1.1.1 irlba_2.3.5.1 httpuv_1.6.8 R6_2.5.1
[121] promises_1.2.0.1 KernSmooth_2.23-20 gridExtra_2.3 vipor_0.4.5
[125] parallelly_1.34.0 codetools_0.2-18 gtools_3.9.4 MASS_7.3-58.2
[129] withr_2.5.0 sctransform_0.3.5 GenomeInfoDbData_1.2.9 parallel_4.2.2
[133] beachmat_2.14.2 grid_4.2.2 tidyr_1.3.0 DelayedMatrixStats_1.20.0 [137] rmarkdown_2.20 Rtsne_0.16 ggforce_0.4.1 spatstat.explore_3.0-6
[141] numDeriv_2016.8-1.1 shiny_1.7.4 ggbeeswarm_0.7.1

MikeDMorgan commented 4 months ago

Hi @tnystul Please see the default behaviour for plotDANhoods, i.e. that alpha=0.05 means it will only plot DA nhoods with spatial FDR ≤ 5%. If you don't have any DA nhoods at this threshold it will plot everything as white points.

tnystul commented 4 months ago

Thank you for the prompt reply!