Merge pull request #86 from LieberInstitute/registration_wrapper_k2

Registration wrapper k2

Author: lcolladotor

R/registration_wrapper.R

@@ -47,7 +47,12 @@
 #' ## Compute all modeling results
 #' example_modeling_results <- registration_wrapper(
 #'     sce,
-#'     "Cell_Cycle", "sample_id", c("age"), "ensembl", "gene_name", "wrapper"
+#'     var_regustration ="Cell_Cycle", 
+#'     var_sample_id ="sample_id", 
+#'     covars = c("age"),
+#'     gene_ensembl = "ensembl", 
+#'     gene_name = "gene_name", 
+#'     suffix = "wrapper"
 #' )
 registration_wrapper <-
     function(sce,
@@ -76,6 +81,12 @@ registration_wrapper <-
 
         block_cor <-
             registration_block_cor(sce_pseudo, registration_model = registration_mod)
+        
+        ## test if registration var has two groups
+        registration_var_k2 <- length(grep("^registration_variable", colnames(registration_mod))) == 2
+        if (registration_var_k2) {
+          warning("You need 'var_registration' to have at least 3 different values to compute an F-statistic, returning Enrichment statistics only", call. = FALSE)
+        }
 
         results_enrichment <-
             registration_stats_enrichment(
@@ -85,29 +96,40 @@ registration_wrapper <-
                 gene_ensembl = gene_ensembl,
                 gene_name = gene_name
             )
-        results_pairwise <-
+        
+        ## with more than 2 groups run ANOVA and pairwise data
+        if(!registration_var_k2){
+          results_pairwise <-
             registration_stats_pairwise(
-                sce_pseudo,
-                registration_model = registration_mod,
-                block_cor = block_cor,
-                gene_ensembl = gene_ensembl,
-                gene_name = gene_name
+              sce_pseudo,
+              registration_model = registration_mod,
+              block_cor = block_cor,
+              gene_ensembl = gene_ensembl,
+              gene_name = gene_name
             )
-        results_anova <-
+          results_anova <-
             registration_stats_anova(
-                sce_pseudo,
-                block_cor = block_cor,
-                covars = covars,
-                gene_ensembl = gene_ensembl,
-                gene_name = gene_name,
-                suffix = suffix
+              sce_pseudo,
+              block_cor = block_cor,
+              covars = covars,
+              gene_ensembl = gene_ensembl,
+              gene_name = gene_name,
+              suffix = suffix
             )
-
-        modeling_results <- list(
+          
+          modeling_results <- list(
             "anova" = results_anova,
             "enrichment" = results_enrichment,
             "pairwise" = results_pairwise
-        )
+          )
+        } else {
+          modeling_results <- list(
+            "anova" = NULL,
+            "enrichment" = results_enrichment,
+            "pairwise" = NULL
+          )
+        }
+        
 
         return(modeling_results)
     }

man/registration_wrapper.Rd

@@ -0,0 +1,35 @@
+## Ensure reproducibility of example data
+set.seed(20220907)
+
+## Generate example data
+sce <- scuttle::mockSCE()
+
+## Add some sample IDs
+sce$sample_id <- sample(LETTERS[1:5], ncol(sce), replace = TRUE)
+
+## Add a sample-level covariate: age
+ages <- rnorm(5, mean = 20, sd = 4)
+names(ages) <- LETTERS[1:5]
+sce$age <- ages[sce$sample_id]
+
+## add variable with one group
+sce$batch <- "batch1"
+
+## Add gene-level information
+rowData(sce)$ensembl <- paste0("ENSG", seq_len(nrow(sce)))
+rowData(sce)$gene_name <- paste0("gene", seq_len(nrow(sce)))
+
+
+test_that("warning for k=2 variable",
+          example_modeling_results <- expect_warning(
+            registration_wrapper(
+              sce,
+              var_registration ="Treatment",
+              var_sample_id ="sample_id",
+              covars = c("age"),
+              gene_ensembl = "ensembl",
+              gene_name = "gene_name",
+              suffix = "wrapper"
+            )
+          )
+)