Introduction

Mutational signature analysis is commonly formulated as a non-negative factorization of a $96 \times J$ matrix of mutation counts, whose columns are tumours and whose rows are trinucleotide substitution types. Aggregating a genome into 96 counts discards the genomic position of each mutation, while mutation density varies by orders of magnitude along the genome with chromatin state, replication timing and sequence composition.

In SignaturePPF, mutations are modelled as an inhomogeneous Poisson point process along the genome whose intensity depends on genomic covariates through a log-linear term, applied to a signature factorization of the 96 substitution channels.

This tutorial illustrates how to use the package on whole-genome sequences mutation data.

Installation

# install.packages("remotes")
remotes::install_github("alessandrozito/SignaturePPF", build_vignettes = TRUE)

library(SignaturePPF)
suppressPackageStartupMessages({
  library(GenomicRanges)
  library(ggplot2)
})

The model

Mutations as a point process

Let $t \in [0,T)$ index position along the genome, $i = 1,\dots,96$ index the substitution channels, and $j = 1,\dots,J$ index the tumours. The mutations of tumour $j$ in channel $i$ are modelled as an inhomogeneous Poisson process with intensity function

\[ \lambda_{ij}(t) \;=\; \tfrac{1}{2}\,c_j(t)\sum_{k=1}^{K} r_{ik}\,\phi_{kj}\,e^{\beta_k^{\top}x(t)} . \]

Here

$r_{\cdot k}\in\Delta^{96}$	signature $k$, a probability vector over the 96 channels
$\phi_{kj}\ge 0$	baseline activity of signature $k$ in tumour $j$
$x(t)\in\mathbb{R}^{p}$	genomic covariates at position $t$, shared across tumours
$\beta_k\in\mathbb{R}^{p}$	covariate coefficients of signature $k$
$c_j(t)$	copy number of tumour $j$ at position $t$

Both $x(t)$ and $c_j(t)$ are treated as fixed and known. The factor $c_j(t)/2$ is the ratio of observed to expected copy number, so it equals one at a diploid locus.

The term $e^{\beta_k^{\top}x(t)}$ modulates the intensity of signature $k$ along the genome, and each signature has its own coefficient vector. If $\beta_k = 0$ for all $k$ and $c_j(t) = 2$ throughout, the integrated intensity reduces to the Poisson non-negative matrix factorization.

The activity prior

Define

\[ q_j(\beta_k) \;=\; \int_0^T \tfrac{1}{2}\,c_j(t)\,e^{\beta_k^{\top}x(t)}\,\mathrm{d}t . \]

Integrating the intensity over the genome shows that the expected number of mutations generated by signature $k$ in tumour $j$ equals $\theta_{kj} = \phi_{kj}\,q_j(\beta_k)$. The prior is specified for $\theta_{kj}$,

\[ \theta_{kj}\mid\mu_k \;\sim\; \mathrm{Ga}(a,\ a/\mu_k), \qquad a = 1.01 \text{ by default}, \]

or equivalently $\phi_{kj}\mid\mu_k,\beta_k \sim \mathrm{Ga}\big(a,\ a\,q_j(\beta_k)/\mu_k\big)$. This is referred to as the activity prior.

Automatic selection of the number of signatures via compressive hyperpriors

The number of signatures $K$ is specified as an upper bound. The relevance weights $\mu_k$ are assigned the prior

\[ \mu_k \sim \mathrm{InvGa}(a_0, b_0), \qquad a_0 = aJ+1,\quad b_0 = \varepsilon\,aJ, \]

with $\varepsilon = 10^{-3}$ by default. This places prior mass near $\varepsilon$, so that for a signature unsupported by the data the posterior for $\mu_k$ concentrates near $\varepsilon$ and $\theta_{kj}$ is shrunk towards zero for every $j$. Such signatures are identified and removed by prune_signatures(). A value of $K$ larger than the number of signatures present is therefore admissible.

The specification is completed by

