-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #156 from jr-leary7/dev
Dev
- Loading branch information
Showing
15 changed files
with
399 additions
and
30 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,79 @@ | ||
| #' Add per-cell module scores for gene programs. | ||
| #' | ||
| #' @name geneProgramScoring | ||
| #' @author Jack Leary | ||
| #' @import magrittr | ||
| #' @importFrom Matrix Matrix | ||
| #' @importFrom ggplot2 ggplot aes geom_point geom_vline geom_ribbon geom_line scale_x_continuous labs | ||
| #' @importFrom scales label_number | ||
| #' @description This function uses \code{\link[UCell]{ScoreSignatures_UCell}} to create a per-cell module score for each of the provided gene clusters. If the | ||
| #' @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 genes A character vector of gene IDs. Defaults to NULL. | ||
| #' @param gene.clusters A factor containing the cluster assignment of each gene in \code{genes}. Defaults to NULL. | ||
| #' @param program.labels (Optional) A character vector specifying a label for each gene cluster. Defaults to NULL. | ||
| #' @param n.cores (Optional) The number of cores used under the hood in \code{\link[UCell]{ScoreSignatures_UCell}}. Defaults to 2. | ||
| #' @return Either a \code{Seurat} or \code{SingleCellExperiment} object if \code{expr.mat} is in either form, or a data.frame containing per-cell program scores if \code{expr.mat} is a matrix. | ||
| #' @export | ||
| #' @examples | ||
| #' \dontrun{ | ||
| #' geneProgramScoring(seu_obj, | ||
| #' genes = gene_embed$gene, | ||
| #' gene.clusters = gene_embed$leiden, | ||
| #' program.labels = c("cell cycle", "organogenesis")) | ||
| #' } | ||
|
|
||
| geneProgramScoring <- function(expr.mat = NULL, | ||
| genes = NULL, | ||
| gene.clusters = NULL, | ||
| program.labels = NULL, | ||
| n.cores = 2) { | ||
| # check inputs | ||
| if (is.null(expr.mat) || is.null(genes) || is.null(gene.clusters)) { stop("Arguments to geneProgramScoring() are missing.") } | ||
| if (!is.factor(gene.clusters)) { | ||
| gene.clusters <- as.factor(gene.clusters) | ||
| } | ||
| # set program labels | ||
| cluster.labels <- unique(gene.clusters) | ||
| if (is.null(program.labels)) { | ||
| program.labels <- paste0("cluster_", cluster.labels) | ||
| } else { | ||
| program.labels <- gsub(" ", "_", program.labels) | ||
| } | ||
| # set up query matrix | ||
| if (inherits(expr.mat, "SingleCellExperiment")) { | ||
| counts_matrix <- BiocGenerics::counts(expr.mat) | ||
| } else if (inherits(expr.mat, "Seurat")) { | ||
| counts_matrix <- Seurat::GetAssayData(expr.mat, | ||
| slot = "counts", | ||
| assay = Seurat::DefaultAssay(expr.mat)) | ||
| } else if (inherits(expr.mat, "matrix") || inherits(expr.mat, "array")) { | ||
| counts_matrix <- Matrix::Matrix(expr.mat, sparse = TRUE) | ||
| } | ||
| # set up feature list | ||
| program_list <- split(genes, gene.clusters) | ||
| names(program_list) <- program.labels | ||
| # run UCell | ||
| program_scores <- UCell::ScoreSignatures_UCell(counts_matrix, | ||
| features = program_list, | ||
| ncores = n.cores) | ||
| # reformat program scores depending on input format | ||
| if (inherits(expr.mat, "matrix") || inherits(expr.mat, "array") || inherits(expr.mat, "dgCMatrix")) { | ||
| colnames(program_scores) <- program.labels | ||
| } else { | ||
| for (g in seq(ncol(program_scores))) { | ||
| if (inherits(expr.mat, "SingleCellExperiment")) { | ||
| SummarizedExperiment::colData(expr.mat)[, program.labels[g]] <- program_scores[, g] | ||
| } else if (inherits(expr.mat, "Seurat")) { | ||
| expr.mat <- Seurat::AddMetaData(expr.mat, | ||
| metadata = program_scores[, g], | ||
| program.labels[g]) | ||
| } | ||
| } | ||
| } | ||
| # return results | ||
| if (inherits(expr.mat, "matrix") || inherits(expr.mat, "array") || inherits(expr.mat, "dgCMatrix")) { | ||
| return(program_scores) | ||
| } else { | ||
| return(expr.mat) | ||
| } | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,137 @@ | ||
| #' Plot gene dynamics with estimated coefficients. | ||
| #' | ||
| #' @name plotModelCoefs | ||
| #' @author Jack Leary | ||
| #' @import magrittr | ||
| #' @importFrom dplyr select mutate lag lead | ||
| #' @importFrom tidyr pivot_longer | ||
| #' @importFrom ggplot2 ggplot aes geom_point geom_vline geom_ribbon geom_line scale_x_continuous labs | ||
| #' @importFrom scales label_number | ||
| #' @description Generate a plot of gene dynamics over a single pseudotime lineage, along with a table of coefficients across pseudotime intervals. | ||
| #' @param test.dyn.res The output from \code{\link{testDynamic}}. Defaults to NULL. | ||
| #' @param gene A character specifying which gene's dynamics should be plotted. Defaults to NULL. | ||
| #' @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 lineage A character vector specifying which lineage should be plotted. Should be letters, i.e. lineage "A" or "B". Defaults to "A". | ||
| #' @return A \code{ggplot2} object displaying a gene dynamics plot & a table of coefficients across pseudotime intervals. | ||
| #' @export | ||
| #' @examples | ||
| #' \dontrun{ | ||
| #' plotModelCoefs(gene_stats, | ||
| #' gene = "BRCA2", | ||
| #' pt = pt_df, | ||
| #' expr.mat = seu_obj, | ||
| #' size.factor.offset = cell_offset) | ||
| #' } | ||
|
|
||
| plotModelCoefs <- function(test.dyn.res = NULL, | ||
| gene = NULL, | ||
| pt = NULL, | ||
| expr.mat = NULL, | ||
| size.factor.offset = NULL, | ||
| lineage = "A", | ||
| log1p.norm = TRUE) { | ||
| # check inputs | ||
| if (is.null(test.dyn.res) || is.null(gene) || is.null(pt) || is.null(expr.mat)) { stop("Arguments to plotModelCoefs() are missing.") } | ||
| # pull fitted values | ||
| all_lineages <- gsub("Lineage_", "", names(test.dyn.res[[1]])) | ||
| if (length(all_lineages) == 1) { | ||
| gfv_filter <- NULL | ||
| } else { | ||
| gfv_filter <- all_lineages[all_lineages != lineage] | ||
| } | ||
| fitted_vals <- getFittedValues(test.dyn.res, | ||
| genes = gene, | ||
| pt = pt, | ||
| expr.mat = expr.mat, | ||
| size.factor.offset = size.factor.offset, | ||
| log1p.norm = log1p.norm, | ||
| filter.lineage = gfv_filter) | ||
| if (log1p.norm) { | ||
| fitted_vals <- dplyr::select(fitted_vals, | ||
| cell, | ||
| lineage, | ||
| pt, | ||
| gene, | ||
| rna = rna_log1p, | ||
| scLANE_pred = scLANE_pred_log1p, | ||
| scLANE_ci_ll = scLANE_ci_ll_log1p, | ||
| scLANE_ci_ul = scLANE_ci_ul_log1p) | ||
| } else { | ||
| fitted_vals <- dplyr::select(fitted_vals, | ||
| cell, | ||
| lineage, | ||
| pt, | ||
| gene, | ||
| rna, | ||
| scLANE_pred, | ||
| scLANE_ci_ll, | ||
| scLANE_ci_ul) | ||
|
|
||
| } | ||
| # generate dynamics plot | ||
| dyn_plot <- ggplot2::ggplot(fitted_vals, ggplot2::aes(x = pt, y = rna)) + | ||
| ggplot2::geom_point(size = 1.5, | ||
| alpha = 0.6, | ||
| stroke = 0, | ||
| color = "grey30") + | ||
| ggplot2::geom_vline(data = data.frame(gene = gene, knot = unique(test.dyn.res[[gene]][[paste0("Lineage_", lineage)]]$MARGE_Slope_Data$Breakpoint)), | ||
| mapping = ggplot2::aes(xintercept = knot), | ||
| linetype = "dashed", | ||
| color = "black", | ||
| linewidth = 0.75) + | ||
| ggplot2::geom_ribbon(ggplot2::aes(ymin = scLANE_ci_ll, ymax = scLANE_ci_ul), | ||
| linewidth = 0, | ||
| fill = "darkgreen", | ||
| alpha = 0.35) + | ||
| ggplot2::geom_line(ggplot2::aes(y = scLANE_pred), | ||
| color = "darkgreen", | ||
| linewidth = 0.75) + | ||
| ggplot2::scale_x_continuous(labels = scales::label_number(accuracy = 0.01)) + | ||
| ggplot2::labs(x = "Pseudotime", | ||
| y = ifelse(log1p.norm, "Normalized Expression", "Expression")) + | ||
| theme_scLANE() | ||
| # generate coefficient summary | ||
| min_pt <- min(pt[, which(LETTERS == lineage)], na.rm = TRUE) | ||
| max_pt <- max(pt[, which(LETTERS == lineage)], na.rm = TRUE) | ||
| coef_sumy <- dplyr::select(test.dyn.res[[gene]][[paste0("Lineage_", lineage)]]$Gene_Dynamics, | ||
| -dplyr::starts_with("Trend")) %>% | ||
| tidyr::pivot_longer(dplyr::starts_with("Slope"), | ||
| names_to = "Segment", | ||
| values_to = "Coef") %>% | ||
| dplyr::mutate(Breakpoint_Lag = dplyr::lag(Breakpoint), | ||
| Breakpoint_Lead = dplyr::lead(Breakpoint), | ||
| Interval = NA_character_, | ||
| .before = 4) %>% | ||
| dplyr::mutate(Breakpoint_Lag = dplyr::if_else(is.na(Breakpoint_Lag), min_pt, Breakpoint_Lag), | ||
| Breakpoint_Lead = dplyr::if_else(is.na(Breakpoint_Lead), max_pt, Breakpoint_Lead), | ||
| Interval = paste0("(", round(Breakpoint_Lag, 3), ", ", round(Breakpoint_Lead, 3), ")")) %>% | ||
| dplyr::select(Interval, Coef) %>% | ||
| dplyr::mutate(Coef = round(Coef, 3)) | ||
| # convert coefficient summary to grob | ||
| coef_sumy_grob <- gridExtra::tableGrob(coef_sumy, | ||
| cols = c("Interval", "Coefficient"), | ||
| theme = gridExtra::ttheme_minimal(base_size = 11, | ||
| core = list(fg_params = list(hjust = 0, x = 0.05)), | ||
| colhead = list(fg_params = list(hjust = 0, x = 0.05))), | ||
| rows = NULL) %>% | ||
| gtable::gtable_add_grob(grobs = grid::rectGrob(gp = grid::gpar(fill = NA, lwd = 3)), | ||
| t = 1, | ||
| b = nrow(.), | ||
| l = 1, | ||
| r = ncol(.)) %>% | ||
| gtable::gtable_add_grob(grobs = grid::rectGrob(gp = grid::gpar(fill = NA, lwd = 3)), | ||
| t = 1, | ||
| l = 1, | ||
| r = ncol(.)) | ||
|
|
||
| # combine objects | ||
| dyn_plot_anno <- ggpubr::ggarrange(dyn_plot, | ||
| coef_sumy_grob, | ||
| ncol = 2, | ||
| nrow = 1, | ||
| widths = c(2, 1)) | ||
| return(dyn_plot_anno) | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.