Skip to content

Commit

Permalink
Merge pull request #114 from jr-leary7/dev
Browse files Browse the repository at this point in the history 
Dev
  • Loading branch information
@jr-leary7
jr-leary7 authored Sep 5, 2023
2 parents 05dbef3 + f64e675 commit 6664fe8
Show file tree
Hide file tree
Showing 6 changed files with 46 additions and 24 deletions.
2 changes: 1 addition & 1 deletion R/GetResultsDE.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
52 changes: 34 additions & 18 deletions R/getFittedValues.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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{
Expand All @@ -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,
Expand All @@ -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)
Expand All @@ -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)
Expand Down
6 changes: 4 additions & 2 deletions R/plotModels.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
}
Expand Down Expand Up @@ -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)
}
Expand Down
5 changes: 4 additions & 1 deletion man/getFittedValues.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion man/getResultsDE.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

3 changes: 2 additions & 1 deletion tests/testthat/test_scLANE.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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", {
Expand Down

0 comments on commit 6664fe8

Please sign in to comment.