brgew
commented
4 years ago Hi, I believe that this problem is fixed in the Monocle3 develop branch on GitHub. You can install it using the command
devtools::install_github('https://github.com/cole-trapnell-lab/monocle3',ref='develop')
Please let me know if the problem persists for you.
hpliner
If this is a question and not a bug report or enhancement request, please post to our google group at https://groups.google.com/forum/#!forum/monocle-3-users
Describe the bug Psuedotime estimates fail to span across partitions, even when trajectory graph does (i.e. when "learn_graph(use.partition = FALSE)"
To Reproduce The code that produced the bug: ######################################## suppressPackageStartupMessages({ library(tidymodels) library(devtools) library(monocle3) library(furrr) library(viridis) library(piano) library(UpSetR) library(snowfall) plan(multicore) })
DelayedArray:::set_verbose_block_processing(TRUE) options(DelayedArray.block.size=1000e7)
load cds and set up colData
sciPlex_cds <- readRDS("sciPlex2_cds.RDS")
CDS can be downloaded here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4150377
sciPlex_cds <- sciPlex_cds[, !is.na(pData(sciPlex_cds)$top_to_second_best_ratio)]
Subset cells that have more than a 10 fold enrichment in hash oligos, and over 30 hash oligo UMIs
sciPlex_cds <- sciPlex_cds[, pData(sciPlex_cds)$top_to_second_best_ratio > 10 & pData(sciPlex_cds)$qval < 0.01 & pData(sciPlex_cds)$hash_umis > 30]
colData(sciPlex_cds)$cell_type = "A549"
meta_data = stringr::str_split_fixed(colData(sciPlex_cds)$topoligo, pattern = "", n = 3)
colData(sciPlex_cds)$cell = colData(sciPlex_cds)$Cell
colData(sciPlex_cds)$dose = meta_data[, 2] %>% as.numeric()
colData(sciPlex_cds)$treatment = meta_data[, 1] %>% as.character()
colData(sciPlex_cds)$culture_well = meta_data[, 3] %>% stringr::str_sub(start = 2, end = 4) %>% as.character()
colData(sciPlex_cds)$well_oligo = meta_data[, 3] %>% as.character()
colData(sciPlex_cds)$culture_plate = meta_data[, 3] %>% stringr::str_sub(start = 1, end = 1) %>% as.character()
colData(sciPlex_cds)$vehicle = colData(sciPlex_cds)$dose == 0
colData(sciPlex_cds)$dose_character = factor(colData(sciPlex_cds)$dose, levels = c("0", "0.1", "0.5", "1", "5", "10", "50", "100"))
colData(sciPlex_cds)$new_treatment_label = sapply(colData(sciPlex_cds)$treatment, function(x) { if (grepl("BMS", x)) return("BMS345541") if (grepl("Dex", x)) return("Dex") if (grepl("Nutlin", x)) return("Nutlin3A") if (grepl("SAHA", x)) return("SAHA") })
colData(sciPlex_cds)$solvent = sapply(pData(sciPlex_cds)$treatment, function(x) { if (grepl("Dex", x)) return("Ethanol") else return("DMSO") })
colData(sciPlex_cds)$replicate <- unlist(sapply(colData(sciPlex_cds)$culture_well, function(x) { if (stringr::str_sub(x, -2) %in% c("01", "04", "07", "10")) return("rep1") if (stringr::str_sub(x, -2) %in% c("02", "05", "08", "11")) return("rep2") if (stringr::str_sub(x, -2) %in% c("03", "06", "09", "12")) return("rep3") else return(NA) })) ##################################
isolate HDAC cells
cds = sciPlex_cds[,colData(sciPlex_cds)$treatment == "SAHA"] cds = detect_genes(cds) cds = estimate_size_factors(cds)
cds <- preprocess_cds(cds, method = 'PCA', num_dim = 40, norm_method = 'log', verbose = T)
################################### cds = reduce_dimension(cds, max_components = 2, reduction_method = 'UMAP', umap.metric = 'cosine', umap.n_neighbors = 10, umap.min_dist = 0.1, verbose = TRUE)
cds <- cluster_cells(cds,reduction_method = "PCA")
colData(cds)$Cluster = clusters(cds, reduction_method="PCA")
cds <- cluster_cells(cds,reduction_method = "UMAP") colData(cds)$louvain_component = cds@clusters[["UMAP"]]$partitions
graph_parameters = list() graph_parameters[["minimal_branch_len"]] = 10 graph_parameters[["ncenter"]] = 750
cds <- learn_graph(cds, learn_graph_control= graph_parameters, use_partition = FALSE, close_loop = FALSE) colData(cds)$UMAP1 = cds@reducedDims[["UMAP"]][,1] colData(cds)$UMAP2 = cds@reducedDims[["UMAP"]][,2]
Get the closest vertex for every cell
colData(cds)$closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex[,1]
ordering_summary = colData(cds) %>% as.data.frame() %>% dplyr::group_by(closest_vertex) %>% dplyr::count(dose) %>% dplyr::mutate(total_cells = sum(n), fraction_dose = n / total_cells)
root_nodes =
ordering_summary %>% filter(dose == 0 & fraction_dose > 0.5)
root_nodes = root_nodes$closest_vertex
colData(cds)$root_node = colData(cds)$closest_vertex %in% root_nodes
root_cells = colData(cds) %>% as.data.frame() %>% filter(root_node) %>% pull(cell) %>% as.character()
cds <- order_cells(cds, root_cells=root_cells)
colData(cds)$Pseudotime = cds@principal_graph_aux[["UMAP"]]$pseudotime
plot_cells(cds, color_cells_by = "Pseudotime")
################################################
traceback() No errors produced
Expected behavior Pseudotime spans the full graph (not constrained to one partition)
Screenshots
sessionInfo():
R version 3.6.1 (2019-07-05) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Catalina 10.15.3
Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
Random number generation: RNG: Mersenne-Twister Normal: Inversion Sample: Rounding
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 parallel stats graphics grDevices utils datasets methods base
other attached packages: [1] snowfall_1.84-6.1 snow_0.4-3 UpSetR_1.4.0
[4] piano_2.0.2 viridis_0.5.1 viridisLite_0.3.0
[7] furrr_0.1.0 future_1.15.1 monocle3_0.2.0
[10] SingleCellExperiment_1.6.0 SummarizedExperiment_1.14.1 DelayedArray_0.10.0
[13] BiocParallel_1.18.0 matrixStats_0.55.0 GenomicRanges_1.36.0
[16] GenomeInfoDb_1.20.0 IRanges_2.18.1 S4Vectors_0.22.1
[19] Biobase_2.44.0 BiocGenerics_0.30.0 devtools_2.2.1
[22] usethis_1.5.1 yardstick_0.0.6 workflows_0.1.1
[25] tune_0.1.0 tibble_2.1.3 rsample_0.0.5
[28] tidyr_0.8.3 recipes_0.1.10 purrr_0.3.3
[31] parsnip_0.1.0 infer_0.5.1 ggplot2_3.3.0
[34] dplyr_0.8.5 dials_0.0.6 scales_1.0.0
[37] broom_0.5.5 tidymodels_0.1.0
loaded via a namespace (and not attached): [1] shinydashboard_0.7.1 tidyselect_0.2.5 lme4_1.1-23
[4] htmlwidgets_1.3 grid_3.6.1 pROC_1.16.2
[7] munsell_0.5.0 codetools_0.2-16 statmod_1.4.34
[10] DT_0.9 miniUI_0.1.1.1 withr_2.1.2
[13] colorspace_1.4-1 knitr_1.23 rstudioapi_0.11
[16] bayesplot_1.7.1 listenv_0.7.0 labeling_0.3
[19] rstan_2.19.3 slam_0.1-45 GenomeInfoDbData_1.2.1
[22] DiceDesign_1.8-1 rprojroot_1.3-2 vctrs_0.2.4
[25] generics_0.0.2 ipred_0.9-9 xfun_0.8
[28] sets_1.0-18 R6_2.4.0 markdown_1.0
[31] rstanarm_2.19.3 bitops_1.0-6 lhs_1.0.1
[34] fgsea_1.10.1 assertthat_0.2.1 promises_1.0.1
[37] nnet_7.3-12 gtable_0.3.0 globals_0.12.5
[40] processx_3.4.1 timeDate_3043.102 rlang_0.4.5
[43] splines_3.6.1 inline_0.3.15 reshape2_1.4.3
[46] tidytext_0.2.3 threejs_0.3.3 crosstalk_1.0.0
[49] backports_1.1.4 httpuv_1.5.1 rsconnect_0.8.16
[52] tokenizers_0.2.1 tools_3.6.1 lava_1.6.7
[55] relations_0.6-9 ellipsis_0.3.0 gplots_3.0.3
[58] proxy_0.4-23 sessioninfo_1.1.1 ggridges_0.5.1
[61] Rcpp_1.0.1 plyr_1.8.4 visNetwork_2.0.9
[64] base64enc_0.1-3 zlibbioc_1.30.0 RCurl_1.95-4.12
[67] ps_1.3.0 prettyunits_1.0.2 rpart_4.1-15
[70] zoo_1.8-6 ggrepel_0.8.1 cluster_2.1.0
[73] fs_1.3.1 magrittr_1.5 RSpectra_0.15-0
[76] data.table_1.12.2 RANN_2.6.1 colourpicker_1.0
[79] GPfit_1.0-8 SnowballC_0.7.0 pkgload_1.0.2
[82] tidyposterior_0.0.2 shinyjs_1.1 mime_0.7
[85] xtable_1.8-4 tidypredict_0.4.5 shinystan_2.5.0
[88] gridExtra_2.3 rstantools_2.0.0 testthat_2.1.1
[91] compiler_3.6.1 KernSmooth_2.23-15 crayon_1.3.4
[94] minqa_1.2.4 StanHeaders_2.21.0-1 htmltools_0.3.6
[97] later_0.8.0 RcppParallel_4.4.3 lubridate_1.7.8
[100] MASS_7.3-51.4 boot_1.3-23 leidenbase_0.1.0
[103] Matrix_1.2-17 cli_2.0.2 marray_1.62.0
[106] gdata_2.18.0 gower_0.2.1 igraph_1.2.4.1
[109] pkgconfig_2.0.3 foreach_1.5.0 dygraphs_1.1.1.6
[112] XVector_0.24.0 prodlim_2019.11.13 janeaustenr_0.1.5
[115] stringr_1.4.0 callr_3.3.1 digest_0.6.20
[118] RcppAnnoy_0.0.12 fastmatch_1.1-0 uwot_0.1.5
[121] DelayedMatrixStats_1.6.1 shiny_1.3.2 gtools_3.8.1
[124] nloptr_1.2.2.1 jsonlite_1.6 lifecycle_0.1.0
[127] nlme_3.1-140 desc_1.2.0 limma_3.40.6
[130] fansi_0.4.0 pillar_1.4.3 lattice_0.20-38
[133] loo_2.2.0 pkgbuild_1.0.3 survival_2.44-1.1
[136] glue_1.3.1 xts_0.12-0 remotes_2.1.0
[139] shinythemes_1.1.2 iterators_1.0.12 class_7.3-15
[142] stringi_1.4.3 caTools_1.18.0 memoise_1.1.0
[145] irlba_2.3.3
Additional context The CDS used can be downloaded here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4150377