Skip to contents

Extracting signatures de novo

Source: vignettes/signature-extraction.Rmd
signature-extraction.Rmd

ms_extract() discovers signatures de novo from a catalog. This guide covers the K-selection evidence, the default consensus pipeline, and swapping extraction engines through the registry. Building the catalog is the subject of Mutation catalogs; using a known reference instead is the subject of Fitting known signatures.

A catalog to work against

A small synthetic cohort seeded with three well-separated bundled COSMIC signatures keeps the guide offline and fast:

ref <- ms_refdb_bundled()$SBS_GRCh38@matrices[[1]]
S   <- ref[, c(1, 5, 9)]
colnames(S) <- paste0("CS", 1:3)
E   <- matrix(rep(c(0.5, 0.3, 0.2), each = 24), 3, 24)
sim <- ms_simulate(S, exposures = E, burden = 4000, seed = 7)
sim
#> <MsCatalog> SBS96: 96 channels x 24 samples

K-selection evidence

ms_select_k() gathers two evidence layers over a rank grid: a cross-validation reconstruction curve and a stability panel (multiple replicates per rank, clustered and matched). The arbitration rule is a hard gate: a rank is only selected if it passes the stability veto, and an undecidable grid raises a structured error carrying the full evidence table rather than silently picking a winner.

err <- tryCatch(
  ms_select_k(sim, k_grid = 2:5, replicates = 4, max_iter = 80,
              seed = 1, k_folds = 3, n_seeds = 3),
  error = function(e) e)
class(err)[1]
#> [1] "msuiter_error_kselect"

The tiny synthetic grid above is undecidable — the veto rejects every rank, and the error says exactly why each candidate failed. The evidence table rides on the condition:

ev <- err$k_evidence
head(ev[, c("k", "cv_test_deviance", "avg_stability", "passes_veto",
            "veto_reason")])
#>   k cv_test_deviance avg_stability passes_veto
#> 1 2         288.2824     0.3198347       FALSE
#> 2 3         312.5883     0.5600331       FALSE
#> 3 4         322.5716     0.3730421       FALSE
#> 4 5         324.2326     0.4335074       FALSE
#>                                                                                                veto_reason
#> 1                                     avg_stability 0.3198 < 0.8; avg + min_cluster_stability 0.6147 < 1.0
#> 2                                       avg_stability 0.56 < 0.8; avg + min_cluster_stability 0.9698 < 1.0
#> 3  avg_stability 0.373 < 0.8; min_cluster_stability 0.1338 < 0.2; avg + min_cluster_stability 0.5069 < 1.0
#> 4 avg_stability 0.4335 < 0.8; min_cluster_stability 0.1826 < 0.2; avg + min_cluster_stability 0.6161 < 1.0

This is deliberate: on a small or noisy cohort, K is genuinely undecidable and the package refuses to pretend otherwise. In production analyses the fix is exactly what the error message says — widen the rank grid, raise the replicate count, and read plot_k_cv_curve() / plot_k_stability() / plot_reconstruction_panel() (see The visualization gallery).

The default pipeline

ms_extract() with no method runs the consensus pipeline: an ensemble of NMF replicates, Hungarian-matched consensus clustering, and a refit of the consensus matrix. Pass the K you adjudicated (here the generative truth, since the grid above was refused):

ext <- ms_extract(sim, k = 3)
ext
#> <MsSignature> 3 signatures x 24 samples | engine: nmf-pipeline | seed: 1

The result is an MsSignature whose slots carry the extracted matrix, the per-sample exposures, the pipeline stability evidence, and a summary of the catalog it was extracted from (channel name and hash — the provenance pin used by export and comparison):

dim(ext@signatures)
#> [1] 96  3
head(ext@exposures[, 1L:4L])
#>            sim1      sim2   sim3     sim4
#> Sig1 0.33698814 0.2077078 0.0000 1.689444
#> Sig2 0.03106471 0.0000000 2.2794 0.000000
#> Sig3 0.20429735 0.4051029 0.0000 0.000000
ext@catalog_summary$channel_name
#> [1] "SBS96"

The reconstruction quality of the extraction can be eyeballed against the input catalog with plot_catalog_profile() (the visualization gallery shows it).

Swapping engines

Extraction methods are engines registered in the package registry. The default pipeline is consensus-based; a single-method extraction passes a spec — for example the NMF engine with its update rule and iteration budget:

ext_nmf <- ms_extract(sim, k = 3, method = ms_nmf(engine = "kl",
                                                  max_iter = 80,
                                                  seed = 1))
dim(ext_nmf@signatures)
#> [1] 96  3

The registry also carries the ARD and sparse engines (ms_ard(), ms_sparse()), and third-party packages can register their own — the contract and the onboarding checklist live in the extending guide (Extending msuiter). Because every engine speaks the same MsSignature output contract, comparing engines is a one-line swap, and the results can be compared against each other or against a reference with ms_compare() (Comparing signatures).

Where to go next