Merge pull request #51 from stemangiola/master

rebase
MangiolaLaboratory · Jun 12, 2022 · 1ba467f · 1ba467f
2 parents 766349d + e112ed3
commit 1ba467f
Show file tree

Hide file tree

Showing 26 changed files with 1,021 additions and 806 deletions.
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -1,11 +1,15 @@
 ^sccomp\.Rproj$
 ^\.Rproj\.user$
 ^\.github
+^\.git
 ^README*
 ^dev$
-^dev*
 ^dev/*
 ^rsconnect
 model_glm_dirichlet_multinomial.rds
 model_glm_dirichlet_multinomial_generate_quantities.rds
 model_glm_dirichlet_multinomial_imputation.rds
+^sccode_env
+^sccomp_draws_files
+^glm_multi*
+^model_glm*
diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -0,0 +1,12 @@
+Thanks for submitting an issue.
+
+Please add the following information to the issue
+
+1. Describe the issue/bug
+2. Print out the input dataset immediately before the bug occurs
+3. Paste the code immediately leading to the bug
+4. Print out of the output, if any
+5. Print out of the complete error/warning message, if any
+6. sessionInfo()
+
+Thanks!
diff --git a/.github/workflows/check-bioc.yml b/.github/workflows/check-bioc.yml
diff --git a/.gitignore b/.gitignore
@@ -7,19 +7,16 @@ vignettes/introduction_cache*
 sccomp.Rproj
 Meta
 doc
-
 sccode_env*
 src/stanExports*
 slurm*
 makeflow*
-
 dmbvs-pid-50624.out
 dmbvs-pid-39940.out
 model_glm_dirichlet_multinomial.rds
 model_glm_dirichlet_multinomial_generate_quantities.rds
 model_glm_dirichlet_multinomial_imputation.rds
 rsconnect
-
 dev/*csv
 dev/*rds
 dev/*rda
@@ -37,3 +34,9 @@ dev/benchmark*
 dev/*.gz
 dev/*csv
 dev/reproducible_example_beta_binomial_truncated
+glm_multi_beta_binomial_cmdstanr
+glm_multi_beta_binomial_cmdstanr.rds
+glm_multi_beta_binomial_cmdstanr.stan
+glm_multi_beta_binomial_generate_cmdstanr
+glm_multi_beta_binomial_generate_cmdstanr.rds
+glm_multi_beta_binomial_generate_cmdstanr.stan
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: sccomp
 Title: Robust Outlier-aware Estimation of Composition and Heterogeneity for Single-cell Data
-Version: 0.99.18
+Version: 0.99.20
 Authors@R: c(person("Stefano", "Mangiola", email = "[email protected]",
                   role = c("aut", "cre"))
                   )
@@ -12,13 +12,13 @@ Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.1.2
 Biarch: true
 Depends: 
-    R (>= 4.1.0),
-    rstan (>= 2.18.1)
+    R (>= 4.1.0)
 Imports: 
     methods,
     Rcpp (>= 0.12.0),
     RcppParallel (>= 5.0.1),
     rstantools (>= 2.1.1),
+    rstan (>= 2.18.1),
     SeuratObject,
     SingleCellExperiment,
     parallel,
@@ -35,14 +35,17 @@ Imports:
     utils,
     ggplot2,
     ggrepel,
-    patchwork
+    patchwork,
+    forcats,
+    readr,
+    scales,
+    stringr
 Suggests:
     BiocStyle,
     testthat (>= 3.0.0),
     markdown,
     ggplot2,
     knitr,
-    forcats,
     tidyseurat,
     tidySingleCellExperiment
 Enhances:

diff --git a/NAMESPACE b/NAMESPACE
@@ -29,6 +29,7 @@ importFrom(dplyr,enquo)
 importFrom(dplyr,left_join)
 importFrom(dplyr,select)
 importFrom(dplyr,with_groups)
+importFrom(forcats,fct_reorder)
 importFrom(ggrepel,geom_text_repel)
 importFrom(magrittr,"%$%")
 importFrom(magrittr,divide_by)
@@ -47,6 +48,7 @@ importFrom(purrr,map_int)
 importFrom(purrr,map_lgl)
 importFrom(purrr,pmap)
 importFrom(purrr,when)
+importFrom(readr,read_file)
 importFrom(rlang,":=")
 importFrom(rlang,enquo)
 importFrom(rlang,ensym)
@@ -61,12 +63,15 @@ importFrom(rstan,sflist2stanfit)
 importFrom(rstan,stan_model)
 importFrom(rstan,summary)
 importFrom(rstan,vb)
+importFrom(scales,trans_new)
 importFrom(stats,C)
 importFrom(stats,model.matrix)
+importFrom(stats,quantile)
 importFrom(stats,rnorm)
 importFrom(stats,sd)
 importFrom(stats,setNames)
 importFrom(stats,terms)
+importFrom(stringr,str_replace)
 importFrom(tibble,as_tibble)
 importFrom(tibble,enframe)
 importFrom(tibble,rowid_to_column)

diff --git a/R/functions_multi_beta_binomial.R b/R/functions_multi_beta_binomial.R
@@ -57,6 +57,8 @@ estimate_multi_beta_binomial_glm = function(.data,
                                             seed = sample(1e5, 1),
                                             verbose = FALSE,
                                             exclude_priors = FALSE,
+                                            bimodal_mean_variability_association = bimodal_mean_variability_association,
+                                            use_data = use_data,
                                             max_sampling_iterations = 20000
 ) {
   # Prepare column same enquo
@@ -86,9 +88,15 @@ estimate_multi_beta_binomial_glm = function(.data,
         variance_association = variance_association,
         truncation_ajustment = 1.1,
         approximate_posterior_inference = approximate_posterior_inference == "all",
-        formula_variability = formula_variability
+        formula_variability = formula_variability,
+        bimodal_mean_variability_association = bimodal_mean_variability_association,
+        use_data = use_data
       )
 
+    # Print design matrix
+    message(sprintf("sccomp says: the composition design matrix has columns: %s", data_for_model$X %>% colnames %>% paste(collapse=", ")))
+    message(sprintf("sccomp says: the variability design matrix has columns: %s", data_for_model$Xa %>% colnames %>% paste(collapse=", ")))
+
     # Pior
     data_for_model$prior_prec_intercept = prior_mean_variable_association$intercept
     data_for_model$prior_prec_slope  = prior_mean_variable_association$slope
@@ -136,7 +144,9 @@ estimate_multi_beta_binomial_glm = function(.data,
         variance_association = FALSE,
         truncation_ajustment = 1.1,
         approximate_posterior_inference = approximate_posterior_inference %in% c("outlier_detection", "all"),
-        formula_variability = ~1
+        formula_variability = ~1,
+        bimodal_mean_variability_association = bimodal_mean_variability_association,
+        use_data = use_data
       )
 
     # Pior
@@ -206,7 +216,9 @@ estimate_multi_beta_binomial_glm = function(.data,
         variance_association = variance_association,
         truncation_ajustment = 1.1,
         approximate_posterior_inference = approximate_posterior_inference %in% c("outlier_detection", "all"),
-        formula_variability = formula_variability
+        formula_variability = formula_variability,
+        bimodal_mean_variability_association = bimodal_mean_variability_association,
+        use_data = use_data
       )
 
     # Pior
@@ -296,6 +308,10 @@ estimate_multi_beta_binomial_glm = function(.data,
 
     message("sccomp says: outlier-free model fitting - step 3/3 [ETA: ~20s]")
 
+    # Print design matrix
+    message(sprintf("sccomp says: the composition design matrix has columns: %s", data_for_model$X %>% colnames %>% paste(collapse=", ")))
+    message(sprintf("sccomp says: the variability design matrix has columns: %s", data_for_model$Xa %>% colnames %>% paste(collapse=", ")))
+
     fit3 =
       data_for_model %>%
       # Run the first discovery phase with permissive false discovery rate
@@ -464,7 +480,14 @@ hypothesis_test_multi_beta_binomial_glm = function( .sample,
       ),
       ~ (.)
     ) %>%
-    select(parameter  = C_name, everything())
+
+    # Add easy to understand covariate labels
+    left_join(
+      data_for_model$covariate_parameter_dictionary %>%
+        select(covariate, design_matrix_col),
+      by = c("C_name" = "design_matrix_col" )
+    ) %>%
+    select(parameter  = C_name, covariate, everything())
 
 }
 
@@ -526,6 +549,8 @@ multi_beta_binomial_glm = function(.data,
                                    seed = sample(1e5, 1),
                                    verbose = FALSE,
                                    exclude_priors = FALSE,
+                                   bimodal_mean_variability_association = FALSE,
+                                   use_data = TRUE,
                                    test_composition_above_logit_fold_change,
                                    max_sampling_iterations = 20000,
                                    pass_fit = TRUE
@@ -554,6 +579,8 @@ multi_beta_binomial_glm = function(.data,
       seed = seed,
       verbose = verbose,
       exclude_priors = exclude_priors,
+      bimodal_mean_variability_association = bimodal_mean_variability_association,
+      use_data = use_data,
       max_sampling_iterations = max_sampling_iterations
     )