SciBet reference training + cell type prediction

Overview

SciBet (scibet, Zhang Lab / Peking University) is a supervised cell type identifier: train a
multinomial model on a labelled reference dataset, then assign a cell type to every cell of a query
dataset. It works directly on the expression matrix (no PCA/clustering embedding is used), it is
order-of-magnitude faster than correlation/nearest-neighbour label transfer, and a trained model is
tiny (roughly <1 MB for 100 cell types x 1000 genes), so it is well suited to
cross-dataset / cross-platform / cross-species annotation.

The pipeline is three steps:

SelectGene(expr, k)                  -> informative gene names   (supervised, entropy reduction)
gindex  <- which(colnames==gene)     -> 1-based positions
model   <- Learn(expr, geneset=gindex-1)   -> 0-based positions, returns a predictor function
model(query_expr)                    -> predicted cell type per cell

Reference: Li et al., Nature Communications 2020, doi:10.1038/s41467-020-15523-2;
source <https://github.com/PaulingLiu/scibet>; guided tutorial <http://scibet.cancer-pku.cn/document.html>.

Inputs

• Reference Seurat .rds: must contain a meta.data column with the cell type labels, and an
  expression matrix (log-normalized data layer is the established convention).
• Query Seurat .rds (optional if you only want to train/save a model): same species and same
  gene symbol convention as the reference.
• Optional a previously trained model .rds (the Learn() return value), or a fixed gene list.
• Both objects are read with readRDS(); anything with @meta.data and an extractable assay works.

Steps

The reusable script is scripts/scibet_predict.R.txt — fully command-line parameterised, no hard-coded
paths:

/opt/conda/bin/Rscript scripts/scibet_predict.R.txt \
  --reference-rds /path/ref.rds \
  --reference-label-col subCluster1 \
  --query-rds /path/query.rds \
  --query-label-col Annotation \
  --output-dir ./scibet_out \
  --k 2000

# train once, reuse the model later (second run skips training)
/opt/conda/bin/Rscript scripts/scibet_predict.R.txt \
  --reference-rds ref.rds --output-dir out --k 2000 --model-rds out/ref.model.rds
/opt/conda/bin/Rscript scripts/scibet_predict.R.txt \
  --reference-rds ref.rds --query-rds query.rds --model-rds out/ref.model.rds \
  --output-dir out

/opt/conda/bin/Rscript scripts/scibet_predict.R.txt --help

The script was exercised end to end on a synthetic 3-cell-type benchmark (train+predict, model
reuse, and --retrain), reaching a 180/180 correct confusion matrix.

What it does, and the same steps by hand:

1. Extract the expression matrix. SciBet wants cells as rows, genes as columns, label as the
   last column. Seurat stores genes x cells, so transpose:

r train_set <- tibble::as_tibble(t(as.matrix(train@assays$RNA@data))) train_set[["label"]] <- as.character(train$subCluster1) # MUST be the last column

2. Select informative genes (supervised, by total entropy reduction across cell types):

r genes <- scibet::SelectGene(train_set, k = 2000) # -> character vector of gene names

3. Convert gene names to 0-based column positions (see Pitfalls — this is mandatory):

r gindex <- sapply(genes, function(x) which(colnames(train_set) == x)) gindex <- unname(gindex)

4. Train — Learn() returns a closure, not an S4 object:

r model <- scibet::Learn(train_set, geneset = gindex - 1) # 0-based, exactly as SciBet expects saveRDS(model, "tnk.model.rds")

5. Predict on the query (the closure applies the log transform itself):

r sc_query <- tibble::as_tibble(t(as.matrix(sce.all@assays$RNA@data))) major_type <- model(sc_query) # character vector, one label per cell prob <- model(sc_query, result = "table") # cells x cell types probability matrix sce.all@meta.data$subtype <- major_type table(sce.all@meta.data$Annotation, sce.all@meta.data$subtype) # original vs predicted

6. Check the result: cross-tabulate predicted against the query's own annotation (agreement is
   the sanity metric), inspect the composition, and confirm the model/query gene overlap.

Verified API (read from the package source R/SciBet.R and src/SciBet.cpp):

Parameters

Learn()'s k and a are ignored whenever geneset is supplied — with an explicit gene set the
only thing that matters is the index vector.

Outputs

Written under --output-dir (default .) in the standard three-directory layout — Img/ for
figures, Stat/ for tables/text, Rdata/ for R objects — with prefix --prefix (default: reference
file basename). The script creates all three sub-directories and prints their contents when it exits.

An explicit --model-rds PATH is used verbatim (it is an input as much as an output); only the
default model location follows the Rdata/ convention.

Pitfalls

• ⚠️ Degenerate reference input segfaults the R session (no error message). scibet::SelectGene()
  calls C++ (DsGene) that indexes without bounds checks, so a reference whose expression matrix is
  empty or has a single label level crashes R (An irrecoverable exception occurred. R is aborting now)
  instead of raising an error. Observed failure modes while validating this skill:
  (a) a reference object that never ran Seurat::NormalizeData() — Seurat's data layer is then 0×0,
  so the script reported reference: 0 cells x 0 genes and died inside SelectGene;
  (b) a reference with a single cell type — the entropy difference SelectGene maximises is undefined.
  The script now validates both cases up front and aborts with an actionable message
  (Reference expression matrix is degenerate … / Reference has only 1 label level(s) …).
  When calling the API directly, do the same: check dim(expr) is non-degenerate and that the label
  column has ≥ 2 levels before touching SelectGene().

