A single-cell and spatial atlas of early human olfactory development.

Code ↔ Paper

26 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 26 matches

[1] § Methods › Data integration and dimensionality reduction ↔ R/figS4.Rmd, lines 81–138 · score 1.00 · IQCJ SCHIP1, COL13A1, COL19A1, KIF18B, OLFML2A, PDE1C
[2] § Methods › Detection and removal of maternal erythroid contamination ↔ R/000-functions.R, lines 345–381 · score 1.00 · canonical hemoglobin genes, erythroid contamination, sex score thresholds, erythroid signal, DDX3Y, RPS4Y1
[3] § Methods › Data integration and dimensionality reduction ↔ R/002-annotation.Rmd, lines 269–312 · score 0.99 · IQCJ SCHIP1, COL17A1, CycBasal, GnRH, TMEM119, ACAN
[4] § Results › Single-nucleus atlas of the human fetal nasal compartment ↔ 05.scmusk.ipynb, lines 245–269 · score 0.83 · skeletal muscle, olfactory ensheathing, Schwann cells, neural crest, respiratory epithelial, pericytes
[5] § Results › Single-nucleus atlas of the human fetal nasal compartment ↔ 05.scmusk.ipynb, lines 147–227 · score 0.77 · NOS1 neurons, cycling basal, Schwann cell, GnRH, precursors, GABA
[6] § Methods › Bulk RNA-seq QC and reference for demultiplexing ↔ Python/figS1.ipynb, lines 61–195 · score 0.77 · minor allele frequency, Bulk RNA, missingness, Python, variants, genotype
[7] § Results › Transcriptional regulators and gene networks shaping human OE development ↔ R/fig4.Rmd, lines 98–137 · score 0.76 · E2F1, ATOH1, EGR2, ID1, ID3, KLF15
[8] § Results › Transcriptional regulators and gene networks shaping human OE development ↔ R/fig4.Rmd, lines 338–363 · score 0.73 · ARID1B, PHOX2A, BARX2, NEUROD2, HOXC8, NFKBIZ
[9] § Results › Transcriptional regulators and gene networks shaping human OE development ↔ R/000-functions.R, lines 5364–5394 · score 0.71 · absolute regulatory, Source TFs, regulatory network, transcription factor, Node, scores
[10] § Results › Olfactory receptor gene expression patterns during human OE development ↔ R/000-functions.R, lines 1781–1819 · score 0.70 · sparse single cell, cell identity, olfactory receptor, Dominance scores, iOSN, bin
[11] § Results › Composition and development of the olfactory system in human fetuses ↔ R/000-functions.R, lines 2426–2489 · score 0.69 · Temporally regulated genes, dynamic range, FDR, Slingshot, smoothed, GAMs
[12] § Methods › Genetic demultiplexing ↔ Python/figS1.ipynb, lines 61–195 · score 0.67 · donor genotype, bulk RNA, Allele, variants, SNPs, VCFs
[13] § Methods › Statistics and reproducibility ↔ R/000-functions.R, lines 3491–3539 · score 0.65 · Kruskal Wallis, pairwise comparisons, Wilcoxon, dynamics, genes, cell
[14] § Methods › Transcription factory activity ↔ R/fig4.Rmd, lines 58–96 · score 0.61 · decoupleR, TF activities, Seurat, assay, heatmaps, gene
[15] § Results › Spatially resolved olfactory epithelium OR genes profiling confirms the 1 OR-1 OSN rule at single-cell resolution ↔ R/000-functions.R, lines 2042–2107 · score 0.59 · adjusted dominance score, expression magnitude, metric, iOSN, gene, cell
[16] § Methods › Single-cell olfactory receptor (OR) dominance ↔ R/000-functions.R, lines 1781–1819 · score 0.59 · dominance scores quantified, zero, Seurat, binned, ORs, receptor
[17] § Results › Single-nucleus atlas of the human fetal nasal compartment ↔ Python/figS1.ipynb, lines 509–584 · score 0.58 · bulk genotypes, Bulk RNA seq, missingness, Souporcell, Vireo, SNP
[18] § Results › Single-nucleus atlas of the human fetal nasal compartment ↔ R/000-functions.R, lines 345–381 · score 0.58 · linked gene expression, discordant, erythroid, classifying, contamination, XIST
[19] § Results › Single-nucleus atlas of the human fetal nasal compartment ↔ R/000-custom-colors.R, lines 2–44 · score 0.56 · early human, GnRH, GABA, GLUT, iOSN, NOS1
[20] § Results › Composition and development of the olfactory system in human fetuses ↔ R/figS4.Rmd, lines 81–138 · score 0.56 · TOP2A, MYBL1, NNAT, NEUROG1, NEUROD1, HES6
[21] § Methods › Transcription factory activity ↔ R/000-functions.R, lines 5364–5394 · score 0.54 · regulatory networks, igraph, ggraph, TFs, transcriptional
[22] § Results › Composition and development of the olfactory system in human fetuses ↔ R/002-annotation.Rmd, lines 269–312 · score 0.54 · gnrh1, MUC5AC, iOSN, SOX2, RHBC, MV
[23] § Results › Single-nucleus atlas of the human fetal nasal compartment ↔ R/000-functions.R, lines 173–206 · score 0.53 · quality control, detected genes, UMIs, mitochondrial, fractions, log
[24] § Results › Composition and development of the olfactory system in human fetuses ↔ R/fig3.Rmd, lines 73–86 · score 0.52 · TOP2A, FOXN4, MYBL1, NNAT, NEUROG1, NEUROD1
[25] § Methods › Trajectory inference and pseudotime ↔ R/000-functions.R, lines 2426–2489 · score 0.52 · highly variable genes, Slingshot, iOSN, lineage, pseudotime, OHBC
[26] § Methods › Statistics and reproducibility ↔ 07.proportion.ipynb, lines 42–62 · score 0.51 · Mann Whitney, Wilcoxon, Kruskal, cell

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R · 6,024 lines · 195 KB · GPL-3.0 · 11 matches