\[ r_{\cdot k} \sim \mathrm{Dir}(\alpha \mathbf{1}_{96}), \qquad \beta_{k\ell}\mid\sigma^2_k \sim N(0,\sigma^2_k), \qquad \sigma^2_k \sim \mathrm{InvGa}(c_0,d_0), \]

with $\alpha = 1.01$. All hyperparameters are set through SignaturePPF_prior().

Data formats in SignaturePPF

Three objects, and one alignment rule between them.

`gr_Mutations`	a `GRanges`, one range per mutation, with `sample` and `channel` in `mcols` plus one numeric column per covariate
`SignalTrack`	$\text{bins}\times p$, the covariates over the genome
`CopyTrack`	$\text{bins}\times J$, copy number times usable sequence in each bin

The numeric columns of mcols(gr_Mutations) must be the same covariates, in the same order, as the columns of SignalTrack, and SignalTrack and CopyTrack must describe the same bins in the same order. If the ordering differs, each coefficient is associated with the wrong covariate.

Simulating a cohort

The cohort consists of eight tumours on a 4 Mb chromosome divided into 1 kb bins, with three covariates. Mutations are generated from three COSMIC signatures: the clock-like signature SBS1 and the APOBEC signatures SBS2 and SBS13. The estimates can then be compared with these.

set.seed(1)
nbins <- 4000; J <- 8; bin_width <- 1000

truth_sigs  <- c("SBS1", "SBS2", "SBS13")
R           <- COSMIC_v3.4_SBS96_GRCh37[, truth_sigs]
channels    <- rownames(R)
samples     <- sprintf("S%02d", seq_len(J))
covs        <- c("H3K9me3", "RepliTime", "GC")

# Covariates. Real tracks are spatially correlated along the genome, so each is
# drawn as an AR(1) sequence over bins and then standardised, as the real tracks
# are after preprocessing.
ar1 <- function(n, rho) {
  as.numeric(stats::filter(rnorm(n, 0, sqrt(1 - rho^2)), rho,
                           method = "recursive", init = rnorm(1)))
}
SignalTrack <- scale(sapply(covs, function(nm) ar1(nbins, rho = 0.97)))
attributes(SignalTrack)[c("scaled:center", "scaled:scale")] <- NULL
colnames(SignalTrack) <- covs

# Copy number x bin width: the exposure c_j(t) the process integrates over.
CopyTrack <- matrix(bin_width * pmax(rnorm(nbins * J, 1, 0.15), 0.1),
                    nbins, J, dimnames = list(NULL, samples))

Each signature is assigned coefficients $\beta_k$ and an activity $\theta_{kj}$ in each tumour. The three coefficient vectors are chosen to be distinct: SBS1 increases with H3K9me3 and decreases with replication timing, SBS2 has the opposite signs, and SBS13 depends only on GC content.

Betas <- matrix(c( 0.8, -0.5,  0.0,     # SBS1
                  -0.6,  0.4,  0.3,     # SBS2
                   0.0,  0.0, -0.7),    # SBS13
                nrow = length(covs), ncol = 3,
                dimnames = list(covs, truth_sigs))
Theta <- matrix(rgamma(3 * J, 8, 8 / 500), 3, J,
                dimnames = list(truth_sigs, samples))
round(Theta[, 1:4])
#>       S01 S02 S03 S04
#> SBS1  476 543 594 669
#> SBS2  516 521 400 385
#> SBS13 680 259 307 395

On a binned genome the integral is a sum, so signature $k$ contributes to bin $b$ of tumour $j$ a Poisson count with mean

