diff --git a/DESCRIPTION b/DESCRIPTION index 121a3ee9..8726b357 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: tidybulk Title: Brings transcriptomics to the tidyverse -Version: 1.17.3 +Version: 1.17.4 Authors@R: c(person("Stefano", "Mangiola", email = "mangiolastefano@gmail.com", role = c("aut", "cre")), person("Maria", "Doyle", email = "Maria.Doyle@petermac.org", diff --git a/R/functions_SE.R b/R/functions_SE.R index 23c9c7e5..14b177e4 100755 --- a/R/functions_SE.R +++ b/R/functions_SE.R @@ -794,7 +794,7 @@ get_differential_transcript_abundance_bulk_SE <- function(.data, "with `___` in the design matrix, in order to work with edgeR.") colnames(design) = design |> colnames() |> str_replace(":", "___") } - + # Print the design column names in case I want contrasts message( sprintf( @@ -1137,6 +1137,7 @@ get_differential_transcript_abundance_bulk_voom_SE <- function(.data, #' "edgeR_likelihood_ratio" (i.e., LRT) #' @param scaling_method A character string. The scaling method passed to the #' backend function (i.e., edgeR::calcNormFactors; "TMM","TMMwsp","RLE","upperquartile") +#' @param .scaling_factor A tidyeval (column name) for the precalculated TMM scaling #' @param omit_contrast_in_colnames If just one contrast is specified you can #' choose to omit the contrast label in the colnames. #' @param ... Additional arguments for glmmSeq @@ -1151,15 +1152,16 @@ get_differential_transcript_abundance_glmmSeq_SE <- function(.data, method, test_above_log2_fold_change = NULL, - - scaling_method = "TMM", - omit_contrast_in_colnames = FALSE, - prefix = "", - .dispersion = NULL, - ...) { + scaling_method = "TMM", + .scaling_factor = NULL, + omit_contrast_in_colnames = FALSE, + prefix = "", + .dispersion = NULL, + ...) { .abundance = enquo(.abundance) .dispersion = enquo(.dispersion) + .scaling_factor = enquo(.scaling_factor) # Check if contrasts are of the same form if( @@ -1233,7 +1235,10 @@ get_differential_transcript_abundance_glmmSeq_SE <- function(.data, dispersion = dispersion[rownames(counts)] # Scaling - sizeFactors <- counts |> edgeR::calcNormFactors(method = scaling_method) + if(.scaling_factor |> quo_is_symbolic()) + sizeFactors = .data |> pivot_sample() |> pull(!!.scaling_factor) + else + sizeFactors <- counts |> edgeR::calcNormFactors(method = scaling_method) glmmSeq_object = diff --git a/R/methods.R b/R/methods.R index 63be8bdc..7508d888 100755 --- a/R/methods.R +++ b/R/methods.R @@ -627,14 +627,14 @@ setMethod("scale_abundance", "tidybulk", .scale_abundance) #' The scaling inference is then applied back to all unfiltered data. #' #' Underlying method -#' +#' #' If `limma_normalize_quantiles` is chosen -#' +#' #' .data |>limma::normalizeQuantiles() -#' +#' #' If `preprocesscore_normalize_quantiles_use_target` is chosen -#' -#' .data |> +#' +#' .data |> #' preprocessCore::normalize.quantiles.use.target( #' target = preprocessCore::normalize.quantiles.determine.target(.data) #' ) @@ -675,6 +675,7 @@ setGeneric("quantile_normalise_abundance", function( target_distribution = NULL, action = "add") { + # Fix NOTEs . = NULL @@ -727,7 +728,7 @@ setGeneric("quantile_normalise_abundance", function( BiocManager::install("preprocessCore", ask = FALSE) } if(is.null(target_distribution)) target_distribution = preprocessCore::normalize.quantiles.determine.target(.data_norm) - + .data_norm_quant = .data_norm |> preprocessCore::normalize.quantiles.use.target( diff --git a/R/methods_SE.R b/R/methods_SE.R index 9eaf4c8c..cba366f0 100755 --- a/R/methods_SE.R +++ b/R/methods_SE.R @@ -186,17 +186,16 @@ setMethod("tidybulk", "RangedSummarizedExperiment", .tidybulk_se) my_counts_scaled = list( - assays(.data) %>% - as.list() %>% - .[[1]] %>% - multiply_by( - rep(multiplier, rep(nrow(.),length(multiplier))) - ) + assay(.data) %*% + diag(multiplier) + ) %>% setNames(value_scaled) + colnames(my_counts_scaled[[1]]) = assay(.data) |> colnames() + # Add the assay - assays(.data) = assays(.data) %>% c(my_counts_scaled) + assays(.data, withDimnames=FALSE) = assays(.data) %>% c(my_counts_scaled) .data %>% @@ -309,7 +308,7 @@ setMethod("scale_abundance", as.matrix() if(is.null(target_distribution)) target_distribution = preprocessCore::normalize.quantiles.determine.target(.data_norm) - + .data_norm = .data_norm |> preprocessCore::normalize.quantiles.use.target( @@ -2922,6 +2921,8 @@ setMethod("describe_transcript", "RangedSummarizedExperiment", .describe_transcr combination_of_factors_of_NON_interest = # Factors se[1,1, drop=FALSE] |> + colData() |> + as_tibble(rownames = ".sample") |> select(...) |> suppressWarnings() |> colnames() |> diff --git a/man/quantile_normalise_abundance-methods.Rd b/man/quantile_normalise_abundance-methods.Rd index e46176ef..b5f04566 100644 --- a/man/quantile_normalise_abundance-methods.Rd +++ b/man/quantile_normalise_abundance-methods.Rd @@ -141,10 +141,10 @@ Underlying method If `limma_normalize_quantiles` is chosen .data |>limma::normalizeQuantiles() - + If `preprocesscore_normalize_quantiles_use_target` is chosen - -.data |> + +.data |> preprocessCore::normalize.quantiles.use.target( target = preprocessCore::normalize.quantiles.determine.target(.data) ) diff --git a/man/test_differential_abundance-methods.Rd b/man/test_differential_abundance-methods.Rd index cdec9119..21bd0246 100755 --- a/man/test_differential_abundance-methods.Rd +++ b/man/test_differential_abundance-methods.Rd @@ -146,13 +146,9 @@ of the method of choice. For edgeR and limma-voom is a character vector. For DESeq2 is a list including a character vector of length three. The first covariate is the one the model is tested against (e.g., ~ factor_of_interest)} -<<<<<<< HEAD -\item{method}{A string character. Either "edgeR_quasi_likelihood" (i.e., QLF), "edgeR_likelihood_ratio" (i.e., LRT), "edger_robust_likelihood_ratio", "DESeq2", "limma_voom", "limma_voom_sample_weights", "glmmseq_lme4", "glmmseq_glmmtmb"} -======= \item{method}{A string character. Either "edgeR_quasi_likelihood" (i.e., QLF), "edgeR_likelihood_ratio" (i.e., LRT), "edger_robust_likelihood_ratio", -"DESeq2", "limma_voom", "limma_voom_sample_weights"} ->>>>>>> chilampoon-improve-documentation +"DESeq2", "limma_voom", "limma_voom_sample_weights", "glmmseq_lme4", "glmmseq_glmmtmb"} \item{test_above_log2_fold_change}{A positive real value. This works for edgeR and limma_voom methods. It uses the `treat` function, which tests that the @@ -173,8 +169,11 @@ result columns. It is useful if you want to compare several methods.} \item{action}{A character string. Whether to join the new information to the input tbl (add), or just get the non-redundant tbl with the new information (get).} -\item{...}{Further arguments passed to some of the internal functions. -Currently, it is needed just for internal debug.} +\item{...}{Further arguments passed to some of the internal experimental functions. +For example for glmmSeq, it is possible to pass .dispersion, and .scaling_factor +column tidyeval to skip the caluclation of dispersion and scaling and use precalculated values. +This is helpful is you want to calculate those quantities on many genes and do DE testing on fewer genes. +.scaling_factor is the TMM value that can be obtained with tidybulk::scale_abundance.} \item{significance_threshold}{DEPRECATED - A real between 0 and 1 (usually 0.05).} diff --git a/tests/testthat/test-bulk_methods_SummarizedExperiment.R b/tests/testthat/test-bulk_methods_SummarizedExperiment.R index 84bc5a8f..3aae6428 100755 --- a/tests/testthat/test-bulk_methods_SummarizedExperiment.R +++ b/tests/testthat/test-bulk_methods_SummarizedExperiment.R @@ -747,7 +747,7 @@ test_that("gene over representation",{ species="Homo sapiens" ) - expect_equal( ncol(res), 10 ) + expect_equal( ncol(res), 13 ) @@ -854,8 +854,8 @@ test_that("Only reduced dimensions UMAP - no object",{ test_that("resolve_complete_confounders_of_non_interest",{ - library(tidySummarizedExperiment) - library(tidybulk) + #library(tidySummarizedExperiment) + library(SummarizedExperiment) # Sample annotations sample_annotations <- data.frame( @@ -890,7 +890,9 @@ test_that("resolve_complete_confounders_of_non_interest",{ se |> resolve_complete_confounders_of_non_interest(A, B, C) |> - distinct(.sample, A, B, C) |> + colData() |> + _[, c("A", "B", "C")] |> + as_tibble(rownames = ".sample") |> expect_identical(expected_tibble )