• ⚠️ The geneset argument is 0-based, so gindex - 1 is mandatory, not cosmetic. SelectGene()
  returns gene names; Learn() does prob <- GenProb(expr, labels_in, geneset, ...) and
  genes <- colnames(expr)[geneset + 1]. The C++ backend (GenProb) casts geneset to
  std::vector<int> and indexes the columns from 0. Passing the 1-based which() result — or the
  gene names — shifts every selected gene by one column or fails outright, silently producing a
  model trained on the wrong genes. Always convert which(colnames(train_set) == x) with - 1.
  Verified on synthetic data (independent cells, 3 cell types): Learn(geneset = gindex - 1) scored
  100% agreement and ExportModel() returned exactly the SelectGene() genes, whereas the 1-based
  vector aborted with REAL() can only be applied to a 'numeric', not a 'character' and then
  segfaulted, and passing the gene names failed with
  Not compatible with requested type: [type=character; target=integer]. The C++ backend indexes
  columns without bounds checking, so the wrong index convention is a crash, not a warning.

• ⚠️ Gene names must match between reference and query. The predictor resolves genes by name
  (have_genes <- which(genes %in% colnames(expr))) and then silently drops whatever is missing —
  there is no error and no warning. A species/case mismatch (human CD3D vs mouse Cd3d) or a
  different annotation build therefore still "runs" and yields meaningless predictions. Always check
  the overlap; the script warns below 50%.
• Train and predict on the same kind of matrix. GenProb sums log(x + 1)/log(2) over the
  reference and the predictor closure applies the identical log1p(x)/log(2) to the query, so SciBet
  transforms internally: do not pre-log your data and do not feed z-scored/ScaleData values.
  The in-house scripts use log-normalized @assays$RNA@data for both reference and query — keep
  that consistent, that is what makes the numbers comparable.
• Seurat v4 vs v5 expression access. train@assays$RNA@data is the classic v4 form of
  GetAssayData(obj, "RNA", slot = "data"). In Seurat v5 the RNA assay is an Assay5 whose data
  lives in @layers$data, so the @data slot access errors. Use
  SeuratObject::LayerData(obj, assay = "RNA", layer = "data") (or GetAssayData(..., layer = "data")).
  The script's get_assay_data() tries the v5 path first and falls back to the v4 slot.
• The label must be the last column and named exactly label. SelectGene/Learn treat
  ncol(expr) as the label (n_gene = LENGTH(expr_r) - 1 in C++), so an extra column anywhere
  corrupts the feature space. Never leave label in the query frame — the query has genes only.
• Labels must be non-missing. factor(expr$label) feeds as.integer(...) directly; NA labels
  break the class numbering. Drop cells with NA/empty labels before training.
• Storage type matters. The C++ reads every column with REAL(), which expects a double
  (REALSXP) vector. A matrix of integers (raw counts kept as integer) is misread. as.matrix()
  on a dgCMatrix already yields doubles; the script forces storage.mode(m) <- "double" as a guard.
• Cell order must line up. as_tibble(t(as.matrix(...))) drops cell names, so assign the label
  vector from the same object in the same order (or index meta.data by the matrix's colnames, as
  the script does).
• t(as.matrix(...)) densifies the matrix. An all-genes dense matrix for a >100k-cell object is
  a large memory allocation; subset to the genes you need, or at least size the job accordingly
  (rm() + gc() after training).
• Do not call scibet::Bet(). It exists only as an empty stub for documentation; the real
  predictor is the closure returned by Learn()/LoadModel().
• Portability / web tool. For a model matrix to load elsewhere (including the JavaScript web
  client) use ExportModel(model) to get the matrix, then LoadModel(x). A plain saveRDS() of the
  closure only works in R.
• Subset the reference if you want a cell-type-specific model; e.g.
  train <- subset(train, subset = maingroups == "Myeloid cells") before SelectGene(). Training on
  the full atlas and querying a subset are two different tasks — keep them straight.

Environment

• R: /opt/conda/bin/Rscript — R 4.5.3. Library search path: /opt/user_R_libs, /opt/public_R_libs,
  /opt/conda/lib/R/library.
• Dependency status in this container (checked with requireNamespace()):

• Install scibet (it is not on CRAN). It is already present in this container; the commands
  below are for a fresh environment. The package compiles C++ (src/Makevars, OpenMP), so build
  tools must be present:

r install.packages(c("Rcpp", "RcppEigen", "ggsci", "viridis", "tidyverse")) remotes::install_github("PaulingLiu/scibet") # default branch: master

For installing R/Bioconductor packages into this container's shared library, see the
abio-r-package-install skill.
- Unofficial alternative implementing the same algorithm: zhangwenjie-tina/scibetR
(SciBet_R(train, test, k = 1000, result = "list")); prefer the original PaulingLiu/scibet.

> Packaging note. Scripts and tables in this skill carry a .txt suffix on their
> original extension (e.g. .R.txt, .tsv.txt) because the SkillHub upload allowlist
> does not accept .R or .tsv. This is cosmetic: Rscript is extension-agnostic, so
> Rscript scripts/<name>.R.txt ... behaves exactly like .R, and tab-separated data
> parses identically from .tsv.txt. Keep the suffix when repackaging.