\[ \theta_{kj}\;\frac{c_j(b)\,e^{\beta_k^{\top}x(b)}}{\sum_{b'} c_j(b')\,e^{\beta_k^{\top}x(b')}} , \]

that is, $\theta_{kj}$ distributed across bins in proportion to the exposure times the covariate effect. Theta determines the number of mutations, Betas their distribution along the genome.

EBS <- exp(SignalTrack %*% Betas)      # bins x K
Q   <- crossprod(EBS, CopyTrack)       # K x J

rows <- list()
for (j in seq_len(J)) for (k in seq_len(3)) {
  n   <- rpois(nbins, Theta[k, j] * CopyTrack[, j] * EBS[, k] / Q[k, j])
  idx <- rep(seq_len(nbins), n)
  if (!length(idx)) next
  rows[[length(rows) + 1L]] <- data.frame(
    bin     = idx,
    sample  = samples[j],
    channel = channels[sample.int(nrow(R), length(idx), TRUE, prob = R[, k])])
}
mut <- do.call(rbind, rows)

The GRanges is assembled with the covariates evaluated at each mutation, taken from the bin containing it and ordered as the columns of SignalTrack.

gr <- GRanges("chr1", IRanges((mut$bin - 1L) * bin_width + 1L, width = 1L))
mcols(gr) <- cbind(
  data.frame(sample  = factor(mut$sample,  levels = samples),
             channel = factor(mut$channel, levels = channels)),
  as.data.frame(SignalTrack[mut$bin, , drop = FALSE]))

data <- list(gr_Mutations = gr, SignalTrack = SignalTrack, CopyTrack = CopyTrack)
length(gr)
#> [1] 12033

What the covariates did

Aggregating the cohort into 100 kb windows shows only part of the structure. Total mutation count is clearly associated with GC content, but shows almost no association with H3K9me3 or with replication timing, despite both entering the intensity with coefficients of magnitude 0.4 to 0.8.

The reason is that the signatures respond in opposite directions. SBS1 increases with H3K9me3 and SBS2 decreases with it, and the two have comparable activities, so the effects largely cancel once the signatures are summed. Only GC content, where SBS13 dominates, survives aggregation. Recovering the individual effects requires a model in which each signature carries its own coefficients.

win  <- 100                                  # 100 bins = 100 kb
grp  <- ceiling(seq_len(nbins) / win)
cnt  <- as.vector(table(factor(ceiling(mut$bin / win), levels = unique(grp))))
Xwin <- apply(SignalTrack, 2, function(z) tapply(z, grp, mean))

ggplot(data.frame(count     = rep(cnt, times = length(covs)),
                  value     = as.vector(Xwin),
                  covariate = rep(covs, each = length(cnt))),
       aes(value, count)) +
  geom_point(alpha = 0.5, size = 1, colour = "#2E6189") +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
              colour = "#B3312C", linewidth = 0.7) +
  facet_wrap(~ covariate) +
  labs(x = "covariate in 100 kb window (standardised)", y = "mutations") +
  theme_bw(base_size = 11)

Checking the data

SignaturePPF_validate() checks the requirements above and returns the data in the form the samplers use. It is called internally by SignaturePPF(), but running it first is inexpensive and reports any inconsistency immediately.

v <- SignaturePPF_validate(data)
str(v[c("I", "J", "p", "N", "nbins")])
#> List of 5
#>  $ I    : int 96
#>  $ J    : int 8
#>  $ p    : int 3
#>  $ N    : int 12033
#>  $ nbins: int 4000

De novo estimation

The model is fitted with $K = 6$, above the number of signatures present.

fit <- SignaturePPF(data, K = 6, method = "map", seed = 1, verbose = FALSE)
fit <- prune_signatures(fit)
fit
#> Poisson process factorization of mutational signatures
#>   fit    : denovo, MAP
#>   data   : 12,033 mutations, 8 samples, 96 channels, 3 covariates, 4,000 bins
#>   K      : 3 of 6 fitted (3 pruned at Mu <= 0.01)
#>   optim  : 240 iterations, relative change 3.01e-07
#>   time   : 00:00:01

Three of the six are retained. Cosine similarity to the signatures used to simulate the data:

cosine <- function(a, b) sum(a * b) / sqrt(sum(a^2) * sum(b^2))
M <- outer(seq_len(3), seq_len(ncol(fit$Signatures)),
           Vectorize(function(i, j) cosine(R[, i], fit$Signatures[, j])))
dimnames(M) <- list(truth_sigs, colnames(fit$Signatures))
round(M, 3)
#>       SigN02 SigN03 SigN04
#> SBS1   0.999  0.011  0.016
#> SBS2   0.020  0.006  1.000
#> SBS13  0.016  1.000  0.018

Each of the three is recovered by a single column, at cosine similarity 0.999 or above.

The estimated coefficients, with columns reordered to match the simulated signatures:

ord <- apply(M, 1, which.max)
round(Betas, 2)
#>           SBS1 SBS2 SBS13
#> H3K9me3    0.8 -0.6   0.0
#> RepliTime -0.5  0.4   0.0
#> GC         0.0  0.3  -0.7
round(fit$Betas[, ord], 2)
#>           SigN02 SigN04 SigN03
#> H3K9me3     0.77  -0.57  -0.01
#> RepliTime  -0.49   0.37   0.00
#> GC          0.00   0.28  -0.74

plot_Signatures(fit)

plot_Betas(fit)

Rows are signatures and columns covariates, corresponding to the matrix printed above.

Fitted mutation rate along the genome

Summing the fitted intensity over signatures and tumours gives the expected number of mutations in each bin,

\[ \widehat{\Lambda}(b) \;=\; \sum_{j=1}^{J} c_j(b) \sum_{k=1}^{K} \widehat{\phi}_{kj}\, e^{\widehat{\beta}_k^{\top} x(b)} , \]

where $\widehat{\phi}_{kj} = \widehat{\theta}_{kj} / q_j(\widehat{\beta}_k)$ is returned in fit$Baseline.

lambda <- rowSums((exp(SignalTrack %*% fit$Betas) %*% fit$Baseline) * CopyTrack)
observed <- tabulate(mut$bin, nbins)

c(observed = sum(observed), fitted = round(sum(lambda), 1))
#> observed   fitted 
#>  12033.0  11978.8

Aggregated into 50 kb windows and plotted along the chromosome, against the observed counts:

w    <- 50
grpw <- ceiling(seq_len(nbins) / w)
track <- data.frame(
  position = as.vector(tapply(seq_len(nbins), grpw, mean)) * bin_width / 1e6,
  observed = as.vector(tapply(observed, grpw, sum)),
  fitted   = as.vector(tapply(lambda,   grpw, sum)))

ggplot(track, aes(position)) +
  geom_point(aes(y = observed), shape = 21, size = 1.4, colour = "grey30") +
  geom_line(aes(y = fitted), colour = "#B3312C", linewidth = 0.7) +
  labs(x = "position on chr1 (Mb)", y = "mutations per 50 kb") +
  theme_bw(base_size = 11)

round(cor(track$observed, track$fitted), 3)
#> [1] 0.924

The fitted curve follows the observed counts because the covariates vary along the genome and the estimated $\beta_k$ reproduce their effect. A model with $\beta_k = 0$ would predict a rate proportional to copy number alone, which here is close to constant.

Posterior inference

Setting method = "mcmc" samples from the posterior. The $\beta$ block is updated by adaptive generalized elliptical slice sampling, which adapts the sampling ellipse using the chain rather than drawing it from the prior. Initialising the chain at the posterior mode shortens the burn-in required.

fitm <- SignaturePPF(data, K = 6, method = "mcmc", seed = 1, verbose = FALSE,
                     init_mcmc_from_map = TRUE, prune_after_map = TRUE,
                     controls = SignaturePPF_control(nsamples = 2000,
                                                     burnin = 500, thin = 2))
fitm
#> Poisson process factorization of mutational signatures
#>   fit    : denovo, MCMC (agess)
#>   data   : 12,033 mutations, 8 samples, 96 channels, 3 covariates, 4,000 bins
#>   K      : 3
#>   chain  : 2000 iterations, thin 2, burn-in 500, 750 draws summarised
#>   time   : 00:00:05

posterior_CI() returns the posterior mean and credible bounds for a given block of parameters:

ci <- posterior_CI(fitm, "Betas")
round(ci$lowCI, 2)
#>           SigN02 SigN03 SigN04
#> H3K9me3     0.75  -0.05  -0.61
#> RepliTime  -0.53  -0.04   0.35
#> GC         -0.04  -0.80   0.25
round(ci$highCI, 2)
#>           SigN02 SigN03 SigN04
#> H3K9me3     0.82   0.03  -0.54
#> RepliTime  -0.47   0.04   0.40
#> GC          0.03  -0.72   0.32

plot_Betas(fitm)

For an MCMC fit the credible intervals are used in the display: coefficients whose interval contains zero are shown in grey.

Refitting a catalogue

Supplying sigs holds the signatures fixed at reference values and estimates only the activities and the coefficients. Five reference signatures are supplied below, three of which were used to generate the data:

sigs <- COSMIC_v3.4_SBS96_GRCh37[, c("SBS1", "SBS2", "SBS3", "SBS5", "SBS13")]
fitr <- SignaturePPF(data, sigs = sigs, method = "map", seed = 1, verbose = FALSE)
fitr
#> Poisson process factorization of mutational signatures
#>   fit    : refit, MAP
#>   data   : 12,033 mutations, 8 samples, 96 channels, 3 covariates, 4,000 bins
#>   K      : 3 of 5 fitted (2 pruned at Mu <= 0.01)
#>   optim  : 30 iterations, relative change 0.00e+00
#>   time   : 00:00:00

The compressive prior applies to reference signatures as well, so SBS3 and SBS5 are shrunk to $\varepsilon$ and pruned:

df_assign(fitr, data)
#>   best_sig    m       mu
#> 1     SBS1 4310 235.2089
#> 2    SBS13 3186 178.4378
#> 3     SBS2 4537 252.5742

This differs from non-negative least squares attribution, under which every supplied reference signature receives a non-zero activity.

plot_Mu(fitr, data)

A semi-supervised specification is obtained with sigs_fixed = FALSE. Each reference signature $s_h$ is then assigned the informative prior $r_{\cdot h} \sim \mathrm{Dir}(b_h s_h)$, centred at $s_h$ with concentration $b_h$, and $K$ additional signatures are estimated de novo. Setting betah = "auto" selects $b_h$ for each signature as the concentration at which a prior draw has median cosine similarity $0.975$ to $s_h$. See ?SignaturePPF and ?tune_betah.

Activities

plot_Theta() displays $\theta_{kj}$, the expected number of mutations per signature and tumour.

plot_Theta(fit)

Where to go next

SignaturePPF_prior() and SignaturePPF_control() hold the hyperparameters and the algorithm settings.
SignaturePPF_checkpoint() makes long chains resumable. A whole-genome fit requires several hours; re-issuing the same call after an interruption resumes the chain exactly.
The preprocessing that turns real chromatin tracks into SignalTrack and CopyTrack, and the full analysis for the paper, are in SignaturePPF-paper.

References

Zito, A., Parmigiani, G. and Miller, J. W. (2025). Poisson process factorization for modeling mutational processes along cancer genomes. arXiv:2510.26090

Marco, N. and Tokdar, S. T. (2026). Adaptive generalized elliptical slice sampling. arXiv:2605.21659

Alexandrov, L. B. et al. (2020). The repertoire of mutational signatures in human cancer. Nature 578, 94–101.

\(r_{\cdot k}\in\Delta^{96}\)	signature \(k\), a probability vector over the 96 channels
\(\phi_{kj}\ge 0\)	baseline activity of signature \(k\) in tumour \(j\)
\(x(t)\in\mathbb{R}^{p}\)	genomic covariates at position \(t\), shared across tumours
\(\beta_k\in\mathbb{R}^{p}\)	covariate coefficients of signature \(k\)
\(c_j(t)\)	copy number of tumour \(j\) at position \(t\)

A tutorial of SignaturePPF

Poisson process factorization for mutational signatures with genomic covariates