A robust mutation spectrum analysis framework for high-depth error-corrected sequencing data

Comparing mutational spectra across experimental groups requires inference at the level of the biological replicate. Error-corrected sequencing (ECS) generates hundreds to thousands of mutations per replicate, yet classical pooled G² tests treat individual mutations as independent observations and ignore replicate-level variability, inflating statistical significance (pseudo-replication). We present a replicate-aware Ridge-Penalized Dirichlet-Multinomial (RP-DMM) framework that models biological replicates directly and accounts for overdispersion. An L2 (ridge) penalty, with strength chosen by leave-one-out cross-validation, stabilizes spectrum estimates by shrinking weakly supported mutation categories toward a common baseline, and inference uses a parametric bootstrap. We also provide a complementary non-parametric rank-permutation test for sparse datasets.

In simulations, the pooled G² test became increasingly anti-conservative as mutation depth increased, with false-positive rates approaching 100% under the null hypothesis. RP-DMM maintained Type I error near the nominal 5% level and had substantially greater power than the rank-permutation test. Regularization was most consequential for high-dimensional spectra. False positive rates approached the nominal 5% level once the lower-mutation group contained approximately 300-500 mutations, depending on the number of biological replicates. A practical target of approximately 400 mutations in the lower-mutation group provided robust Type I error control across all scenarios examined; below this, the rank-permutation test performed robustly and is preferred. In Duplex Sequencing data from benzo[a]pyrene-exposed mice, RP-DMM distinguished exposure-associated mutational signatures from comparisons with higher mutation counts but null effect size, whereas the pooled G² test declared all comparisons significant. The cross-validated penalty decreased as mutation depth increased, indicating that regularization strength should be tuned to the data.

The framework is implemented in the MutSeqR R package.