# Customized R functions
# Author= Yvon Mbouamboua ([email hidden])
# Utility: %||%
`%||%` <- function(a, b) if(!is.null(a)) a else b
#' Remove duplicated cell barcodes within and/or across Seurat objects
#'
#' This function cleans duplicated cell barcodes in a list of Seurat objects.
#' It can remove duplicates occurring **within each object**, **across objects**,
#' or **both**. A detailed log is returned reporting how many cells were removed
#' per object at each step.
#'
#' @param objs Named list of \code{Seurat} objects.
#' @param mode Character string specifying which duplicates to remove:
#' \itemize{
#' \item \code{"within"} – remove duplicated barcodes within each object
#' \item \code{"across"} – remove duplicated barcodes shared across objects
#' \item \code{"both"} – remove both within- and across-object duplicates (default)
#' }
#' @param save Logical; whether to save cleaned Seurat objects to disk (default: \code{TRUE}).
#' @param outdir Output directory where cleaned objects will be written.
#' @param suffix Filename suffix for saved objects (default: \code{"_clean.rds"}).
#' @param verbose Logical; whether to print a summary table to the console (default: \code{TRUE}).
#'
#' @return A list with two elements:
#' \itemize{
#' \item \code{objs} – cleaned list of Seurat objects
#' \item \code{log} – data.frame summarizing barcode removal per object
#' }
#'
#' @details
#' Duplicate removal strategy:
#' \enumerate{
#' \item Remove duplicated cell barcodes *within* each Seurat object using
#' \code{duplicated(colnames(object))}.
#' \item Identify barcodes shared across multiple objects and remove them
#' from all objects.
#' }
#'
#' The function does not rename barcodes; duplicates are removed entirely.
#' This behavior is appropriate for demultiplexed or merged datasets where
#' barcode reuse indicates technical duplication.
#'
#' @examples
#' \dontrun{
#' res <- clean_seurat_barcodes(
#' objs,
#' mode = "both",
#' outdir = "cleaned_objects"
#' )
#'
#' objs_clean <- res$objs
#' removal_log <- res$log
#' }
#'
#' @seealso \code{\link[Seurat]{Seurat}}, \code{\link[base]{duplicated}}
#'
#' @export
#'
#' Clean duplicated cell barcodes in Seurat objects
#'
#' Removes duplicated cell barcodes within individual Seurat objects,
#' across multiple Seurat objects, or both.
#'
#' @param objs A Seurat object or a named list of Seurat objects
#' @param mode One of "within", "across", or "both"
#' @param outdir Optional output directory to save cleaned objects
#' @param verbose Logical; print summary messages
#' @param save Logical; save data
#' @return A list with cleaned Seurat objects and a log table
#' @export
#'
clean_seurat_barcodes <- function(
objs,
mode = c("within", "across", "both"),
outdir = NULL,
save = FALSE,
verbose = TRUE
) {
mode <- match.arg(mode)
# CRITICAL FIX — MUST BE FIRST
if (inherits(objs, "Seurat")) {
objs <- list(object = objs)
}
if (!is.list(objs)) {
stop("`objs` must be a Seurat object or a list of Seurat objects.")
}
# Initialize log
log <- data.frame(
object = names(objs),
removed_within = 0L,
removed_across = 0L,
kept_cells = 0L,
stringsAsFactors = FALSE
)
# Remove duplicates WITHIN each object
if (mode %in% c("within", "both")) {
for (i in seq_along(objs)) {
obj <- objs[[i]]
cn <- colnames(obj)
dup_within <- duplicated(cn)
n_removed <- sum(dup_within)
if (n_removed > 0) {
obj <- obj[, !dup_within]
}
objs[[i]] <- obj
log$removed_within[i] <- n_removed
}
}
# Remove duplicates ACROSS objects
if (mode %in% c("across", "both") && length(objs) > 1) {
all_cells <- unlist(lapply(objs, colnames))
dup_cells <- all_cells[duplicated(all_cells)]
if (length(dup_cells) > 0) {
dup_cells <- unique(dup_cells)
for (i in seq_along(objs)) {
obj <- objs[[i]]
cn <- colnames(obj)
to_remove <- cn %in% dup_cells
n_removed <- sum(to_remove)
if (n_removed > 0) {
obj <- obj[, !to_remove]
}
objs[[i]] <- obj
log$removed_across[i] <- n_removed
}
}
}
# Final counts
for (i in seq_along(objs)) {
log$kept_cells[i] <- ncol(objs[[i]])
}
# Save cleaned objects
if (save) {
if (!is.null(outdir)) {
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
for (nm in names(objs)) {
saveRDS(objs[[nm]], file.path(outdir, paste0(nm, ".rds")))
}
}
}
if (verbose) {
message("✔ Barcode cleaning complete")
print(log)
}
return(list(
objects = objs,
log = log
))
}
#' Run Quality Control on Seurat Objects
#'
#' Performs automated cell filtering based on gene/UMI counts, mitochondrial
#' and ribosomal content, dropout rate, and optional doublet detection.
#' Works on single or multiple Seurat objects.
#'
#' @param x Seurat object or list of Seurat objects.
#' @param min_feat Minimum number of detected genes.
#' @param min_umi Minimum number of UMIs.
#' @param mad_n Number of MADs for upper cutoffs (if method = "MAD").
#' @param max_mito Maximum mitochondrial percentage.
#' @param calc_ribo Logical; compute ribosomal gene percentage.
#' @param max_ribo Maximum ribosomal percentage.
#' @param rm_dbl Logical; perform doublet removal using scDblFinder.
#' @param calc_drop Logical; compute dropout fraction.
#' @param max_drop Maximum dropout fraction.
#' @param method Thresholding method: "MAD", "fixed", or "none".
#' @param fixed_thr Named list of fixed upper cutoffs.
#' @param mito_pat Regex for mitochondrial genes.
#' @param ribo_pat Regex for ribosomal genes.
#' @param species "human" or "mouse" (sets default gene patterns).
#' @param log_g2u Logical; filter by log10(genes/UMI).
#' @param min_g2u Minimum genes/UMI ratio.
#' @param sample_col Metadata column for sample ID.
#' @param outdir Output directory.
#' @param merge Logical; merge outputs into a single Seurat object.
#'
#' @return A filtered Seurat object, list, or merged object.
#' @export
#'
#' @examples
#' qc_data <- run_qc(seurat_list, calc_ribo = TRUE, rm_dbl = TRUE)
#'
run_qc <- function(
x,
min_feat = 200,
min_umi = 500,
mad_n = 5,
max_mito = 5,
calc_ribo = FALSE,
max_ribo = 3,
calc_drop = FALSE,
max_drop = 0.95,
rm_dbl = FALSE,
method = c("MAD", "fixed", "none"),
fixed_thr = list(max_feat = 6000, max_umi = 20000),
mito_pat = NULL,
ribo_pat = NULL,
species = c("human", "mouse"),
sample_col = "orig.ident",
outdir = "QC",
merge_output = TRUE,
verbose = TRUE
) {
library(Seurat)
library(dplyr)
library(Matrix)
method <- match.arg(method)
species <- match.arg(species)
# Default gene patterns
if (is.null(mito_pat)) mito_pat <- if (species == "human") "^MT-" else "^mt-"
if (is.null(ribo_pat)) ribo_pat <- if (species == "human") "^RP[LS]" else "^Rp[ls]"
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)
# Input validation
if (inherits(x, "Seurat")) {
obj_list <- list(x)
single <- TRUE
} else if (is.list(x) && all(sapply(x, inherits, "Seurat"))) {
obj_list <- x
single <- FALSE
} else stop("Input must be a Seurat object or list of Seurat objects.")
qc_sum <- data.frame()
dbl_sum <- data.frame()
if (verbose) cat("\n>>> Starting quality control for", length(obj_list), "sample(s)...\n")
for (i in seq_along(obj_list)) {
obj <- obj_list[[i]]
sample <- if (sample_col %in% colnames([email hidden]))
unique([email hidden][[sample_col]]) else paste0("Sample_", i)
pre_cells <- ncol(obj)
if (verbose) cat("\n────────────────────────────\n Sample:", sample, "| Cells:", pre_cells, "\n")
# Compute QC metrics
obj$percent_mito <- PercentageFeatureSet(obj, pattern = mito_pat)
if (calc_ribo) obj$percent_ribo <- PercentageFeatureSet(obj, pattern = ribo_pat)
if (calc_drop) {
counts <- GetAssayData(obj, slot = "counts")
obj$dropout <- Matrix::colSums(counts == 0) / nrow(counts)
}
# Compute thresholds
if (method == "MAD") {
max_feat <- median(obj$nFeature_RNA) + mad_n * mad(obj$nFeature_RNA)
max_umi <- median(obj$nCount_RNA) + mad_n * mad(obj$nCount_RNA)
} else if (method == "fixed") {
max_feat <- fixed_thr$max_feat
max_umi <- fixed_thr$max_umi
} else {
max_feat <- Inf; max_umi <- Inf
}
# Apply permissive filters
filt <- subset(
obj,
nFeature_RNA > min_feat &
nFeature_RNA < max_feat &
nCount_RNA > min_umi &
percent_mito < max_mito &
(if (calc_ribo) percent_ribo < max_ribo else TRUE) &
(if (calc_drop) dropout < max_drop else TRUE)
# No max UMI cutoff — preserves high-depth nuclei
)
post_cells <- ncol(filt)
pct_rm <- round((1 - post_cells / pre_cells) * 100, 2)
if (verbose) cat(">>> Retained:", post_cells, "cells (", pct_rm, "% removed)\n")
# Optional doublet filtering
if (rm_dbl) {
if (verbose) cat(">>> Running scDblFinder...\n")
sce <- scDblFinder::scDblFinder(as.SingleCellExperiment(filt))
filt$scDblFinder <- sce$scDblFinder.class
filt <- subset(filt, scDblFinder == "singlet")
dbl_sum <- rbind(dbl_sum, data.frame(Sample = sample, Doublets = sum(sce$scDblFinder.class == "doublet")))
}
# Summarize
qc_sum <- rbind(qc_sum, data.frame(
Sample = sample,
Pre_Cells = pre_cells,
Post_Cells = post_cells,
Removed_Pct = pct_rm,
min_feat, max_feat, min_umi,
max_mito,
Method = method, MADs = mad_n
))
obj_list[[i]] <- filt
}
# Save summary
write.csv(qc_sum, file.path(outdir, "QC_summary.csv"), row.names = FALSE)
if (rm_dbl && nrow(dbl_sum) > 0)
write.csv(dbl_sum, file.path(outdir, "Doublet_summary.csv"), row.names = FALSE)
# Return object
if (merge_output && length(obj_list) > 1) {
merged <- merge(x = obj_list[[1]], y = obj_list[-1], merge.data = TRUE)
if (verbose) cat("\n>>> Merged object returned.\n")
if (verbose) cat("\n>>> Join leyers.\n")
#merged <- JoinLayers(merged)
#merged <- fix_seurat_matrix_names(merged)
return(merged)
} else if (single) {
if (verbose) cat("\n>>> Returning single filtered object.\n")
#obj_list[[1]] <- fix_seurat_matrix_names(obj_list[[1]])
return(obj_list[[1]])
} else {
if (verbose) cat("\n>>> Returning list of filtered objects.\n")
return(obj_list)
}
}
#' Sex and Erythroid Contamination QC for Seurat Objects
#'
#' This function computes per-cell sex scores based on XIST and Y-linked gene expression,
#' flags sex-discordant cells, calculates erythroid contamination using canonical hemoglobin genes,
#' and filters cells that are both sex-discordant and above a specified erythroid threshold. Optional QC plots are produced.
#'
#' @param obj A Seurat object containing single-cell RNA-seq data.
#' @param sample_col Character. Column name in `[email hidden]` corresponding to sample IDs. Default is "sample".
#' @param y_genes Character vector. List of Y-chromosome genes to use for sex score. Default is c("UTY","RPS4Y1","ZFY","DDX3Y","KDM5D").
#' @param xist_gene Character. Gene name for XIST. Default is "XIST".
#' @param eps Numeric. Small pseudocount added to avoid division by zero. Default is 1e-6.
#' @param female_thresh Numeric. Sex score threshold above which a cell is classified as female-like. Default is 1.
#' @param male_thresh Numeric. Sex score threshold below which a cell is classified as male-like. Default is -1.
#' @param eryth_genes Character vector. Canonical hemoglobin genes for erythroid contamination assessment. Default is c("HBB","HBA1","HBA2","HBE1","HBG1","HBG2","HBM").
#' @param eryth_percentile Numeric. Percentile to define high erythroid signal for filtering (0-1). Default is 0.95.
#' @param make_plots Logical. Whether to produce QC plots. Default is TRUE.
#' @param plot_prefix Character. Prefix for saved QC plot filenames. Default is "qc_".
#'
#' @return A list with the following components:
#' \item{object}{The input Seurat object with added metadata columns: sex_score, sex_call, eryth_sum, inferred_sex, sex_discordant.}
#' \item{filtered}{A Seurat object filtered to remove cells that are both sex-discordant and above the erythroid threshold.}
#' \item{summary}{A per-sample summary table including cell counts, fractions female/male/ambiguous, median sex score, fraction of cells with erythroid expression, and median erythroid sum.}
#' \item{eryth_threshold}{Numeric value corresponding to the threshold for high erythroid signal (computed from `eryth_percentile`).}
#' \item{n_removed}{Number of cells removed by filtering.}
#'
#' @examples
#' \dontrun{
#' qc <- sex_contamination_qc(obj, sample_col = "sample")
#' qc$summary
#' filtered_obj <- qc$filtered
#' }
#'
#' @import Seurat
#' @import dplyr
#' @import ggplot2
#' @export
sex_contamination_qc <- function(
obj,
sample_col = "sample",
y_genes = c("UTY","RPS4Y1","ZFY","DDX3Y","KDM5D"),
xist_gene = "XIST",
eps = 1e-6,
female_thresh = 1,
male_thresh = -1,
eryth_genes = c("HBB","HBA1","HBA2","HBE1","HBG1","HBG2","HBM"),
majority_threshold = 0.5,
eryth_percentile = 0.95,
remove_ambiguous = TRUE,
min_sex_score_diff = NULL,
out_file = NULL
){
cat("Extracting expression matrix (slot = data)...")
Seurat::DefaultAssay(obj) <- "RNA"
expr <- Seurat::GetAssayData(obj, slot = "data")
meta_df <- [email hidden]
# Sex score
cat("Computing sex scores...")
meta_df$XIST <- if (xist_gene %in% rownames(expr)) expr[xist_gene, ] else 0
y_present <- intersect(y_genes, rownames(expr))
meta_df$Ysum <- if (length(y_present) > 0)
Matrix::colSums(expr[y_present, , drop = FALSE])
else 0
meta_df$Ysum[is.na(meta_df$Ysum)] <- 0
meta_df$XIST[is.na(meta_df$XIST)] <- 0
meta_df$sex_score <- log2((meta_df$XIST + eps)/(meta_df$Ysum + eps))
meta_df$sex_call <- ifelse(
meta_df$sex_score > female_thresh, "Female-like",
ifelse(meta_df$sex_score < male_thresh, "Male-like", "Ambiguous")
)
# Erythroid score
cat("Computing erythroid scores...")
eryth_present <- intersect(eryth_genes, rownames(expr))
meta_df$eryth_sum <- if (length(eryth_present) > 0)
Matrix::colSums(expr[eryth_present, , drop = FALSE])
else 0
meta_df$eryth_sum[is.na(meta_df$eryth_sum)] <- 0
eryth_threshold <- as.numeric(
stats::quantile(meta_df$eryth_sum, eryth_percentile, na.rm = TRUE)
)
# Summary per sample
cat("Summarizing per-sample sex composition...")
tbl_summary <- meta_df %>%
dplyr::group_by(.data[[sample_col]]) %>%
dplyr::summarise(
n_cells = dplyr::n(),
frac_female_like = mean(sex_call == "Female-like"),
frac_male_like = mean(sex_call == "Male-like"),
frac_ambiguous = mean(sex_call == "Ambiguous"),
median_sex_score = stats::median(sex_score),
median_eryth_sum = stats::median(eryth_sum),
.groups = "drop"
)
# Infer sample sex
inferred <- tbl_summary %>%
dplyr::mutate(inferred_sex =
ifelse(frac_female_like > majority_threshold, "Female", "Male")) %>%
dplyr::select(.data[[sample_col]], inferred_sex) %>%
tibble::deframe()
meta_df$inferred_sex <- unname(inferred[meta_df[[sample_col]]])
# Discordance
meta_df$sex_discordant <-
(meta_df$inferred_sex == "Female" & meta_df$sex_call == "Male-like") |
(meta_df$inferred_sex == "Male" & meta_df$sex_call == "Female-like")
if (!is.null(min_sex_score_diff)) {
global_median <- stats::median(meta_df$sex_score, na.rm = TRUE)
meta_df$strong_discordance <- meta_df$sex_discordant &
abs(meta_df$sex_score - global_median) >= min_sex_score_diff
} else {
meta_df$strong_discordance <- meta_df$sex_discordant
}
# Flag cells for removal
cat("Flagging cells for removal...")
remove_cells <- meta_df$strong_discordance & (meta_df$eryth_sum > eryth_threshold)
if (remove_ambiguous)
remove_cells <- remove_cells | meta_df$sex_call == "Ambiguous"
cat(sum(remove_cells), " cells flagged.")
# Write metadat back
obj <- Seurat::AddMetaData(obj, metadata = meta_df)
# Subset filtered object
cat("Creating filtered object...")
keep_cells <- rownames(meta_df)[!remove_cells]
filt_obj <- subset(obj, cells = keep_cells)
# Summary for output
removal_summary <- data.frame(
sample = meta_df[[sample_col]],
removed = remove_cells
) %>%
dplyr::group_by(sample) %>%
dplyr::summarise(
n_removed = sum(removed),
frac_removed = mean(removed),
.groups = "drop"
)
tbl_summary <- tbl_summary %>%
dplyr::left_join(removal_summary, by = "sample") %>%
dplyr::mutate(eryth_threshold = eryth_threshold)
# Output file
if (!is.null(out_file)) {
utils::write.table(tbl_summary, file = out_file, sep = "\t",
quote = FALSE, row.names = FALSE)
}
# Fix seurat matrix names
obj <- fix_seurat_matrix_names(obj)
filt_obj <- fix_seurat_matrix_names(filt_obj)
# Return output
return(list(
object = obj,
filtered = filt_obj,
tbl_summary = tbl_summary,
eryth_threshold = eryth_threshold,
n_removed = sum(remove_cells)
))
}
genefilter <- function(
obj,
mode = c("B","A","C"),
keep.genes = c("MEG3"),
excl.genes = NULL,
assay = "RNA",
normalize = TRUE,
norm.method = "LogNormalize",
scale.fac = 10000,
find.hvg = TRUE,
n.hvg = 2000,
scale.data = TRUE,
verbose = TRUE
){
stopifnot(inherits(obj, "Seurat"))
mode <- match.arg(mode)
DefaultAssay(obj) <- assay
genes <- rownames(obj[[assay]])
# robust meta.features accessor (works for Seurat v5)
get_gene_biotype <- function(obj, assay) {
assay_obj <- obj[[assay]]
# Assay5: meta.features only exists if user added it
if (!"meta.features" %in% slotNames(assay_obj)) {
return(NULL)
}
mf <- [email hidden]
if (is.null(mf) || nrow(mf) == 0) return(NULL)
mf <- as.data.frame(mf)
if (!"gene_biotype" %in% colnames(mf)) return(NULL)
x <- mf$gene_biotype
names(x) <- rownames(mf)
return(x)
}
gb <- get_gene_biotype(obj, assay)
# filtering rules
filters <- list(
biotypes_rm = c(
"lincRNA","lncRNA","antisense_RNA","processed_transcript",
"pseudogene","processed_pseudogene","unprocessed_pseudogene",
"transcribed_unprocessed_pseudogene","transcribed_processed_pseudogene",
"unitary_pseudogene","Mt_rRNA","Mt_tRNA","rRNA"
),
regex_rm = c(
"^MT-",
"^(RPL|RPS)[0-9]+",
"^HB[ABDGEMQZ]",
"^LINC[0-9]+",
"^(AC|AL|AP|RP)[0-9]+",
"^(CTD|CTC)-",
"^LOC[0-9]+",
"^C[0-9]+orf", "-AS[0-9]+$",
"^FAM[0-9]+","^DAZ[0-9]+$",
"^XXBAC", "^RPS4Y[12]$",
"\\.[0-9]+$", "^TTTY[0-9]+$",
"^ENSG[0-9]+$","^USP9Y$",
"^XIST$", "^TSIX$",
"^SRY$", "^ZFY$", "^UTY$"
)
)
genes_up <- toupper(genes)
names(genes_up) <- genes
rm_biotype <- c()
if (!is.null(gb)) {
for (bt in filters$biotypes_rm)
rm_biotype <- c(rm_biotype, names(gb)[gb == bt])
}
rm_regex <- c()
for (rx in filters$regex_rm)
rm_regex <- c(rm_regex, genes[grepl(rx, genes_up)])
must_remove <- unique(c(rm_biotype, rm_regex))
# keep sets
if (!is.null(gb)) {
prot_ok <- names(gb)[gb == "protein_coding"]
} else {
prot_ok <- setdiff(genes, must_remove)
}
if (mode == "A") {
keep <- prot_ok
} else if (mode == "C") {
keep <- setdiff(genes, must_remove)
} else {
keep <- union(prot_ok, keep.genes)
}
if (!is.null(excl.genes))
keep <- setdiff(keep, excl.genes)
keep <- setdiff(keep, must_remove)
removed <- setdiff(genes, keep)
if (verbose) {
message("=== filter_genes (Seurat v5-safe) ===")
message("Genes:", length(genes))
message("Kept :", length(keep))
message("Removed:", length(removed))
}
# subset + normalization
obj_f <- subset(obj, features = keep)
if (normalize)
obj_f <- NormalizeData(obj_f, normalization.method = norm.method, scale.factor = scale.fac)
if (find.hvg)
obj_f <- FindVariableFeatures(obj_f, nfeatures = n.hvg)
if (scale.data)
obj_f <- ScaleData(obj_f, features = VariableFeatures(obj_f))
return(obj_f)
}
#' Filter Outliers in Seurat Clusters Using Local Density
#'
#' @description
#' Removes low-density (outlier) cells within specified clusters of a Seurat object,
#' based on k-nearest neighbor (KNN) distances in a given dimensional reduction (e.g., UMAP or PCA).
#' Optionally visualizes before/after filtering for each cluster and summarizes retained/removed cell counts.
#'
#' @param obj A Seurat object.
#' @param group.by Character string. Metadata column used to define clusters (default: `"seurat_clusters"`).
#' @param reduction Character string. Dimensional reduction to use (e.g., `"umap"`, `"pca"`).
#' @param density.threshold Numeric in (0,1). Quantile of local density defining cutoff (default: `0.95`).
#' @param k Integer. Number of neighbors for KNN-based local density (default: `30`).
#' @param relax Numeric multiplier (>1) to loosen the density threshold (default: `2`).
#' @param clean.all Logical. If `TRUE`, filter all clusters; otherwise specify `target.clusters`.
#' @param target.clusters Vector of cluster identities to clean. Ignored if `clean.all = TRUE`.
#' @param verbose Logical. Print progress and summary (default: `TRUE`).
#' @param color Color for highlighted cells in plots (default: `"red"`).
#' @param pt.size Point size for plotting (default: `1`).
#' @param base.size Base text size for plot titles (default: `10`).
#' @param ncol Integer. Number of columns for the combined plot layout (default: `2`).
#' @param return.plot Logical. Return combined patchwork plot (default: `FALSE`).
#'
#' @return A list with:
#' \itemize{
#' \item \code{obj}: The filtered Seurat object.
#' \item \code{summary}: A data.frame summarizing before/after cell counts per cluster.
#' \item \code{plot}: (optional) Combined before/after UMAP visualization if `return.plot = TRUE`.
#' }
#'
#' @details
#' The function estimates local cell density for each cluster using mean KNN distances.
#' Cells below the density cutoff (low-density regions) are considered outliers and removed.
#' Clusters not listed in `target.clusters` remain untouched.
#'
#' @examples
#' \dontrun{
#' res <- filter_outliers(
#' obj = obj,
#' group.by = "seurat_clusters",
#' target.clusters = c(0,1,2,3),
#' reduction = "umap",
#' density.threshold = 0.95,
#' k = 30,
#' relax = 2,
#' verbose = TRUE,
#' return.plot = TRUE
#' )
#' res$summary
#' print(res$plot)
#' }
#'
#' @import Seurat
#' @import FNN
#' @import patchwork
#' @export
filter_outliers <- function(obj,
group.by = "seurat_clusters",
reduction = "umap",
density.threshold = 0.95,
k = 30,
relax = 2,
clean.all = FALSE,
target.clusters = NULL,
verbose = TRUE,
color = "red",
pt.size = 1,
base.size = 10,
ncol = 2,
return.plot = FALSE) {
require(Seurat)
require(FNN)
require(patchwork)
#Checks
if (!group.by %in% colnames([email hidden]))
stop("Column ", group.by, " not found in metadata.")
if (!reduction %in% Reductions(obj))
stop("Reduction ", reduction, " not found in object.")
obj <- SetIdent(obj, value = group.by)
uniq.clust <- unique(na.omit(as.character([email hidden][[group.by]])))
if (clean.all) {
clusts <- uniq.clust
} else if (!is.null(target.clusters)) {
bad <- setdiff(target.clusters, uniq.clust)
if (length(bad) > 0) stop("Invalid clusters: ", paste(bad, collapse = ", "))
clusts <- target.clusters
} else {
stop("Specify either clean.all = TRUE or target.clusters.")
}
if (verbose) cat("Processing clusters:", paste(clusts, collapse = ", "), "\n")
#Loop
all.keep <- c()
plots <- list()
summ <- data.frame()
emb <- Embeddings(obj, reduction)
for (cl in clusts) {
cells <- rownames([email hidden])[[email hidden][[group.by]] == cl]
if (length(cells) < 2) next
coord <- emb[cells, , drop = FALSE]
kk <- min(k, length(cells) - 1)
dens <- rowMeans(knn.dist(coord, k = kk), na.rm = TRUE)
[email hidden][cells, "dens"] <- dens
thr <- quantile(dens, density.threshold, na.rm = TRUE)
keep <- cells[dens <= thr * relax]
all.keep <- c(all.keep, keep)
summ <- rbind(summ, data.frame(
cluster = cl,
before = length(cells),
after = length(keep),
removed = length(cells) - length(keep),
pct.rm = round((1 - length(keep)/length(cells)) * 100, 2)
))
p1 <- DimPlot(obj, reduction = reduction, group.by = group.by,
cells.highlight = cells, cols.highlight = color,
sizes.highlight = pt.size) +
NoAxes() + NoLegend() +
ggtitle(paste("Cluster", cl, "(n=", length(cells), ")")) +
theme(plot.title = element_text(hjust = 0.5, size = base.size))
p2 <- DimPlot(obj, reduction = reduction, group.by = group.by,
cells.highlight = keep, cols.highlight = color,
sizes.highlight = pt.size) +
NoAxes() + NoLegend() +
ggtitle(paste("Filtered", cl, "(n=", length(keep),
")\nCutoff=", round(thr, 2), " K=", kk)) +
theme(plot.title = element_text(hjust = 0.5, size = base.size))
plots[[as.character(cl)]] <- p1 + p2 + patchwork::plot_layout(ncol = 2)
}
#Keep filtered + untouched clusters
all.target <- rownames([email hidden])[[email hidden][[group.by]] %in% clusts]
other.cells <- setdiff(colnames(obj), all.target)
keep.cells <- unique(c(all.keep, other.cells))
filt <- subset(obj, cells = keep.cells)
if (verbose) {
cat("\n=== FILTER SUMMARY ===\n")
print(summ, row.names = FALSE)
cat("Total before:", ncol(obj), "\n")
cat("Total after:", ncol(filt), "\n")
cat("Removed:", ncol(obj) - ncol(filt), "\n")
}
if (return.plot) {
combo <- wrap_plots(plots, ncol = ncol)
return(list(obj = filt, summary = summ, plot = combo))
} else {
return(list(obj = filt, summary = summ))
}
}
#' Compute cluster statistics from Seurat object (multiple group.by)
#'
#' @param object Seurat object
#' @param ident Metadata column for cluster identities (default: "seurat_clusters")
#' @param group.by One or more metadata columns to group by (default: "orig.ident")
#' @param rm.sum Logical; remove total row (default: FALSE)
#' @param rm.pct Logical; hide percentage columns (default: FALSE)
#' @return Data frame of cluster counts and percentages
#' @export
cluster_stats <- function(object, ident = "seurat_clusters",
group.by = "orig.ident", rm.sum = FALSE, rm.pct = FALSE) {
if (!inherits(object, "Seurat")) stop("`object` must be a Seurat object.")
if (!missing(ident)) object <- SetIdent(object, value = ident)
group.by <- as.character(group.by)
missing.cols <- setdiff(group.by, colnames([email hidden]))
if(length(missing.cols)) stop(sprintf("'%s' not in meta.data.", paste(missing.cols, collapse = ", ")))
# Base cluster totals
tbl <- table([email hidden])
stats <- data.frame(cluster = names(tbl), count = as.numeric(tbl), pct = as.numeric(prop.table(tbl)*100))
# Function to compute counts and percent per group.by
compute_group <- function(g) {
tab <- table([email hidden], [email hidden][[g]])
df <- as.data.frame.matrix(tab) %>% tibble::rownames_to_column("cluster")
pct <- prop.table(tab, margin = 2) * 100
pct.df <- as.data.frame.matrix(pct) %>% tibble::rownames_to_column("cluster")
colnames(pct.df)[-1] <- paste0(colnames(pct.df)[-1], ".", g, ".pct")
cbind(df[,-1], pct.df[,-1])
}
if(length(group.by) > 0) {
stats <- cbind(stats, purrr::map_dfc(group.by, compute_group))
}
stats <- janitor::adorn_totals(stats, "row")
if(rm.sum) stats <- stats[1:(nrow(stats)-1), ]
if(rm.pct) stats <- stats[, !grepl("\\.pct$", colnames(stats))]
stats
}
#' Automatically Select the Optimal Number of PCs
#'
#' This function computes an optimal number of principal components (PCs)
#' to use for downstream analysis based on variance explained. It supports
#' any dimensional reduction stored in a Seurat object (e.g., PCA, Harmony).
#'
#' @param object A Seurat object.
#' @param reduction.name Name of the reduction slot to inspect
#' (default: \code{"pca"}). Must exist in \code{object[[reduction.name]]}.
#'
#' @details
#' The selection heuristic uses two criteria:
#'
#' \itemize{
#' \item \strong{co1}: First PC where cumulative variance exceeds 90\% AND
#' the PC-specific variance drops below 5\%.
#'
#' \item \strong{co2}: Last PC where the drop in variance from the previous PC
#' is greater than 0.1\%.
#' }
#'
#' The function returns the larger of the two values, ensuring sufficient
#' dimensionality for clustering and UMAP.
#'
#' @return An integer giving the recommended number of dimensions.
#'
#' @examples
#' \dontrun{
#' pcs <- get_pcs(object, reduction.name = "pca")
#' }
#'
#' @export
#'
get_pcs <- function(object, reduction.name="pca") {
# Check seurat object
if(class(object) != "Seurat") {
message("WARNING: this rds file does not contain a Seurat object! STOP RUNNING THIS SCRIPT")
message("Check the data type by running:")
message("class(obj)")
stop()
}
# Determine percent of variation associated with each PC
pct <- object[[reduction.name]]@stdev / sum(object[[reduction.name]]@stdev)*100
# Calculate cumulative percents for each PC
cumu <- cumsum(pct)
# Determine which PC exhibits cumulative percent greater than 90% and %
# variation associated with the PC as less than 5
co1 <- which(cumu > 90 & pct < 5)[1]
co1
# Determine the difference between variation of PC and subsequent PC and
# selecting last point where change of % of variation is more than 0.1%.
co2 <- sort(which((pct[1:length(pct) - 1] - pct[2:length(pct)]) > 0.1), decreasing = T)[1] + 1
# Minimum of the two calculation
#pcs <- min(co1, co2)
c(co1, co2)
pcs <- max(co1, co2)
#pcs <- max(co1, co2)
return(pcs)
}
#' Compute Unique Gene Expression Statistics
#'
#' This function evaluates a set of genes in a Seurat object and returns:
#' \enumerate{
#' \item The number of selected genes expressed by each cell.
#' \item The number of cells expressing each gene.
#' }
#'
#' @param object A Seurat object containing an RNA assay.
#' @param gene.list Character vector of gene names to evaluate.
#'
#' @return A list containing:
#' \describe{
#' \item{\code{per.cell}}{A data frame with columns \code{Cell} and \code{Unique}.}
#' \item{\code{per.gene}}{A data frame with columns \code{gene} and \code{cells.expressing}.}
#' }
#'
#' @export
get_unique_gene_table <- function(object, gene.list, assay = "RNA", layer = "counts") {
if (!inherits(object, "Seurat"))
stop("`object` must be a Seurat object.")
if (!assay %in% names(object@assays))
stop("Assay not found in object.")
assay_obj <- object[[assay]]
# Access the matrix via @layers (Seurat v5)
if (!layer %in% names(assay_obj@layers))
stop(paste0("Layer '", layer, "' not found in the assay."))
mat <- assay_obj@layers[[layer]]
# Ensure dimnames are set
if (is.null(rownames(mat))) rownames(mat) <- rownames(assay_obj)
if (is.null(colnames(mat))) colnames(mat) <- colnames(object)
genes.present <- intersect(gene.list, rownames(mat))
if (length(genes.present) == 0)
stop("None of the provided genes exist in the object.")
mat <- mat[genes.present, , drop = FALSE]
mat <- mat[Matrix::rowSums(mat > 0) > 0, , drop = FALSE]
per.cell <- data.frame(
Cell = colnames(mat),
Unique = Matrix::colSums(mat > 0)
)
per.gene <- data.frame(
gene = rownames(mat),
cells.expressing = Matrix::rowSums(mat > 0)
)
per.gene <- per.gene[order(per.gene$cells.expressing), ]
list(per.cell = per.cell, per.gene = per.gene)
}
run_decontx <- function(
input,
assay_name = "RNA",
verbose = TRUE,
...
) {
library(Seurat)
library(SingleCellExperiment)
library(scater)
library(celda)
library(Matrix)
# Load input
if (inherits(input, "Seurat")) {
if (verbose) cat("Input is a Seurat object\n")
sce <- as.SingleCellExperiment(input)
# Extract original UMAP from Seurat
if ("umap" %in% names(input@reductions)) {
reducedDims(sce)$UMAP_before <- Embeddings(input, "umap")
if (verbose) cat("UMAP_before extracted from Seurat\n")
} else {
reducedDims(sce)$UMAP_before <- NULL
if (verbose) cat("No UMAP found in Seurat\n")
}
} else if (inherits(input, "SingleCellExperiment")) {
sce <- input
if (verbose) cat("Input is SingleCellExperiment\n")
# Preserve existing UMAP
if ("UMAP" %in% names(reducedDims(sce))) {
reducedDims(sce)$UMAP_before <- reducedDims(sce)$UMAP
if (verbose) cat("UMAP_before extracted from SCE\n")
} else {
reducedDims(sce)$UMAP_before <- NULL
if (verbose) cat("No UMAP present in SCE\n")
}
} else if (is.matrix(input) || inherits(input, "dgCMatrix")) {
if (verbose) cat("Input is counts matrix\n")
sce <- SingleCellExperiment(assays=list(counts=input))
reducedDims(sce)$UMAP_before <- NULL
} else {
stop("Input must be Seurat, SCE, or matrix.")
}
# Run decontx
if (verbose) cat("Running DecontX...\n")
sce <- decontX(sce, ...)
cleaned <- decontXcounts(sce)
cleaned <- round(cleaned)
# Compute umap after
if (verbose) cat("Computing UMAP_after from DecontX-cleaned counts...\n")
sce <- runUMAP(sce, exprs_values = "decontXcounts")
reducedDims(sce)$UMAP_after <- reducedDims(sce)$UMAP
reducedDims(sce)$UMAP <- NULL # avoid confusion
# Return both
return(list(
sce = sce,
UMAP_before = reducedDims(sce)$UMAP_before,
UMAP_after = reducedDims(sce)$UMAP_after,
cleaned_counts = cleaned
))
}
save_heatmap <- function(
ht,
filename,
width = 8,
height = 6,
formats = c("pdf"),
dpi = 600
) {
stopifnot(!missing(ht))
stopifnot(is.character(filename))
stopifnot(is.numeric(width), is.numeric(height))
for (fmt in formats) {
fmt <- tolower(fmt)
file_out <- paste0(filename, ".", fmt)
if (fmt == "pdf") {
pdf(
file_out,
width = width,
height = height,
useDingbats = FALSE
)
} else if (fmt == "png") {
png(
file_out,
width = width,
height = height,
units = "in",
res = dpi
)
} else if (fmt == "tiff") {
tiff(
file_out,
width = width,
height = height,
units = "in",
res = dpi,
compression = "lzw"
)
} else if (fmt %in% c("jpg", "jpeg")) {
jpeg(
file_out,
width = width,
height = height,
units = "in",
res = dpi,
quality = 100
)
} else {
stop("Unsupported format: ", fmt)
}
# Draw AFTER opening device
ComplexHeatmap::draw(ht)
dev.off()
}
invisible(TRUE)
}
#' Fix matrix names in Seurat v5 objects (robust)
#'
#' Ensures that all matrices in a Seurat RNA assay (`layers`, `data`, `scale.data`)
#' have proper row and column names. Skips slots that do not exist or are empty.
#'
#' @param sobj A Seurat object
#' @return A Seurat object with all matrix names fixed
fix_seurat_matrix_names <- function(sobj) {
assay <- sobj[["RNA"]]
for (slot.name in c("layers", "data", "scale.data")) {
# Handle layers specially
if (slot.name == "layers" && !is.null(assay@layers$counts)) {
mat <- assay@layers$counts
if (!is.null(mat) && length(dim(mat)) == 2) {
if (is.null(rownames(mat))) rownames(mat) <- rownames(assay)
if (is.null(colnames(mat))) colnames(mat) <- colnames(sobj)
assay@layers$counts <- mat
}
}
# Other slots
else if (slot.name %in% slotNames(assay) && !is.null(slot(assay, slot.name))) {
mat <- slot(assay, slot.name)
if (!is.null(mat) && length(dim(mat)) == 2) {
if (is.null(rownames(mat))) rownames(mat) <- rownames(assay)
if (is.null(colnames(mat))) colnames(mat) <- colnames(sobj)
slot(assay, slot.name) <- mat
}
}
}
sobj[["RNA"]] <- assay
return(sobj)
}
savefig <- function(
filename,
fig = NULL,
layout = c("single", "double"),
type = c("ggplot", "heatmap", "baseplot"),
format = c("pdf", "png", "tiff"),
width = NULL,
height = NULL,
dpi = 600
){
layout <- match.arg(layout)
type <- match.arg(type)
formats <- format # allow vector of formats
# Nature Communications sizing defaults
default_width <- ifelse(layout == "single", 3.5, 7.2)
default_height <- ifelse(type %in% c("heatmap", "baseplot"), 8, 5)
width <- width %||% default_width
height <- height %||% default_height
for(fmt in formats) {
out_file <- filename
if(!grepl(paste0("\\.", fmt, "$"), out_file))
out_file <- paste0(out_file, ".", fmt)
cat("Saving ", out_file,
"\n layout=", layout,
"\n type=", type,
"\n size=", width, "×", height, " inches @ DPI=", dpi,
"\n format=", fmt, "\n\n")
# PDF
if(fmt == "pdf") {
pdf(out_file, width = width, height = height, useDingbats = FALSE)
if(type == "ggplot") print(fig)
if(type %in% c("heatmap", "baseplot")) eval(fig)
dev.off()
next
}
# PNG / TIFF
if(fmt %in% c("png", "tiff")) {
if(type == "ggplot") {
ggsave(out_file, plot = fig,
width = width, height = height,
units = "in", dpi = dpi)
} else if(type %in% c("heatmap", "baseplot")) {
if(fmt == "png")
png(out_file, width = width, height = height, units = "in", res = dpi)
if(fmt == "tiff")
tiff(out_file, width = width, height = height, units = "in",
res = dpi, compression = "lzw")
eval(fig)
dev.off()
}
}
}
invisible(filename)
}
#' Colorize Values Using Various Color Palettes
#'
#' This function assigns colors to unique values in a vector using a variety of distinct and continuous color palettes.
#' It supports both discrete and continuous color schemes, including RColorBrewer, Viridis, and custom palettes.
#'
#' @param x A vector of values to be colorized.
#' @param theme A string specifying the color theme. Available themes:
#' * **Distinct palettes**: `"alphabet"`, `"tube"`, `"tol"`, `"glasbey"`, `"tableau"`, `"bright"`, `"pastel"`, `"cList"``
#' * **RColorBrewer palettes**:
#' - *Sequential*: `"Blues"`, `"BuGn"`, `"BuPu"`, `"GnBu"`, `"Greens"`, `"Greys"`, `"Oranges"`, `"OrRd"`, `"PuBu"`,
#' `"PuBuGn"`, `"PuRd"`, `"Purples"`, `"RdPu"`, `"Reds"`, `"YlGn"`, `"YlGnBu"`, `"YlOrBr"`, `"YlOrRd"`
#' - *Qualitative*: `"Accent"`, `"Dark2"`, `"Paired"`, `"Pastel1"`, `"Pastel2"`, `"Set1"`, `"Set2"`, `"Set3"`
#' - *Divergent*: `"BrBG"`, `"PiYG"`, `"PRGn"`, `"PuOr"`, `"RdBu"`, `"RdGy"`, `"RdYlBu"`, `"RdYlGn"`, `"Spectral"`
#' * **Viridis palettes**: `"viridis"`, `"magma"`, `"plasma"`, `"cividis"`, `"inferno"`
#' * **Random colors**: `"random"`
#' @param bgval A value in `x` that should be assigned a background color (`"#CCCCCC"` by default). Default is `NULL` (no special background value).
#' @param n_colors Number of colors to generate for continuous palettes. Default is `100`.
#'
#' @return A named vector of hex color codes, where names correspond to unique values in `x`.
#'
#' @examples
#' # Apply Viridis color palette
#' colorize(1:10, theme = "viridis")
#'
#' # Apply RColorBrewer "Blues" palette
#' colorize(letters[1:10], theme = "Blues")
#'
#' # Use a distinct color palette (Glasbey)
#' colorize(1:15, theme = "glasbey")
#'
#' @import viridis
#' @import RColorBrewer
#' @export
colorize <- function(x, theme = "alphabet", bgval = NULL, n_colors = 100, reverse = FALSE) {
require(viridis)
require(RColorBrewer)
# Define distinct color palettes
color_alphabet <- matrix(c(
240, 163, 255, 0, 117, 220, 153, 63, 0, 76, 0, 92, 0, 92, 49, 43, 206, 72, 255, 204, 153,
128, 128, 128, 148, 255, 181, 143, 124, 0, 157, 204, 0, 194, 0, 136, 0, 51, 128, 255, 164, 5,
255, 168, 187, 66, 102, 0, 255, 0, 16, 94, 241, 242, 0, 153, 143, 224, 255, 102, 116, 10, 255,
153, 0, 0, 255, 255, 128, 255, 255, 0, 255, 80, 5
), ncol = 3, byrow = TRUE) / 256
tube_colors <- c(
"#B36305", "#E32017", "#FFD300", "#00782A", "#F3A9BB",
"#A0A5A9", "#9B0056", "#000000", "#003688", "#0098D4",
"#95CDBA", "#00A4A7", "#EE7C0E", "#84B817", "#E21836",
"#7156A5"
)
tol_colors <- c(
"#332288", "#88CCEE", "#44AA99", "#117733", "#999933",
"#DDCC77", "#CC6677", "#882255", "#AA4499", "#DDDDDD",
"#E69F00", "#56B4E9"
)
glasbey_colors <- c(
"#FF0000", "#00FF00", "#0000FF", "#FFFF00", "#FF00FF",
"#00FFFF", "#800000", "#808000", "#008000", "#800080",
"#808080", "#C0C0C0", "#008080", "#000080", "#FFA500"
)
tableau_colors <- c(
"#1F77B4", "#FF7F0E", "#2CA02C", "#D62728", "#9467BD",
"#8C564B", "#E377C2", "#BCBD22", "#17BECF",
"#AEC7E8", "#FFBB78", "#98DF8A", "#FF9896", "#C5B0D5",
"#C49C94", "#F7B6D2", "#DBDB8D", "#9EDAE5"
)
bright_colors <- c(
"#E6194B", "#3CB44B", "#FFE119", "#4363D8", "#F58231",
"#911EB4", "#42D4F4", "#F032E6", "#BFEF45", "#FABEBE"
)
pastel_colors <- c(
"#FFB3BA", "#FFDFBA", "#FFFFBA", "#BAFFC9", "#BAE1FF",
"#E6C0E9", "#D9C3A1", "#C4E5F5", "#F6D7A7", "#D1E8E2"
)
cList = list(c("grey85","#FFF7EC","#FEE8C8","#FDD49E","#FDBB84",
"#FC8D59","#EF6548","#D7301F","#B30000","#7F0000"),
c("#4575B4","#74ADD1","#ABD9E9","#E0F3F8","#FFFFBF",
"#FEE090","#FDAE61","#F46D43","#D73027")[c(1,1:9,9)],
c("#FDE725","#AADC32","#5DC863","#27AD81","#21908C",
"#2C728E","#3B528B","#472D7B","#440154"))
# cList = c("grey85","#FFF7EC","#FEE8C8","#FDD49E","#FDBB84", "#FC8D59","#EF6548","#D7301F","#B30000","#7F0000", "#4575B4","#74ADD1","#ABD9E9","#E0F3F8","#FFFFBF",
# "#FEE090","#FDAE61","#F46D43","#D73027", "#FDE725","#AADC32","#5DC863","#27AD81","#21908C", "#2C728E","#3B528B","#472D7B","#440154")
viridis_colors <- viridis(n_colors, option = "viridis")
magma_colors <- viridis(n_colors, option = "magma")
plasma_colors <- viridis(n_colors, option = "plasma")
inferno_colors <- viridis(n_colors, option = "inferno")
cividis_colors <- viridis(n_colors, option = "cividis")
brewer_palettes <- c('Blues', 'BuGn', 'BuPu', 'GnBu', 'Greens', 'Greys', 'Oranges', 'OrRd', 'PuBu',
'PuBuGn', 'PuRd', 'Purples', 'RdPu', 'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd',
"Accent", "Dark2", "Paired", "Pastel1", "Pastel2", "Set1", "Set2", "Set3",
"BrBG", "PiYG", "PRGn", "PuOr", "RdBu", "RdGy", "RdYlBu", "RdYlGn", "Spectral")
set.seed(42)
random_colors <- rgb(runif(n_colors), runif(n_colors), runif(n_colors))
# Choose color theme
if (theme == "alphabet") {
colors <- apply(rbind(color_alphabet, 1 - (1 - color_alphabet) / 2, color_alphabet / 2), 1,
function(row) rgb(row[1], row[2], row[3]))
} else if (theme == "tube") {
colors <- tube_colors
} else if (theme == "cList") {
colors <- cList
} else if (theme == "tol") {
colors <- tol_colors
} else if (theme == "glasbey") {
colors <- glasbey_colors
} else if (theme == "tableau") {
colors <- tableau_colors
} else if (theme == "bright") {
colors <- bright_colors
} else if (theme == "pastel") {
colors <- pastel_colors
} else if (theme == "random") {
colors <- random_colors
} else if (theme %in% brewer_palettes) {
colors <- colorRampPalette(brewer.pal(8, theme))(n_colors)
} else if (theme == "viridis") {
colors <- viridis_colors
} else if (theme == "magma") {
colors <- magma_colors
} else if (theme == "plasma") {
colors <- plasma_colors
} else if (theme == "cividis") {
colors <- cividis_colors
} else if (theme == "inferno") {
colors <- inferno_colors
} else {
stop("Invalid theme.")
}
# Reverse colors if needed
if (reverse) {
colors <- rev(colors)
}
unique_vals <- unique(x)
encoded_vals <- match(x, unique_vals) - 1
colors <- colors[(encoded_vals %% length(colors)) + 1]
if (!is.null(bgval)) {
colors[x == bgval] <- "#CCCCCC"
}
names(colors) <- unique_vals
return(colors)
}
#' Safe subset of Seurat object
#'
#' Wrapper around Seurat::subset that fixes matrix names automatically.
#'
#' @param sobj Seurat object
#' @param ... Arguments passed to Seurat::subset
#' @return Subsetted Seurat object
get_subset <- function(sobj, group.by = NULL, invert = FALSE,...) {
if (!is.null(group.by)) sobj <- SetIdent(sobj, value = group.by)
sobj2 <- subset(sobj, invert = invert,...)
sobj2 <- fix_seurat_matrix_names(sobj2)
return(sobj2)
}
#' Safe merge of Seurat objects
#'
#' Wrapper around Seurat::merge that fixes matrix names automatically.
#'
#' @param x Seurat object
#' @param y Seurat object
#' @param ... Additional arguments to Seurat::merge
#' @return Merged Seurat object
get_merge <- function(obj.list) {
if (length(obj.list) == 0) stop("No Seurat objects provided")
merged <- Reduce(merge, obj.list)
merged <- JoinLayers(merged)
assay <- merged[["RNA"]]
assay@layers$counts <- GetAssayData(merged, slot = "counts")
merged[["RNA"]] <- assay
merged <- fix_seurat_matrix_names(merged)
return(merged)
}
#' Load multiple 10X/Seurat samples with optional merge and save
#'
#' @param sample.dirs Character vector of sample folder names
#' @param parent.dir Parent directory containing sample folders
#' @param use.filtered Logical, whether to use filtered matrices (default TRUE)
#' @param save.dir Directory to save individual Seurat objects (default NULL)
#' @param merge Logical, whether to merge all Seurat objects (default FALSE)
#' @return List of Seurat objects or a single merged Seurat object
#' @examples
#' objs <- load_seurat(c("FN_S1256","FN_S3478"), parent.dir = "data/processed")
load_seurat <- function(sample.dirs,
parent.dir,
use.filtered = TRUE,
save.dir = NULL,
merge = FALSE,
verbose = TRUE) {
suppressPackageStartupMessages({
library(Seurat)
library(Matrix)
})
load_one <- function(sample.name) {
log <- function(...) if (verbose) cat(..., "\n")
log("\n[INFO] Loading sample:", sample.name)
samp.path <- file.path(parent.dir, sample.name)
if (!dir.exists(samp.path)) {
stop("[ERROR] Sample folder not found: ", samp.path)
}
# Patterns
if (use.filtered) {
h5.pattern <- "filtered_feature_bc_matrix\\.h5$"
dir.pattern <- "filtered_feature_bc_matrix$"
} else {
h5.pattern <- "raw_feature_bc_matrix\\.h5$"
dir.pattern <- "raw_feature_bc_matrix$"
}
# Find H5 files
h5.files <- list.files(
samp.path,
pattern = h5.pattern,
recursive = TRUE,
full.names = TRUE
)
# Find matrix directories
dirs <- list.dirs(
samp.path,
recursive = TRUE,
full.names = TRUE
)
dirs <- dirs[grepl(dir.pattern, dirs)]
log(" Found", length(h5.files), "H5 file(s)")
log(" Found", length(dirs), "matrix folder(s)")
if (length(h5.files) == 0 && length(dirs) == 0) {
stop("[ERROR] No valid 10X input found for sample: ", sample.name)
}
# Read counts
counts <- if (length(h5.files) > 0) {
log(" → Reading H5:", basename(h5.files[1]))
Read10X_h5(h5.files[1])
} else {
log(" → Reading MTX folder:", dirs[1])
Read10X(dirs[1])
}
# Create Seurat object
log(" Creating Seurat object...")
sobj <- CreateSeuratObject(
counts = counts,
project = sample.name
)
# Optional matrix-name fix
if (exists("fix_seurat_matrix_names")) {
log(" Fixing matrix names...")
sobj <- fix_seurat_matrix_names(sobj)
}
log(" Done:",
ncol(sobj), "cells x",
nrow(sobj), "features")
return(sobj)
}
cat("[INFO] Loading samples...\n")
objs <- setNames(lapply(sample.dirs, load_one), sample.dirs)
cat("[INFO] All samples loaded\n")
# Save individual objects
if (!is.null(save.dir)) {
dir.create(save.dir, recursive = TRUE, showWarnings = FALSE)
cat("[INFO] Saving Seurat objects to:", save.dir, "\n")
for (n in names(objs)) {
saveRDS(objs[[n]], file.path(save.dir, paste0(n, ".rds")))
}
}
# Merge if requested
if (merge && length(objs) > 1) {
cat("[INFO] Merging objects...\n")
merged <- Reduce(function(x, y) merge(x, y), objs)
if (exists("fix_seurat_matrix_names")) {
merged <- fix_seurat_matrix_names(merged)
}
cat("[INFO] Merged object:",
ncol(merged), "cells x",
nrow(merged), "features\n")
return(merged)
}
return(objs)
}
#' Split a Seurat object by metadata and fix matrix names
#'
#' @param object A Seurat object
#' @param split.by Metadata column name to split by
#' @return A named list of Seurat objects with @layers$counts fixed
#' @examples
#' seurat.list <- split_objectects(merged, split.by = "BulkSample")
split_objectects <- function(object, split.by) {
require(Seurat)
objs <- Seurat::SplitObject(object, split.by = split.by)
objs <- lapply(objs, function(sobj) {
# Ensure @layers$counts exist
assay <- sobj[["RNA"]]
assay@layers$counts <- Seurat::GetAssayData(sobj, slot = "counts")
sobj[["RNA"]] <- assay
# Fix row/colnames
sobj <- fix_seurat_matrix_names(sobj)
return(sobj)
})
return(objs)
}
subset_seurat <- function(object,
group.by = NULL,
use.harmony = TRUE,
harmony.group = "orig.ident",
ndims = 50,
resolution = 0.5,
min.dist = 0.3,
spread = 1,
preprocess = TRUE,
...) {
stopifnot(inherits(object, "Seurat"))
# Set identities
if (!is.null(group.by)) {
object <- Seurat::SetIdent(object, value = object[[group.by]][,1])
}
# Subset
result <- tryCatch(subset(object, ...), error = function(e) {
stop("Subsetting failed: ", e$message)
})
if (preprocess) {
# Normalize + variable features
result <- Seurat::NormalizeData(result)
result <- Seurat::FindVariableFeatures(result)
if (length(unique(Seurat::Idents(result))) > 1) {
result <- Seurat::ScaleData(result)
}
# PCA / Harmony
if (use.harmony && length(unique([email hidden][[harmony.group]])) > 1) {
result <- harmony::RunHarmony(result, group.by.vars = harmony.group)
reduction <- "harmony"
n_reduc_dims <- ncol(result@reductions[[reduction]]@cell.embeddings)
} else {
result <- Seurat::RunPCA(result, npcs = min(ndims, 50))
reduction <- "pca"
n_reduc_dims <- ncol(result@reductions[[reduction]]@cell.embeddings)
}
dims_use <- seq_len(min(ndims, n_reduc_dims))
# Neighbors + Clusters + UMAP
graph_name <- "RNA_nn"
result <- Seurat::FindNeighbors(result, dims = dims_use, reduction = reduction, graph.name = graph_name)
result <- Seurat::FindClusters(result, resolution = resolution, graph.name = graph_name)
result <- Seurat::RunUMAP(result, dims = dims_use, reduction = reduction, min.dist = min.dist, spread = spread)
}
return(result)
}
#' Select cells from a Seurat object based on gene expression, clusters, and metadata
#'
#' @description
#' `cellpick()` selects cells from a Seurat object using positive and/or negative
#' gene expression criteria, optional cluster restrictions, and flexible metadata
#' filtering (categorical or numeric ranges).
#'
#' @param obj A Seurat object.
#' @param pos.genes Character vector of genes that should be expressed.
#' @param neg.genes Character vector of genes that must NOT be expressed (default NULL).
#' @param slot Assay slot to use for expression values (default "data").
#' @param expr.thresh Numeric expr.thresh above which a gene is considered expressed (default 0).
#' @param min.genes Minimum number of `pos.genes` that must be expressed
#' (default 1; use length(pos.genes) to require all).
#' @param clusters Optional vector of cluster IDs to retain (default NULL).
#' @param cluster.col Metadata column containing cluster identities
#' (default NULL; uses `Idents(obj)`).
#' @param meta.col Metadata column to filter on (default NULL).
#' @param meta.vals Values or numeric range (length 2) used to filter `meta.col`.
#' @param return.obj Logical; if TRUE, return a Seurat object subset
#' instead of cell names (default FALSE).
#' @param verbose Logical; print filtering diagnostics (default FALSE).
#'
#' @return
#' Character vector of selected cell names, or a Seurat object if `return.obj = TRUE`.
#'
#' @export
#'
#' @examples
#' # Olfactory HBC selection
#' cells <- cellpick(
#' obj,
#' pos.genes = c("TP63", "KRT5", "KRT14"),
#' neg.genes = c("KRT13"),
#' min.genes = 2,
#' cluster.col = "seurat_clusters",
#' clusters = 16,
#' meta.col = "group",
#' meta.vals = c("PCW10", "PCW12"),
#' verbose = TRUE
#' )
cellpick <- function(
obj,
pos.genes,
neg.genes = NULL,
slot = "data",
expr.thresh = 0,
min.genes = 1,
clusters = NULL,
cluster.col = NULL,
meta.col = NULL,
meta.vals = NULL,
return.obj = FALSE,
verbose = FALSE
) {
stopifnot(inherits(obj, "Seurat"))
stopifnot(length(pos.genes) >= 1)
# Expression matrix
mat <- Seurat::GetAssayData(obj, slot = slot)
missing <- setdiff(pos.genes, rownames(mat))
if (length(missing) > 0) {
stop("Positive gene(s) not found: ", paste(missing, collapse = ", "))
}
# 2. Positive gene filtering
expr.pos <- mat[pos.genes, , drop = FALSE] > expr.thresh
keep.pos <- Matrix::colSums(expr.pos) >= min.genes
cells <- colnames(mat)[keep.pos]
if (verbose) message("After positive genes: ", length(cells))
# Negative gene filtering
if (!is.null(neg.genes)) {
missing.neg <- setdiff(neg.genes, rownames(mat))
if (length(missing.neg) > 0) {
stop("Negative gene(s) not found: ", paste(missing.neg, collapse = ", "))
}
expr.neg <- mat[neg.genes, , drop = FALSE] > expr.thresh
keep.neg <- Matrix::colSums(expr.neg) == 0
cells <- intersect(cells, colnames(mat)[keep.neg])
if (verbose) message("After negative genes: ", length(cells))
}
# Cluster restriction
if (!is.null(clusters)) {
clust <- if (is.null(cluster.col)) {
as.character(Seurat::Idents(obj))
} else {
if (!cluster.col %in% colnames([email hidden])) {
stop("cluster.col not found in [email hidden]")
}
as.character([email hidden][[cluster.col]])
}
names(clust) <- colnames(obj)
cells <- intersect(cells, names(clust)[clust %in% clusters])
if (verbose) message("After cluster filter: ", length(cells))
}
# Metadata filtering (generalized)
if (!is.null(meta.col)) {
if (!meta.col %in% colnames([email hidden])) {
stop("meta.col not found in [email hidden]")
}
if (is.null(meta.vals)) {
stop("meta.vals must be provided when meta.col is used")
}
meta <- [email hidden][[meta.col]]
names(meta) <- colnames(obj)
# Attempt numeric coercion for character metadata
meta.num <- suppressWarnings(as.numeric(gsub("[^0-9.-]", "", meta)))
if (all(!is.na(meta.num)) && length(meta.vals) == 2) {
keep.meta <- meta.num >= meta.vals[1] & meta.num <= meta.vals[2]
} else {
keep.meta <- meta %in% meta.vals
}
cells <- intersect(cells, names(meta)[keep.meta])
if (verbose) message("After metadata filter: ", length(cells))
}
# Return
if (return.obj) {
return(subset(obj, cells = cells))
}
return(cells)
}
#'#' Compute OR Gene Dominance Scores in Single-Cell RNA-seq Data
#'
#' This function computes a dominance score for olfactory receptor (OR) genes in each cell
#' of a Seurat object. The dominance score quantifies how strongly the top-expressed OR
#' dominates over the second-most expressed OR, adjusted for sparse single-cell expression.
#'
#' @param object A Seurat object containing single-cell RNA-seq data.
#' @param or_genes A character vector of OR gene names to consider.
#' @param assay Character; the assay name in Seurat object (default: "RNA").
#' @param slot Character; the assay slot to use (default: "data").
#' @param epsilon Numeric; a small number to prevent division by zero (default: 1e-9).
#' @param scale_factor Numeric; multiplier for top OR expression to scale dominance (default: 10).
#' @param exponent Numeric; exponent to raise the scaled log expression (default: 1).
#'
#' @return A data frame with one row per cell containing:
#' \describe{
#' \item{top_expr}{Expression of the top-expressed OR gene.}
#' \item{second_expr}{Expression of the second-highest OR gene.}
#' \item{dominance_score}{Computed dominance score.}
#' \item{cell}{Cell barcode or identifier.}
#' \item{n_OR_expressed}{Number of OR genes with nonzero expression in the cell.}
#' \item{PCW}{PCW metadata from the Seurat object.}
#' \item{cluster}{Cell identity from Seurat object.}
#' \item{cell_group}{Coarse cell grouping (mOSN, iOSN, Others).}
#' \item{dominance_bin}{Binned dominance score (0–1, 1–2, …, >6).}
#' }
#'
#' @examples
#' \dontrun{
#' Idents(sub) <- "ann_level_2"
#' or_genes <- intersect(ORs$Gene_name, rownames(sub))
#' df <- compute_OR_dominance(sub, or_genes)
#' head(df)
#' table(df$dominance_bin)
#' }
#'
#' @export
#'
compute_OR_dominance <- function(
object,
or_genes,
assay = "RNA",
slot = "data",
epsilon = 1e-9,
scale_factor = 10,
exponent = 1
){
# --- Load required packages ---
require(dplyr)
require(Seurat)
require(rlang)
# --- Filter OR genes present in Seurat object ---
or_genes <- intersect(or_genes, rownames(object))
if(length(or_genes) == 0) stop("No OR genes found in the Seurat object.")
# --- Extract expression matrix ---
or_mat <- as.matrix(GetAssayData(object, assay = assay, slot = slot)[or_genes, , drop = FALSE])
# --- Function to compute scaled dominance per cell ---
compute_scaled_dominance <- function(x) {
x_sorted <- sort(x, decreasing = TRUE)
x1 <- x_sorted[1]
x2 <- ifelse(length(x_sorted) >= 2, x_sorted[2], 0)
if (x1 == 0) return(c(top_expr = 0, second_expr = 0, dominance_score = 0))
# --- Scaled for sparse single-cell data ---
score <- ((x1 - x2) / (x1 + x2 + epsilon)) * (log1p(scale_factor * x1) ^ exponent)
return(c(top_expr = x1, second_expr = x2, dominance_score = score))
}
# --- Apply dominance calculation to all cells ---
dom_vals <- t(apply(or_mat, 2, compute_scaled_dominance))
dom_vals <- as.data.frame(dom_vals)
# --- Ensure correct column names ---
colnames(dom_vals) <- c("top_expr", "second_expr", "dominance_score")
# --- Construct output dataframe ---
df <- dom_vals %>%
mutate(
cell = colnames(or_mat),
n_OR_expressed = colSums(or_mat > 0),
PCW = object$PCW,
cluster = Idents(object)
) %>%
mutate(
cell_group = case_when(
cluster %in% c("GBC") ~ "GBC",
cluster %in% c("INP") ~ "INP",
cluster %in% c("iOSN") ~ "iOSN",
TRUE ~ "Others"
),
dominance_bin = cut(
dominance_score,
breaks = c(0, 1, 2, 3, 4, 5, 6, Inf),
labels = c("0–1", "1–2", "2–3", "3–4", "4–5", "5–6", ">6"),
right = FALSE,
include.lowest = TRUE
)
)
return(df)
}
run_slingshot <- function(
obj,
cluster_col = "ann_level_2",
reduction = "pca",
n_pcs = NULL,
start_clust = "OHBC",
end_clusts = NULL,
approx_points = 100,
verbose = TRUE
) {
suppressPackageStartupMessages({
library(slingshot)
library(Seurat)
library(dplyr)
library(mgcv)
library(scales)
})
if (verbose) cat("Preparing clustering info...\n")
Idents(obj) <- obj[[cluster_col, drop = TRUE]]
clusters <- factor(as.character(Idents(obj)))
if (verbose) cat(paste0("Using ", reduction, " reduction...\n"))
emb <- Embeddings(obj, reduction)
if (is.null(n_pcs)) {
if (verbose) cat("Auto-detecting PCs (all available)...\n")
pca_mat <- emb
} else {
pca_mat <- emb[, seq_len(n_pcs), drop = FALSE]
}
if (verbose) cat("Running Slingshot...\n")
sds <- slingshot(
pca_mat,
clusterLabels = clusters,
start.clus = start_clust,
end.clus = end_clusts,
approx_points = approx_points
)
if (verbose) cat("Extracting pseudotime...\n")
pt <- slingPseudotime(sds)
colnames(pt) <- paste0("Lineage", seq_len(ncol(pt)))
# add pseudotime back into metadata
pt_df <- as.data.frame(pt)
obj <- AddMetaData(obj, metadata = pt_df)
if (verbose) {
cat("Done!\n")
cat("Lineages detected: ", paste(colnames(pt), collapse = ", "),"\n")
cat("Cells with pseudotime per lineage:\n")
print(colSums(!is.na(pt)))
}
return(list(
seurat = obj,
sds = sds,
pseudotime = pt
))
}
#' Split a Seurat object by metadata and fix matrix names
#'
#' @param object A Seurat object
#' @param split.by Metadata column name to split by
#' @return A named list of Seurat objects with @layers$counts fixed
#' @examples
#' seurat.list <- split_seurat_objects(merged, split.by = "BulkSample")
split_seurat <- function(object, split.by) {
require(Seurat)
objs <- Seurat::SplitObject(object, split.by = split.by)
objs <- lapply(objs, function(sobj) {
# Ensure @layers$counts exist
assay <- sobj[["RNA"]]
assay@layers$counts <- Seurat::GetAssayData(sobj, slot = "counts")
sobj[["RNA"]] <- assay
# Fix row/colnames
sobj <- fix_seurat_matrix_names(sobj)
return(sobj)
})
return(objs)
}
#' Add combined bulk/demux mapping to a Seurat object
#'
#' This function merges a combined cell-level mapping table (e.g., from souporcell, vireo, or scds)
#' into a Seurat object's metadata. Optionally, it can remove cells labeled as "doublet" or "unassigned".
#'
#' @param object A Seurat object.
#' @param combined_tsv Path to the combined mapping TSV file or a data.frame.
#' @param barcode_col Name of the barcode column in the TSV (default: "Barcode").
#' @param remove_doublets Logical, if TRUE removes cells labeled as "doublet" or "unassigned" in `BulkSample` (default: TRUE).
#' @return A Seurat object with updated metadata.
#' @examples
#' FN_S1256 <- add_bulk_mapping_to_seurat(FN_S1256, combined_tsv = "combined_results_with_bulk_mapping.tsv")
#' @export
add_bulk_mapping_to_seurat <- function(
object,
combined_tsv,
barcode_col = "Barcode",
remove_doublets = TRUE
) {
library(Seurat)
library(dplyr)
# Load combined mapping
if (is.character(combined_tsv)) {
df <- read.delim(combined_tsv, stringsAsFactors = FALSE)
} else if (is.data.frame(combined_tsv)) {
df <- combined_tsv
} else {
stop("combined_tsv must be a file path or a data.frame")
}
# Ensure barcode column exists
if (!barcode_col %in% colnames(df)) stop(paste("Column", barcode_col, "not found in combined mapping"))
# Clean barcodes to match Seurat
df$Barcode_clean <- gsub("_1$", "", df[[barcode_col]])
rownames(df) <- df$Barcode_clean
# Keep only barcodes present in Seurat
common_cells <- intersect(colnames(object), rownames(df))
if (length(common_cells) == 0) stop("No overlapping barcodes between Seurat object and mapping")
df_sub <- df[common_cells, , drop = FALSE]
# Add metadata
object <- AddMetaData(object, df_sub)
# Remove doublets/unassigned if requested
if (remove_doublets && "BulkSample" %in% colnames([email hidden])) {
bad_labels <- c("doublet", "doublets", "Doublet", "Doublets",
"unassigned", "Unassigned")
current_idents <- Idents(object)
to_remove <- intersect(bad_labels, levels(current_idents))
if (length(to_remove) > 0) {
object <- subset(object, idents = to_remove, invert = TRUE)
}
}
return(object)
}
#' Compute gene-set dominance score per cell
#'
#' Computes a dominance score for a gene set in each cell based on the
#' relative expression of the top-expressed gene versus the second-highest
#' expressed gene. Optionally scales dominance by expression magnitude,
#' bins scores, assigns QC flags, and appends selected metadata for
#' downstream analysis.
#'
#' @param obj A \code{Seurat} object.
#' @param geneset Character vector of gene names defining the gene set.
#' @param assay Assay to use (default: \code{"RNA"}).
#' @param slot Expression slot to use (default: \code{"data"}).
#' @param epsilon Small constant to avoid division by zero.
#' @param scale.log Logical; whether to scale dominance by log-transformed
#' top gene expression.
#' @param exponent Numeric exponent applied to log-scaled expression.
#' @param bins Numeric vector defining bins for adjusted dominance scores.
#' @param bin.labels Optional character labels for dominance bins.
#' @param strong.dom.thresh Threshold above which a cell is flagged as
#' strongly dominant.
#' @param group.by Metadata column used for cell grouping
#' (e.g. \code{"ann2"}).
#' @param idents Character vector of identity values to retain;
#' all others are labeled \code{"Others"}.
#' @param add.meta Character vector of additional metadata columns
#' (e.g. \code{c("sample","stage","sex")}) to append to output.
#' @param verbose Logical; whether to print informative messages.
#'
#' @details
#' For each cell, the dominance score is calculated as:
#'
#' \deqn{(Top1 - Top2) / (Top1 + Top2 + epsilon)}
#'
#' If \code{scale.log = TRUE}, the score is multiplied by
#' \code{log1p(Top1)^exponent}.
#'
#' Only the metadata columns specified in \code{group.by} and
#' \code{add.meta} are joined, ensuring minimal memory overhead.
#'
#' @return A list with:
#' \itemize{
#' \item \code{df}: Data frame with per-cell dominance metrics,
#' bin assignments, QC flags, grouping identity, and selected metadata.
#' \item \code{genes_detected}: Character vector of detected genes from
#' \code{geneset}.
#' }
#'
#' @seealso \code{\link[Seurat]{GetAssayData}}
#'
#' @examples
#' \dontrun{
#' res <- domscore(
#' obj = object,
#' geneset = OR_genes,
#' group.by = "ann2",
#' idents = c("INP", "iOSN"),
#' add.meta = c("sample", "stage", "sex"),
#' bins = c(-Inf, 0.5, 1, 1.5, 2, Inf)
#' )
#'
#' head(res$df)
#' }
#'
#' @export
#'
domscore <- function(
obj,
geneset = NULL,
assay = "RNA",
slot = "data",
epsilon = 1e-9,
scale.log = TRUE,
exponent = 1,
bins = c(-Inf, 0.1, 0.5, 1, 2, Inf),
bin.labels = NULL,
strong.dom.thresh = 1.5,
group.by = NULL,
idents = NULL,
add.meta = NULL,
verbose = TRUE
) {
require(Seurat)
require(Matrix)
require(dplyr)
stopifnot(inherits(obj, "Seurat"))
## Gene detection
all_genes <- rownames(obj[[assay]])
detected_genes <- intersect(unique(geneset), all_genes)
if (length(detected_genes) == 0)
stop("No genes from geneset found in Seurat object.")
expr_mat <- as.matrix(
GetAssayData(obj, assay = assay, leyer = slot)[
detected_genes, , drop = FALSE
]
)
## Top1 / Top2 per cell
compute_top2_cell <- function(v) {
nz <- which(v != 0)
if (length(nz) == 0) return(c(top1 = 0, top2 = 0))
if (length(nz) == 1) return(c(top1 = v[nz], top2 = 0))
ord <- order(v, decreasing = TRUE)
c(top1 = v[ord[1]], top2 = v[ord[2]])
}
top2_mat <- apply(expr_mat, 2, compute_top2_cell)
if (is.null(dim(top2_mat)))
top2_mat <- matrix(top2_mat, nrow = 2)
top1_vec <- top2_mat[1, ]
top2_vec <- top2_mat[2, ]
dominance <- (top1_vec - top2_vec) /
(top1_vec + top2_vec + epsilon)
adjust_score <- if (scale.log) {
dominance * (log1p(top1_vec) ^ exponent)
} else {
dominance
}
## Top gene & counts
get_top_gene <- function(v) {
if (all(v == 0)) return(NA_character_)
rownames(expr_mat)[which.max(v)]
}
top_gene_vec <- apply(expr_mat, 2, get_top_gene)
n_genes_expr <- Matrix::colSums(expr_mat > 0)
## Base dataframe
df <- data.frame(
cell = colnames(expr_mat),
top_gene = top_gene_vec,
top_expr = as.numeric(top1_vec),
second_expr = as.numeric(top2_vec),
dominance = as.numeric(dominance),
adjust_score = as.numeric(adjust_score),
n_genes_expr = as.integer(n_genes_expr),
stringsAsFactors = FALSE
)
## Selective metadata join
meta <- [email hidden]
meta$cell <- rownames(meta)
cols_to_add <- unique(c(group.by, add.meta))
cols_to_add <- intersect(cols_to_add, colnames(meta))
if (length(cols_to_add) > 0) {
df <- left_join(
df,
meta[, c("cell", cols_to_add), drop = FALSE],
by = "cell"
)
} else if (verbose) {
message("No metadata columns added.")
}
## Cell grouping
## Cell grouping (SAFE)
if (!is.null(group.by) && group.by %in% colnames(df)) {
group_vals <- as.character(df[[group.by]])
if (!is.null(idents)) {
df$ident <- ifelse(
group_vals %in% idents,
group_vals,
"Others"
)
} else {
df$ident <- group_vals
}
} else {
df$ident <- "Others"
if (verbose)
warning("group.by not found; all cells set to 'Others'.")
}
## Binning
if (is.null(bin.labels)) {
bin.labels <- sapply(seq_len(length(bins) - 1), function(i) {
lb <- bins[i]; ub <- bins[i + 1]
if (is.finite(lb) && is.finite(ub)) paste0(lb, "-", ub)
else if (!is.finite(lb)) paste0("<", ub)
else paste0(">=", lb)
})
}
df$adjusted_bin <- cut(
df$adjust_score,
breaks = bins,
labels = bin.labels,
include.lowest = TRUE,
right = FALSE
)
## QC flags
df$qc_flag <- "ok"
df$qc_flag[df$n_genes_expressed == 0] <- "no_gene"
df$qc_flag[df$top_expr == 0] <- "no_top_expr"
df$qc_flag[df$adjust_score >= strong.dom.thresh] <- "strong_dominant"
## Output
list(df = df, genes_detected = detected_genes)
}
#' Find Marker Genes in a Seurat Object
#' This function identifies marker genes for specified clusters in a Seurat object
#' using differential expression testing within each cluster.
#' Run Seurat::FindMarkers across clusters with multiple tests + consensus support + visualization
#'
#' @param object Seurat object
#' @param group.by Metadata column for clustering/identity
#' @param test.use Differential test(s) ("wilcox", "bimod", "roc", "t",
#' "negbinom", "poisson", "LR", "MAST" or "all")
#' @param only.pos Keep only positive markers (default TRUE)
#' @param min.pct Minimum fraction of cells expressing a gene (default 0.25)
#' @param min.diff.pct Minimum difference in pct (default -Inf)
#' @param man.logfc.threshold LogFC threshold (default 0.25)
#' @param clusters.to.exclude Clusters to skip (default none)
#' @param max.cells.per.ident Downsample max cells per cluster (default Inf)
#' @param consensus Logical, if TRUE build consensus markers across tests
#' @param consensus_min_tests Integer, minimum number of tests a gene must be significant in
#' @param alpha Adjusted p-value threshold for significance (default 0.05)
#' @param return_both Logical, return both raw + consensus results (default FALSE)
#' @param plot_type "none", "bar", or "upset" (default "none")
#' @param plot_cluster Cluster name for upset plot (default first cluster)
#' @param ... Passed to Seurat::FindMarkers
#'
#' @return Data frame(s) of markers, optionally with ggplot object
#' @export
cellmarker <- function(object,
group.by,
assay = "RNA",
features = NULL,
test.use = 'wilcox',
only.pos = TRUE,
min.pct = 0.01,
min.diff.pct = -Inf,
logfc.threshold = 0.1,
clusters.to.exclude = c(),
max.cells.per.ident = Inf,
latent.vars = NULL,
consensus = FALSE,
consensus_min_tests = 2,
alpha = 0.05,
return_both = FALSE,
plot_type = c("none","bar","upset"),
plot_cluster = NULL,
...) {
requireNamespace("ggplot2")
if ("upset" %in% plot_type) requireNamespace("UpSetR")
if (!inherits(object, "Seurat")) stop("The provided object is not a Seurat object.")
if (!missing(group.by)) object <- Seurat::SetIdent(object = object, value = group.by)
valid_tests <- c("wilcox", "bimod", "roc", "t",
"negbinom", "poisson", "LR", "MAST")
if (identical(test.use, "all")) {
tests_to_run <- valid_tests
} else {
if (!all(test.use %in% valid_tests)) stop(paste("Invalid test.use. Choose from:", paste(valid_tests, collapse = ", ")))
tests_to_run <- test.use
}
clusters.to.test <- sort(unique([email hidden]))
clusters.to.test <- setdiff(clusters.to.test, clusters.to.exclude)
results_all <- list()
for (test in tests_to_run) {
message("=== Running test: ", test, " ===")
joined <- data.frame()
for (i in seq_along(clusters.to.test)) {
cluster <- clusters.to.test[i]
message(sprintf("[ %d / %d ] %s - cluster: %s ...",
i, length(clusters.to.test), test, cluster))
markers <- Seurat::FindMarkers(
object,
ident.1 = cluster,
ident.2 = NULL,
features = features,
only.pos = only.pos,
assay = assay,
slot = "data",
test.use = test,
min.pct = min.pct,
min.cells.group = 3,
min.diff.pct = min.diff.pct,
logfc.threshold = logfc.threshold,
max.cells.per.ident = max.cells.per.ident,
latent.vars = latent.vars,
...
)
if (!is.null(markers) && nrow(markers) > 0) {
markers$cluster <- cluster
markers$gene <- rownames(markers)
markers$test.use <- test
rownames(markers) <- NULL
joined <- rbind(joined, markers)
} else {
message(sprintf("Skipping cluster '%s': no markers found.", cluster))
}
}
results_all[[test]] <- joined
}
all_results <- dplyr::bind_rows(results_all)
#- Single test shortcut-
if (length(tests_to_run) == 1 && !consensus && !return_both) {
return(all_results)
}
#- Consensus / visualization-
consensus_df <- NULL
if (consensus || return_both) {
consensus_df <- all_results %>%
dplyr::filter(p_val_adj <= alpha) %>%
dplyr::group_by(cluster, gene) %>%
dplyr::summarise(
n_tests = dplyr::n_distinct(test.use),
mean_logFC = mean(avg_log2FC, na.rm = TRUE),
min_p_val_adj = min(p_val_adj, na.rm = TRUE),
.groups = "drop"
) %>%
dplyr::filter(n_tests >= consensus_min_tests)
}
plot_type <- match.arg(plot_type)
plt <- NULL
if (plot_type == "bar") {
marker_counts <- all_results %>%
dplyr::filter(p_val_adj <= alpha) %>%
dplyr::group_by(cluster, test.use) %>%
dplyr::summarise(n_genes = dplyr::n(), .groups = "drop")
plt <- ggplot2::ggplot(marker_counts,
ggplot2::aes(x = cluster, y = n_genes, fill = test.use)) +
ggplot2::geom_bar(stat="identity", position="dodge") +
plot_theme(theme.type = "classic", x.angle = 45) +
ggplot2::labs(title="Significant markers per test per cluster",
y="Number of markers", x="Cluster")
}
if (plot_type == "upset") {
if (is.null(plot_cluster)) plot_cluster <- clusters.to.test[1]
df_upset <- all_results %>%
dplyr::filter(cluster == plot_cluster, p_val_adj <= alpha) %>%
dplyr::select(gene, test.use) %>%
dplyr::distinct()
mat <- table(df_upset$gene, df_upset$test.use) > 0
plt <- UpSetR::upset(UpSetR::fromMatrix(mat),
mainbar.y.label = paste("Overlap of markers -", plot_cluster))
}
if (return_both) return(list(raw = all_results, consensus = consensus_df, plot = plt))
if (consensus) return(list(consensus = consensus_df, plot = plt))
return(list(raw = all_results, plot = plt))
}
#' Identify Temporally Regulated Genes Along Pseudotime
#'
#' This function identifies genes whose expression varies significantly
#' along pseudotime within selected cell types. It supports Slingshot
#' pseudotime objects (`sds_obj`) or a user-provided pseudotime matrix (`pt_matrix`).
#' Genes are filtered by variability, minimum expression, pseudotime range,
#' and minimal dynamic range (logFC). Candidate genes are tested via GAM
#' (`mgcv::gam`) for significant smooth changes along pseudotime.
#'
#' @param object Seurat object containing expression data.
#' @param sds_obj Optional \code{SlingshotDataSet} object with pseudotime estimates.
#' @param pt_matrix Optional numeric matrix of pseudotime values
#' (cells × lineages). Used if \code{sds_obj} is not provided.
#' @param annotation_level Metadata column name specifying cell type labels
#' (default: `"ann_level_2"`).
#' @param selected_celltypes Vector of cell types to include (default:
#' `c("OHBC","GBC","INP","iOSN")`).
#' @param assay Assay name in Seurat object (default: `"RNA"`).
#' @param expr_slot Expression slot to use (default: `"data"`).
#' @param n_var_genes Number of highly variable genes to preselect for testing.
#' @param min_cells_expressed Minimum number of cells a gene must be expressed in.
#' @param min_pct_cells Minimum fraction of cells expressing a gene.
#' @param min_pt_range Minimum pseudotime range over which a gene is expressed.
#' @param min_logFC Minimum dynamic range (max − min expression) required.
#' @param gam_k Number of basis functions used in GAM smoothing (default: 6).
#' @param top_n Optional integer. If provided, extracts the top-n significant genes for heatmap use.
#' @param qval_cutoff Adjusted p-value (FDR) threshold for significance.
#' @param verbose Logical; print progress messages.
#'
#' @return A list containing:
#' \describe{
#' \item{lineage_results}{A list of data frames per lineage,
#' containing p-values, q-values, and significant genes.}
#' \item{heatmap_matrices}{A list of matrices per lineage:
#' raw expression, scaled expression (all significant genes),
#' scaled expression for top N genes, and ordered cell order.}
#' \item{pt_used}{The pseudotime matrix used for analysis.}
#' }
#'
#' @details
#' Filtering steps:
#' \enumerate{
#' \item Highest-variable genes preselected.
#' \item Genes must be expressed in sufficient cells.
#' \item Genes must span a minimum pseudotime window.
#' \item Genes must show sufficient dynamic range (logFC).
#' }
#'
#' GAM testing is performed using: \code{y ~ s(pt, k = gam_k)}.
#'
#' @import mgcv Seurat dplyr matrixStats tibble
#'
#' @examples
#' \dontrun{
#' result <- identify_temporal_genes(
#' object = object,
#' sds_obj = slingshot_obj,
#' selected_celltypes = c("iOSN")
#' )
#' }
#'
#' @export
#'
identify_temporal_genes <- function(object,
sds_obj = NULL,
pt_matrix = NULL,
annotation_level = "ann_level_2",
selected_celltypes = c("OHBC","GBC","INP","iOSN"),
assay = "RNA",
expr_slot = "data",
n_var_genes = 500,
min_cells_expressed = 10,
min_pct_cells = 0.05,
min_pt_range = 0.15,
min_logFC = 0.5,
gam_k = 6,
top_n = NULL, # optional
qval_cutoff = 0.05,
verbose = TRUE) {
require(mgcv); require(Seurat); require(dplyr); require(matrixStats); require(tibble)
if (is.null(sds_obj) && is.null(pt_matrix)) stop("Provide sds_obj or pt_matrix")
cells_sel <- WhichCells(object, idents = selected_celltypes)
pt_full <- if(is.null(pt_matrix)) slingPseudotime(sds_obj) else pt_matrix
common_cells <- intersect(rownames(pt_full), colnames(object))
cells_use <- intersect(common_cells, cells_sel)
pt_full <- pt_full[cells_use,, drop=FALSE]
expr_all <- as.matrix(GetAssayData(object, assay=assay, slot=expr_slot))[,cells_use, drop = FALSE]
var_feats <- rownames(expr_all)[order(rowVars(expr_all), decreasing=TRUE)][1:min(n_var_genes, nrow(expr_all))]
lineage_results <- list()
heatmap_matrices <- list()
for(li in seq_len(ncol(pt_full))){
lineage_name <- colnames(pt_full)[li]
if(verbose) message("Processing ", lineage_name, " ...")
pts <- pt_full[, li]
valid_cells <- names(pts)[!is.na(pts)]
pts <- pts[valid_cells]
expr <- expr_all[var_feats, valid_cells, drop = FALSE]
# filtering
expressed_counts <- rowSums(expr > 0)
cells_min <- max(min_cells_expressed, ceiling(min_pct_cells * length(valid_cells)))
keep1 <- names(expressed_counts)[expressed_counts >= cells_min]
if(length(keep1) == 0) {
if(verbose) message(" No genes after expression filter; skipping")
lineage_results[[lineage_name]] <- tibble()
heatmap_matrices[[lineage_name]] <- list(expr_raw = NULL, expr_scaled_top = NULL, ordered_cells = names(sort(pts)), top_genes = NULL)
next
}
expr_pt_ranges <- t(apply(expr[keep1, , drop=FALSE] > 0, 1, function(x) {
r <- range(pts[x], na.rm=TRUE)
if (any(is.na(r))) return(c(NA, NA))
r
}))
expr_pt_range <- expr_pt_ranges[,2] - expr_pt_ranges[,1]
keep2 <- names(expr_pt_range)[!is.na(expr_pt_range) & expr_pt_range >= min_pt_range]
if(length(keep2) == 0) {
if(verbose) message(" No genes after pseudotime-range filter; skipping")
lineage_results[[lineage_name]] <- tibble()
heatmap_matrices[[lineage_name]] <- list(expr_raw = NULL, expr_scaled_top = NULL, ordered_cells = names(sort(pts)), top_genes = NULL)
next
}
logfc <- matrixStats::rowMaxs(expr[keep2, , drop=FALSE]) - matrixStats::rowMins(expr[keep2, , drop=FALSE])
keep3 <- names(logfc)[!is.na(logfc) & logfc >= min_logFC]
if(length(keep3) == 0) {
if(verbose) message(" No genes after logFC filter; skipping")
lineage_results[[lineage_name]] <- tibble()
heatmap_matrices[[lineage_name]] <- list(expr_raw = NULL, expr_scaled_top = NULL, ordered_cells = names(sort(pts)), top_genes = NULL)
next
}
if(verbose) message(" Candidate genes: ", length(keep3))
# GAM testing
pvals <- setNames(rep(NA_real_, length(keep3)), keep3)
for(g in keep3){
d <- data.frame(y = as.numeric(expr[g, ]), pt = as.numeric(pts))
if(sd(d$y, na.rm = TRUE) == 0) next
m <- tryCatch(mgcv::gam(y ~ s(pt, k = gam_k), data = d), error = function(e) NULL)
if(!is.null(m)){
sst <- summary(m)$s.table
if(!is.null(sst)) {
pvals[g] <- if("p-value" %in% colnames(sst)) sst[1, "p-value"] else sst[1, ncol(sst)]
}
}
}
res_df <- tibble(id = names(pvals), raw_p = as.numeric(pvals)) %>%
filter(!is.na(raw_p)) %>%
mutate(qval = p.adjust(raw_p, method = "fdr")) %>%
arrange(qval)
# all significant genes (full table)
sig_df <- res_df %>% filter(qval <= qval_cutoff)
if(verbose) message(" Significant genes: ", nrow(sig_df))
ord_cells <- names(sort(pts))
# Always store expr_raw for the significant genes (could be zero rows)
if(nrow(sig_df) > 0) {
expr_raw_sig <- as.matrix(expr[sig_df$id, ord_cells, drop = FALSE])
# scale for inspection (not required for plotting; plot function may re-scale or smooth)
expr_scaled_all_sig <- t(scale(t(expr_raw_sig)))
expr_scaled_all_sig[is.na(expr_scaled_all_sig)] <- 0
} else {
expr_raw_sig <- NULL
expr_scaled_all_sig <- NULL
}
# If user requested top_n, compute expr_scaled_top and top_genes
if(!is.null(top_n) && nrow(sig_df) > 0) {
top_genes <- head(sig_df$id, top_n)
expr_sub <- as.matrix(expr[top_genes, ord_cells, drop = FALSE])
expr_scaled_top <- t(scale(t(expr_sub))); expr_scaled_top[is.na(expr_scaled_top)] <- 0
} else {
top_genes <- NULL
expr_scaled_top <- NULL
}
lineage_results[[lineage_name]] <- sig_df
heatmap_matrices[[lineage_name]] <- list(
expr_raw = expr_raw_sig, # raw expr for *all* significant genes
expr_scaled_top = expr_scaled_top, # scaled expr matrix for top_n (NULL if not requested)
expr_scaled_all = expr_scaled_all_sig, # scaled expr matrix for all sig genes (may be NULL)
ordered_cells = ord_cells,
top_genes = top_genes
)
}
return(list(
lineage_results = lineage_results,
heatmap_matrices = heatmap_matrices,
pt_used = pt_full
))
}
plot_theme <- function(
theme.style = c("minimal","classic","bw","test","void","dirty","gray"),
font.size = 8,
xy.val = TRUE,
x.angle = 0,
hjust = NULL,
vjust = NULL,
xlab = TRUE,
ylab = TRUE,
xy.lab = TRUE,
facet.face = "bold",
ttl.face = "bold",
txt.face = c("plain","italic","bold"),
ttl.pos = c("center","left","right"),
x.ttl = TRUE,
y.ttl = TRUE,
ticks = NULL,
line = NULL,
border = NULL,
grid.major = NULL,
grid.minor = NULL,
panel.fill = "white",
facet.bg = TRUE,
mode = c("light","dark"),
leg.pos = "right",
leg.dir = "vertical",
leg.size = 8,
leg.ttl = 8,
leg.ttl.size = 8,
leg.just = "center",
leg.ttl.text = NULL,
...
) {
require(ggplot2)
theme.style <- match.arg(theme.style)
ttl.pos <- match.arg(ttl.pos)
txt.face <- match.arg(txt.face)
mode <- match.arg(mode)
# Canonical line width (THIS FIXES YOUR PROBLEM)
lw <- 0.3
if (is.null(line)) {
line <- theme.style == "classic"
}
# Colors
if (mode == "light") {
col.txt <- "#1A1A1A"
col.grid <- "#D9D9D9"
col.panel <- panel.fill
col.strip <- "#EFEFEF"
} else {
col.txt <- "#DDDDDD"
col.grid <- "#444444"
col.panel <- "#1E1E1E"
col.strip <- "#383838"
}
# Automatic x-label alignment
if (is.null(hjust) || is.null(vjust)) {
if (x.angle == 0) { hjust <- .5; vjust <- 1 }
else if (x.angle == 45) { hjust <- 1; vjust <- 1 }
else if (x.angle == 90) { hjust <- 1; vjust <- .5 }
else if (x.angle == 270){ hjust <- 0; vjust <- .5 }
else { hjust <- 1; vjust <- 1 }
}
ttl.pos <- switch(ttl.pos, left = 0, center = .5, right = 1)
# Base theme components
base <- theme(
text = element_text(color = col.txt, size = font.size, family = "Helvetica"),
axis.text.x = element_text(color = col.txt, size = font.size),
axis.text.y = element_text(color = col.txt, size = font.size),
axis.title = element_text(size = font.size),
plot.title = element_text(
hjust = ttl.pos, face = ttl.face,
size = font.size + 1, color = col.txt
),
strip.text = element_text(face = facet.face, color = col.txt),
legend.title = element_text(size = leg.ttl.size, face = "bold"),
legend.text = element_text(size = leg.size),
legend.position = leg.pos,
legend.key.height = unit(.4, "cm"),
legend.key.width = unit(.4, "cm"),
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.key = element_blank(),
legend.box = "vertical",
legend.spacing.y = unit(0.05, "cm"),
legend.margin = margin(1,1,1,1),
...
)
# Preset themes
preset <- switch(
theme.style,
minimal = theme_minimal(base_size = font.size),
classic = theme_classic(base_size = font.size),
bw = theme_bw(base_size = font.size),
test = theme_test(base_size = font.size),
void = theme_void(base_size = font.size),
dirty = theme_minimal(base_size = font.size) +
theme(
panel.grid = element_blank(),
panel.border = element_blank(),
axis.ticks = element_blank()
),
gray = theme_gray(base_size = font.size) +
theme(
panel.background = element_rect(fill = "#EDEDED", color = NA),
panel.grid.major = element_line(color = "#CCCCCC", linewidth = lw),
panel.grid.minor = element_line(color = "#DDDDDD", linewidth = lw/2)
)
)
th <- preset + base
# Normalize panel borders for border-based themes
if (theme.style %in% c("bw", "test", "gray")) {
th <- th + theme(
panel.border = element_rect(
linewidth = lw,
color = col.txt,
fill = NA
),
axis.line = element_blank() # ⬅️ critical
)
}
# Axis lines (classic-style)
if (line && theme.style == "classic") {
th <- th + theme(
axis.line.x = element_line(color = col.txt, linewidth = lw),
axis.line.y = element_line(color = col.txt, linewidth = lw)
)
} else {
th <- th + theme(axis.line = element_blank())
}
# X-axis angle
if (xy.val && xlab) {
th <- th + theme(
axis.text.x = element_text(angle = x.angle, hjust = hjust, vjust = vjust)
)
}
# Label / tick / grid overrides
if (!xy.lab) th <- th + theme(axis.text = element_blank(), axis.ticks = element_blank())
if (!xlab) th <- th + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
if (!ylab) th <- th + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
if (!x.ttl) th <- th + theme(axis.title.x = element_blank())
if (!y.ttl) th <- th + theme(axis.title.y = element_blank())
if (!is.null(ticks)) {
th <- th + if (ticks)
theme(axis.ticks = element_line(color = col.txt, linewidth = lw))
else
theme(axis.ticks = element_blank())
}
if (!is.null(border)) {
th <- th + if (border)
theme(panel.border = element_rect(color = col.grid, fill = NA, linewidth = lw))
else
theme(panel.border = element_blank())
}
if (!is.null(grid.major)) {
th <- th + if (grid.major)
theme(panel.grid.major = element_line(color = col.grid, linewidth = lw))
else
theme(panel.grid.major = element_blank())
}
if (!is.null(grid.minor)) {
th <- th + if (grid.minor)
theme(panel.grid.minor = element_line(color = col.grid, linewidth = lw/2))
else
theme(panel.grid.minor = element_blank())
}
if (!facet.bg) {
th <- th + theme(strip.background = element_blank())
}
if (!is.null(leg.ttl.text)) {
th <- th + labs(color = leg.ttl.text)
}
# Void cleanup
# Void cleanup
if (theme.style == "void") {
th <- th + theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.line = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank(),
strip.text = element_blank(),
strip.background = element_blank()
)
line <- FALSE
ticks <- FALSE
border <- FALSE
grid.major <- FALSE
grid.minor <- FALSE
facet.bg <- FALSE
}
# Dirty theme
if (theme.style == "dirty") {
line <- TRUE
ticks <- FALSE
border <- FALSE
grid.major <- FALSE
grid.minor <- FALSE
facet.bg <- FALSE
}
th
}
#' Generate a Color Palette with Presets and Interpolation
#'
#' This function generates a vector of colors based on predefined palette presets
#' or a user-supplied set of base colors. Colors can be interpolated in multiple
#' color spaces and adjusted for saturation, lightness, and vividness.
#'
#' @param n Integer. Number of colors to generate.
#' @param preset Character. One of \code{"bright"}, \code{"pastel"}, \code{"warm"},
#' \code{"cool"}, \code{"contrast"}, \code{"earth"}, \code{"base"}, or \code{"custom"}.
#' If \code{"custom"} is chosen, \code{base_colors} must be supplied.
#' @param base_colors Character vector of HEX color codes. Required if
#' \code{preset = "custom"}.
#' @param space Character. Color interpolation space. One of \code{"Lab"},
#' \code{"rgb"}, or \code{"HCL"}.
#' @param oversample_factor Numeric. Factor by which to oversample colors before filtering.
#' Larger values create smoother gradients. Default is \code{1.3}.
#' @param remove_gray Logical. If \code{TRUE} (default), removes grayish colors
#' (low saturation).
#' @param reverse Logical. If \code{TRUE}, reverses the order of colors.
#' @param adjust_saturation Numeric multiplier for color saturation (chroma) in HCL space.
#' Default is \code{1} (no change).
#' @param adjust_lightness Numeric multiplier for lightness in HCL space.
#' Default is \code{1} (no change).
#'
#' @details
#' This function provides flexibility for both discrete and continuous color needs.
#' If the number of requested colors (\code{n}) is greater than the number of colors
#' in the base palette, the palette is interpolated in the chosen color space.
#'
#' The \code{remove_gray} option filters out low-chroma colors that appear grayish.
#'
#' @return A character vector of HEX color codes.
#'
#' @examples
#' # Get 5 bright colors
#' custom_palette(5, preset = "bright", space = "Lab")
#'
#' # Get 15 colors from the custom "base" preset
#' custom_palette(15, preset = "base", space = "HCL")
#'
#' # Use a fully custom palette
#' my_cols <- c("#123456", "#abcdef", "#ff0000")
#' custom_palette(10, preset = "custom", base_colors = my_cols)
#'
#' @import colorspace
#' @export
custom_palette <- function(n,
preset = c("base","bright", "pastel", "warm", "cool", "contrast", "earth", "custom"),
base_colors = NULL,
space = c("Lab", "rgb", "HCL"),
oversample_factor = 1.3,
remove_gray = TRUE,
reverse = FALSE,
adjust_saturation = 1,
adjust_lightness = 1) {
library(colorspace)
# Presets
presets <- list(
bright = c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00",
"#FFFF33", "#A65628", "#F781BF", "#999999"),
pastel = c("#FDB462", "#B3DE69", "#BC80BD", "#CCEBC5", "#FFED6F",
"#FB9A99", "#B2DF8A", "#CAB2D6", "#FFFFB3"),
warm = c("#8C510A", "#BF812D", "#DFC27D", "#F6E8C3", "#FDD49E",
"#F4A582", "#D6604D", "#B2182B"),
cool = c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0",
"#0571B0", "#74ADD1", "#ABD9E9", "#E0F3F8"),
contrast = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A",
"#66A61E", "#E6AB02", "#A6761D"),
earth = c("#A6611A", "#DFC27D", "#80CDC1", "#018571",
"#E5E5E5", "#F5F5F5", "#B2182B", "#D6604D"),
# base = c(
# "#E41A1C", "#57A156", "#06A5FF", "#7289da", "#8D5B96", "#CB7647", "#F38E38",
# "#F781BE", "#CC95C8", "#B27E85", "#9A6242", "#5FA3C9", "#3766A4", "#204E75",
# "#1B9D77", "#86CC84", "#D3EC90", "#FBF583", "#E7C715", "#F0A957", "#F57994",
# "#E7298A", "#A90D55", "#52587E", "#17CDD3", "#8ECDE0", "#BC6298", "#AE2373",
# "#5E4EA1", "#7E8D86", "#507C51", "#1F5917", "#BEC603", "#C5DD3B", "#A8DA83",
# "#8DD3C7"
# ),
base = c(
"#E41A1C", "#68618B", "#409388", "#57A156", "#8D5B96", "#CB7647", "#F38E38",
"#F781BE", "#CC95C8", "#B27E85", "#9A6242", "#5FA3C9", "#3766A4", "#204E75",
"#1B9D77", "#86CC84", "#D3EC90", "#FBF583", "#E7C715", "#F0A957", "#F57994",
"#E7298A", "#A90D55", "#52587E", "#17CDD3", "#8ECDE0", "#BC6298", "#AE2373",
"#5E4EA1", "#7E8D86", "#507C51", "#1F5917", "#BEC603", "#C5DD3B", "#A8DA83",
"#8DD3C7"
),
custom = NULL
)
# Pick preset
preset <- match.arg(preset)
if (preset != "custom") {
base_colors <- presets[[preset]]
} else if (is.null(base_colors)) {
stop("For preset = 'custom', you must provide base_colors.")
}
# Match space argument
space <- match.arg(space)
# Short-circuit if enough colors
if (n <= length(base_colors)) {
cols <- base_colors[1:n]
if (reverse) cols <- rev(cols)
return(cols)
}
# Oversample
extra_colors <- ceiling(n * oversample_factor)
# Interpolation
if (space %in% c("rgb", "Lab")) {
extended_colors <- grDevices::colorRampPalette(base_colors, space = space)(extra_colors)
} else if (space == "HCL") {
extended_colors <- grDevices::colorRampPalette(base_colors, space = "hcl")(extra_colors)
}
# Remove grayish
if (remove_gray) {
is_grayish <- function(col) {
rgb <- grDevices::col2rgb(col)
sd(rgb) < 15
}
extended_colors <- Filter(function(c) !is_grayish(c), extended_colors)
}
# Adjust saturation / lightness
if (adjust_saturation != 1 || adjust_lightness != 1) {
hcl_vals <- colorspace::coords(
methods::as(colorspace::hex2RGB(extended_colors), "polarLUV")
)
hcl_vals[, "C"] <- pmax(0, hcl_vals[, "C"] * adjust_saturation)
hcl_vals[, "L"] <- pmax(0, pmin(100, hcl_vals[, "L"] * adjust_lightness))
extended_colors <- colorspace::hex(colorspace::polarLUV(hcl_vals))
}
# Take first n
final_colors <- extended_colors[1:n]
if (reverse) final_colors <- rev(final_colors)
return(final_colors)
}
#' Cell map visualization for Seurat objects
#'
#' @param object Seurat object
#' @param group.by Column(s) in metadata to group cells
#' @param reduction Dimensionality reduction to use (default "umap")
#' @param dims Dimensions to plot (1,2) or 3 for 3D
#' @param shuffle Logical, shuffle cells
#' @param raster Logical, rasterize plot
#' @param alpha Point transparency
#' @param repel Logical, repel labels
#' @param n.cells Logical, show number of cells in labels
#' @param label Logical, show cluster labels
#' @param label.size Cluster label size
#' @param label.face Cluster label font face
#' @param colors Named vector of colors
#' @param figplot Logical, minimal figure plot for figure panels
#' @param no.axes Logical, hide axes
#' @param plot.ttl Plot title
#' @param legend Logical, show legend
#' @param leg.ttl Legend title
#' @param item.size Legend item size
#' @param leg.pos Legend position
#' @param leg.just Legend justification
#' @param leg.dir Legend direction
#' @param leg.size Legend font size
#' @param leg.ncol Legend number of columns
#' @param item.border Logical, border around legend items
#' @param font.size Base font size
#' @param pt.size Point size
#' @param dark Logical, dark theme
#' @param total.cells Logical, include total cells in title
#' @param threeD Logical, 3D plot
#' @param theme Theme name
#' @param facet.bg Logical, add facet background
#' @param ... Additional arguments passed to DimPlot
#' @return ggplot or plotly object
#' @export
#'
cellmap <- function(
object,
group.by = NULL,
reduction = "umap",
dims = c(1,2),
shuffle = FALSE,
raster = NULL,
raster.dpi = c(512, 512),
alpha = 1,
repel = FALSE,
n.cells = TRUE,
label = FALSE,
label.size = 3.5,
label.face = "plain",
cols = NULL,
figplot = FALSE,
no.axes = FALSE,
plot.ttl = NULL,
legend = TRUE,
leg.ttl = NULL,
leg.ttl.size = font.size,
item.size = 3.5,
leg.pos = "right",
leg.just = "center",
leg.dir = "vertical",
leg.size = 10,
leg.ncol = NULL,
item.border = TRUE,
font.size = 10,
pt.size = 0.5,
dark = FALSE,
total.cells = FALSE,
threeD = FALSE,
theme.style = "classic",
facet.bg = FALSE,
...
) {
if (!is.null(list(...)$theme)) {
theme.style <- list(...)$theme
}
# Helper: prepare object and colors
.prepare_object <- function(obj, group, cols=NULL){
stopifnot(inherits(obj,"Seurat"))
if(is.null(group)) group <- "ident"
if(group=="ident"){
[email hidden]$ident <- Idents(obj)
} else if(!group %in% colnames([email hidden])){
stop(paste("Grouping column", group, "not found."))
}
# Drop unused levels
[email hidden][[group]] <- droplevels(factor([email hidden][[group]]))
values <- as.character([email hidden][[group]])
values[is.na(values)] <- "Unknown"
[email hidden][[group]] <- factor(values)
Idents(obj) <- [email hidden][[group]]
levels_group <- levels([email hidden][[group]])
if(is.null(cols)){
cols <- custom_palette(length(levels_group))
names(cols) <- levels_group
if("Unknown" %in% levels_group) cols["Unknown"] <- "gray70"
} else {
if(is.null(names(cols))) cols <- setNames(cols[seq_along(levels_group)], levels_group)
missing <- setdiff(levels_group, names(cols))
if(length(missing)) cols[missing] <- "gray70"
cols <- cols[levels_group]
}
list(obj=obj, cols=cols, levels_group=levels_group)
}
# 3D plotting helper
.plot_3D <- function(obj, dims, cols, alpha, pt.size, label, label.size, label.face, n.cells){
emb <- obj@reductions[[reduction]]@cell.embeddings
df <- data.frame(x=emb[, dims[1]], y=emb[, dims[2]], z=emb[, dims[3]], cluster=Idents(obj))
hover_labels <- if(n.cells){
tbl <- table(df$cluster)
paste0(df$cluster, " (", tbl[as.character(df$cluster)], ")")
} else as.character(df$cluster)
p3d <- plotly::plot_ly(df, x=~x, y=~y, z=~z, color=~cluster, colors=cols,
type="scatter3d", mode="markers",
marker=list(size=pt.size, opacity=alpha, line=list(width=0)),
text=hover_labels, hoverinfo="text")
if(label){
centers <- df %>% dplyr::group_by(cluster) %>% dplyr::summarise(x=median(x), y=median(y), z=median(z))
p3d <- p3d %>% plotly::add_text(data=centers, x=~x, y=~y, z=~z, text=~cluster, textposition="top center")
}
p3d
}
if(length(group.by) > 1){
plots <- lapply(group.by, function(g){
cellmap(
object = object,
group.by = g,
shuffle = shuffle,
raster = raster,
alpha = alpha,
repel = repel,
reduction = reduction,
dims = dims,
n.cells = n.cells,
label = label,
label.size = label.size,
label.face = label.face,
cols = cols,
figplot = figplot,
plot.ttl = g,
legend = legend,
leg.ttl = g,
item.size = item.size,
leg.pos = leg.pos,
leg.just = leg.just,
leg.dir = leg.dir,
leg.ncol = leg.ncol,
font.size = font.size,
item.border = item.border,
pt.size = pt.size,
dark = dark,
total.cells = total.cells,
threeD = threeD,
theme.style = theme.style,
facet.bg = facet.bg,
...
)
})
return(patchwork::wrap_plots(plots))
}
# Prepare object & colors
prep <- .prepare_object(object, group.by, cols)
object <- prep$obj
cols <- prep$cols
levels_group <- prep$levels_group
if(is.null(leg.ncol)) leg.ncol <- if(length(levels_group) > 18) 2 else 1
# Validate reduction
if(!(reduction %in% names(object@reductions))){
stop(paste0("Reduction '", reduction, "' not found. Available: ", paste(names(object@reductions), collapse=", ")))
}
emb <- object@reductions[[reduction]]@cell.embeddings
if(max(dims) > ncol(emb)) stop("Selected dims exceed available dimensions in reduction.")
# 3D plotting
if(threeD || length(dims) == 3) return(.plot_3D(object, dims, cols, alpha, pt.size, label, label.size, label.face, n.cells))
# 2D plotting
plt <- Seurat::DimPlot(object, group.by=group.by, shuffle=shuffle, raster=raster, pt.size=pt.size,
repel=repel, alpha=alpha, reduction=reduction, dims=dims, raster.dpi=raster.dpi,...)
# reset Seurat's forced theme_classic
plt <- plt + ggplot2::theme_void()
present_levels <- levels(droplevels([email hidden]))
if (n.cells) {
cell.nb <- table([email hidden])[present_levels]
clust.lab <- paste0(present_levels, " (", cell.nb, ")")
} else {
clust.lab <- present_levels
}
cols_use <- cols[present_levels]
leg.ttl <- if(is.null(leg.ttl)) group.by else leg.ttl
plt <- plt + ggplot2::scale_color_manual(
breaks = present_levels,
labels = clust.lab,
values = cols_use
)
if (legend) {
plt <- plt & ggplot2::guides(
color = ggplot2::guide_legend(
override.aes = if(item.border)
list(size = item.size, shape = 21, color = "black",
stroke = 0.2, fill = unname(cols))
else
list(size = item.size),
ncol = leg.ncol,
title = leg.ttl,
keyheight = grid::unit(0.25,"cm"),
keywidth = grid::unit(0.25,"cm")
)
)
} else {
plt <- plt & ggplot2::guides(color = "none")
}
# Plot title with total cells
if(total.cells){
plot.ttl <- paste0(plot.ttl, " (n=", format(ncol(object), big.mark=","), ")")
}
#if(!is.null(plot.ttl)) plt <- plt + labs(title = plot.ttl)
if(!is.null(plot.ttl)) plt <- plt + labs(title = plot.ttl) else plt <- plt + labs(title = NULL)
# Set default legend title size
if (is.null(leg.ttl.size)) leg.ttl.size <- font.size
# Apply plot_theme using do.call
theme_args <- list(
theme.style = theme.style,
font.size = font.size,
leg.size = leg.size,
leg.pos = leg.pos,
leg.dir = leg.dir,
leg.ttl = leg.ttl,
leg.ttl.size = leg.ttl.size,
facet.bg = facet.bg,
mode = if(dark) "dark" else "light"
)
if(figplot){
# Warn if the user specified a non-classic theme
if(!missing(theme.style) && theme.style != "classic"){
warning(sprintf(
"figplot = TRUE ignores custom themes (theme.style = '%s'). Use figplot = FALSE for full theming.",
theme.style
))
}
# Apply plot_theme for figplot figure
plt <- plt &
do.call(plot_theme, c(
theme_args,
list(
x.ttl = FALSE,
ticks = FALSE,
line = FALSE,
border = FALSE,
grid.major = FALSE,
grid.minor = FALSE,
panel.fill = "white"
),
list(...)
))
text_col <- "black"
} else {
# Apply plot_theme normally
plt <- plt & do.call(plot_theme, c(theme_args, list(...)))
}
# Add cluster labels if requested
if(label){
umap_data <- dplyr::tibble(x=emb[, dims[1]], y=emb[, dims[2]], cluster=as.character([email hidden])) %>%
dplyr::group_by(cluster) %>% dplyr::summarise(x=median(x), y=median(y), .groups="drop")
plt <- plt + ggrepel::geom_text_repel(
data=umap_data, aes(x, y, label=cluster),
color = if(dark) "white" else "black",
fontface = label.face,
bg.color = if(dark) "#3A3A3A" else "grey95",
bg.r = 0.1, size = label.size, seed = 42
)
}
# figplot arrow axes (minimal figure)
if(figplot){
x.lab.reduc <- plt$labels$x %||% paste0(toupper(reduction), dims[1])
y.lab.reduc <- plt$labels$y %||% paste0(toupper(reduction), dims[2])
plt <- plt & Seurat::NoAxes()
L <- 0.12
axis.df <- data.frame(x0=c(0,0), y0=c(0,0), x1=c(L,0), y1=c(0,L))
axis.plot <- ggplot2::ggplot(axis.df) +
ggplot2::geom_segment(ggplot2::aes(x=x0, y=y0, xend=x1, yend=y1), linewidth=0.4, lineend="round") +
ggplot2::xlab(x.lab.reduc) + ggplot2::ylab(y.lab.reduc) +
ggplot2::coord_fixed() + ggplot2::theme_classic(base_size=font.size) +
ggplot2::theme(plot.background=ggplot2::element_rect(fill="transparent", colour=NA),
panel.background=ggplot2::element_rect(fill="transparent", colour=NA),
axis.text=ggplot2::element_blank(),
axis.ticks=ggplot2::element_blank(),
axis.line=ggplot2::element_blank(),
panel.border=ggplot2::element_blank(),
axis.title=ggplot2::element_text(size=font.size, face="plain"),
plot.margin=ggplot2::margin(0,0,0,0))
figure.layout <- c(patchwork::area(t=1,l=1,b=11,r=11), patchwork::area(t=10,l=1,b=11,r=2))
return(plt + axis.plot + patchwork::plot_layout(design=figure.layout))
}
if(!legend) plt <- plt & Seurat::NoLegend()
if(no.axes) plt <- plt & Seurat::NoAxes()
plt
}
#' Cell Dot Plot (Enhanced Seurat DotPlot)
#'
#' A cleaner, more customizable wrapper around **Seurat::DotPlot**, providing
#' improved color handling, optional dot outlines, flexible axis formatting,
#' legend placement, and theme control. Useful for visualizing gene expression
#' patterns across clusters or metadata-defined groups.
#'
#' @param object A Seurat object.
#' @param features Character vector of features (genes or metadata fields) to plot.
#' @param group.by Column in `[email hidden]` used to group cells.
#' Default: `"seurat_clusters"`.
#'
#' @param th.cols Color palette name from **RColorBrewer** used for the
#' expression gradient. Default: `"Reds"`.
#' @param rev.th.cols Logical; reverse the gradient palette. Default: FALSE.
#'
#' @param dot.scale Numeric scale factor controlling the dot size range.
#' Passed to `Seurat::DotPlot`. Default: 4.5.
#' @param dot.outline Logical; draw outlines around dots. Default: FALSE.
#'
#' @param x.angle Angle for x-axis labels (degrees). Default: 90.
#' @param vjust.x,hjust.x Vertical and horizontal justification for x labels.
#'
#' @param flip Logical; swap x and y axes using `coord_flip()`. Default: FALSE.
#'
#' @param font.size Base font size passed to internal theme helper. Default: 8.
#'
#' @param plot.title Optional plot title.
#'
#' @param leg.size Legend text size. Default: 8.
#' @param leg.pos Position of the legend (e.g., `"right"`, `"bottom"`).
#' @param leg.just Legend justification.
#' @param leg.hjust Logical; if TRUE, use a horizontal legend layout when possible.
#'
#' @param x.axis.pos Position of the x-axis (`"top"` or `"bottom"`).
#' @param theme ggplot2 theme name used by the internal theme helper.
#'
#' @param x.face,y.face Logical; italic styling for x and/or y-axis labels.
#' @param x.ttl,y.ttl Logical; italic styling for x and/or y-axis titles.
#'
#' @param ... Additional parameters passed to `Seurat::DotPlot()`.
#'
#' @details
#' This function enhances the standard Seurat dot plot by providing:
#' * Customizable Brewer color gradients
#' * Optional dot outlines
#' * Flexible axis label styling
#' * Improved legend customization and ordering
#' * Optional axis flipping
#'
#' It retains all functionality of `Seurat::DotPlot` while adding cleaner,
#' publication-ready defaults.
#'
#' @return A ggplot object.
#'
#' @examples
#' \dontrun{
#' celldot(pbmc, features = c("MS4A1","CD3D"))
#'
#' celldot(pbmc, features = c("MS4A1","CD14"), th.cols = "Blues",
#' dot.outline = TRUE, flip = TRUE)
#' }
#'
#' @export
#'
celldot <- function(object, features, group.by="seurat_clusters", th.cols="Reds",
rev.th.cols=FALSE, dot.scale=4.5, x.angle=90, vjust.x=NULL,
hjust.x=NULL, flip=FALSE, font.size=8, plot.title=NULL,
leg.size=10, leg.pos="right", leg.just="bottom", leg.hjust=FALSE,
x.axis.pos="bottom", theme="classic", x.face=FALSE, y.face=FALSE,
x.ttl=FALSE, y.ttl=FALSE, dot.outline=FALSE, ...) {
stopifnot(inherits(object,"Seurat"))
object <- Seurat::SetIdent(object, value=group.by)
features <- unique(features)
pal <- RColorBrewer::brewer.pal(9, th.cols)
if (rev.th.cols) pal <- rev(pal)
outline_col <- if (dot.outline) "gray60" else NA
outline_stroke <- if (dot.outline) 0.5 else 0
plt <- suppressWarnings({
suppressMessages({
Seurat::DotPlot(object, features=features, dot.scale=dot.scale, ...)
})
}) +
scale_color_gradientn(colors=pal, oob=scales::squish) +
geom_point(aes(size=pct.exp), shape=21, colour=outline_col, stroke=outline_stroke) +
labs(title=plot.title, color="Average\nExpression", size="Percent\nExpressed") +
plot_theme(theme=theme, font.size=font.size, x.angle=x.angle,
x.hjust=hjust.x, x.vjust=vjust.x, xy.val=TRUE, x.lab=TRUE, y.lab=TRUE, ...) +
theme(
axis.text.x = if (x.face || (flip && y.face)) element_text(face="italic") else element_text(),
axis.text.y = if (y.face || (flip && x.face)) element_text(face="italic") else element_text(),
axis.title = element_blank(),
legend.spacing.y = unit(0.05, "cm"),
legend.spacing.x = unit(0.05, "cm"),
legend.box.spacing = unit(0.05, "cm"),
legend.margin = margin(2,2,2,2)
)
if (flip) plt <- plt + coord_flip()
if (is.list(features)) plt <- plt + theme(strip.text.x=element_text(angle=45))
# Legend positioning
# Legend positioning outside plot, bottom-right
if (!is.null(leg.pos)) {
if (leg.pos == "right") {
plt <- plt + theme(
legend.position = "right",
legend.justification = c("right","bottom"),
legend.box.just = "right",
legend.box.margin = margin(0,0,0,0)
)
} else if (leg.pos == "left") {
plt <- plt + theme(
legend.position = "left",
legend.justification = c("left","bottom"),
legend.box.just = "left",
legend.box.margin = margin(0,0,0,0)
)
} else if (leg.pos == "top") {
plt <- plt + theme(
legend.position = "top",
legend.justification = c("right","top"),
legend.box.just = "right"
)
} else if (leg.pos == "bottom") {
plt <- plt + theme(
legend.position = "bottom",
legend.justification = c("right","bottom"),
legend.box.just = "right"
)
}
}
# Keep original guides for color and size
guide_color <- guide_colorbar(frame.colour="black", ticks.colour="black")
guide_size <- guide_legend(override.aes=list(shape=21, colour=outline_col, fill="black"))
guide_color$order <- 1
guide_size$order <- 2
plt <- plt + guides(color=guide_color, size=guide_size)
plt
}
#' Cell Feature Violin Plot with Statistics
#'
#' Plots expression or metadata features as violin plots for a Seurat object,
#' optionally adding median points, shared y-axis scaling, flipped axes, and
#' statistical comparisons (Wilcoxon for 2 groups, Kruskal-Wallis for >2 groups).
#'
#' @param obj A Seurat object.
#' @param features Character vector of feature names (genes or metadata columns) to plot.
#' @param ncol Number of columns in the output patchwork plot. Defaults to sqrt(#features / 1.5).
#' @param stack Logical; if TRUE, plots are stacked in a single column.
#' @param shared.y Logical; if TRUE, all violins share the same y-axis.
#' @param ttl.pos Position of subplot titles: "center", "left", or "right".
#' @param group.by Metadata column to group by. Defaults to "seurat_clusters".
#' @param split.by Optional metadata column to split violins by.
#' @param assay Assay to pull data from. Default is "RNA".
#' @param slot Slot to use for expression values. One of "data", "counts", or "scale.data".
#' @param log Logical; if TRUE, log-transform the expression values.
#' @param cols Optional named vector of colors for each group. If NULL, defaults are used.
#' @param med Logical; if TRUE, overlay median points on each violin.
#' @param med.size Size of median points if med = TRUE.
#' @param pt.size Size of jittered points. Set to 0 to hide points.
#' @param border.size Size of the violin border lines.
#' @param txtsize Base font size for titles and labels.
#' @param theme ggplot2 theme to use: "classic", "minimal", etc.
#' @param x.ang Rotation angle of x-axis labels.
#' @param leg.pos Position of legend: "none", "right", "left", etc.
#' @param title Optional overall title for the patchwork plot.
#' @param rm.subtitles Logical; if TRUE, removes individual subplot titles.
#' @param flip Logical; if TRUE, flips x and y axes.
#' @param auto.resize Logical; if TRUE, sets dynamic width/height attributes.
#' @param ylab.global Global y-axis label. Defaults to expression level.
#' @param xlab.global Global x-axis label. Defaults to blank.
#' @param pairwise Logical; if TRUE, perform pairwise comparisons between groups.
#' @param add.stats Logical; if TRUE, add p-values to plots.
#' @param show.pval Logical; if TRUE, show p-values above violins.
#' @param pval.label Character; label type for p-values, e.g., "p.signif" or "p.format".
#' @param ... Additional arguments passed to ggplot2 layers.
#'
#' @return A patchwork object containing the violin plots.
#' @examples
#' \dontrun{
#' cellvio(sub, features = c("MEG3","TP63","HES6"),
#' group.by = "ann_level_2",
#' pt.size = 0.1, ncol = 3, pairwise = TRUE,
#' txtsize = 10, show.pval = TRUE)
#' }
#' @export
#'
cellvio <- function(
obj, features,
ncol = NULL,
shared.y = FALSE,
ttl.pos = c("center", "left", "right"),
group.by = "seurat_clusters",
split.by = NULL,
stack = FALSE,
assay = "RNA",
slot = "data",
log = FALSE,
cols = NULL,
med = FALSE,
med.size = 1,
pt.size = 0,
border.size = 0.1,
style = "classic",
leg.pos = "none",
x.ang = 45,
title = NULL,
rm.subttl = FALSE,
flip = FALSE,
auto.resize = TRUE,
ylab.global = NULL,
xlab.global = NULL,
add.stats = FALSE,
show.pval = FALSE,
pairwise = FALSE,
pval.label = "p.signif",
txtsize = 10,
...
) {
stopifnot(inherits(obj, "Seurat"))
if (length(features) == 0) stop("features must be provided.")
ttl.pos <- match.arg(ttl.pos)
# determine ncol
if (is.null(ncol)) {
ncol <- if (stack) 1 else max(1, ceiling(sqrt(length(features) / 1.5)))
}
# Set identities if grouping
if (!is.null(group.by)) {
if (!group.by %in% colnames([email hidden]))
stop(paste(group.by, "not found in metadata."))
Idents(obj) <- group.by
}
# Determine which features exist
f.expr <- intersect(features, rownames(obj[[assay]]))
f.meta <- intersect(features, colnames([email hidden]))
features <- unique(c(f.expr, f.meta))
if (!length(features)) stop("No features found in assay or metadata.")
# Shared y-scale
ymax <- NULL
if (shared.y) {
vals <- c(
if (length(f.expr)) as.numeric(Seurat::GetAssayData(obj, assay, slot)[f.expr, ]),
if (length(f.meta)) as.numeric(as.matrix([email hidden][, f.meta, drop = FALSE]))
)
vals <- vals[is.finite(vals)]
if (length(vals)) ymax <- max(vals)
}
# Colors
if (is.null(cols)) {
g <- tryCatch(unique(obj[[group.by]][, 1]), error = \(e) NULL)
cols <- custom_palette(length(g))
#cols <- scales::hue_pal()(if (is.null(g)) 8 else length(g))
#cols <- ggpubr::get_palette("npg", length(g))
names(cols) <- g
}
# Helper to build a single violin with stats
vln <- function(f) {
if (log && slot == "counts") {
obj[[assay]]@data[f, ] <- log1p(obj[[assay]]@counts[f, ])
slot_use <- "data"
} else {
slot_use <- slot
}
df <- Seurat::FetchData(obj, vars = c(group.by, f))
names(df)[2] <- "value"
df[[group.by]] <- factor(df[[group.by]]) # ensure factor
# Base plot
# p <- ggplot(df, aes_string(group.by, "value", fill = group.by)) +
# geom_violin(scale = "width", color = "black", size = border.size) +
# scale_fill_manual(values = cols)
p <- suppressWarnings({
suppressMessages({Seurat::VlnPlot(
obj,
features = f,
group.by = group.by,
split.by = split.by,
assay = assay,
slot = slot,
pt.size = pt.size,
cols = cols,
...
) + scale_y_continuous(
expand = expansion(mult = c(0.05, 0.25))
)
})
})
# Points
#if (pt.size > 0) p <- p + geom_jitter(width = 0.1, size = pt.size, alpha = 0.6)
# Add statistics
if (add.stats && show.pval) {
if (nlevels(df[[group.by]]) > 1) {
df$.grp <- df[[group.by]]
if (!is.null(split.by)) {
df$.grp <- interaction(df[[group.by]], [email hidden][[split.by]], drop = TRUE)
}
y_max <- max(df$value, na.rm = TRUE)
y_step <- (ymax %||% y_max) * 0.08
if (pairwise && nlevels(df$.grp) > 1) {
cmp <- utils::combn(levels(df$.grp), 2, simplify = FALSE)
stat_df <- ggpubr::compare_means(value ~ .grp, data = df, method = "wilcox.test", comparisons = cmp)
stat_df$y.position <- y_max + seq_len(nrow(stat_df)) * y_step
p <- p + ggpubr::stat_pvalue_manual(
stat_df,
label = pval.label,
y.position = "y.position",
tip.length = 0.02,
size = txtsize * 0.25
)
} else {
# Kruskal test if not pairwise
stat_df <- ggpubr::compare_means(value ~ .grp, data = df, method = "kruskal.test")
p <- p + annotate(
"text",
x = 1,
y = y_max + y_step,
label = paste0(signif(stat_df$p, 3)),
hjust = 0,
size = txtsize * 0.25
)
}
}
}
# Titles
p <- if (!rm.subttl) p + labs(title = f) else p + labs(title = NULL)
# Make feature titles bold + italic
p <- p + theme(
plot.title = element_text(face = "bold.italic")
)
.style_layers <- function(
p,
violin_lw = 0.15,
point_size = NULL,
jitter_width = NULL
) {
for (i in seq_along(p$layers)) {
layer <- p$layers[[i]]
# Violin outline
if (inherits(layer$geom, "GeomViolin")) {
layer$aes_params$linewidth <- violin_lw
}
# Points (Seurat uses GeomPoint + position_jitterdodge)
if (inherits(layer$geom, "GeomPoint")) {
if (!is.null(point_size)) {
layer$aes_params$size <- point_size
layer$aes_params$alpha <- 0.6
}
if (!is.null(jitter_width) &&
inherits(layer$position, "PositionJitterdodge")) {
layer$position$width <- jitter_width
}
}
p$layers[[i]] <- layer
}
p
}
p <- .style_layers(
p,
violin_lw = border.size,
point_size = if (pt.size > 0) pt.size else NULL,
jitter_width = 0.08
)
# Theme & formatting
p <- p + plot_theme(style = style, txtsize = txtsize, x.ang = x.ang,
leg.pos = leg.pos, x.ttl = FALSE, ttl.pos = ttl.pos,...) +
theme(
plot.title = element_text(face = "bold.italic")
)
# Force final legend position
if (!is.null(leg.pos)) {
p <- p + theme(legend.position = leg.pos)
}
if (med) p <- p + stat_summary(fun = median, geom = "point", shape = 3, size = med.size)
if (!is.null(ymax)) p <- p + ylim(0, ymax)
if (flip) p <- p + coord_flip()
p + ylab(NULL)
}
# number of cols
if (is.null(ncol))
ncol <- max(1, ceiling(sqrt(length(features) / 1.5)))
plist <- lapply(features, vln)
total <- length(plist)
# Only show x-axis on bottom plots
bottom <- sapply(1:ncol, \(i) max(seq(i, total, by = ncol)))
for (i in seq_along(plist)) {
if (!(i %in% bottom)) {
plist[[i]] <- plist[[i]] +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
}
}
# Layout
combo <- patchwork::wrap_plots(plist, ncol = ncol) +
patchwork::plot_layout(guides = "collect")
if (!is.null(title)) {
combo <- combo +
patchwork::plot_annotation(
title = title,
theme = theme(title = element_text(face = "bold"))
)
}
# Global labels
auto_y <- switch(slot,
data = "Expression level",
counts = "Raw counts",
scale.data = "Scaled expression",
"Expression level")
ylab <- if (is.null(ylab.global)) auto_y else ylab.global
xlab <- if (is.null(xlab.global)) "" else xlab.global
plt <- cowplot::ggdraw(combo) +
cowplot::draw_label(ylab, x = -0.01, y = 0.55, angle = 90, size = txtsize) +
cowplot::draw_label(xlab, x = 0.5, y = 0.02, size = txtsize) +
theme(plot.margin = margin(15, 15, 15, 15))
# Auto resize attributes
if (auto.resize) {
ng <- length(unique(obj[[group.by]][, 1]))
attr(plt, "dynamic_width") <- 6 + ng * 0.3
attr(plt, "dynamic_height") <- 4 + length(features) * 0.25
}
plt
}
#' Plot Unique Gene Counts Per Cell (Robust Version)
#'
#' This function visualizes the number of selected genes expressed per cell.
#' It accepts either a Seurat object with a gene list or a precomputed per-cell table.
#' Missing or unexpressed genes are automatically handled with warnings.
#'
#' @param object Either a Seurat object or a data.frame with columns `Cell` and `Unique`.
#' @param gene.list Character vector of genes to evaluate (required if `object` is a Seurat object).
#' @param plot.type One of "bar", "hist", or "violin".
#' @param font.size Numeric font size.
#' @param theme Theme type passed to `plot_theme()`.
#' @param x.lab X-axis title.
#' @param y.lab Y-axis title.
#' @param ... Additional arguments passed to `plot_theme()`.
#'
#' @return A ggplot2 object.
#' @export
plot_unique_gene_counts <- function(
object,
gene.list = NULL,
plot.type = c("bar", "hist", "violin"),
font.size = 8,
theme = "classic",
x.lab = "Number of cells",
y.lab = "Number of ORs",
color = "#20679B",
...
) {
# --- Get per-cell table ---
if (inherits(object, "Seurat")) {
if (is.null(gene.list)) stop("If `object` is a Seurat object, you must supply `gene.list`.")
counts <- rownames(object[["RNA"]]@counts)
genes.present <- intersect(gene.list, counts)
if (length(genes.present) == 0) stop("None of the genes in `gene.list` exist in the object.")
if (length(genes.present) < length(gene.list)) {
warning(sprintf("Only %d/%d genes found in the object. Proceeding with available genes.",
length(genes.present), length(gene.list)))
}
tbl <- get_unique_gene_table(object, genes.present)
per.cell <- as.data.frame(tbl$per.cell)
} else {
per.cell <- as.data.frame(object)
if (!all(c("Cell", "Unique") %in% colnames(per.cell))) {
stop("`object` must be a Seurat object OR a data.frame with columns: Cell, Unique")
}
}
# --- Filter non-expressing cells ---
df <- per.cell[per.cell$Unique > 0, , drop = FALSE]
if (nrow(df) == 0) stop("No cells express the selected genes.")
df$Unique.factor <- factor(df$Unique)
plot.type <- match.arg(plot.type)
# --- Generate plot ---
plt <- switch(
plot.type,
"bar" = ggplot(df, aes(x = Cell, y = Unique.factor)) +
geom_bar(stat = "identity", colour = color) +
labs(x = x.lab, y = y.lab) +
plot_theme(theme = theme, legend.position = "none", font.size = font.size) +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()),
"hist" = ggplot(df, aes(x = Unique)) +
geom_histogram(binwidth = 1, fill = color, colour = "black") +
labs(x = y.lab, y = "Number of cells") +
plot_theme(theme = theme, font.size = font.size),
"violin" = ggplot(df, aes(x = "", y = Unique)) +
geom_violin(trim = FALSE, fill = color) +
geom_jitter(width = 0.1, alpha = 0.4) +
labs(x = "", y = y.lab) +
plot_theme(theme = theme, font.size = font.size, ...)
)
return(plt)
}
cellpct <- function(
object,
cell.col = "ann2",
group.col = "sex",
donor.col = "sample",
xttl = NULL,
yttl = "Cell type fraction per donor (log-shifted)",
pseudo = 1e-3,
jitter.size = 1.5,
alpha = 0.85,
plot_theme_fn = plot_theme,
group.colors = NULL,
donor.colors = NULL,
font.size = 10,
leg.pos = "right",
x.angle = 45,
theme.style = "classic",
leg.ncol = 1,
show.pval = TRUE,
...
) {
library(dplyr)
library(ggplot2)
library(rlang)
# 1. Compute donor-level fractions
frac_df <- [email hidden] %>%
group_by(!!sym(donor.col), !!sym(group.col), !!sym(cell.col)) %>%
summarise(n_cells = n(), .groups = "drop") %>%
group_by(!!sym(donor.col)) %>%
mutate(
total_cells = sum(n_cells),
fraction = n_cells / total_cells,
frac_log = log10(fraction + pseudo)
) %>%
ungroup()
# 2. Shift log values
min_val <- min(frac_df$frac_log, na.rm = TRUE)
frac_df <- frac_df %>%
mutate(frac_shift = frac_log - min_val)
# 3. Statistics per cell type
stat_df <- frac_df %>%
distinct(!!sym(donor.col), !!sym(group.col), !!sym(cell.col), fraction, frac_shift) %>%
group_by(!!sym(cell.col)) %>%
summarise(
n_group1 = sum(!!sym(group.col) == levels(factor(!!sym(group.col)))[1]),
n_group2 = sum(!!sym(group.col) == levels(factor(!!sym(group.col)))[2]),
p_value = if (n_group1 >= 2 && n_group2 >= 2)
wilcox.test(fraction ~ !!sym(group.col), exact = FALSE)$p.value
else NA_real_,
max_val = max(frac_shift),
.groups = "drop"
) %>%
mutate(
sig = case_when(
is.na(p_value) ~ "ns",
p_value < 0.001 ~ "***",
p_value < 0.01 ~ "**",
p_value < 0.05 ~ "*",
TRUE ~ "ns"
),
y_pos = max_val * 1.08
)
# 3b. Merge stats into source data
source_df <- frac_df %>%
left_join(stat_df %>% dplyr::select(!!sym(cell.col), p_value, sig),by = cell.col)
# 4. Donor colors
if (!is.null(donor.colors)) {
frac_df[[donor.col]] <- factor(frac_df[[donor.col]], levels = names(donor.colors))
} else {
donors <- unique(frac_df[[donor.col]])
donor.colors <- setNames(scales::hue_pal()(length(donors)), donors)
}
# 5. Base plot
p <- ggplot(frac_df, aes(x = !!sym(cell.col), y = frac_shift, fill = !!sym(group.col))) +
stat_boxplot(geom = "errorbar", width = 0.8, position = position_dodge(0.8)) +
geom_boxplot(outlier.shape = NA, position = position_dodge(0.8)) +
geom_jitter( aes(color = !!sym(donor.col)),size = jitter.size,alpha = alpha,
position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.8)) +
labs(x = xttl, y = yttl, fill = group.col) +
plot_theme_fn(font.size = font.size,x.angle = x.angle,theme.style = theme.style,leg.pos = leg.pos,...)
# 6. Add significance labels
if (show.pval) {
p <- p + geom_text(data = stat_df,aes(x = !!sym(cell.col), y = y_pos, label = sig),
inherit.aes = FALSE,size = 4)
}
# 7. Custom colors
if (!is.null(group.colors)) p <- p + scale_fill_manual(values = group.colors)
if (!is.null(donor.colors)) p <- p + scale_color_manual(values = donor.colors)
# 8. Guides (number of columns)
if (!is.null(leg.ncol)) {
p <- p + guides(
color = guide_legend(ncol = leg.ncol, override.aes = list(size = 3)),
fill = guide_legend(ncol = leg.ncol)
)
}
# 9. Force legend position after all guides/scales
if (!is.null(leg.pos)) {p <- p + theme(legend.position = leg.pos)}
return(list(plot = p, source_data = source_df))
}
#' Plot Counts or Proportions by Group
#'
#' This function generates a ggplot2 bar, point, or box plot showing counts or
#' proportions of a variable (y) across groups (x). It supports Seurat objects
#' (uses `meta.data`), custom color palettes, stacking, coordinate flipping, and facetting.
#'
#' @param data Data frame or Seurat object (if Seurat, uses meta.data)
#' @param x Grouping variable (bare name)
#' @param y Identity variable (bare name)
#' @param plot.type Plot type: "bar", "point", or "box" (default "count")
#' @param prop Logical; if TRUE, plot proportions instead of counts
#' @param prop.multi Logical; compute proportions per group for multiple facets
#' @param stack Logical; stack bars (for bar plot) instead of dodging
#' @param coord.flip Logical; flip x and y axes
#' @param colors Vector of colors; if NULL, defaults to hue palette or RColorBrewer
#' @param use.brewer Logical; use RColorBrewer palette
#' @param brew.pal Brewer palette name (default "Set1")
#' @param raster Logical; not implemented (reserved)
#' @param theme ggplot2 theme type (default "classic")
#' @param font.size Base font size for text
#' @param x.angle Rotation angle for x-axis labels
#' @param ncol Number of columns for facet_wrap (if prop.multi)
#' @param legend Logical; whether to show legend
#' @param legend.title Legend title; if NULL, uses y variable
#' @param legend.text.size Legend text size
#' @param legend.position Legend position: "right", "bottom", etc.
#' @param legend.ncol Number of columns for legend
#' @param x.title Custom x-axis title
#' @param y.title Custom y-axis title
#' @param x.lab Logical; show x-axis labels
#' @param xy.lab Logical; show both x and y axis labels
#' @param show.contour Logical; add black border around bars (bar plot)
#'
#' @return ggplot object
#' @examples
#' \dontrun{
#' cellprop(df, x = dominance.bin, y = cluster, plot.type = "bar", prop = TRUE)
#' cellprop(sub, x = cell_group, y = dominance_bin, prop.multi = TRUE)
#' }
#' @export
cellprop <- function(
data,
x, y,
plot.type = "bar",
prop = FALSE,
prop.multi = FALSE,
percent.stack = FALSE,
stack = FALSE,
coord.flip = FALSE,
x.reverse = FALSE,
y.reverse = FALSE,
colors = NULL,
use.brewer = FALSE,
brew.pal = "Set1",
raster = FALSE,
theme = "classic",
font.size = 8,
x.angle = 90,
ncol = NULL,
legend = TRUE,
legend.title = NULL,
legend.text.size = 8,
legend.position = "right",
legend.ncol = 1,
x.title = NULL,
y.title = NULL,
x.lab = TRUE,
xy.lab = TRUE,
show.contour = TRUE,
add.pval = FALSE,
pval.test = "chisq",
pval.size = 3,
...
) {
require(ggplot2)
require(dplyr)
require(RColorBrewer)
# Handle Seurat
if (inherits(data, "Seurat")) data <- [email hidden]
x_var <- rlang::enquo(x)
y_var <- rlang::enquo(y)
# Count table
stat <- data %>%
dplyr::select(!!x_var, !!y_var) %>%
dplyr::rename(.group = !!x_var, .ident = !!y_var) %>%
dplyr::group_by(.group, .ident) %>%
dplyr::summarise(.n = n(), .groups = "drop")
# proportions (per x-group)
if (prop || prop.multi || percent.stack) {
stat <- stat %>%
dplyr::group_by(.group) %>%
dplyr::mutate(.value = .n / sum(.n) * 100) %>%
dplyr::ungroup()
} else {
stat$.value <- stat$.n
}
# reverse x-axis
if (x.reverse) {
stat$.group <- factor(stat$.group, levels = rev(sort(unique(stat$.group))))
} else {
stat$.group <- factor(stat$.group, levels = sort(unique(stat$.group)))
}
stat$.ident <- factor(stat$.ident)
# colors
n.colors <- length(unique(stat$.ident))
if (is.null(colors)) {
if (use.brewer) {
colors <- colorRampPalette(brewer.pal(min(9, n.colors), brew.pal))(n.colors)
} else {
colors <- scales::hue_pal()(n.colors)
}
}
# Build main plot
p <- ggplot(stat, aes(x = .group, y = .value, fill = .ident))
# stacked or dodged bars
if (plot.type == "bar") {
p <- p + geom_bar(
stat = "identity",
position = if (stack || percent.stack) "stack" else "dodge",
color = if (show.contour) "black" else NA,
linewidth = 0.2
)
}
# point/box optional
if (plot.type == "point") p <- p + geom_point(size = 2)
if (plot.type == "box") p <- p + geom_boxplot()
# percent stacked bar → fix y-scale to 100%
if (percent.stack) {
p <- p + scale_y_continuous(limits = c(0,100), expand = expansion(mult = c(0, 0.05)))
y.title <- "Percent (%)"
}
# prop.multi facet mode
if (prop.multi) {
p <- p +
facet_wrap(~.ident, scales = "free_y", ncol = ncol) +
scale_y_continuous(labels = function(x) paste0(x, "%"))
}
# axis reversing
if (coord.flip) p <- p + coord_flip()
if (y.reverse && !percent.stack && !prop.multi) p <- p + scale_y_reverse()
# labels + theme
p <- p +
scale_fill_manual(values = colors) +
plot_theme(theme = theme, font.size = font.size, x.angle = x.angle, x.lab = x.lab, xy.lab = xy.lab, ...) +
labs(
x = if (!is.null(x.title)) x.title else rlang::as_name(x_var),
y = if (!is.null(y.title)) y.title else if (prop || prop.multi || percent.stack) "Percent (%)" else "Count",
fill = if (!is.null(legend.title)) legend.title else rlang::as_name(y_var)
) +
theme(
legend.title = element_text(size = legend.text.size),
legend.text = element_text(size = legend.text.size),
legend.position = if (legend) legend.position else "none"
) +
guides(fill = guide_legend(ncol = legend.ncol))
# ADD P-VALUES
if (add.pval) {
pval.df <- stat %>%
tidyr::pivot_wider(names_from = .ident, values_from = .n, values_fill = 0) %>%
dplyr::rowwise() %>%
mutate(
pval =
if (pval.test == "chisq")
chisq.test(c_across(!.group))$p.value
else if (pval.test == "fisher")
fisher.test(matrix(c_across(!.group), nrow = 1))$p.value
else NA_real_
) %>%
ungroup() %>%
mutate(
label = paste0("p=", signif(pval, 2)),
y.position = max(stat$.value) * 1.05
)
p <- p +
geom_text(
data = pval.df,
aes(x = .group, y = y.position, label = label),
inherit.aes = FALSE,
size = pval.size
)
}
return(p)
}
##' Pseudobulk cell-type proportion plotting with statistics
#'
#' Generate pseudobulk cell-type proportion plots from single-cell data,
#' supporting Seurat objects or precomputed data frames. The function
#' computes donor-level proportions, summarizes data as mean ± SEM,
#' performs statistical testing, and returns both the plot and
#' Nature-style source data.
#'
#' @param input A \code{Seurat} object or a data frame containing
#' donor-, stage-, and cell-level information.
#' @param donor.col Column name identifying biological replicates
#' (e.g. patient, sample, donor).
#' @param stage.col Column name identifying experimental groups or stages.
#' @param cell.col Column name identifying cell types or clusters.
#' @param prop.col Optional column name containing precomputed proportions.
#' If NULL, proportions are calculated automatically.
#' @param plot.type Plot type: \code{"bar"} (mean ± SEM with points)
#' or \code{"box"} (distribution with jittered points).
#' @param min.lmm Minimum number of donors per group required to apply
#' linear mixed models. Below this threshold, Kruskal–Wallis is used.
#' @param facet.by Variable used for faceting (typically \code{"cell"}).
#' @param cell.order Optional character vector specifying cell-type order.
#' @param palette Brewer palette name used when \code{stage.cols} is NULL.
#' @param donor.cols Optional named vector of colors for donors.
#' @param stage.cols Optional named vector of colors for stages/groups.
#' @param show.p Logical; whether to display p-values on the plot.
#' @param title Plot title.
#' @param x.lab X-axis label.
#' @param y.lab Y-axis label.
#' @param leg.pos Legend position passed to \code{plot_theme()}.
#' @param pt.size Point size for jittered donor-level points.
#' @param ncol Number of columns for facet wrapping.
#' @param ... Additional arguments passed directly to \code{plot_theme()}.
#'
#' @details
#' For Seurat objects, cell-type proportions are computed per donor
#' (pseudobulk). Summary statistics are calculated as mean ± SEM
#' across donors. Statistical testing is performed independently
#' for each cell type:
#'
#' \itemize{
#' \item Linear mixed model (LMM): \code{prop ~ stage + (1 | donor)}
#' \item Kruskal–Wallis test when donor numbers are insufficient
#' }
#'
#' The returned \code{source_data} table is suitable for direct inclusion
#' as Nature-style Source Data and includes raw proportions, donor counts,
#' summary statistics, and p-values.
#'
#' @return A list with two elements:
#' \describe{
#' \item{plot}{A \code{ggplot2} object.}
#' \item{source_data}{A data frame containing raw and summarized data,
#' statistical results, and donor counts.}
#' }
#'
#' @examples
#' res <- cellpbulk(
#' input = seurat_obj,
#' donor.col = "sample",
#' stage.col = "condition",
#' cell.col = "celltype",
#' plot.type = "bar",
#' show.p = TRUE
#' )
#'
#' res$plot
#' head(res$source_data)
#'
#' @export
cellpbulk <- function(
input,
donor.col = "sample",
stage.col = "stage",
cell.col = "ann2",
prop.col = NULL,
plot.type = c("bar", "box"),
min.lmm = 3,
facet.by = "cell",
cell.order = NULL,
palette = "Set3",
donor.cols = NULL,
stage.cols = NULL,
show.p = TRUE,
title = NULL,
x.lab = NULL,
y.lab = "Proportion (%)",
leg.pos = "right",
pt.size = 2,
label.size = 3,
ncol = NULL,
theme.style = "classic",
facet.title.size = 10,
...
) {
library(dplyr)
library(ggplot2)
library(lme4)
plot.type <- match.arg(plot.type)
# 1. Build pseudobulk dataframe
if (inherits(input, "Seurat")) {
meta <- [email hidden] %>% as.data.frame()
df <- meta %>%
mutate(
donor = .data[[donor.col]],
stage = .data[[stage.col]],
cell = .data[[cell.col]]
) %>%
group_by(donor, stage, cell) %>%
summarise(n_cells = n(), .groups = "drop") %>%
group_by(donor, stage) %>%
mutate(prop = n_cells / sum(n_cells)) %>%
ungroup()
} else {
df <- input %>%
rename(
donor = all_of(donor.col),
stage = all_of(stage.col),
cell = all_of(cell.col)
)
if (!is.null(prop.col)) {
df <- df %>% rename(prop = all_of(prop.col))
} else if (!"prop" %in% colnames(df)) {
stop("prop.col must be provided or 'prop' column must exist")
}
}
# 2. Factor ordering
if (!is.null(cell.order)) {
df$cell <- factor(df$cell, levels = intersect(cell.order, unique(df$cell)))
} else {
df$cell <- factor(df$cell)
}
# 3. Donor counts + Mean ± SEM
donor_n <- df %>%
group_by(cell, stage) %>%
summarise(n_donor = n_distinct(donor), .groups = "drop")
summary_stats <- df %>%
group_by(cell, stage) %>%
summarise(
n = n_distinct(donor),
mean = mean(prop, na.rm = TRUE),
sd = sd(prop, na.rm = TRUE),
sem = sd / sqrt(n),
.groups = "drop"
)
# 4. Statistics (LMM / KW)
stats <- lapply(levels(df$cell), function(ct) {
tmp <- df %>% filter(cell == ct)
n_ds <- donor_n %>% filter(cell == ct) %>% pull(n_donor)
if (all(n_ds >= min.lmm)) {
fit <- try(lmer(prop ~ stage + (1 | donor), data = tmp), silent = TRUE)
p <- if (!inherits(fit, "try-error")) anova(fit)$`Pr(>F)`[1] else NA
data.frame(cell = ct, method = "LMM", p.value = p)
} else if (length(unique(tmp$stage)) > 1) {
kw <- kruskal.test(prop ~ stage, data = tmp)
data.frame(cell = ct, method = "KW", p.value = kw$p.value)
} else {
data.frame(cell = ct, method = NA, p.value = NA)
}
}) %>% bind_rows()
# 5. Plot
fill.scale <- if (!is.null(stage.cols)) {
scale_fill_manual(values = stage.cols)
} else {
scale_fill_brewer(palette = palette)
}
if (plot.type == "bar") {
p <- ggplot(df, aes(stage, prop, fill = stage)) +
stat_summary(fun = mean, geom = "bar",
color = "black", linewidth = 0.25, width = 0.7) +
stat_summary(fun.data = mean_se, geom = "errorbar",
width = 0.25, linewidth = 0.3) +
geom_jitter(aes(color = donor),
color = "black",
stroke = 0.3,
shape = 21, size = pt.size,
width = 0.1, alpha = 0.8)
} else {
p <- ggplot(df, aes(stage, prop, fill = stage)) +
geom_boxplot(alpha = 0.5) +
geom_jitter(aes(color = donor),
color = "black",
shape = 21, size = pt.size,
width = 0.1, alpha = 0.8)
}
if (!is.null(donor.cols)) {
p <- p + scale_color_manual(values = donor.cols)
}
p <- p +
fill.scale +
facet_wrap(as.formula(paste("~", facet.by)),scales = "free_y",ncol = ncol) +
plot_theme(theme.style = theme.style, leg.pos = leg.pos, ...) +
labs(title = title, x = x.lab, y = y.lab) +
theme(legend.position = leg.pos, strip.text = element_text(size = facet.title.size))
# 6. P-value annotation
if (show.p) {
annot <- df %>%
group_by(cell) %>%
summarise(y = max(prop, na.rm = TRUE), .groups = "drop") %>%
left_join(stats, by = "cell") %>%
mutate(label = ifelse(is.na(p.value), "",
paste0("p=", signif(p.value, 2))))
annot$cell <- factor(annot$cell, levels = levels(df$cell))
p <- p + geom_text(
data = annot,
aes(x = 2, y = y * 1.05, label = label),
inherit.aes = FALSE,
size = label.size
)
}
# 7. Source data (Nature-ready)
source_data <- df %>%
left_join(donor_n, by = c("cell", "stage")) %>%
left_join(summary_stats, by = c("cell", "stage")) %>%
left_join(stats, by = "cell") %>%
arrange(cell, stage, donor)
return(list(
plot = p,
source_data = source_data
))
}
#' Summarize and Plot Gene Set Expression
#'
#' @description
#' Computes the total expression of a given gene set across cells, optionally summarized
#' in pseudobulk per group. Produces a bar plot showing summed expression per group/fill variable.
#'
#' @param object A Seurat object containing gene expression data.
#' @param gene.list A character vector of gene names to include.
#' @param group.var Metadata column to group cells by (x-axis). Default: "ann2".
#' @param fill.var Metadata column to use for fill colors (stacked or dodged bars). Default: "PCW".
#' @param sample.var Metadata column for pseudobulk (optional). If provided, aggregates by sample.
#' @param pseudobulk.mode Character. If "sum", sums expression per sample; if "cpm", normalizes counts per million. Default: NULL (no pseudobulk).
#' @param group.levels Optional factor levels for group.var.
#' @param fill.levels Optional factor levels for fill.var.
#' @param fill.colors Optional named vector of colors for fill levels.
#' @param expr.threshold Minimum expression to retain a gene. Default: 0.
#' @param theme ggplot2 theme. Default: "minimal".
#' @param leg.pos Legend position. Default: "right".
#' @param leg.dir Legend direction. Default: "vertical".
#' @param leg.ttl Legend title. Default: "".
#' @param font.size Base font size. Default: 10.
#' @param x.angle X-axis text angle. Default: 45.
#' @param x.lab X-axis label. Default: NULL.
#' @param y.lab Y-axis label. Default: "Total transcripts".
#' @param plot.ttl Plot title. Default: NULL.
#' @param flip Logical. If TRUE, flips coordinates. Default: FALSE.
#' @param return.data Logical. If TRUE, returns a list with plot and summarized data. Default: FALSE.
#' @param ... Additional arguments passed to `plot_theme`.
#'
#' @return A ggplot object (or list with plot and summarized data if `return.data = TRUE`).
#' @export
#'
#' @examples
#' genesum(
#' object = seurat_obj,
#' gene.list = c("OR1", "OR2"),
#' group.var = "ann2",
#' fill.var = "stage",
#' sample.var = "sample",
#' pseudobulk.mode = "cpm",
#' plot.ttl = "OR pseudobulk across development"
#' )
#'
genesum <- function(
object,
gene.list,
group.var = "ann2",
fill.var = "stage",
sample.var = NULL,
mode = c("cell", "pseudobulk"),
pb.norm = c("none", "cpm"),
group.levels = NULL,
fill.levels = NULL,
fill.colors = NULL,
expr.threshold = 0,
theme = "classic",
leg.pos = "right",
leg.dir = "vertical",
leg.ttl = "",
font.size = 10,
x.angle = 45,
x.lab = NULL,
y.lab = "Total transcripts",
plot.ttl = NULL,
flip = FALSE,
return.data = FALSE,
...
) {
require(ggplot2)
require(dplyr)
require(reshape2)
mode <- match.arg(mode)
pb.norm <- match.arg(pb.norm)
## 1. Filter genes (FAST, sparse-safe)
counts <- Seurat::GetAssayData(object, assay = "RNA", slot = "counts")
valid.genes <- intersect(gene.list, rownames(counts))
if (!length(valid.genes)) stop("No valid genes found")
counts <- counts[valid.genes, , drop = FALSE]
# Remove low-expression genes
keep.genes <- Matrix::rowSums(counts) > expr.threshold
counts <- counts[keep.genes, , drop = FALSE]
## 2. Total expression per cell (FAST)
cell_totals <- Matrix::colSums(counts)
meta <- [email hidden]
meta$total_expr <- cell_totals[colnames(object)]
## 3. Aggregation
if (mode == "cell") {
summary.df <- meta %>%
group_by(.data[[group.var]], .data[[fill.var]]) %>%
summarise(
total = sum(total_expr),
n = n(),
sd = sd(total_expr),
sem = sd / sqrt(n),
.groups = "drop"
)
}
if (mode == "pseudobulk") {
if (is.null(sample.var))
stop("sample.var must be provided for pseudobulk mode")
pb <- meta %>%
group_by(.data[[sample.var]],
.data[[group.var]],
.data[[fill.var]]) %>%
summarise(total_expr = sum(total_expr), .groups = "drop")
if (pb.norm == "cpm") {
libsize <- pb %>%
group_by(.data[[sample.var]]) %>%
summarise(lib = sum(total_expr), .groups = "drop")
pb <- pb %>%
left_join(libsize, by = sample.var) %>%
mutate(total_expr = total_expr / lib * 1e6)
}
summary.df <- pb %>%
group_by(.data[[group.var]], .data[[fill.var]]) %>%
summarise(
total = sum(total_expr),
n = n(),
sd = sd(total_expr),
sem = sd / sqrt(n),
.groups = "drop"
)
}
## 4. Factors
if (!is.null(group.levels))
summary.df[[group.var]] <- factor(summary.df[[group.var]], levels = group.levels)
if (!is.null(fill.levels))
summary.df[[fill.var]] <- factor(summary.df[[fill.var]], levels = fill.levels)
## 5. Plot
pd <- position_dodge(width = 0.7) # increase width to create space between bars within fill.var
bar_width <- 0.5 # set bar width
min_frac <- 0.05
p <- ggplot(summary.df, aes(x = .data[[group.var]], y = total, fill = .data[[fill.var]])) +
geom_bar(stat = "identity", position = pd, color = "black", linewidth = 0.25, width = bar_width) +
geom_errorbar(aes(ymin = total - pmax(sem, min_frac * total),
ymax = total + pmax(sem, min_frac * total)),
position = pd, width = 0.25, linewidth = 0.25) +
guides(fill = guide_legend(title = leg.ttl)) +
labs(x = x.lab, y = y.lab, title = plot.ttl) +
plot_theme(theme = theme, leg.pos = leg.pos, leg.dir = leg.dir,
font.size = font.size, x.angle = x.angle, leg.size = font.size, ...)
if (!is.null(fill.colors))
p <- p + scale_fill_manual(values = fill.colors)
if (flip)
p <- p + coord_flip()
if (return.data)
return(list(plot = p, data = summary.df))
return(p)
}
#'#' Highlight selected clusters on UMAP/tSNE plots
#'
#' @description
#' Plots a dimensionality reduction (UMAP/tSNE) and highlights specified clusters.
#' Can highlight clusters individually (mode = "single") or
#' simultaneously (mode = "multi").
#'
#' @param object Seurat object
#' @param cluster.names Character vector of clusters to highlight.
#' If NULL, all identities are used.
#' @param group.by Metadata column used for cluster identity.
#' @param reduction Dimensional reduction name (default: "umap")
#' @param mode "single" = one plot per cluster, "multi" = combined highlight
#' @param ncol Number of columns for "single" mode layout
#' @param custom.colors Named vector of colors corresponding to highlighted clusters.
#' @param pt.size Baseline point size for non-highlight cells
#' @param highlight.size Point size for highlighted cells
#' @param fontsize Base text size
#' @param show.cell.counts Logical, adds cell numbers to legend
#' @param background.color Color for all non-highlighted cells
#' @param label Whether to label clusters (multi mode only)
#' @param label.size Cluster label font size
#' @param plot.title Custom plot title (NULL removes the default title)
#' @param no_axes Remove axes (recommended for UMAP/tSNE)
#' @param legend.position Legend position ("none" by default)
#' @param item.size Size of legend point markers
#' @param legend.ncol Number of legend columns
#' @param legend.title Title for legend (NULL hides it)
#' @param legend.text.size Legend text size
#' @param ... Additional arguments to DimPlot
#'
#' @return A ggplot2 or patchwork object
#' @export
#'
cellmark <- function(
object,
cluster.names = NULL,
group.by = "seurat_clusters",
reduction = "umap",
mode = c("single", "multi"),
ncol = 3,
custom.colors = NULL,
pt.size = 0.1,
highlight.size = 0.5,
fontsize = 10,
show.cell.counts = FALSE,
background.color = "lightgray",
label = FALSE,
label.size = 4,
plot.title = NULL,
no.axes = TRUE,
legend.position = "none",
item.size = 3,
legend.ncol = 1,
legend.title = NULL,
legend.text.size = 10,
...
) {
require(Seurat)
require(ggplot2)
require(patchwork)
require(dplyr)
require(ggrepel)
mode <- match.arg(mode)
object <- SetIdent(object, value = group.by)
clusters <- levels(Idents(object))
if (is.null(cluster.names)) cluster.names <- clusters
valid.clusters <- cluster.names[cluster.names %in% clusters]
if (length(valid.clusters) == 0)
stop("No valid clusters found in metadata.")
if (is.null(custom.colors)) {
pal <- custom_palette(n = length(valid.clusters))
names(pal) <- valid.clusters
custom.colors <- pal
}
# MODE = SINGLE
if (mode == "single") {
plots <- lapply(valid.clusters, function(cl) {
cells <- WhichCells(object, idents = cl)
cnt <- length(cells)
title <- if (show.cell.counts) paste0(cl, " (n=", cnt, ")") else cl
p <- suppressWarnings({
suppressMessages({DimPlot(
object,
reduction = reduction,
group.by = group.by,
cells.highlight = cells,
cols = background.color,
cols.highlight = custom.colors[cl],
pt.size = pt.size,
sizes.highlight = highlight.size,
order = TRUE,
...
) })})+
ggtitle(title) +
theme(
plot.title = element_text(hjust = 0.5, size = fontsize, face = "bold"),
text = element_text(size = fontsize)
) &
NoLegend()
if (no_axes) p <- p & NoAxes()
return(p)
})
return(wrap_plots(plots, ncol = ncol))
}
# MODE = MULTI
cells.to.highlight <- CellsByIdentities(object, idents = valid.clusters)
plt <- suppressWarnings({
suppressMessages({DimPlot(
object,
reduction = reduction,
group.by = group.by,
cells.highlight = cells.to.highlight,
cols.highlight = custom.colors[valid.clusters],
cols = background.color,
pt.size = pt.size,
sizes.highlight = highlight.size,
order = TRUE,
...
) }) })
# Ensure consistent color scale ONLY once
plt <- plt +
scale_color_manual(
values = custom.colors[valid.clusters],
breaks = valid.clusters,
na.value = background.color
) +
guides(
color = guide_legend(
override.aes = list(
size = item.size,
shape = 21,
fill = unname(custom.colors[valid.clusters]),
stroke = 0.3,
color = "black"
),
ncol = legend.ncol,
title = legend.title
)
) #+
# theme(
# legend.position = legend.position,
# legend.text = element_text(size = legend.text.size),
# text = element_text(size = fontsize)
# )
plt <- plt +
theme(
legend.position = legend.position,
legend.text = element_text(size = legend.text.size),
legend.title = element_text(size = legend.text.size, face = "bold"),
legend.spacing.y = unit(1, "pt"),
legend.key.height = unit(4, "pt"),
legend.key.width = unit(4, "pt"),
legend.margin = margin(1, 1, 1, 1),
text = element_text(size = fontsize)
)
# Plot title rule
if (is.null(plot.title)) {
plt <- plt + theme(plot.title = element_blank())
} else {
plt <- plt + ggtitle(plot.title) +
theme(plot.title = element_text(
hjust = 0.5, size = fontsize + 2, face = "bold"
))
}
# Labels on highlighted clusters
if (label) {
emb <- Embeddings(object[[reduction]])
df <- data.frame(emb, cluster = Idents(object)) %>%
filter(cluster %in% valid.clusters) %>%
group_by(cluster) %>%
summarise(x = median(umap_1), y = median(umap_2), .groups = "drop")
plt <- plt +
geom_text_repel(
data = df,
aes(x = x, y = y, label = cluster),
size = label.size,
fontface = "bold"
)
}
# Cell count labels in legend
if (show.cell.counts) {
cell.nb <- table(Idents(object))
lbl <- paste0(valid.clusters, " (",cell.nb[valid.clusters], ")")
names(lbl) <- valid.clusters
plt <- plt +
scale_color_manual(
values = custom.colors[valid.clusters],
breaks = valid.clusters,
labels = lbl,
na.value = background.color
)
}
if (no.axes) plt <- plt & NoAxes()
return(plt)
}
#' Cell Feature Plot for Seurat Objects
#'
#' A flexible wrapper around Seurat's `FeaturePlot` to visualize gene expression or metadata
#' features in a Seurat object. Supports custom color palettes, viridis, and hotspot/rainbow palettes.
#'
#' @param object A `Seurat` object.
#' @param features Character vector of features (genes or metadata columns) to plot.
#' @param cols Optional character vector of colors for plotting.
#' @param theme.cols Character. Predefined theme color palette (default: `"Reds"`). Options include `"Reds"`, `"Blues"`, etc., or custom list palettes `"hotspot"` and `"rainbow"`.
#' @param viridis Logical. If TRUE, use viridis palette instead of RColorBrewer or custom palettes.
#' @param viridis.opt Character. Viridis palette option (default: `"D"`).
#' @param rev.cols Logical. Reverse the color palette (default: `FALSE`).
#' @param na.col Color for NA or below-cutoff expression values (default: `"lightgray"`).
#' @param order Logical. If TRUE, plot high-expression cells on top (default: `FALSE`).
#' @param pt.size Numeric. Point size. If NULL, automatically calculated based on number of cells.
#' @param font.size Numeric. Base font size for plot titles and axis labels (default: 10).
#' @param reduction Character. Dimensional reduction to use (default: first available in Seurat object).
#' @param na.cutoff Numeric. Minimum expression value for coloring; below this will be NA if palette requires (default: 1e-9).
#' @param raster Logical. If TRUE, rasterize points for faster plotting of large datasets.
#' @param raster.dpi Numeric vector of length 2. DPI for rasterization (default: c(512,512)).
#' @param split.by Character. Metadata column to split the plot.
#' @param ncol Numeric. Number of columns when combining multiple plots.
#' @param layer Character. Seurat assay slot to fetch data from (default: `"data"`).
#' @param label Logical. Whether to label clusters (default: FALSE).
#' @param axes Logical. Whether to show axes (default: TRUE).
#' @param combine Logical. Whether to return a single combined plot (default: TRUE).
#' @param blend Logical. Whether to blend exactly two features (default: FALSE).
#' @param merge.leg Logical. Whether to merge multiple legends into one (default: FALSE).
#' @param theme.style Character. ggplot2 theme to apply (default: `"classic"`).
#' @param ... Additional arguments passed to `Seurat::FeaturePlot`.
#'
#' @return A `ggplot` object (or `patchwork` object if multiple features).
#' @export
#'
#' @examples
#' # Single feature with default Reds palette
#' cellfeat(seurat_obj, features = "POMC")
#'
#' # Multiple features with viridis palette and merged legend
#' cellfeat(seurat_obj, features = c("POMC", "NPY"), viridis = TRUE, merge.leg = TRUE)
#'
#' # Blend two features
#' cellfeat(seurat_obj, features = c("POMC", "NPY"), blend = TRUE)
#'
cellfeat <- function(
object,
features,
cols = NULL,
theme.cols = "Reds",
viridis = FALSE,
viridis.opt = "D",
rev.cols = FALSE,
na.col = "lightgray",
order = FALSE,
pt.size = NULL,
font.size = 10,
reduction = NULL,
na.cutoff = 1e-9,
raster = NULL,
raster.dpi = c(512, 512),
split.by = NULL,
ncol = NULL,
layer = "data",
label = FALSE,
axes = TRUE,
combine = TRUE,
blend = FALSE,
merge.leg = FALSE,
theme.style = "classic",
...
) {
suppressWarnings(suppressMessages({
# Checks
if (!inherits(object, "Seurat"))
stop("object must be a Seurat object")
features <- intersect(
unique(as.character(features)),
c(rownames(object), colnames([email hidden]))
)
if (length(features) == 0)
stop("No valid features found")
reduction <- reduction %||% SeuratObject::DefaultDimReduc(object)
# Colors
# if (is.null(cols)) {
# cols <- RColorBrewer::brewer.pal(9, theme.cols)
# }
# if (rev.cols) cols <- rev(cols)
# Colors
if (!is.null(cols)) {
if (length(cols) < 2)
stop("colors_use must contain at least 2 colors")
cols <- cols
if (rev.cols) cols <- rev(cols)
} else if (isTRUE(viridis)) {
if (!requireNamespace("viridis", quietly = TRUE))
stop("Package 'viridis' is required when viridis = TRUE")
cols <- viridis::viridis(9, option = viridis.opt)
if (rev.cols) cols <- rev(cols)
} else {
pal <- list(
hotspot = c("navy", "lightblue", "yellow", "orange", "red"),
rainbow = c("blue", "green", "yellow", "orange", "red")
)
if (theme.cols %in% names(pal)) {
cols <- pal[[theme.cols]]
} else {
cols <- RColorBrewer::brewer.pal(9, theme.cols)
}
if (rev.cols) cols <- rev(cols)
}
# Point size
raster <- raster %||% (ncol(object) > 2e5)
if (is.null(pt.size)) {
pt.size <- if (raster) 1 else min(1583 / ncol(object), 1)
}
## Global max for merged legend
global_max <- NA
if (merge.leg && is.null(split.by) && length(features) > 1 && !blend) {
expr_data <- Seurat::FetchData(object, vars = features, layer = layer)
global_max <- max(expr_data, na.rm = TRUE)
}
# Plot
plt <- Seurat::FeaturePlot(
object = object,
features = features,
reduction = reduction,
order = order,
pt.size = pt.size,
raster = raster,
raster.dpi = raster.dpi,
split.by = split.by,
combine = combine,
blend = blend,
label = label,
ncol = ncol,
...
)
# Apply color scale (always)
# if (!blend) {
# if (merge.leg && is.null(split.by) && length(features) > 1) {
# # Shared global scale
# plt <- plt & scale_color_gradientn(
# colors = cols,
# limits = c(na.cutoff, global_max),
# na.value = na.col
# )
# } else {
# # Per-feature scale
# plt <- plt & scale_color_gradientn(
# colors = cols,
# limits = c(na.cutoff, NA),
# na.value = na.col,
# name = NULL
# )
# }
# }
# Determine palette behavior
like <- isTRUE(viridis) || theme.cols %in% c("hotspot", "rainbow")
# Brewer palettes should gray-out low expression
if (!like && !blend) {
plt <- plt & scale_color_gradientn(
colors = cols,
limits = c(na.cutoff, NA),
oob = scales::censor, # force < cutoff → NA
na.value = na.col, # show gray
name = NULL
)} else if (!blend) {
# No gray background
plt <- plt & scale_color_gradientn(
colors = cols,
limits = c(na.cutoff, NA),
oob = scales::squish, # clip instead of NA
na.value = NA
)
}
# Theme
plt <- plt & plot_theme(theme.style = theme.style, font.size = font.size, ...) &
theme(plot.title = element_text(
hjust = 0.5, size = font.size + 2, face = "bold.italic"
))
# ONLY add colorbar guide for non-blend plots
if (!blend) {
plt <- plt & guides(
color = guide_colorbar(
frame.colour = "black",
ticks.colour = "black"
)
)
}
# Manage axes
if (!axes) {
plt <- plt & Seurat::NoAxes()
}
## Merge legend
if (merge.leg && is.null(split.by) && length(features) > 1) {
plt <- Seurat::CombinePlots(plots = plt, legend = "right", ncol = ncol)
}
return(plt)
}))
}
#' Generate Expression Heatmaps for Stage or Celltype–Sample Modes
#'
#' This function generates heatmaps from a Seurat object using either:
#' 1) **Stage mode**: aggregated expression across developmental or experimental stages; or
#' 2) **Celltype–Sample mode**: paired annotation of cell types and samples.
#'
#' It supports z-score normalization, gene filtering, palette selection, optional
#' flipping of heatmaps, row/column label position control, and integrated saving
#' through a custom `save.fig()` function.
#'
#' @param object A Seurat object.
#' @param features Character vector of gene names to include.
#' @param mode Plot mode: `"stage"` (default) or `"celltype-sample"`.
#' @param assay Assay to use (default `"RNA"`).
#' @param slot Data slot to extract (default `"counts"`).
#' @param zlim Numeric vector defining z-score color scaling boundaries
#' (default `c(-2, 0, 2)`).
#' @param fontsize Base font size for labels.
#' @param palette Color palette. Either `"viridis"` or manual colors.
#' @param manual.colors Colors used when `palette != "viridis"`.
#' @param flip.heatmap Logical. If TRUE, transpose the heatmap and switch label sides.
#' @param group.by Metadata column used for averaging in `"stage"` mode.
#' @param column.title Optional column title for the heatmap.
#' @param celltype.col Metadata column defining cell types (for `"celltype-sample"` mode).
#' @param sample.col Metadata column defining samples (for `"celltype-sample"` mode).
#' @param celltype.palette Named color vector for cell types.
#' @param sample.palette Named color vector for samples.
#' @param show.annotation.legend Logical; show annotation legends (default TRUE).
#' @param row.names.side Side for row names: `"left"` or `"right"` (default `"left"`).
#' @param column.names.side Side for column names: `"top"` or `"bottom"` (default `"top"`).
#' @param row.names.rot Rotation angle of row labels.
#' @param column.names.rot Rotation angle of column labels.
#' @param output.pdf File path to save heatmap as PDF. If NULL, no file is saved.
#' @param pdf.width PDF width (inches).
#' @param pdf.height PDF height (inches).
#' @param heatmap.legend.side Side for heatmap legend.
#' @param annotation.legend.side Side for annotation legend.
#'
#' @details
#' The function uses Seurat's `AverageExpression()` to compute aggregated values,
#' applies per-gene z-scoring, and uses `ComplexHeatmap` for rendering.
#'
#' If `flip.heatmap = TRUE`, the matrix is transposed and label sides are swapped
#' automatically to maintain readability.
#'
#' If `output.pdf` is provided, the function saves the heatmap using the user’s
#' custom `save.fig()` helper, ensuring directory creation and figure consistency.
#'
#' @return Invisibly returns the ComplexHeatmap object used for plotting.
#'
#' @examples
#' \dontrun{
#' expression_heatmap(
#' object,
#' features = c("GATA3", "SOX2"),
#' mode = "stage",
#' group.by = "PCW.group",
#' output.pdf = "plots/stage.heatmap.pdf"
#' )
#' }
#'
#' @import Seurat
#' @import dplyr
#' @import ComplexHeatmap
#' @import viridis
#' @import circlize
#' @import stringr
#' @import grid
#'
#' @export
#'
expression_heatmap <- function(
object,
features,
mode = c("stage", "celltype-sample"),
assay = "RNA",
slot = "counts",
zlim = c(-2, 0, 2),
fontsize = 12,
palette = "viridis",
manual.colors = c("purple","black","yellow"),
flip.heatmap = FALSE,
group.by = "PCW_group",
column.title = NULL,
celltype.col = "ann_level_2",
sample.col = "sample",
celltype.palette = NULL,
sample.palette = NULL,
show.annotation.legend = TRUE,
row.names.side = "left",
column.names.side = "top",
row.names.rot = 0,
column.names.rot = 90,
output.pdf = NULL,
pdf.width = 8,
pdf.height = 6,
heatmap.legend.side = "right",
annotation.legend.side = "right"
){
require(Seurat)
require(ComplexHeatmap)
require(circlize)
require(viridis)
require(stringr)
require(grid)
mode <- match.arg(mode)
## helpers
italic_labels <- function(x) {
as.expression(lapply(x, function(g) bquote(italic(.(g)))))
}
make_labels <- function(mat, flip) {
if (!flip) {
list(row = italic_labels(rownames(mat)), col = NULL)
} else {
list(row = NULL, col = italic_labels(colnames(mat)))
}
}
make_col_fun <- function() {
if (palette == "viridis") {
colorRamp2(seq(zlim[1], zlim[3], length.out = 256), viridis(256))
} else {
colorRamp2(zlim, manual.colors)
}
}
## gene filtering
genes <- intersect(features, rownames(object))
if (length(genes) == 0) stop("None of the genes exist in Seurat object.")
sub <- subset(object, features = genes)
DefaultAssay(sub) <- assay
## compute average expression
if (mode == "stage") {
avg <- AverageExpression(
sub, group.by = group.by, assays = assay, slot = slot
)[[assay]]
} else {
avg <- AverageExpression(
sub, group.by = c(celltype.col, sample.col),
assays = assay, slot = slot
)[[assay]]
}
## z-score & ordering
avg <- avg[!rowSums(is.na(avg)), , drop = FALSE]
avg <- t(scale(t(avg)))
avg[is.na(avg)] <- 0
avg <- avg[order(apply(avg, 1, which.max)), , drop = FALSE]
## annotations (celltype × sample)
top.anno <- NULL
if (mode == "celltype-sample") {
parts <- str_split_fixed(colnames(avg), "_", 2)
anno.df <- data.frame(
CellTypes = factor(parts[,1], levels = levels(object[[celltype.col]][,1])),
Samples = factor(parts[,2], levels = levels(object[[sample.col]][,1]))
)
rownames(anno.df) <- colnames(avg)
if (is.null(celltype.palette)) {
celltype.palette <- setNames(
viridis(length(levels(anno.df$CellTypes))),
levels(anno.df$CellTypes)
)
}
if (is.null(sample.palette)) {
sample.palette <- setNames(
viridis(length(levels(anno.df$Samples))),
levels(anno.df$Samples)
)
}
top.anno <- HeatmapAnnotation(
df = anno.df,
col = list(CellTypes = celltype.palette, Samples = sample.palette),
show_legend = show.annotation.legend,
annotation_name_gp = gpar(fontsize = fontsize, fontface = "bold")
)
}
## flip heatmap
if (flip.heatmap) {
avg <- t(avg)
row.names.side <- ifelse(row.names.side == "left", "right", "left")
column.names.side <- ifelse(column.names.side == "top", "bottom", "top")
}
labels <- make_labels(avg, flip.heatmap)
## build heatmap
ht_args <- list(
matrix = avg,
name = "Z-score",
top_annotation = top.anno,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = make_col_fun(),
border = TRUE,
row_names_side = row.names.side,
column_names_side = column.names.side,
row_names_rot = row.names.rot,
column_names_rot = column.names.rot,
row_names_gp = gpar(fontsize = fontsize),
column_names_gp = gpar(fontsize = fontsize),
column_title = column.title,
heatmap_legend_param = list(
title = "Z-score",
title_gp = gpar(fontface = "bold"),
border = "black"
)
)
if (!flip.heatmap) {
ht_args$row_labels <- italic_labels(rownames(avg))
} else {
ht_args$column_labels <- italic_labels(colnames(avg))
}
ht <- do.call(Heatmap, ht_args)
## draw / save
if (!is.null(output.pdf)) {
out.dir <- dirname(output.pdf)
if (!dir.exists(out.dir)) dir.create(out.dir, recursive = TRUE)
ComplexHeatmap::draw(
ht,
heatmap_legend_side = heatmap.legend.side,
annotation_legend_side = annotation.legend.side,
merge_legend = TRUE
)
save_heatmap(
ht = ComplexHeatmap::draw(
ht,
heatmap_legend_side = heatmap.legend.side,
annotation_legend_side = annotation.legend.side,
merge_legend = TRUE
),
filename = output.pdf,
formats = c("pdf", "png", "tiff"),
width = pdf.width,
height = pdf.height,
dpi = 600
)
cat("Saved heatmap to:", output.pdf, "\n")
} else {
ComplexHeatmap::draw(
ht,
heatmap_legend_side = heatmap.legend.side,
annotation_legend_side = annotation.legend.side,
merge_legend = TRUE
)
}
}
#'#' Plot a Transcription Factor Regulatory Network
#'
#' Visualizes a transcription factor (TF) regulatory network using selected source TFs and their top N targets based on absolute regulatory strength (mor score).
#' Supports visual differentiation of positive and negative regulation, with shape and color coding for source vs. target nodes.
#'
#' @param net A data frame representing the TF regulatory network with columns: \code{source}, \code{target}, and \code{mor} (mode of regulation).
#' @param selected_sources A character vector of TFs to be used as source nodes.
#' @param n_targets Integer specifying the number of top targets to display per source (default: 5).
#' @param label_size Numeric; font size of node labels (default: 2).
#' @param node_size Numeric; size of node points (default: 15).
#' @param repel Logical; whether to repel node labels to avoid overlap using \code{geom_node_text} (default: FALSE).
#' @param layout_type Layout type for the network graph; passed to \code{ggraph::ggraph()} (default: "fr").
#' @param positive_color Color for positively regulated targets (default: "darkgreen").
#' @param negative_color Color for negatively regulated targets (default: "#e59866").
#' @param plot_title Title of the plot (default: "TF Regulatory Network").
#'
#' @return A ggplot2 object representing the regulatory network.
#'
#' @import ggraph
#' @import igraph
#' @import dplyr
#' @export
#'
#' @examples
#' \dontrun{
#' plot_tf_network(net = tf_network_df,
#' selected_sources = c("SOX2", "FOXA1"),
#' n_targets = 10,
#' repel = TRUE)
#' }
plot_tf_network <- function(net,
selected_sources,
n_targets = 5,
label_size = 2,
node_size = 15,
linewidth = 1.2,
repel = FALSE,
layout_type = "graphopt", # deterministic
max_edges = 500,
positive_color = "darkred",
negative_color = "navy",
plot_title = "TF Regulatory Network") {
require(ggraph)
require(igraph)
require(dplyr)
# Filter network
filtered_net <- net %>%
filter(source %in% selected_sources) %>%
group_by(source) %>%
top_n(n = n_targets, wt = abs(mor)) %>%
ungroup()
if (nrow(filtered_net) > max_edges) {
filtered_net <- filtered_net %>%
arrange(desc(abs(mor))) %>%
slice(1:max_edges)
warning(paste("Edge count exceeds", max_edges, "— truncating to top", max_edges, "edges."))
}
# Build graph
g <- graph_from_data_frame(filtered_net, directed = TRUE)
filtered_net$edge_color <- ifelse(filtered_net$mor > 0, "Positive", "Negative")
node_df <- data.frame(name = V(g)$name) %>%
mutate(node_type = ifelse(name %in% selected_sources, "Source", "Target"),
shape = ifelse(node_type == "Source", 22, 21),
edge_color = ifelse(node_type == "Source", "Source",
ifelse(name %in% filtered_net$target[filtered_net$mor > 0], "Positive", "Negative")),
label_color = ifelse(node_type == "Source", "black", "white"))
node_color_map <- c("Positive" = positive_color,
"Negative" = negative_color,
"Source" = "white")
# FIXED LAYOUT → stable, reproducible!
set.seed(42) # <-- makes layout identical every run
# Plot
p <- ggraph(g, layout = layout_type) +
geom_edge_link0(aes(edge_alpha = abs(mor), edge_colour = filtered_net$edge_color),
show.legend = TRUE, linewidth = linewidth) +
geom_node_point(aes(fill = factor(name, levels = node_df$name, labels = node_df$edge_color)),
size = node_size, shape = node_df$shape, show.legend = FALSE) +
geom_node_text(aes(label = name, color = factor(name, levels = node_df$name, labels = node_df$label_color)),
size = label_size, fontface = "bold.italic", repel = repel) +
# Legend
scale_edge_color_manual(
name = "Regulation",
values = c("Positive" = positive_color, "Negative" = negative_color),
guide = guide_legend(
override.aes = list(size = 6),
title.position = "top",
title.theme = element_text(size = 18, face = "bold"),
label.theme = element_text(size = 16)
)
) +
scale_fill_manual(values = node_color_map, guide = "none") +
scale_color_manual(values = c("black", "white"), guide = "none") +
guides(edge_alpha = "none") +
theme_void() +
ggtitle(plot_title) +
theme(
plot.margin = unit(c(1,1,1,1), "cm"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18, face = "bold"),
legend.key.size = unit(1.4, "cm"),
legend.key.height = unit(1.4, "cm"),
legend.key.width = unit(1.4, "cm"),
legend.position = "right"
) +
scale_x_discrete(expand = expansion(mult = c(0.1, 0.1))) +
scale_y_discrete(expand = expansion(mult = c(0.1, 0.1)))
return(p)
}
plot_tf_activity <- function(
seurat_object,
tf_list,
assay_activity = "tfsulm",
group.by = "seurat_clusters",
low_color = "navy",
high_color = "darkred",
zscore = TRUE,
zscore_mode = c("per_tf", "global"),
fontsize = 10,
legend.text.size = 10,
theme = "classic",
ncol = 1,
leg.pos = "bottom",
return_data = TRUE
) {
require(Seurat)
require(ggplot2)
require(dplyr)
require(patchwork)
zscore_mode <- match.arg(zscore_mode)
DefaultAssay(seurat_object) <- assay_activity
found_activity <- intersect(tf_list, rownames(GetAssayData(seurat_object, slot = "data")))
if (length(found_activity) == 0)
stop("None of the TFs were found in the assay.")
# Compute TF activity data (ACCUMULATE CORRECTLY)
all_tf_data <- list()
for (tf in found_activity) {
df <- FetchData(seurat_object, vars = c(tf, group.by)) %>%
group_by(.data[[group.by]]) %>%
summarise(
mean_val = mean(.data[[tf]], na.rm = TRUE),
n_cells = n(),
.groups = "drop"
) %>%
mutate(TF = tf)
df$zval <- if (zscore) scale(df$mean_val)[, 1] else df$mean_val
all_tf_data[[tf]] <- df
}
tf_activity_df <- bind_rows(all_tf_data)
# Global z-score limits (if requested)
if (zscore && zscore_mode == "global") {
global_abs_lim <- max(abs(tf_activity_df$zval), na.rm = TRUE)
} else {
global_abs_lim <- NULL
}
# Plotting
plot_list <- list()
total_tf <- length(tf_list)
for (i in seq_along(tf_list)) {
tf <- tf_list[i]
df <- tf_activity_df %>% filter(TF == tf)
if (nrow(df) == 0) {
p <- ggplot() + theme_void() + ggtitle(paste0(tf, " (not found)"))
} else {
show_x <- i == total_tf
p <- ggplot(df, aes(x = .data[[group.by]], y = zval, fill = zval)) +
geom_col(color = "black", linewidth = 0.2, width = 0.55) +
scale_fill_gradient2(
low = low_color,
mid = "white",
high = high_color,
midpoint = 0,
limits = if (!is.null(global_abs_lim)) c(-global_abs_lim, global_abs_lim) else NULL,
name = if (zscore) "z-score" else "Mean activity"
) +
plot_theme(font.size = fontsize, theme = theme, axis.title.x.show = FALSE) +
theme(
axis.text.x = if (show_x) element_text(angle = 45, hjust = 1) else element_blank(),
axis.ticks.x = if (show_x) element_line() else element_blank(),
axis.title.x = if (show_x) element_text() else element_blank(),
plot.title = element_text(face = "bold.italic", hjust = 0.5)
) +
ggtitle(tf) +
ylab(if (zscore) "z-score" else "Mean TF activity") +
xlab(if (show_x) group.by else "")
}
plot_list[[tf]] <- p
}
final_plot <- wrap_plots(
plot_list,
ncol = ncol,
guides = "collect",
axis_titles = "collect"
) &
theme(
legend.position = leg.pos,
legend.title = element_text(size = legend.text.size),
legend.text = element_text(size = legend.text.size)
)
# Return
if (return_data) {
return(list(
plot = final_plot,
data = tf_activity_df
))
} else {
return(final_plot)
}
}
#' Plot dominance score results
#'
#' @param df Data frame returned by `dominance_score()`$df
#' @param mode Plot mode: "bins", "high_cells", "top_genes"
#' @param ident.col Grouping column for "bins" and "high_cells" (default "ident")
#' @param stage.col Optional column to facet/stratify plots by, e.g., "stage" or "sample"
#' @param bin.col Column name for dominance bins (default "adjusted_bin")
#' @param gene.col Column name for genes (default "top_gene")
#' @param min.adjust Minimum adjusted dominance score (used in high_cells & top_genes)
#' @param min.top.expr Minimum top gene expression
#' @param max.n.genes Maximum number of expressed genes
#' @param top.n Number of top genes to show (mode = "top_genes")
#' @param colors Optional named vector of colors
#' @param font.size Base font size
#' @param theme ggplot theme name
#' @param leg.pos Legend position
#' @param x.angle X-axis text angle
#'
#' @return ggplot object
#' @export
plot_domscore <- function(
df,
mode = c("bins", "high_cells", "top_genes"),
ident.col = "ident",
stage.col = NULL,
bin.col = "adjusted_bin",
gene.col = "top_gene",
min.adjust = 1,
min.top.expr = 1,
max.n.genes = 3,
top.n = 20,
colors = NULL,
font.size = 10,
theme = "classic",
plot.ttl = NULL,
leg.pos = "top",
leg.dir = "horizontal",
x.angle = 45,
ncol = 2,
flip = TRUE,
facet.bg = FALSE,
...
) {
require(dplyr)
require(ggplot2)
mode <- match.arg(mode)
if (is.null(colors)) {
groups <- unique(df[[ident.col]])
colors <- setNames(
grDevices::colorRampPalette(RColorBrewer::brewer.pal(8, "Set2"))(length(groups)),
groups
)
}
# MODE: dominance bins
if (mode == "bins") {
plot.df <- df %>%
dplyr::filter(!is.na(.data[[bin.col]])) %>%
{
if (!is.null(stage.col) && stage.col %in% colnames(.)) {
group_by(., .data[[stage.col]], .data[[bin.col]], .data[[ident.col]])
} else {
group_by(., .data[[bin.col]], .data[[ident.col]])
}
} %>%
summarise(n = n(), .groups = "drop") %>%
{
if (!is.null(stage.col) && stage.col %in% colnames(.)) {
group_by(., .data[[stage.col]], .data[[bin.col]])
} else group_by(., .data[[bin.col]])
} %>%
mutate(perc = 100 * n / sum(n)) %>%
ungroup()
p <- ggplot(plot.df,
aes(x = .data[[bin.col]], y = perc, fill = .data[[ident.col]])) +
geom_col(color = "black", linewidth = 0.2, width=0.5) +
scale_fill_manual(values = colors) +
labs(x = "Dominance score", y = "Cell type proportion (%)", fill = "", title = plot.ttl) +
plot_theme(theme = theme, font.size = font.size, leg.pos = leg.pos,leg.dir = leg.dir, x.angle = 0)
}
# MODE: high-dominance cells
else if (mode == "high_cells") {
plot.df <- df %>%
filter(
adjust_score > min.adjust,
top_expr >= min.top.expr,
n_genes_expr <= max.n.genes
) %>%
{
if (!is.null(stage.col) && stage.col %in% colnames(.)) {
group_by(., .data[[stage.col]], .data[[ident.col]])
} else group_by(., .data[[ident.col]])
} %>%
summarise(n_cells = n(), .groups = "drop")
x_col <- if (!is.null(stage.col) && stage.col %in% colnames(df)) stage.col else ident.col
p <- ggplot(plot.df,
aes(x = factor(.data[[x_col]]),
y = n_cells,
fill = .data[[ident.col]])) +
geom_col(color = "black", linewidth = 0.2, width = 0.5) +
scale_fill_manual(values = colors) +
labs(x = "", y = paste0("High-dominance cells\n(adjusted > ", min.adjust, ")"), fill = "", title = plot.ttl) +
plot_theme(theme = theme, font.size = font.size, leg.pos = leg.pos, leg.dir = leg.dir, x.angle = x.angle)
}
# MODE: top dominant genes
if (mode == "top_genes") {
high.df <- df %>%
filter(
adjust_score > min.adjust,
top_expr >= min.top.expr,
n_genes_expr <= max.n.genes,
!is.na(.data[[gene.col]])
)
high.df[[gene.col]] <- as.character(high.df[[gene.col]])
if (!is.null(stage.col) && stage.col %in% colnames(high.df)) {
plot.df <- high.df %>%
group_by(.data[[stage.col]], .data[[gene.col]]) %>%
summarise(n = n(), .groups = "drop") %>%
arrange(.data[[stage.col]], desc(n)) %>%
group_by(.data[[stage.col]]) %>%
slice_head(n = top.n) %>%
ungroup()
plot.df$n <- as.integer(plot.df$n)
p <- ggplot(plot.df, aes(x = reorder(.data[[gene.col]], -n), y = n)) + # <- descending
geom_col(fill = "steelblue", color = "black", linewidth = 0.2, width = 0.5) +
facet_wrap(as.formula(paste0("~", stage.col)), ncol = ncol) +
scale_y_continuous(breaks = scales::pretty_breaks()) +
labs(x = "Top dominant OR genes", y = "Number of high-dominance cells", title = plot.ttl) +
plot_theme(font.size = font.size, theme.style = theme, leg.pos = "none", x.angle = x.angle, facet.bg = facet.bg,...)
if (flip) p <- p + coord_flip()
} else {
plot.df <- high.df %>%
dplyr::count(!!rlang::sym(gene.col), sort = TRUE, name = "n") %>%
slice_head(n = top.n)
plot.df$n <- as.integer(plot.df$n)
p <- ggplot(plot.df, aes(x = reorder(!!rlang::sym(gene.col), -n), y = n)) + # <- descending
geom_col(fill = "steelblue", color = "black", linewidth = 0.2, width = 0.5) +
scale_y_continuous(breaks = scales::pretty_breaks()) +
labs(x = "", y = "Number of high-dominance cells", title = plot.ttl) +
plot_theme(font.size = font.size, theme.style = theme, leg.pos = "none", x.angle = x.angle,facet.bg = facet.bg,...)
if (flip) p <- p + coord_flip()
}
}
if (!is.null(colors)) p <- p + scale_fill_manual(values = colors)
return(p)
}
#'------------------------------------------------------------------------------
#'------------------------------------------------------------------------------
#'
#' Plot Heatmap Along Pseudotime
#'
#' Generates a heatmap of gene expression along pseudotime for selected lineages, with
#' flexible annotation and legend options. Supports automatic legend type detection
#' (discrete vs continuous), smooth scaling, z-score transformation, and flexible
#' legend title positioning.
#'
#' @param heatmap_matrices List of heatmap matrices produced by `identify_temporal_genes` or similar. Each element should have `expr_raw`, `ordered_cells`, `top_genes`.
#' @param seurat_obj Seurat object containing cell metadata.
#' @param lineage_name Character. Name of the lineage to plot. Defaults to first element of `heatmap_matrices`.
#' @param pseudotime_vec Numeric vector of pseudotime values. Defaults to sequential ordering of `ordered_cells`.
#' @param gene_list Character vector of genes to plot. Defaults to `hm$top_genes`.
#' @param annotation_col Column in `[email hidden]` to use for annotation. Default `"ann_level_2"`.
#' @param celltype_colors Named vector of colors for annotation clusters. If `NULL`, colors are auto-generated.
#' @param zscore_limits Numeric vector of length 3, c(min, mid, max), used for color scaling. Default c(-2,0,2).
#' @param pseudotime_palette Character. Palette name for pseudotime continuous scale. Default `"YlGnBu"`.
#' @param row_cluster_colors Not currently used. Placeholder for future row annotations.
#' @param font_size Font size for heatmap labels and legends. Default 10.
#' @param output_pdf Path to save heatmap. If `NULL`, heatmap is drawn on screen.
#' @param fig_name Filename for saving heatmap. Default uses `plot_title`.
#' @param pdf_width Width of saved PDF/PNG. Default 6.
#' @param pdf_height Height of saved PDF/PNG. Default 10.
#' @param plot_title Title of heatmap. Default `NULL`.
#' @param show_annotation_legend Logical. Show legend for top annotations. Default `TRUE`.
#' @param smooth_df Degrees of freedom for smoothing spline. Default 3.
#' @param legend_dir `"vertical"` or `"horizontal"`. Orientation of legend. Default `"vertical"`.
#' @param ht_legend_side Side of heatmap legend (`"right"` or `"left"`). Default `"right"`.
#' @param ann_legend_side Side of top annotation legend (`"top"` or `"right"`). Default `"right"`.
#' @param legend_width_cm Width of legends in cm. Default 3.
#' @param cont_legend_cm Height of continuous legends in cm. Default 3.
#' @param merge_legends Logical. Merge annotation and heatmap legends. Default `FALSE`.
#' @param center_legend_titles Logical. Center all heatmap legend titles. Default `TRUE`.
#' @param center_cells_title Logical. Center "Cells" annotation legend title. Default `TRUE`.
#' @param center_pt_title Logical. Center "Pseudotime" annotation legend title. Default `TRUE`.
#' @param cells_title_pos Custom position for "Cells" title (overrides centering). Default `NULL`.
#' @param pt_title_pos Custom position for "Pseudotime" title (overrides centering). Default `NULL`.
#' @param hm_title_pos Custom position for heatmap legend title. Default `NULL`.
#'
#' @return `ComplexHeatmap` object invisibly. If `output_pdf` is provided, the heatmap is saved as PNG.
#' @export
#'
#' @examples
#' plot_heatmap_pseudotime(
#' heatmap_matrices = l1$heatmap_matrices,
#' seurat_obj = sub,
#' lineage_name = "Lineage1",
#' pseudotime_vec = pt_vec,
#' gene_list = gene_list,
#' celltype_colors = cols,
#' font_size = 8,
#' pdf_width = 4,
#' pdf_height = 6,
#' output_pdf = file.path(dir.results),
#' plot_title = "Neuronal lineage",
#' legend_dir = "vertical",
#' ht_legend_side = "right",
#' ann_legend_side = "right",
#' annotation_col = "ann2",
#' merge_legends = TRUE,
#' center_legend_titles = TRUE,
#' center_cells_title = TRUE,
#' center_pt_title = FALSE,
#' cells_title_pos = "topcenter",
#' pt_title_pos = "leftcenter",
#' hm_title_pos = "topcenter"
#' )
plot_heatmap_pseudotime <- function(
heatmap_matrices,
seurat_obj,
lineage_name = NULL,
pseudotime_vec = NULL,
gene_list = NULL,
annotation_col = "ann_level_2",
celltype_colors = NULL,
zscore_limits = c(-2, 0, 2),
pseudotime_palette = "YlGnBu",
row_cluster_colors = NULL,
font_size = 10,
output_pdf = NULL,
fig_name = paste(plot_title),
pdf_width = 6,
pdf_height = 10,
plot_title = NULL,
show_annotation_legend = TRUE,
smooth_df = 3,
legend_dir = "vertical",
ht_legend_side = "right",
ann_legend_side = "right",
legend_width_cm = 3,
cont_legend_cm = 3,
merge_legends = FALSE,
center_legend_titles = TRUE,
center_cells_title = TRUE,
center_pt_title = TRUE,
cells_title_pos = NULL,
pt_title_pos = NULL,
hm_title_pos = NULL,
file_path = file.path(outdir, fig_name)
){
require(ComplexHeatmap)
require(circlize)
require(RColorBrewer)
# --- Setup lineage and genes
if(is.null(lineage_name)) lineage_name <- names(heatmap_matrices)[1]
if(!lineage_name %in% names(heatmap_matrices)) stop("Lineage not found")
hm <- heatmap_matrices[[lineage_name]]
if(is.null(gene_list)) gene_list <- hm$top_genes
cells_hm <- hm$ordered_cells
expr_mat <- as.matrix(hm$expr_raw[gene_list, cells_hm, drop = FALSE])
# --- Smooth
if(!is.null(smooth_df) && smooth_df > 0){
expr_mat <- t(apply(expr_mat, 1, function(x){
if(sd(x) == 0) return(x)
tryCatch(smooth.spline(x, df = smooth_df)$y, error = function(e) x)
}))
}
# --- Z-score scaling
row_means <- rowMeans(expr_mat, na.rm = TRUE)
row_sds <- apply(expr_mat, 1, sd, na.rm = TRUE)
row_sds[row_sds == 0] <- 1
expr_mat_scaled <- (expr_mat - row_means) / row_sds
# --- Order genes by pseudotime peak
gene_peak <- apply(expr_mat_scaled, 1, which.max)
expr_mat_scaled <- expr_mat_scaled[order(gene_peak), ]
# --- Clusters and annotation
clusters <- [email hidden][cells_hm, annotation_col]
if(!is.factor(clusters)){
cluster_order <- unique(clusters)
clusters <- factor(clusters, levels = cluster_order)
} else {
cluster_order <- levels(clusters)
}
clusters[is.na(clusters)] <- "Unknown"
if(is.null(celltype_colors)){
celltype_colors <- setNames(
colorRampPalette(brewer.pal(8,"Dark2"))(length(cluster_order)),
cluster_order
)
}
if(is.null(pseudotime_vec)) pseudotime_vec <- seq_along(cells_hm)
pseudotime <- pseudotime_vec[cells_hm]
pt_col_fun <- colorRamp2(
seq(min(pseudotime), max(pseudotime), length.out = 100),
viridis::viridis(100, option = "C")
)
# --- Helper to resolve title positions
resolve_title_pos <- function(title_pos, center, is_discrete, legend_dir){
if(is_discrete){
if(!is.null(title_pos)) return(title_pos)
if(legend_dir == "vertical") return("left") else return("top")
} else {
if(!is.null(title_pos)) return(title_pos)
if(center){
if(legend_dir == "vertical") return("leftcenter-rot") else return("topcenter")
} else {
if(legend_dir == "vertical") return("lefttop-rot") else return("topcenter")
}
}
}
is_cells_discrete <- is.factor(clusters)
is_pt_discrete <- !is.numeric(pseudotime)
cells_pos <- resolve_title_pos(cells_title_pos, center_cells_title, is_cells_discrete, legend_dir)
pt_pos <- resolve_title_pos(pt_title_pos, center_pt_title, is_pt_discrete, legend_dir)
top_annot <- HeatmapAnnotation(
"Cells" = factor(clusters, levels = cluster_order),
"Pseudotime" = pseudotime,
col = list("Cells" = celltype_colors, "Pseudotime" = pt_col_fun),
show_legend = show_annotation_legend,
annotation_legend_param = list(
"Cells" = list(
border = "black",
labels_gp = gpar(fontsize = font_size),
title_gp = gpar(
fontsize = font_size,
fontface = "bold",
just = if(center_cells_title) "center" else "left"
),
title_position = cells_pos,
legend_direction = legend_dir,
legend_width = unit(legend_width_cm, "cm"),
legend_height = unit(cont_legend_cm, "cm"),
nrow = if(legend_dir == "horizontal") 1 else NULL
),
"Pseudotime" = list(
border = "black",
labels_gp = gpar(fontsize = font_size),
title_gp = gpar(
fontsize = font_size,
fontface = "bold",
just = if(center_pt_title) "center" else "left"
),
title_position = pt_pos,
legend_direction = legend_dir,
legend_width = unit(legend_width_cm, "cm"),
legend_height = unit(cont_legend_cm, "cm")
)
),
annotation_name_gp = gpar(fontsize = font_size, fontface = "bold")
)
ht <- Heatmap(
expr_mat_scaled,
name = "Z-score",
col = colorRamp2(seq(-2, 2, length = 11), rev(brewer.pal(11, "Spectral"))),
cluster_rows = FALSE,
cluster_columns = FALSE,
show_column_names = FALSE,
show_row_names = TRUE,
row_names_gp = gpar(fontsize = font_size, fontface="italic"),
top_annotation = top_annot,
row_title_rot = 0,
cluster_row_slices = FALSE,
show_row_dend = FALSE,
column_title = plot_title,
column_title_gp = gpar(fontsize = font_size+2, fontface="bold"),
heatmap_legend_param = list(
title = "Z-score",
legend_direction = legend_dir,
title_position = resolve_title_pos(hm_title_pos, center_legend_titles, FALSE, legend_dir),
legend_width = unit(legend_width_cm, "cm"),
legend_height = unit(cont_legend_cm, "cm"),
labels_gp = gpar(fontsize = font_size),
title_gp = gpar(
fontsize = font_size,
fontface = "bold",
just = if(center_legend_titles) "center" else "left"
),
border = "black"
)
)
if(!is.null(output_pdf)){
ComplexHeatmap::draw(ht,
heatmap_legen

000-functions.R at commit 2bd43a8, under GPL-3.0 · at the source

Overview

Authors: Yvon Mbouamboua^1,2, Kevin Lebrigand², Sreekala Nampoothiri¹, Marie Couralet², Marie-Jeanne Arguel², Ludovica Cotellessa¹, Cécile Allet¹, Vincent Prevot¹, Pascal Barbry², Paolo Giacobini¹

ORCID iDs: Kevin Lebrigand, Sreekala Nampoothiri, Marie Couralet, Marie-Jeanne Arguel, Ludovica Cotellessa, Cécile Allet, Vincent Prevot, Pascal Barbry, Paolo Giacobini

Institutions: Centre Hospitalier Universitaire de Lille (France); Inserm (France)

Journal: Nature communications, volume 17, issue 1, article 3537

Dates: received 15 September 2025; accepted 16 March 2026; published online 17 April 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41467-026-71595-6 · PMID 41997969 · PMCID PMC13090377 · OpenAlex W7154687751

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), histology / microscopy (modality), human (organism)

Methods: Statistics, Smoothing, state filtering, decompositions

Keywords: Morphogenesis, Neurogenesis

MeSH: Olfactory Mucosa*, Cell Differentiation, Female, Fetus, Gene Expression Regulation, Developmental, Humans, In Situ Hybridization, Fluorescence, Male, Neurodevelopment, Neurogenesis, Olfactory Receptor Neurons, Receptors, Odorant, Single-Cell Analysis, Single-Cell Gene Expression Analysis, Spatial Transcriptomics, Transcriptome (* major topic)

Topic: Olfactory and Sensory Function Studies (Sensory Systems, Neuroscience), according to OpenAlex

Funding: Agence Nationale de la Recherche (French National Research Agency) (ANR-24-CHBS-0002, ANR-19-CE14–0027, ANR-23-IAHU-0007)

Citations: cited by 1 paper (Europe PMC); 67 references in the paper

Abstract

The human nasal region arises from neural crest and placodal lineages, yet its early development remains poorly understood owing to limited fetal tissue access and structural complexity. Here we present an integrated single-nucleus and spatial transcriptomic atlas of the human fetal nasal region, generated from male and female fetuses between 7 and 12 post-conceptional weeks. Single-nucleus RNA sequencing [snRNA-seq] resolved 32 distinct cell types, while integration with multiplexed error-robust fluorescence in situ hybridization (MERFISH) enabled spatial and temporal mapping of gene expression dynamics across the olfactory epithelium (OE) and adjacent tissues. We identify markers of olfactory sensory neuron differentiation and pathways governing epithelial patterning and OE morphogenesis. Notably, spatially resolved snRNA-seq profiles of 169 olfactory receptor genes reveal molecular support for the “one neuron-one receptor” principle already in the first trimester. Together, this work establishes a molecular and spatial framework of early human olfactory development and provides a resource for studies of sensory neurogenesis and congenital disorders.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 26 matches between paragraphs and lines of code.

ymbouamboua/HuDeCa

License: GPL-3.0

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: 2bd43a819a472b3982789e8d5d4e0bd672530a1e, 26 February 2026

Languages: R (18), Shell (3), Jupyter (2)

Size: 28 files, 23 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, license file, 18 notebooks

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: tidyverse (9 files), Seurat (6 files), patchwork (3 files), circlize (2 files), ComplexHeatmap (2 files), ggplot2 (2 files), Matplotlib (2 files), mgcv (2 files), NumPy (2 files), pandas (2 files), reshape2 (2 files), BCFtools (1 file), caret (1 file), cowplot (1 file), ggpubr (1 file), Harmony (1 file), igraph (1 file), lme4 (1 file), Plotly (1 file), seaborn (1 file), SingleCellExperiment (1 file)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

25 files

Python/fig6.f-g.ipynb, Jupyter, 227 lines
Python/figS1.ipynb, Jupyter, 1,268 lines, 3 matches
R/000-custom-colors.R, R, 176 lines, 1 match
R/000-functions.R, R, 6,024 lines, 11 matches
R/000-qc.Rmd, R, 168 lines
R/001-integration.Rmd, R, 113 lines
R/002-annotation.Rmd, R, 638 lines, 2 matches
R/fig1.Rmd, R, 195 lines
R/fig2.Rmd, R, 113 lines
R/fig3.Rmd, R, 638 lines, 1 match
R/fig4.Rmd, R, 371 lines, 3 matches
R/fig6.Rmd, R, 339 lines
R/figS1.Rmd, R, 114 lines
R/figS10.Rmd, R, 203 lines
R/figS2.Rmd, R, 337 lines
R/figS3.Rmd, R, 212 lines
R/figS4.Rmd, R, 197 lines, 2 matches
R/figS5.Rmd, R, 117 lines
R/figS6.Rmd, R, 143 lines
R/figS9.Rmd, R, 106 lines
bash/000-make_bulk_vcf.sh, Shell, 166 lines
bash/000-run_demuxafy.sh, Shell, 68 lines
bash/001-run-combined.sh, Shell, 27 lines
LICENSE, License, 674 lines
README.md, Text, 6 lines

At the source: github.com/ymbouamboua/HuDeCa

cobioda/human_fetal_olfactory_system

License: none: the authors keep all their rights

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: eb07faa5e0c6bd200b885eb8ae4f5646d89d9035, 4 December 2025

Languages: Jupyter (16)

Size: 18 files, 16 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, 16 notebooks

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: pandas (11 files), Scanpy (11 files), anndata (10 files), Matplotlib (9 files), NumPy (8 files), seaborn (8 files), CuPy (1 file), statannotations (1 file)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

12 files

00.load.save.ipynb, Jupyter, 44 lines
01.focus.ipynb, Jupyter, 112 lines
02.harmony.ipynb, Jupyter, 99 lines
04.novae.ipynb, Jupyter, 48 lines
05.scmusk.ipynb, Jupyter, 273 lines, 2 matches
06.xenium.exporter.ipynb, Jupyter, 58 lines
07.exploration.ipynb, Jupyter, 155 lines
07.proportion.ipynb, Jupyter, 87 lines, 1 match
07.pseudobulk.ipynb, Jupyter, 75 lines
08.FIGURES.5bcdei.ipynb, Jupyter, 113 lines
08.FIGURES.5jk.ipynb, Jupyter, 162 lines
08.FIGURES.7ab.ipynb, Jupyter, 148 lines
repository limit reached (2,000 files or 30 MB): the rest is at the source (5 files)

At the source: github.com/cobioda/human_fetal_olfactory_system

Code availability

Custom scripts and codes used for snRNA-seq analysis can be found at https://github.com/ymbouamboua/HuDeCa. Python scripts for re-analysis and figures production of MERFISH experiments can be found at https://github.com/cobioda/human_fetal_olfactory_system.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
35 scripts, each with its path and the digest of its content;
26 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

geo:GSE303809, at NCBI GEO; found in “Data availability”
zenodo:18245692, at Zenodo; found in “Data availability”

Data availability

All data generated or analyzed in the current study are included in the article and its Supplementary Figs. and Supplementary Data. The raw snRNA-seq datasets generated during the current study have been deposited in the European Genome/Phenome Archive (https://ega-archive.org/datasets/EGAD50000001712) and are available through controlled access due to privacy and consent restrictions related to human genomic data. Access to EGA archive datasets can be obtained by formal application to the Data Access Committee (DAC). Each DAC requires users/applicants to sign a Data Access Agreement (DAA), which details the terms and conditions of use for each dataset. An interactive Shiny application associated with processed human snRNA-seq data has been deposited on Zenodo (10.5281/zenodo.18245692). MERFISH data have been deposited in Gene Expression Omnibus under accession code GSE303809 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE303809). Source data are provided with this paper.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 2 keywords, 16 MeSH terms, 1 funder, 67 references.

Cite

This paper

Mbouamboua, Y., Lebrigand, K., Nampoothiri, S., Couralet, M., Arguel, M.-J., Cotellessa, L., Allet, C., Prevot, V., Barbry, P., & Giacobini, P. (2026). A single-cell and spatial atlas of early human olfactory development. Nature communications, 17(1), 3537. https://doi.org/10.1038/s41467-026-71595-6

BibTeX

@article{mbouamboua2026single, author = {Mbouamboua, Yvon and Lebrigand, Kevin and Nampoothiri, Sreekala and Couralet, Marie and Arguel, Marie-Jeanne and Cotellessa, Ludovica and Allet, Cécile and Prevot, Vincent and Barbry, Pascal and Giacobini, Paolo}, title = {{A single-cell and spatial atlas of early human olfactory development}}, journal = {Nature communications}, year = {2026}, month = apr, volume = {17}, number = {1}, pages = {3537}, publisher = {Nature Publishing Group}, issn = {2041-1723}, doi = {10.1038/s41467-026-71595-6}, url = {https://doi.org/10.1038/s41467-026-71595-6}, pmid = {41997969}, pmcid = {PMC13090377} }

RIS

TY - JOUR AU - Mbouamboua, Yvon AU - Lebrigand, Kevin AU - Nampoothiri, Sreekala AU - Couralet, Marie AU - Arguel, Marie-Jeanne AU - Cotellessa, Ludovica AU - Allet, Cécile AU - Prevot, Vincent AU - Barbry, Pascal AU - Giacobini, Paolo TI - A single-cell and spatial atlas of early human olfactory development T2 - Nature communications J2 - Nat Commun PY - 2026 DA - 2026/04/17 VL - 17 IS - 1 SP - 3537 SN - 2041-1723 PB - Nature Publishing Group DO - 10.1038/s41467-026-71595-6 UR - https://doi.org/10.1038/s41467-026-71595-6 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41467-026-71595-6", "type": "article-journal", "title": "A single-cell and spatial atlas of early human olfactory development", "container-title": "Nature communications", "author": [ { "family": "Mbouamboua", "given": "Yvon" }, { "family": "Lebrigand", "given": "Kevin" }, { "family": "Nampoothiri", "given": "Sreekala" }, { "family": "Couralet", "given": "Marie" }, { "family": "Arguel", "given": "Marie-Jeanne" }, { "family": "Cotellessa", "given": "Ludovica" }, { "family": "Allet", "given": "Cécile" }, { "family": "Prevot", "given": "Vincent" }, { "family": "Barbry", "given": "Pascal" }, { "family": "Giacobini", "given": "Paolo" } ], "container-title-short": "Nat Commun", "volume": "17", "issue": "1", "page": "3537", "DOI": "10.1038/s41467-026-71595-6", "PMID": "41997969", "PMCID": "PMC13090377", "ISSN": "2041-1723", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41467-026-71595-6", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 17 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Harmony, SingleCellExperiment, anndata, 17 other tools, genetics / omics, 3 references
[2] doi:10.1038/s41467-026-75722-1 [code]: Single-nucleus analysis of the adult human olfactory epithelium uncovers shared neurogenesis programs with the brain.
Journal: Nature communications
In common: SingleCellExperiment, anndata, circlize, 12 other tools, genetics / omics, 6 references
[3] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: Harmony, SingleCellExperiment, anndata, 16 other tools, genetics / omics, 1 reference
[4] doi:10.1038/s41467-026-74598-5 [code]: A synaptoid connectome differentiates tanycytic subpopulations and underlies neuroglial communication and neuroendocrine regulation.
Journal: Nature communications
In common: Seurat, cowplot, ggplot2, 1 other tool, 4 authors
[5] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: BCFtools, CuPy, caret, 15 other tools
[6] doi:10.3390/ijms27104466 [code]: Uncovering the Key Circuit FOSL2/FOS/EGR3/EGR1, Contributing to the Hyperexcitability of Excitatory Neurons in the Epileptic Temporal Cortex and Hippocampus.
Journal: International journal of molecular sciences
In common: caret, anndata, igraph, 13 other tools, genetics / omics, 2 references
[7] doi:10.1016/j.xcrm.2026.102651 [code]: Integrative CSF profiling identifies disease-specific immune responses in leptomeningeal disease.
Journal: Cell reports. Medicine
In common: Harmony, SingleCellExperiment, anndata, 15 other tools, genetics / omics
[8] doi:10.1038/s41586-026-10214-2 [code]: Multidimensional profiling of heterogeneity in supratentorial ependymomas.
Journal: Nature
In common: Harmony, SingleCellExperiment, anndata, 15 other tools, genetics / omics
[9] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: BCFtools, mgcv, igraph, 14 other tools, genetics / omics, 1 reference
[10] doi:10.1038/s41586-026-10629-x [code]: Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: Harmony, anndata, igraph, 14 other tools, genetics / omics, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.

A single-cell and spatial atlas of early human olfactory development.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

ymbouamboua/HuDeCa

cobioda/human_fetal_olfactory_system

Code availability

Tracing map

Data

Datasets cited

Data availability

Versions

Version 1, 29 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity