SCutils

An R package providing single-cell RNA-seq utility functions to complement Seurat-based workflows.

# install.packages("remotes")
remotes::install_github("jmgs7/SCutils")

Batch-opens .h5 single-cell matrices, writes/loads BPCells matrix directories (*_BP), supports both 10X and AnnData HDF5 inputs, and can optionally convert ENSEMBL IDs to gene symbols; it can also generate per-cell sample provenance metadata. Multiprocessing is supported via future framework. Future options must be set in the global environment before calling BatchOpenH5().

BatchOpenH5(
  files,
  relative = TRUE,
  BP.data.dir = NULL,
  platform = "10X",
  use.names = TRUE,
  ensembl.to.symbol = FALSE,
  species = "human",
  generate.metadata = FALSE
)

Calculates the Cellular Detection Rate (fraction of detected features per cell), z-scales it, and adds a CDR column to SeuratObject@meta.data. Handles both Seurat v3 (@assays$RNA@counts) and v5/split-layer (@assays$RNA@layers$counts) object formats.

SeuratObject <- CalculateCDR(SeuratObject)

Computes common single-cell QC metrics and appends them directly to SeuratObject@meta.data.

Always-computed columns (from raw counts):

• percent.mt: The percentage of mitochondrial gene counts per cell, computed from the ^MT- gene pattern (human) or ^mt- (mouse).
• percent.ribo: The percentage of ribosomal gene counts per cell, computed from the ^RP[SL] gene pattern (human) or ^Rp[sl] (mouse).
• percent.hb: The percentage of hemoglobin gene counts per cell, computed from the ^HB[^(P)]" gene pattern (human) or ^Hb^(p)]" (mouse).
• percent.ig: The percentage of immunoglobulin gene counts per cell, computed from the ^IG gene pattern (human) or ^Ig (mouse).
• percent.plat: The percentage of platelet gene counts per cell, computed from the ^PPBP|^PF4 gene pattern (human) or ^Ppbp|^Pf4 (mouse).
• percent.MALAT1: The percentage of MALAT1 gene counts per cell, computed from the MALAT1 gene pattern (human) or Malat1 (mouse).
• percent.S100A9: The percentage of S100A9 gene counts per cell, computed from the S100A9 gene pattern (human) or S100a9 (mouse).
• percent.S100A8: The percentage of S100A8 gene counts per cell, computed from the S100A8 gene pattern (human) or S100a8 (mouse).
• percent.FCGR3B: The percentage of FCGR3B gene counts per cell, computed from the FCGR3B gene pattern (human) or Fcgr3b (mouse).
• log10_nFeature_RNA: The log10 of the number of features per cell, computed from the nFeature_RNA column.
• log10_nCount_RNA: The log10 of the number of counts per cell, computed from the nCount_RNA column.
• complexity: The complexity of the library per cell, computed as log10(nFeature_RNA) / log10(nCount_RNA).

Conditionally computed (requires normalized counts in data layer):

• nCount_logRNA, nFeature_logRNA: log-normalized count and detected feature totals per cell.
• Cell cycle scoring (S and G2M scores, Phase) — when perform.cell.cycle.scoring = TRUE.
• MALAT1 thresholding — when perform.MALAT1.test = TRUE; adds MALAT1.threshold (numeric) and MALAT1.pass (logical) columns. See the feature_threshold repository for details.

The function supports split Seurat v5 objects: if normalized counts are stored in standard data.* layers or a concrete layer vector is passed via layers, metrics are computed per-layer.

SeuratObject <- CalculateQC(
  SeuratObject,
  assay = "RNA",
  layers = NULL,
  species = "human" , # Also 'mouse' available
  perform.cell.cycle.scoring = TRUE,
  perform.MALAT1.test = TRUE
)

Creates a bar plot showing the number of cells per group, with bars filled by a viridis gradient based on cell count.

CellsHistoGradient(
  SeuratObject,
  group.by = NULL,
  scale.colors = "viridis",
  breaks = scales::extended_breaks()
)

Extends Seurat::FeatureScatter() by coloring points with a continuous viridis gradient defined by a third feature. Supports optional per-group panels with independent per-group correlation subtitles, configurable correlation methods, per-feature layer selection for all three axes, and a single shared legend across grouped panels. The gradient scale is computed globally across all panels.

FeatureScatterGradient(
  SeuratObject,
  feature1,
  feature2,
  gradient,
  group.by = "ident",
  scale.colors = "viridis",
  lower.limit = NULL,
  upper.limit = NULL,
  corr.method = "pearson",
  layer1 = NULL,
  layer2 = NULL,
  layer.gradient = NULL,
  plot.title = NULL,
  pt.size = 0.5,
  common.scales = TRUE,
  collect.axes = FALSE
)

