From 479e0837d7731d615d919a6a294de07cee5c05be Mon Sep 17 00:00:00 2001 From: jr-leary7 Date: Fri, 1 Sep 2023 17:33:25 -0400 Subject: [PATCH 1/2] fixed facet strip widths in plotModels --- R/plotModels.R | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/R/plotModels.R b/R/plotModels.R index b2bf246..a4a667d 100644 --- a/R/plotModels.R +++ b/R/plotModels.R @@ -277,7 +277,8 @@ plotModels <- function(test.dyn.res = NULL, if (requireNamespace("ggh4x", quietly = TRUE)) { p <- p + ggh4x::facet_nested_wrap(~paste0("Lineage ", LINEAGE) + MODEL + ID, nrow = length(levels(counts_df$MODEL)), - strip = ggh4x::strip_nested(background_x = list(ggplot2::element_rect(linewidth = gg.theme$line$linewidth)))) + strip = ggh4x::strip_nested(clip = "off", + background_x = list(ggplot2::element_rect(linewidth = gg.theme$line$linewidth)))) } else { p <- p + ggplot2::facet_wrap(~paste0("Lineage ", LINEAGE) + MODEL + ID) } @@ -310,7 +311,8 @@ plotModels <- function(test.dyn.res = NULL, show.legend = ifelse(ncol(pt) > 1, TRUE, FALSE)) if (requireNamespace("ggh4x", quietly = TRUE)) { p <- p + ggh4x::facet_nested_wrap(~paste0("Lineage ", LINEAGE) + MODEL, - strip = ggh4x::strip_nested(background_x = list(ggplot2::element_rect(linewidth = gg.theme$line$linewidth)))) + strip = ggh4x::strip_nested(clip = "off", + background_x = list(ggplot2::element_rect(linewidth = gg.theme$line$linewidth)))) } else { p <- p + ggplot2::facet_wrap(~paste0("Lineage ", LINEAGE) + MODEL) } From f64e675557efff8b20386c3a1d2c4f7b8d5446ee Mon Sep 17 00:00:00 2001 From: jr-leary7 Date: Tue, 5 Sep 2023 14:23:54 -0400 Subject: [PATCH 2/2] critical changes to getFittedValues function for pub figures --- R/GetResultsDE.R | 2 +- R/getFittedValues.R | 52 +++++++++++++++++++++++------------- man/getFittedValues.Rd | 5 +++- man/getResultsDE.Rd | 2 +- tests/testthat/test_scLANE.R | 3 ++- 5 files changed, 42 insertions(+), 22 deletions(-) diff --git a/R/GetResultsDE.R b/R/GetResultsDE.R index 63f5314..c23265b 100644 --- a/R/GetResultsDE.R +++ b/R/GetResultsDE.R @@ -38,7 +38,7 @@ getResultsDE <- function(test.dyn.res = NULL, dplyr::mutate(dplyr::across(tidyselect::everything(), \(z) unname(unlist(z)))) }) }) %>% - dplyr::arrange(P_Val) %>% + dplyr::arrange(P_Val, Test_Stat) %>% dplyr::mutate(P_Val_Adj = stats::p.adjust(P_Val, method = p.adj.method), Gene_Dynamic_Lineage = dplyr::if_else(P_Val_Adj < fdr.cutoff, 1, 0, missing = 0)) %>% dplyr::with_groups(Gene, diff --git a/R/getFittedValues.R b/R/getFittedValues.R index 915b257..7d59c91 100644 --- a/R/getFittedValues.R +++ b/R/getFittedValues.R @@ -12,11 +12,12 @@ #' @param pt A data.frame of pseudotime values for each cell. Defaults to NULL. #' @param expr.mat Either a \code{SingleCellExperiment} or \code{Seurat} object from which counts can be extracted, or a matrix of integer-valued counts with genes as rows & cells as columns. Defaults to NULL. #' @param size.factor.offset (Optional) An offset to be used to rescale the fitted values. Can be generated easily with \code{\link{createCellOffset}}. No need to provide if the GEE backend was used. Defaults to NULL. +#' @param log1p.norm (Optional) Should log1p-normalized versions of expression & model predictions be returned as well? Defaults to TRUE. #' @param cell.meta.data (Optional) A data.frame of metadata values for each cell (celltype label, subject characteristics, tissue type, etc.) that will be included in the result table. Defaults to NULL. #' @param id.vec (Optional) A vector of subject IDs used in fitting GEE or GLMM models. Defaults to NULL. #' @param ci.alpha (Optional) The pre-specified Type I Error rate used in generating (\eqn{1 - \alpha})\% CIs. Defaults to good old 0.05. #' @param filter.lineage (Optional) A character vector of lineages to filter out before generating the final plot. Should be letters, i.e. lineage "A" or "B". Defaults to NULL. -#' @return A data.frame containing expression, fitted values, and metadata. +#' @return A data.frame containing depth- and log1p-normalized expression, model predictions, and cell-level metadata. #' @export #' @examples #' \dontrun{ @@ -33,6 +34,7 @@ getFittedValues <- function(test.dyn.res = NULL, pt = NULL, expr.mat = NULL, size.factor.offset = NULL, + log1p.norm = TRUE, cell.meta.data = NULL, id.vec = NULL, ci.alpha = 0.05, @@ -41,22 +43,20 @@ getFittedValues <- function(test.dyn.res = NULL, if (is.null(expr.mat) || is.null(pt) || is.null(genes) || is.null(test.dyn.res)) { stop("You forgot one or more of the arguments to getFittedValues().") } # get raw counts from SingleCellExperiment or Seurat object & transpose to cell x gene dense matrix if (inherits(expr.mat, "SingleCellExperiment")) { - expr.mat <- BiocGenerics::counts(expr.mat)[genes, ] + expr.mat <- BiocGenerics::counts(expr.mat)[genes, , drop = FALSE] expr.mat <- as.matrix(expr.mat) } else if (inherits(expr.mat, "Seurat")) { expr.mat <- Seurat::GetAssayData(expr.mat, slot = "counts", assay = Seurat::DefaultAssay(expr.mat)) - expr.mat <- as.matrix(expr.mat[genes, ]) - } else if (inherits(expr.mat, "dgCMatrix")) { - expr.mat <- as.matrix(expr.mat) + expr.mat <- as.matrix(expr.mat[genes, , drop = FALSE]) + } else if (inherits(expr.mat, "dgCMatrix") || inherits(expr.mat, "dgTMatrix")) { + expr.mat <- as.matrix(expr.mat[genes, , drop = FALSE]) } if (!(inherits(expr.mat, "matrix") || inherits(expr.mat, "array"))) { stop("Input expr.mat must be coerceable to a matrix of integer counts.") } expr.mat <- t(expr.mat) # transpose to cell x gene matrix if (is.null(genes)) { genes <- colnames(expr.mat) - } else { - expr.mat <- expr.mat[, genes] } # generate parameters for CIs Z <- stats::qnorm(ci.alpha / 2, lower.tail = FALSE) @@ -69,31 +69,47 @@ getFittedValues <- function(test.dyn.res = NULL, mod_df_list <- purrr::map2(pt, lineages, \(x, y) { mod_df <- purrr::map(genes, \(g) { - gene_df <- data.frame(cell = rownames(pt)[!is.na(x)], + gene_df <- data.frame(cell = rownames(expr.mat)[!is.na(x)], lineage = y, pt = x[!is.na(x)], gene = g, - expression = expr.mat[!is.na(x), g]) + rna = expr.mat[!is.na(x), g]) + if (!is.null(size.factor.offset)) { + gene_df <- dplyr::mutate(gene_df, + size_factor = unname(size.factor.offset)[!is.na(x)], + model_offset = log(1 / unname(size.factor.offset)[!is.na(x)])) + } pred_df <- try({ - data.frame(scLANE_fitted_link = test.dyn.res[[g]][[paste0("Lineage_", y)]]$MARGE_Preds$marge_link_fit, + data.frame(scLANE_pred_link = test.dyn.res[[g]][[paste0("Lineage_", y)]]$MARGE_Preds$marge_link_fit, scLANE_se_link = test.dyn.res[[g]][[paste0("Lineage_", y)]]$MARGE_Preds$marge_link_se) }, silent = TRUE) if (inherits(pred_df, "try-error")) { gene_df <- dplyr::mutate(gene_df, - scLANE_fitted_link = NA_real_, - scLANE_se_link = NA_real_) + scLANE_pred_link = NA_real_, + scLANE_se_link = NA_real_, + scLANE_ci_ll_link = NA_real_, + scLANE_ci_ul_link = NA_real_) } else { gene_df <- dplyr::mutate(gene_df, - scLANE_fitted_link = pred_df$scLANE_fitted_link, - scLANE_se_link = pred_df$scLANE_se_link) + scLANE_pred_link = pred_df$scLANE_pred_link, + scLANE_se_link = pred_df$scLANE_se_link, + scLANE_ci_ll_link = scLANE_pred_link - Z * scLANE_se_link, + scLANE_ci_ul_link = scLANE_pred_link + Z * scLANE_se_link) } gene_df <- dplyr::mutate(gene_df, - scLANE_fitted = exp(scLANE_fitted_link), - scLANE_ci_ll = exp(scLANE_fitted_link - Z * scLANE_se_link), - scLANE_ci_ul = exp(scLANE_fitted_link + Z * scLANE_se_link)) + scLANE_pred = exp(scLANE_pred_link), + scLANE_ci_ll = exp(scLANE_pred_link - Z * scLANE_se_link), + scLANE_ci_ul = exp(scLANE_pred_link + Z * scLANE_se_link)) if (!is.null(size.factor.offset)) { gene_df <- dplyr::mutate(gene_df, - dplyr::across(c(scLANE_fitted, scLANE_ci_ll, scLANE_ci_ul), \(m) m * unname(size.factor.offset)[!is.na(x)])) + dplyr::across(c(rna, scLANE_pred, scLANE_ci_ll, scLANE_ci_ul), \(m) m * size_factor)) + } + if (log1p.norm) { + gene_df <- dplyr::mutate(gene_df, + expression_log1p = log1p(rna), + scLANE_pred_log1p = log1p(scLANE_pred), + scLANE_ci_ll_log1p = log1p(scLANE_ci_ll), + scLANE_ci_ul_log1p = log1p(scLANE_ci_ul)) } gene_df <- dplyr::bind_cols(gene_df, cell.meta.data[!is.na(x), , drop = FALSE]) return(gene_df) diff --git a/man/getFittedValues.Rd b/man/getFittedValues.Rd index 6726a1a..f1aba91 100644 --- a/man/getFittedValues.Rd +++ b/man/getFittedValues.Rd @@ -10,6 +10,7 @@ getFittedValues( pt = NULL, expr.mat = NULL, size.factor.offset = NULL, + log1p.norm = TRUE, cell.meta.data = NULL, id.vec = NULL, ci.alpha = 0.05, @@ -27,6 +28,8 @@ getFittedValues( \item{size.factor.offset}{(Optional) An offset to be used to rescale the fitted values. Can be generated easily with \code{\link{createCellOffset}}. No need to provide if the GEE backend was used. Defaults to NULL.} +\item{log1p.norm}{(Optional) Should log1p-normalized versions of expression & model predictions be returned as well? Defaults to TRUE.} + \item{cell.meta.data}{(Optional) A data.frame of metadata values for each cell (celltype label, subject characteristics, tissue type, etc.) that will be included in the result table. Defaults to NULL.} \item{id.vec}{(Optional) A vector of subject IDs used in fitting GEE or GLMM models. Defaults to NULL.} @@ -36,7 +39,7 @@ getFittedValues( \item{filter.lineage}{(Optional) A character vector of lineages to filter out before generating the final plot. Should be letters, i.e. lineage "A" or "B". Defaults to NULL.} } \value{ -A data.frame containing expression, fitted values, and metadata. +A data.frame containing depth- and log1p-normalized expression, model predictions, and cell-level metadata. } \description{ Generate a table of expression counts, model fitted values, celltype metadata, etc. in order to create custom plots of gene dynamics. diff --git a/man/getResultsDE.Rd b/man/getResultsDE.Rd index a3d12e4..64b6d98 100644 --- a/man/getResultsDE.Rd +++ b/man/getResultsDE.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/GetResultsDE.R +% Please edit documentation in R/getResultsDE.R \name{getResultsDE} \alias{getResultsDE} \title{Tidy the results of \code{\link{testDynamic}}.} diff --git a/tests/testthat/test_scLANE.R b/tests/testthat/test_scLANE.R index 9cc02ff..2f06a44 100644 --- a/tests/testthat/test_scLANE.R +++ b/tests/testthat/test_scLANE.R @@ -168,6 +168,7 @@ withr::with_output_sink(tempfile(), { fitted_values_table <- getFittedValues(test.dyn.res = glm_gene_stats, genes = names(glm_gene_stats), pt = pt_test, + size.factor.offset = cell_offset, expr.mat = sim_data, cell.meta.data = as.data.frame(SummarizedExperiment::colData(sim_data)), id.vec = sim_data$subject) @@ -312,7 +313,7 @@ test_that("sortGenesHeatmap() output", { test_that("getFittedValues() output", { expect_s3_class(fitted_values_table, "data.frame") - expect_equal(ncol(fitted_values_table), 17) + expect_equal(ncol(fitted_values_table), 25) }) test_that("enrichDynamicGenes() output", {