Every downstream claim — an extracted signature, a fitted exposure, a presence call — inherits the quality of the catalog it came from. The QC faces inspect that inheritance before the analysis: mutation burden, strand artifacts, flat “pseudosignature” contributions from technical processes, and the hypermutant samples that can dominate an ensemble. Build catalogs with Mutation catalogs.
The QC report
ms_qc_report() computes the per-sample QC table on a
catalog. With a signature set attached it also scores the
pseudosignature faces; without one it reports the burden and strand
faces alone:
ref <- ms_refdb_bundled()$SBS_GRCh38@matrices[[1]]
S <- ref[, c(1, 5, 9)]
colnames(S) <- paste0("CS", 1:3)
E <- matrix(c(0.5, 0.3, 0.2), 3, 16,
dimnames = list(colnames(S), paste0("S", 1:16)))
sim <- ms_simulate(S, exposures = E, burden = 2000, seed = 2)
qc <- ms_qc_report(sim)
head(qc)
#> sample total log10_total burden_class gt_share ct_share strand_bias_p
#> S1 S1 2000 3.30103 standard 0.0500 0.5860 NA
#> S2 S2 2000 3.30103 standard 0.0465 0.5765 NA
#> S3 S3 2000 3.30103 standard 0.0560 0.5920 NA
#> S4 S4 2000 3.30103 standard 0.0485 0.5800 NA
#> S5 S5 2000 3.30103 standard 0.0565 0.5795 NA
#> S6 S6 2000 3.30103 standard 0.0560 0.5710 NA
#> artifact_share artifact_sigs qc_flag
#> S1 NA <NA> watch
#> S2 NA <NA> watch
#> S3 NA <NA> watch
#> S4 NA <NA> watch
#> S5 NA <NA> watch
#> S6 NA <NA> watchEach row is a sample: the burden columns (total,
log10_total, burden_class), the strand-bias
faces, and (with signatures attached) the artifact share with the
specific technical signatures contributing. qc_flag
condenses the table into the per-sample verdict your pipeline can filter
on.
The sibling sentries
The report is the composite; the sentries are the individual faces, exposed so a pipeline can consume exactly one claim without recomputing the rest:
qb <- ms_qc_burden(sim)
qb
#> sample total log10_total burden_class
#> S1 S1 2000 3.30103 standard
#> S2 S2 2000 3.30103 standard
#> S3 S3 2000 3.30103 standard
#> S4 S4 2000 3.30103 standard
#> S5 S5 2000 3.30103 standard
#> S6 S6 2000 3.30103 standard
#> S7 S7 2000 3.30103 standard
#> S8 S8 2000 3.30103 standard
#> S9 S9 2000 3.30103 standard
#> S10 S10 2000 3.30103 standard
#> S11 S11 2000 3.30103 standard
#> S12 S12 2000 3.30103 standard
#> S13 S13 2000 3.30103 standard
#> S14 S14 2000 3.30103 standard
#> S15 S15 2000 3.30103 standard
#> S16 S16 2000 3.30103 standard
qs <- ms_qc_sentries(sim)
qs
#> sample gt_share ct_share gt_verdict ct_verdict qc_flag
#> S1 S1 0.0500 0.5860 ok watch watch
#> S2 S2 0.0465 0.5765 ok watch watch
#> S3 S3 0.0560 0.5920 ok watch watch
#> S4 S4 0.0485 0.5800 ok watch watch
#> S5 S5 0.0565 0.5795 ok watch watch
#> S6 S6 0.0560 0.5710 ok watch watch
#> S7 S7 0.0510 0.5780 ok watch watch
#> S8 S8 0.0405 0.6015 ok watch watch
#> S9 S9 0.0465 0.5880 ok watch watch
#> S10 S10 0.0440 0.5885 ok watch watch
#> S11 S11 0.0600 0.5565 ok watch watch
#> S12 S12 0.0455 0.5710 ok watch watch
#> S13 S13 0.0505 0.5680 ok watch watch
#> S14 S14 0.0490 0.5940 ok watch watch
#> S15 S15 0.0465 0.5960 ok watch watch
#> S16 S16 0.0560 0.5770 ok watch watchms_qc_mechanism_channels() exposes the channel-level
mechanism roster — which channels each certified mechanism (APOBEC,
tobacco, UV, and the rest of the roster) is read from. It is the
reference table behind the mechanism columns above, and is also useful
standalone when you want to reason about channel subsets:
mc <- ms_qc_mechanism_channels()
head(as.data.frame(mc))
#> label mechanism
#> 1 A[C>A]A oxoG
#> 2 A[C>A]C oxoG
#> 3 A[C>A]G oxoG
#> 4 A[C>A]T oxoG
#> 5 A[C>T]A ffpe
#> 6 A[C>T]C ffpe
#> citation
#> 1 Costello et al., NAR 2013 (doi:10.1093/nar/gks1443)
#> 2 Costello et al., NAR 2013 (doi:10.1093/nar/gks1443)
#> 3 Costello et al., NAR 2013 (doi:10.1093/nar/gks1443)
#> 4 Costello et al., NAR 2013 (doi:10.1093/nar/gks1443)
#> 5 Guo et al., Nat Commun 2022 (doi:10.1038/s41467-022-32041-5)
#> 6 Guo et al., Nat Commun 2022 (doi:10.1038/s41467-022-32041-5)plot_qc_dashboard() and
plot_pseudosig_contribution() render the report’s faces
(shown in the visualization
gallery).
Hypermutant stratification
A handful of extreme-burden samples can drive an entire ensemble. The stratifier splits the cohort into a hypermutant class and a regular class — by a data-driven cutoff or a manual burden threshold — so the two arms can be extracted, fitted, and reported separately:
strat <- ms_stratify_hypermutants(sim, manual_cutoff = 3000, seed = 1)
strat
#> <MsStratification> GMM cutoff 3000: 0/16 hypermutated; excluded from de novo, refit pending (M2)The result is an MsStratification carrying the class
assignment and the cutoff used; plot_hypermutator_burden()
draws the burden distribution with the cutoff (the visualization gallery).
Downsampling
Burden drives most QC verdicts, so a useful question is: how does the
analysis behave at lower sequencing depth? ms_downsample()
thins a catalog to a target per-sample total — the primitive behind the
calibration sanity protocol of Calibrated statistics, and a
useful robustness check in its own right:
ds <- ms_downsample(sim, depth = 800, seed = 1)
range(colSums(ds@counts))
#> [1] 800 800Where to go next
- Extracting signatures: the de novo pipeline the QC faces gate.
- Calibrated statistics: the sanity protocol built on downsampling.
- The visualization gallery: the QC figure types.