Creates violin plots colored by a per-identity aggregated gradient ("nCells" or the mean of a chosen feature), ordered by gradient value. Supports per-feature assay layer selection via layer and optional per-feature custom titles via plot.title. When layer is specified for assay-backed features, the default panel title is feature_layer; metadata-backed features always use bare names.

VlnPlotGradient(
  SeuratObject,
  features,
  gradient = "nCells",
  group.by = "ident",
  scale.colors = "viridis",
  lower.limit = 0,
  upper.limit = NULL,
  pt.size = 0,
  ncol = NULL,
  layer = NULL,
  layer.gradient = NULL,
  plot.title = NULL
)

Creates kernel-density plots for one or more features extracted from a Seurat object via FetchData(), supporting metadata columns, assay features from any layer, and reduction variables. Features grouped by group.by can be shown as overlaid densities or split into per-group panels. Supports dashed red reference lines (vline), independent black median lines (plot.median), custom titles, and per-feature layer selection. When multiple features are requested, returns a named list of plots; names follow feature_layer for assay-backed features and bare feature for metadata-backed features (suffixes _NULL and _NA are stripped).

FeatureDensityPlot(
  SeuratObject,
  features,
  group.by = "ident",
  split.plot = TRUE,
  scale.colors = "viridis",
  ncol = NULL,
  vline = NULL,
  layer = NULL,
  plot.median = TRUE,
  plot.title = NULL,
  nmad = 2,
  alpha = 0.3,
  pt.size = 0,
  common.scales = TRUE,
  collect.axes = FALSE
)

Creates boxplots to visualize the distribution of QC metrics (default: nFeature_RNA, nCount_RNA, percent.mt) grouped by an entity metadata column. Individual cells are overlaid as jitter points, optionally colored by a gradient feature. Multiple metrics are combined into a single patchwork figure.

QCMetricsBoxplot(
  SeuratObject,
  entity_type,
  entity_name = NULL,
  qc_metrics = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
  gradient_col = NULL,
  scale.colors = "viridis",
  lower.limit = 0,
  upper.limit = NULL,
  pt.size = 1,
  pt.alpha = 0.6,
  fill_color = "lightblue",
  outlier.size = 1,
  ncol = NULL,
  gradient_legend = FALSE
)

Runs fgsea per cluster using the output of Seurat::FindAllMarkers() and returns a list of enriched pathways for each cluster. Genes are ranked by avg_log2FC; only genes passing the padj.threshold are included in the ranked list.

results <- scGSEAmarkers(
  cluster_markers,
  reference_markers,
  padj.threshold = 1e-6,
  only.pos = TRUE,
  workers = 4
)

Extracts the per-cell results of .CalculateFeatureThresholdSeurat() for a given feature from a Seurat object's metadata and returns a tidy data.frame with one row per unique batch level, including the computed threshold (optional) and pass/fail statistics.

ExtractFeatureTestResults(
  SeuratObject,
  batch.col,
  feature = "MALAT1",
  extract.threshold = TRUE
)

The following functions are internal helpers not exported from the namespace. They are documented here for reference.

Computes an expression threshold for a single feature from a normalized counts vector using kernel density estimation and smoothing splines. The threshold corresponds to the left quadratic intercept of the density curve, with fallback to a conservative default value.

Applies .ComputeFeatureThreshold() across all layers of a Seurat assay and stores per-cell results as <feature>.threshold (numeric) and <feature>.pass (logical) in meta.data.

Converts ENSEMBL gene IDs to gene symbols for a raw count matrix using biomaRt. Returns the filtered matrix with symbol rownames. Requires biomaRt and dplyr; not exported.

Sets common x and y axis limits across a patchwork joined ggplot object. Used internally by FeatureScatterGradient() and FeatureDensityPlot() when common.scales = TRUE.

Collapses multiple boolean columns into a single logical column, returning TRUE if all of the input columns are TRUE for a given row, and FALSE if any is FALSE. Used internally by FeatureTest()

• CalculateCDR adapted from conquer_comparison by Charlotte Soneson.

• MALAT1 thresholding logic adapted from [BaderLab's MALAT1_threshold] (https://github.com/BaderLab/MALAT1_threshold)

GPLv3 © José Manuel Gómez Silva

If you find bugs or have suggestions, please open an issue or pull request on GitHub.

When contributing code:

• Follow the existing style (based partially on Google's R style guide: CamelCase for exported functions, dot.case for arguments/variables, explicit package::function calls).
• Include clear roxygen documentation and tests where appropriate.
• Respect licensing and attribution, particularly for any upstream code you incorporate.