SuppReferences

# Compressive Bayesian non-negative matrix factorization for mutational
signatures analysis

Alessandro Zito  Department of Biostatistics, Harvard University, Boston, MA,
02115, U.S.A.  Jeffrey W. Miller  Department of Biostatistics, Harvard
University, Boston, MA, 02115, U.S.A.

###### Abstract

Non-negative matrix factorization (NMF) is widely used in many applications
for dimensionality reduction. Inferring an appropriate number of factors for
NMF is a challenging problem, and several approaches based on information
criteria or sparsity-inducing priors have been proposed. However, inference in
these models is often complicated and computationally challenging. In this
paper, we introduce a novel methodology for overfitted Bayesian NMF models
using “compressive hyperpriors” that force unneeded factors down to negligible
values while only imposing mild shrinkage on needed factors. The method is
based on using simple semi-conjugate priors to facilitate inference, while
setting the strength of the hyperprior in a data-dependent way to achieve this
compressive property. We apply our method to mutational signatures analysis in
cancer genomics, where we find that it outperforms state-of-the-art
alternatives. In particular, we illustrate how our compressive hyperprior
enables the use of biologically informed priors on the signatures, yielding
significantly improved accuracy. We provide theoretical results establishing
the compressive property, and we demonstrate the method in simulations and on
real data from a breast cancer application.

##  1 Introduction

Non-negative matrix factorization (NMF) is a dimensionality reduction
technique that decomposes a non-negative matrix into the product of two lower-
dimensional non-negative matrices of a desired rank by minimizing a given loss
function, such as the squared-error loss or Kullback–Leibler divergence (Lee
and Seung, 2000). Several Bayesian NMF methods have been studied, using both
parametric (Cemgil, 2009; Zhou, 2018) and nonparametric models (Hoffman et
al., 2010; Gopalan et al., 2014; Ayed and Caron, 2021). See the general
overview of Wang and Zhang (2013), and refer to Zhou and Carin (2015) for a
detailed summary of Bayesian factorization methods for discrete outcomes.

In cancer genomics, non-negative matrix factorization has been successfully
used to discover a wide range of mutational signatures corresponding to
distinct processes, such as damaged DNA repair mechanisms and environmental
mutagens like tobacco smoking or metabolic byproducts (Alexandrov et al.,
2013). These signatures, defined as vectors of the frequencies with which
different types of point mutations occur, are inferred from mutation counts in
whole-genome or whole-exome sequencing using NMF algorithms (Nik-Zainal et
al., 2012, 2016; Alexandrov et al., 2013, 2020). Identifying these patterns in
the DNA of cancer patients is a significant advance toward understanding the
etiology of cancer (Koh et al., 2021), improving the effectiveness of
precision therapies (Aguirre et al., 2018; Gulhan et al., 2019), and
characterizing the evolution of the disease over time (Rubanova et al., 2020).

A number of NMF-based methods have been proposed for mutational signatures
analysis; see Islam et al. (2022) and references therein. However, a difficult
aspect of NMF is choosing an appropriate number of latent factors, which
corresponds to the number of mutational signatures present in the data.
Choosing too many factors can lead to the discovery of spurious signatures,
while choosing too few factors can lead to incorrectly merging distinct
signatures. Existing techniques for selecting the number of factors include
the Bayesian information criterion (Rosales et al., 2016; Fischer et al.,
2013), cross-validation (Lal et al., 2021), or even neural networks (Nebgen et
al., 2021; Islam et al., 2022). The disadvantage of such approaches is that
they require estimating a separate model for each choice of rank or
regularization parameter, making them computationally intensive. Another
popular approach is to use automatic relevance determination (Tan and Févotte,
2013; Kim et al., 2016), which is extremely fast but only provides point
estimates without any uncertainty quantification. Bayesian nonparametric
models with sparsity-inducing priors have been used to address this problem in
Gaussian factor models (Bhattacharya and Dunson, 2011; Legramanti et al.,
2020), however, these have not yet been applied to NMF, to our knowledge.

In this article, we introduce a novel Bayesian NMF method that yields accurate
and reliable inference in a computationally simple way by using _compressive
hyperpriors_ to drive the weights of unneeded factors to zero. Specifically,
we use a Poisson factorization model with semi-conjugate Dirichlet and gamma
priors over the signatures and the loadings, respectively, and we induce
sparsity by using a shrinkage hyperprior that strengthens with the amount of
data in such a way that (a) unneeded factors are given negligible weight with
probability tending to one, while (b) needed factors are given weights that
are only mildly shrunk. This compressive property enables the method to select
an appropriate number of latent factors in a continuous way, without the need
to fit multiple models or discretely jump between models in Markov chain Monte
Carlo samplers. Further, posterior inference can be carried out using a simple
auxiliary variable Gibbs sampling algorithm as in Dunson and Herring (2005)
and Cemgil (2009), making the method easy to implement. A crucial benefit of
our use of semi-conjugate continuous priors and hyperpriors is that they yield
a posterior surface that is easy to explore, in contrast with spike-and-slab
constructions that could lead to stronger multi-modalilty.

Additionally, the compressive hyperprior makes it straightforward to take
advantage of prior information about the latent factors. This is especially
useful for mutational signatures analysis, since the _Catalog Of Somatic
Mutations In Cancer_ (cosmic) database provides a curated set of mutational
signatures and their putative etiologies (Alexandrov et al., 2020). By using
an informative prior centered at the cosmic signatures, we find that the model
obtains improved sensitivity to detect the presence of signatures and infers
unambiguous matches to the original database.

The paper is organized as follows. Section 2 defines the model, introduces our
compressive hyperprior, and provides a Gibbs sampler algorithm for posterior
inference. In Section 3, we provide theoretical results establishing the
compressive property of the model. Section 4 provides background on mutational
signatures analysis. Section 5 contains a simulation study comparing to
leading methods, and Section 6 presents an application to mutational
signatures analysis using a benchmark breast cancer dataset. We conclude with
a brief discussion in Section 7.

##  2 Methodology

###  2.1 Poisson non-negative matrix factorization model

We describe the model in the context of mutational signatures, our primary
application of interest. Let $X_{ij}$ represent the number of mutations for
channel $i$ in sample $j$, where $i=1,\ldots,I$ and $j=1,\ldots,J$, and let
$X\in\mathds{R}^{I\times J}$ denote the matrix with entries $X_{ij}$.
Typically, the channels would consist of the 96 single-base substitution (SBS)
types; see Sections 4 and S4 for details. Non-negative matrix factorization
(NMF) consists of finding two non-negative matrices
$R\in\mathds{R}_{+}^{I\times K}$ and $\Theta\in\mathds{R}_{+}^{K\times J}$
such that $X\approx R\,\Theta$, with the _rank_ $K$ typically chosen so that
$K\leq\min\\{I,J\\}$. The $k$th column of $R$, denoted
$r_{k}=(r_{1k},\ldots,r_{Ik})$, is referred to as the $k$th _mutational
signature_. The $k$th row of $\Theta$, denoted
$\theta_{k}=(\theta_{k1},\ldots,\theta_{kJ})$, is the vector of weights
representing the _loading_ of signature $k$ in each of the $J$ samples,
sometimes referred as the _signature activity_.

From a probabilistic perspective, it is natural to model the mutation counts
as

$X_{ij}\sim\mathrm{Poisson}\bigg{(}\sum_{k=1}^{K}r_{ik}\theta_{kj}\bigg{)}$
(1)

independently for $i=1,\ldots,I$ and $j=1,\ldots,J$, where
$\mathrm{Poisson}(\lambda)$ denotes the Poisson distribution with mean
$\lambda$. In Section S4, we show that Equation 1 can be derived from first
principles by modeling the occurrences of mutations as continuous-time Markov
processes across the genome. Maximum likelihood estimation for $R$ and
$\Theta$ in Equation 1 is equivalent to minimizing the Kullback–Leibler
divergence in the original formulation of Lee and Seung (2000); see Zhou and
Carin (2015) for a discussion.

In mutational signatures analysis, it is common to impose the constraint that
$\sum_{i=1}^{I}r_{ik}=1$ for all $k=1,\ldots,K$. This avoids scaling
ambiguities in both the signature vectors $r_{k}=(r_{1k},\ldots,r_{Ik})$ and
their loadings $\theta_{k}=(\theta_{k1},\ldots,\theta_{kJ})$. Most methods do
not impose such a constraint during inference, opting to enforce it as a post-
processing step (Tan and Févotte, 2013; Drummond et al., 2023). However, we
find that building $\sum_{i=1}^{I}r_{ik}=1$ into our model has the additional
benefits of simplifying the inference algorithm and enabling direct use of
cosmic signatures for constructing informative priors.

###  2.2 Prior

For the prior distribution on the signatures $r_{k}$ and loadings
$\theta_{k}$, we take

$\displaystyle r_{k}=(r_{1k},\ldots,r_{Ik})$
$\displaystyle\sim\mathrm{Dirichlet}(\alpha,\ldots,\alpha),$ (2)
$\displaystyle\theta_{k1},\ldots,\theta_{kJ}\mid\mu_{k}\;$
$\displaystyle\sim\;\mathrm{Gamma}(a,\,a/\mu_{k}),$ (3) $\displaystyle\mu_{k}$
$\displaystyle\sim\pi(\mu_{k}),$ (4)

independently for $k=1,\ldots,K$, where $\alpha>0$, $a>0$, $\mu_{k}>0$, and
$\pi(\mu_{k})$ is a hyperprior. The Dirichlet prior in Equation 2
automatically enforces the constraint that $\sum_{i=1}^{I}r_{ik}=1$. Here,
$\mathrm{Gamma}(a,b)$ denotes the gamma distribution with mean $a/b$ and
variance $a/b^{2}$. Thus, the prior mean of the loadings is
$\mathds{E}(\theta_{kj}\mid\mu_{k})=\mu_{k}$, implying that $\mu_{k}$ controls
the overall contribution of signature $r_{k}$ to the factorization and, in
turn, to the total number of mutations generated by process $k$.

We refer to $\mu_{1},\ldots,\mu_{K}$ as _relevance weights_ , following the
usage of this type of prior structure in automatic relevance determination
(ARD) for shrinking the weights of unneeded factors to near-zero values (Tan
and Févotte, 2013). However, unlike Tan and Févotte (2013), we take a fully
Bayesian approach, quantifying uncertainty in $r_{k}$ and $\theta_{k}$ rather
than just optimizing them. Also, unlike typical uses of ARD in Bayesian neural
networks and Gaussian processes (Neal, 1996; Bishop, 2006), the marginal
likelihood in the Poisson NMF model does not have a closed-form expression and
thus is not amenable to direct optimization of $\mu_{k}$. Nonetheless, it
turns out that with a certain choice of data-dependent hyperprior on
$\mu_{k}$, we can obtain appealing computational properties similar to ARD,
while performing inference with simple Gibbs sampling updates; we discuss this
next.

###  2.3 Compressive hyperprior

The hyperprior on the relevance weights in Equation 4 plays a crucial role in
inferring the number of factors. Several approaches have been developed to
provide sparsity in Gaussian factorization models, such as spike-and-slab
priors that introduce exact zeros in the loadings (Carvalho et al., 2008;
Ročková and George, 2016) or induce cumulative shrinkage to near-zero values
for redundant factors (Legramanti et al., 2020; Frühwirth-Schnatter, 2023);
also see Liu et al. (2019) for an extension to NMF settings. However, in our
model, we found that spike-and-slab priors over $\mu_{k}$ tend to make
posterior inference difficult, likely due to the strong multimodal nature of
the resulting posterior.

Instead, we propose a simpler alternative based on continuous shrinkage
priors. Specifically, we let

$\mu_{k}\sim\mathrm{InvGamma}(aJ+1,\,\varepsilon aJ)$ (5)

independently for $k=1,\ldots,K$, where $\mathrm{InvGamma}(a_{0},b_{0})$
denotes the inverse-gamma distribution with mean $b_{0}/(a_{0}-1)$ when
$a_{0}>1$. Here, $a$ is the shape parameter from the prior on $\theta_{kj}$ in
Equation 3, and we set $\varepsilon>0$ to be a small constant, such as
$\varepsilon=0.001$. Note that this makes $\mu_{k}$ small _a priori_ , since
$\mathds{E}(\mu_{k})=\varepsilon$. Further, the full conditional mean is

$\mathds{E}(\mu_{k}\mid-)=\frac{1}{2}\,\varepsilon+\frac{1}{2}\,\bar{\theta}_{k},$
(6)

where $\bar{\theta}_{k}=\frac{1}{J}\sum_{j=1}^{J}\theta_{kj}$, implying that
the prior mean $\varepsilon$ and the average loading $\bar{\theta}_{k}$ for
signature $k$ have equal influence on the posterior for $\mu_{k}$. This
_strength-matching_ property of the hyperprior remains stable as the sample
size $J$ increases. It is important to note that this hyperprior should not be
interpreted as a meaningful representation of prior uncertainty regarding
$\mu_{k}$. Instead, it is designed to yield a posterior with good properties
in terms of computation and accuracy.

This deceptively simple choice of hyperprior has some key features. First, it
favors sparse solutions, since $\mathds{E}(\mu_{k})=\varepsilon$. This makes
it so that for any extra unneeded signatures, the hyperprior encourages the
corresponding relevance weights $\mu_{k}$ to be small, on the order of
$\varepsilon$. This leads to the compressive property of the model, which we
analyze in more detail in Section 3.2.

Second, despite its growing strength with $J$, this hyperprior does not overly
shrink the loadings $\theta_{kj}$ for factors that are needed to fit the data.
To see why, suppose we instead used a fixed-strength hyperprior, say,
$\mu_{k}\sim\mathrm{InvGamma}(a_{0},b_{0})$. Then instead of
$\varepsilon/2+\bar{\theta}_{k}/2$, the full conditional mean of $\mu_{k}$
would be $(b_{0}+aJ\bar{\theta}_{k})/(a_{0}+aJ)$, which is approximately
$\bar{\theta}_{k}$ when $J$ is sufficiently large. Consequently, when $J$ is
large and $\bar{\theta}_{k}\geq\varepsilon$, the prior mean of $\theta_{kj}$
differs by less than a multiple of 2 under the fixed-strength and strength-
matching hyperpriors since $\mathds{E}(\theta_{kj}\mid\mu_{k})=\mu_{k}$, and
this has only a small effect on the full conditional of $\theta_{kj}$ due to
the overwhelming contribution from the likelihood.

Another important feature of the strength-matching hyperprior in Equation 5 is
that small departures from the assumed Poisson NMF model do not strongly
affect the number of factors used by the model. As $J$ grows, a fixed-strength
hyperprior on $\mu_{k}$ would be overwhelmed by the likelihood since the
number of parameters $\theta_{kj}$ grows with $J$, leading to the inclusion of
spurious extra signatures when the model is slightly misspecified. Meanwhile,
the strength-matching drives out spurious extra signatures by maintaining a
balance between the contribution from the loadings and from the hyperprior, as
is evident from Equation 6.

Finally, the resulting posterior density is highly tractable due to the use of
semi-conjugate distributions. This improves the performance of the sampling
algorithm that we employ for inference, providing expeditious convergence to
NMF solutions along with Bayesian uncertainty quantification.

###  2.4 Posterior inference

Posterior inference for the hierarchical model defined by Equations 1, 2, 3, 4
and 5 can be efficiently performed via Gibbs sampling. Since the sum of
independent Poisson random variables is Poisson, we can equivalently write the
hierarchical model as

$\begin{split}X_{ij}&=\textstyle\sum_{k=1}^{K}Y_{ijk},\\\
Y_{ijk}\mid\mu_{k},r_{k},\theta_{k}\;&\sim\;\mathrm{Poisson}(r_{ik}\theta_{kj}),\\\
(r_{1k},\ldots,r_{Ik})\;&\sim\;\mathrm{Dirichlet}(\alpha,\ldots,\alpha),\\\
\theta_{k1},\ldots,\theta_{kJ}\mid\mu_{k}\;&\sim\;\mathrm{Gamma}(a,\,a/\mu_{k}),\\\
\mu_{k}\;&\sim\;\mathrm{InvGamma}(aJ+1,\,\varepsilon aJ).\end{split}$ (7)

Each auxiliary variable $Y_{ijk}$ can be interpreted as the number of
mutations due to signature $k$ in channel $i$ for sample $j$. Defining the
vector $Y_{ij}=(Y_{ij1},\ldots,Y_{ijK})$, it turns out that $Y_{ij}\mid
X_{ij},R,\Theta$ follows a
$\mathrm{Multinomial}\big{(}X_{ij},(q_{ij1},\ldots,q_{ijK})\big{)}$
distribution, where $q_{ijk}=r_{ik}\theta_{kj}/Q_{ij}$ and
$Q_{ij}=\sum_{k=1}^{K}r_{ik}\theta_{kj}$. This auxiliary variable
decomposition has been used in several previous methods (Dunson and Herring,
2005; Cemgil, 2009; Rosales et al., 2016; Zhou and Carin, 2015). The rest of
the sampler relies on standard semi-conjugate updates, which are
straightforward to derive.

##### Gibbs sampler algorithm

  1. 1.

For $i=1,\ldots,I$ and $j=1,\ldots,J$, update the latent mutation counts by
drawing

$(Y_{ij}\mid-)\sim\mathrm{Multinomial}\big{(}X_{ij},(q_{ij1},\ldots,q_{ijK})\big{)}$

where $q_{ijk}=r_{ik}\theta_{kj}/Q_{ij}$ and
$Q_{ij}=\sum_{k=1}^{K}r_{ik}\theta_{kj}$.

  2. 2.

For $k=1,\ldots,K$, update the signatures by drawing

$(r_{k}\mid-)\sim\mathrm{Dirichlet}\bigg{(}\alpha+\sum_{j=1}^{J}Y_{1jk},\,\ldots,\,\alpha+\sum_{j=1}^{J}Y_{Ijk}\bigg{)}.$

  3. 3.

For $k=1,\ldots,K$ and $j=1,\ldots,J$, update the loadings by drawing

$(\theta_{kj}\mid-)\sim\mathrm{Gamma}\bigg{(}a+\sum_{i=1}^{I}Y_{ijk},\;\frac{a}{\mu_{k}}+1\bigg{)}.$

  4. 4.

For $k=1\ldots,K$, update the relevance weights by drawing

$(\mu_{k}\mid-)\sim\mathrm{InvGamma}\bigg{(}2aJ+1,\;\varepsilon
aJ+a\sum_{j=1}^{J}\theta_{kj}\bigg{)}.$

The model is symmetric with respect to the order of the factors, in the sense
that the priors and likelihood are invariant to permutations of
$k=1,\ldots,K$. While attractive from a modeling standpoint, this symmetry
could potentially lead to label switching when running the Gibbs sampler,
complicating the calculation of posterior expectations. However, we have not
encountered label switching on either simulated or real data, so this has not
been an issue in practice.

###  2.5 Informative priors based on known signatures

A favorable aspect of mutational signatures analysis is the abundance of
historical data on signatures across many cancer types. The cosmic database
contains a curated collection of signatures, annotated with associated cancer
types and inferred etiologies (Alexandrov et al., 2020). It is natural to
leverage such prior information as follows.

Suppose $s_{k}=(s_{1k},\ldots,s_{Ik})$, for $k=1,\ldots,K_{\mathrm{pre}}$, are
pre-defined mutational signatures known to occur in cancer. To allow for
variation in signatures across studies, we let
$\rho_{k}=(\rho_{1k},\ldots,\rho_{Ik})$ denote a study-specific version of
$s_{k}$. We then generalize Equation 1 by modeling

$X_{ij}\sim\mathrm{Poisson}\bigg{(}\sum_{k=1}^{K_{\mathrm{pre}}}\rho_{ik}\omega_{kj}+\sum_{k=1}^{K_{\mathrm{new}}}r_{ik}\theta_{kj}\bigg{)}$
(8)

independently, where $r_{ik}$ and $\theta_{ik}$ are given the priors in
Equations 2 and 3, respectively, and

$\displaystyle\begin{split}\rho_{k}\;&\sim\;\mathrm{Dirichlet}(\beta_{k}s_{1k},\ldots,\beta_{k}s_{1k}),\\\
\omega_{kj}\mid\tau_{k}\;&\sim\;\mathrm{Gamma}(b,\,b/\tau_{k}),\\\
\tau_{k}\;&\sim\;\mathrm{InvGamma}(bJ+1,\,\varepsilon bJ).\end{split}$ (9)

Thus, the prior on $\rho_{k}$ is centered at $s_{k}$, with concentration
parameter $\beta_{k}$. The loadings $\omega_{kj}$ and corresponding relevance
weights $\tau_{k}$ are given the same prior and compressive hyperprior as
$\theta_{kj}$ and $\mu_{k}$, respectively, but with $b$ in place of $a$.

The model in Equation 8 is reminiscent of the recovery-discovery model
discussed in Grabski et al. (2023), when only a single study is taken into
consideration. Notice that in such a framework, the prior rank for the lower-
dimensional matrices, $K_{\mathrm{pre}}+K_{\mathrm{new}}$, is often greater
than $J$. This is at odds with the classic approach to NMF, where the
factorization rank is typically smaller than the rank of $X$. However, the
compressive mechanism behind our priors in Equation 9 still ensures a
parsimonious representation in the posterior, such that only the active
signatures have a nonnegligible relevance weight $\tau_{k}$. Posterior
inference can be performed using the same steps as in Section 2.4, with minor
adjustments for handling $\rho_{k}$, $\omega_{kj}$, and $\tau_{k}$; see
Section S3 for details.

##  3 Theory

In this section, we establish certain theoretical properties of the
hierarchical Poisson NMF model in Section 2. In particular, we study the
relationship between the latent mutation counts $Y_{ijk}$ and the relevance
weights $\mu_{k}$, and we prove a concentration result for $\mu_{k}$ given the
$Y_{ijk}$ values.

###  3.1 Conditional distribution of the relevance weights

We show that the distribution of $\mu_{k}\mid Y$ under the model in Equation
7, where $Y$ is the tensor $(Y_{ijk})\in\mathds{R}^{I\times J\times K}$, has a
closed-form expression in terms of confluent hypergeometric functions. We
refer to the resulting family of distributions, which appears to be novel, as
_inverse Kummer_.

######  Definition 1.

The _inverse Kummer_ distribution with parameters $\lambda>0$, $\beta>0$,
$\delta>0$, and $\gamma\in\mathds{R}$ is a continuous distribution on
$(0,\infty)$ with probability density function

$\displaystyle\pi(\mu)=\frac{\mu^{-(\lambda-\gamma)-1}(1+\mu/\delta)^{-\gamma}e^{-\beta/\mu}}{\delta^{\gamma-\lambda}\,\Gamma(\lambda)\,U(\lambda,\lambda+1-\gamma,\beta/\delta)}.$
(10)

We write $\mu\sim\mathrm{InvKummer}(\lambda,\beta,\gamma,\delta)$ to denote
that $\mu$ has the density in Equation 10.

Here, $U(a,b,z)$ denotes the confluent hypergeometric function of the second
kind,

$U(a,b,z)=\frac{1}{\Gamma(a)}\int_{0}^{\infty}t^{a-1}(1+t)^{b-a-1}e^{-zt}\mathrm{d}t,$

with $z>0$ (Abramowitz and Stegun, 1972). We call this an inverse Kummer
distribution since if $\mu\sim\mathrm{InvKummer}(\lambda,\beta,\gamma,\delta)$
then $1/\mu$ follows a Kummer distribution, which was introduced by Armero and
Bayarri (1997) when studying a $\mathrm{M/M/}\infty$ queuing problem. The
moments have closed-form expressions in terms of the hypergeometric function,
following Equation (6.7) in Armero and Bayarri (1997): for $m<\lambda$, the
$m$th moment of $\mu\sim\mathrm{InvKummer}(\lambda,\beta,\gamma,\delta)$ is

$\displaystyle\mathds{E}(\mu^{m})=\delta^{m}\frac{\Gamma(\lambda-m)}{\Gamma(\lambda)}\frac{U(\lambda-m,\,\lambda-m+1-\gamma,\,\beta/\delta)}{U(\lambda,\,\lambda+1-\gamma,\,\beta/\delta)}.$
(11)

The inverse Kummer is a generalization of the inverse gamma distribution,
since
$\mathrm{InvKummer}(\lambda,\beta,0,\delta)=\mathrm{InvGamma}(\lambda,\beta)$.
To elucidate the relationship between the inverse Kummer and inverse gamma, in
Section S2.1 we show that when $\lambda>2$, the mean $\mathds{E}(\mu)$ is
monotonically increasing as a function of $\gamma$.

In the compressive NMF model in Equation 7, the inverse Kummer arises as the
posterior distribution of the relevance weights given the latent counts,
integrating out the signatures and loadings.

######  Theorem 1.

Let $Y=(Y_{ijk})\in\mathds{R}^{I\times J\times K}$ denote the tensor of latent
counts. Under the hierarchical model in Equation 7, it holds that

$(\mu_{k}\mid Y)\sim\mathrm{InvKummer}\big{(}2aJ+1,\,\varepsilon
aJ,\,J\bar{Y}_{k}+aJ,\,a\big{)},$

where $\bar{Y}_{k}=\frac{1}{J}\sum_{i=1}^{I}\sum_{j=1}^{J}Y_{ijk}$ is the
average number of mutations assigned to signature $k$.

See Section S1 for the proof. Figure 1(A) shows the density of $\mu_{k}\mid Y$
for various $\bar{Y}_{k}$ and $J$ values.

Figure 1: (A) Density of $\mu_{k}\mid Y$ as the average number of mutations
$\bar{Y}_{k}$ and the sample size $J$ varies. Here, $a=1$ and
$\varepsilon=0.001$. (B) Mean of $\mu_{k}\mid Y$ as a function of
$\bar{Y}_{k}$, for varying values of $J$. Dashed lines indicate the 10th and
the 90th percentiles of $\mu_{k}\mid Y$. The black line shows where
$\mu_{k}=\bar{Y}_{k}$.

###  3.2 Compressive property

The essence of the compressive hyperprior is that for unneeded factors, the
relevance weights are shrunk to $\approx\varepsilon$, while for needed factors
they are only partly shrunk towards $\varepsilon$. Precisely this behavior is
seen in Figure 1(B), which shows the relationship between $\bar{Y}_{k}$ and
the mean of $\mu_{k}\mid Y$, based on the result in Theorem 1. The average
number of mutations due to the $k$th signature, $\bar{Y}_{k}$, is a minimal
sufficient statistic for $\mu_{k}$ and plays a fundamental role in determining
its posterior distribution. When $\bar{Y}_{k}\leq a$, we see that $\mu_{k}\mid
Y$ becomes concentrated near $\varepsilon$. Meanwhile, when $\bar{Y}_{k}\geq
a$, $\mathds{E}(\mu_{k}\mid Y)$ grows approximately linearly as a function of
$\bar{Y}_{k}$. Thus, given the latent counts $Y$, the model induces shrinkage
by a fixed factor when $\bar{Y}_{k}$ is large, while it has a sparsity
inducing effect when $\bar{Y}_{k}$ is small. Further, both Figure 1(A) and (B)
illustrate that $\mu_{k}\mid Y$ concentrates as $J$ increases.

In the following result, we make this rigorous by providing a concentration
result for the posterior distribution of $\mu_{k}$ in Theorem 1 as
$J\to\infty$.

######  Theorem 2.

Consider the hierarchical model in Equation 7. If $\bar{Y}_{k}\to
y\in[0,\infty)$ as $J\to\infty$, then for all $d>0$,

$\mathds{P}\big{(}|\mu_{k}-\mu_{*}|\leq
d\;\big{|}\;Y\big{)}\xrightarrow[J\to\infty]{}1,$

where
$\mu_{*}=2a\varepsilon/\big{(}\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon}-(y-a+\varepsilon)\big{)}$.

In Section S2.2, we show that $\mu_{*}$ is monotone increasing as a function
of $y$. Further, by a first-order Taylor approximation to the denominator of
$\mu_{*}$, when $\varepsilon\ll|y-a|$,

$\mu_{*}\approx\begin{cases}\displaystyle\frac{y-a}{2}&\text{ if $y>a$}\\\
\displaystyle\frac{\varepsilon a(a-y)}{(a-y)^{2}+(a+y)\varepsilon}&\text{ if
$0\leq y<a$}\end{cases}$ (12)

This explains the shape of $\mathds{E}(\mu_{k}\mid Y)$ seen in Figure 1; in
particular, $\mu_{*}\approx\varepsilon a/(a+\varepsilon)$ as $y\approx 0$. The
following corollary of Theorem 2 provides a criterion for thresholding the
relevance weights.

######  Corollary 1.

Under the assumptions of Theorem 2, if $\bar{Y}_{k}\to 0$ as $J\to\infty$,
then for all $C>1$,

$\mathds{P}(\mu_{k}>C\varepsilon\mid Y)\xrightarrow[J\to\infty]{}0.$

Hence, when signature $k$ is not being used by the model, Corollary 1 shows
that the posterior for $\mu_{k}$ concentrates on the interval
$(0,C\varepsilon)$, for any $C>1$. We refer to this as the _compressive
property_ of the model. This provides a natural criterion for selecting
signatures for inclusion in the model, by using a threshold of
$\mu_{k}>C\varepsilon$ to decide which signatures to keep and which to
discard. Inspection of the proofs of Theorems 1 and 2 shows that they hold for
any prior on signatures $r_{k}$ such that $\sum_{i}r_{ik}=1$. Consequently,
the concentration results in Theorem 2 and Corollary 1 also hold for the
relevance weights $\tau_{k}$ of the augmented model in Equation 9, which
employs informative priors on $r_{k}$.

This compressive property parallels the automatic selection of the number of
active components in overfitted mixture models described by Rousseau and
Mengersen (2011), which has also been used in factor analysis (e.g., Ferrari
and Dunson, 2021).

##  4 Background on mutational signatures analysis

Cancer development in humans is connected to the accumulation of mutations in
the DNA of somatic cells. When considering single-base substitutions,
mutations are classified according to which of the four nucleotide bases was
present before and after the mutation, on the strand containing the pyrimidine
before the mutation occurred. Recalling that adenine (A) and guanine (G) are
_purines_ while cytosine (C) and thymine (T) are _pyrimidines_ and that C
always binds with G and T with A, there are six possible types of
substitutions, namely, C$>$T, C$>$G, C$>$A, T$>$A, T$>$C and T$>$G. To account
for context-specific variability due to adjacent bases, mutations are further
classified according to which bases (A, G, C, or T) occur on the $5$’ and the
$3$’ sides on the strand containing the pre-substitution pyrimidine. This
makes for a total of $6\times 4\times 4=96$ types of single-base
substitutions, referred to as _mutational channels_ (Alexandrov et al., 2013);
see Section S4 for details.

It has been observed that many mutation-causing processes consistently produce
each type of mutation at a particular rate: for instance, ultraviolet
radiation has been observed to produce a large number of C$>$T substitutions
in melanoma and glioma (Greenman et al., 2007). Due to this, the mutational
processes acting on somatic cells can be characterized according to their
_mutational signatures_ (Nik-Zainal et al., 2012; Alexandrov et al., 2013),
where the term “signature” refers to a vector containing the probability of
each type of mutation under consideration.

A curated set of known signatures is maintained in the cosmic
database111https://cancer.sanger.ac.uk/signatures/ (Alexandrov et al., 2020),
which currently lists 86 single-base substitution (SBS) signatures. Many
signatures can be attributed to a specific etiology that has been
experimentally validated. For example, signatures SBS7a, b, c, and d are all
linked to ultraviolet light exposure. Other signatures, such as SBS60, require
further investigation to understand whether they arise from true biological
processes or are due to technical artifacts. See Koh et al. (2021) for an
overview.

The Poisson NMF model commonly used in mutational signatures analysis can be
derived from first principles by modeling the occurrence of nucleotide
substitutions at each base in the genome as a continuous-time Markov process
(Section S4). Aggregating the resulting substitution counts across the entire
genome, and modeling each mutational process as acting independently, it turns
out that the counts are well approximated by the Poisson NMF model in Equation
1. This provides a compelling biological justification for the use of this
model. We provide a detailed derivation in Section S4.

###  4.1 Existing methods

Several methods have been developed for mutational signature analysis; see the
review by Islam et al. (2022). Currently, the most prominent method is
SigProfiler (Alexandrov et al., 2020) and its successor, SigProfilerExtractor
(Islam et al., 2022). SigProfilerExtractor uses a bootstrap-like procedure to
resample the mutation count matrix from a Poisson model, and fits an NMF model
to each resampled data matrix by minimizing the Poisson Kullback–Leibler
divergence. The number of signatures is selected by applying the algorithm for
a range of $K$ values and using a neural network to choose $K$ (Nebgen et al.,
2021).

A leading Bayesian method is signeR (Rosales et al., 2016; Drummond et al.,
2023) which is based on the Poisson NMF model in Equation 1, but employs
independent hierarchical gamma priors over both the signatures and the
loadings. Selection of $K$ is performed using the Bayesian information
criteria (BIC) after running a separate model for each $K$. Hence, both
SigProfilerExtractor and signeR are particularly slow when the range of
possible $K$ values is moderate to large.

A faster alternative is offered by SignatureAnalyzer (Kim et al., 2016), which
fits the Poisson NMF model in Equation 1 using a _maximum a posteriori_
estimation algorithm (Tan and Févotte, 2013). Similarly to our approach,
SignatureAnalyzer uses ARD with inverse-gamma hyperpriors on the relevance
weights to determine the number of signatures. While the method is efficient
and flexible in terms of the choice of the objective function
(Kullback–Leibler or squared error) and prior (exponential or half-normal), it
does not provide any uncertainty quantification.

A natural Bayesian approach to selecting the number of latent factors is to
use spike-and-slab priors, where the relevance weight of each of the $K$
signatures has some probability of being sampled from a spike close to zero.
This is the approach taken in the elegant nonparametric factorization model
proposed by Legramanti et al. (2020), which infers the number of factors using
a cumulative shrinkage spike-and-slab process prior (cusp). We adapt cusp to
the Poisson NMF model in Equation 1 by defining the loadings to be
$\theta_{kj}=\mu_{k}\vartheta_{kj}$, where
$\vartheta_{kj}\sim\mathrm{Gamma}(a,a)$, and
$\mu_{k}\sim(1-\pi_{k})\mathrm{Gamma}(a_{0},b_{0})+\pi_{k}\delta_{\mu_{\infty}}$,
with $\delta_{x}$ denoting the point mass at $x$. The spike location
$\mu_{\infty}$ is typically chosen to be very small but nonzero, such as
$\mu_{\infty}=0.01$, while $\mathrm{Gamma}(a_{0},b_{0})$ represents the slab.
The spike probability $\pi_{k}$ is constructed via a stick-breaking process,
specifically, $\pi_{k}=\sum_{\ell=1}^{k}v_{\ell}\prod_{j=1}^{\ell-1}(1-v_{j})$
where
$v_{\ell}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mathrm{Beta}(1,\alpha)$
and $\alpha>0$. Hence, this model assumes an infinite number of factors _a
priori_ , but the probability that $\mu_{k}$ comes from the spike—effectively
removing that factor from the model—increases with $k$. Posterior inference is
performed via an adaptive Metropolis algorithm; see Section S5.1 for details.

It is worth emphasizing that, unlike our model in Equation 8, none of these
methods allow one to simultaneously use informative priors based on the cosmic
signatures and vague priors for _de novo_ analysis. Some of these methods
could potentially be adapted to Equation 8, but others cannot – for instance,
the cusp prior would result in an unnatural asymmetry across the known
signatures.

##  5 Simulations

In this section, we conduct a simulation study to evaluate the performance of
our compressive Poisson NMF method in terms of (a) detecting the true number
of signatures that are active in the data and (b) accurately recovering the
true signatures and their associated loadings.

###  5.1 Setup of the simulations

We simulate data and true parameters as follows. The mutation counts are
generated as
$X_{ij}\sim\mathrm{NegBin}\big{(}1/\tau,\,1/(1+\tau\lambda_{ij}^{0})\big{)}$,
where
$\lambda^{0}_{ij}=\sum_{k=1}^{K_{\mathrm{pre}}^{0}}\rho^{0}_{ik}\omega^{0}_{kj}+\sum_{k=1}^{K_{\mathrm{new}}^{0}}r^{0}_{ik}\theta_{kj}^{0}$.
We parametrize the negative binomial such that the mean and variance of
$X_{ij}$ are $\lambda_{ij}^{0}$ and
$\lambda_{ij}^{0}(1+\tau\lambda_{ij}^{0})$, respectively, where $\tau>0$ is a
parameter controlling overdispersion. We set $K_{\mathrm{pre}}^{0}=4$, and for
$k=1,\ldots,4$, we define $\rho_{k}^{0}=(\rho_{1k}^{0},\ldots,\rho_{Ik}^{0})$
to be cosmic signatures SBS1, SBS2, SBS3, and SBS13, respectively. SBS1 is a
sparse signature present in every cancer type, arising from the spontaneous
deamination of 5-methylcytosine. SBS3 is a rather flat signature due to
homologous recombination deficiencies in breast and pancreatic cancers. SBS2
and SBS13 are commonly occurring signatures associated with APOBEC activity.
Meanwhile, we randomly generate $r_{k}^{0}=(r_{1k}^{0},\ldots,r_{Ik}^{0})$ as
$r_{k}^{0}\sim\mathrm{Dirichlet}(0.25,\ldots,0.25)$, independently for
$k=1,\ldots,K_{\mathrm{new}}^{0}$. We generate loadings by setting
$\omega^{0}_{kj}=w_{k}\xi_{kj}$ where $w_{k}\sim\mathrm{Gamma}(100,1)$ and
$\xi_{kj}\sim\mathrm{Gamma}(0.5,0.5)$ independently, and $\theta^{0}_{kj}$ in
the same way as $\omega^{0}_{kj}$.

For the overdispersion $\tau$, we consider two settings: $\tau=0$, in which
case the negative binomial reduces to a Poisson (and thus, the assumed model
is correct), and $\tau=0.15$, resulting in mild misspecification. For the
number of non-cosmic signatures $K_{\mathrm{new}}^{0}$, we consider
$K_{\mathrm{new}}^{0}=2$ and $K_{\mathrm{new}}^{0}=6$, so that the total
number $K^{0}=K_{\mathrm{pre}}^{0}+K_{\mathrm{new}}^{0}$ is either $6$ or
$10$. We consider a range of sample sizes $J\in\\{50,100,200\\}$. For each
combination of $\tau$, $J$, and $K^{0}$, we generate $20$ replicate sets of
parameters and data matrices as described above. On each simulated data
matrix, we run six methods:

  1. (i)

CompNMF: our compressive NMF model in Equation 7 with $K=20$ signatures and
$\varepsilon=0.001$,

  2. (ii)

CompNMF+cosmic: our enhanced model in Equation 8 with $K_{\mathrm{new}}=15$
_de novo_ signatures and the $K_{\mathrm{pre}}=67$ cosmic v3.4 signatures that
are not regarded as “possible sequencing artifacts”, and $\varepsilon=0.001$,

  3. (iii)

PoissonCUSP: the cusp model with spike location $\mu_{\infty}=0.01$, slab
parameters $a_{0}=b_{0}=1$, starting at $K=20$ and using Algorithm 2 of
Legramanti et al. (2020) to adaptively tune $K$,

  4. (iv)

signeR: the signeR model (Rosales et al., 2016; Drummond et al., 2023) with
default parameters (estimate_hyper = FALSE) and with $K$ ranging from $2$ to
$20$,

  5. (v)

SignatureAnalyzer: as implemented in the sig_auto_extract function of the
sigminer package (Wang et al., 2020), with selection method set to L1KL and
$K=20$,

  6. (vi)

SigProfiler: SigProfilerExtractor v1.1.23 (Islam et al., 2022) with random
initialization and 5 replicates for each $K\in\\{2,\ldots,20\\}$, using the
sigprofiler_extract wrapper in sigminer.

For methods (i), (ii), and (iii), we set $a=1$ and $\alpha=0.5$, and run the
sampler for $5000$ iterations, discarding the first $4000$ as burn-in. In
method (iii), adaptation of $K$ was started after 500 iterations. In method
(ii), we set the parameter $\beta_{k}$ such that under the Dirichlet prior in
Equation 9, the median cosine similarity between the prior mean and a sample
from the prior is approximately 0.975. This makes $\beta_{k}$ depend on the
sparsity of the signature $s_{k}$. For example, we set $\beta_{k}=17.29$ for
the sparse signature SBS2, while $\beta_{k}=1337.26$ for the rather flat SBS3.

Figure 2: (A) Estimated number of signatures by each model, for varying
choices of $J$, overdispersion, and true number of signatures, across 20
replicated data set in each scenario. The horizontal grey line indicates the
true number ($K^{0}=6$ or $K^{0}=10$). (B) Average precision and sensitivity
across 20 replicate data sets in each scenario, with a $0.9$ cutoff for the
cosine similarity. The dashed contour lines in the background indicate the
$F_{1}$ score. (C) Average $F_{1}$ score as a function of the cosine
similarity cutoff, across 20 replicates in each scenario, when the
overdispersion is set to 0.15.

###  5.2 Simulation results

Figure 2 shows the main results of the simulation study. Figure 2(A) reports
the estimated number of signatures for each method, for each combination of
$\tau$, $J$, and $K^{0}$; the boxplots summarize the distribution of estimated
values across the 20 replicates. The estimated number of signatures $K^{*}$ is
defined as follows for each method: For (i) and (ii), $K^{*}$ is the number of
signatures for which the posterior mean of $\mu_{k}$ is greater than
$5\varepsilon=0.005$; for (iii), $K^{*}$ is the number of signatures for which
the posterior probability of being assigned to the spike is less than $0.05$;
and for (iv), (v), and (vi), $K^{*}$ is the suggested solution returned by the
corresponding package.

As expected, all methods accurately estimate $K^{0}$ when $\tau=0$ (zero
overdispersion), since the Poisson NMF model is correct in this case.
Meanwhile, when $\tau=0.15$ (mild overdispersion), we observe a noticeably
different pattern of results. Our CompNMF+cosmic method, which leverages
informative priors based on all $67$ cosmic signatures, correctly recovers
$K^{0}$ more often than all other methods, both when $K^{0}=6$ and $K^{0}=10$.
CompNMF, which does not rely on cosmic, tends to slightly overestimate
$K^{0}$. This is likely due to the introduction of spurious signatures used by
the model to accommodate the overdispersion. Interestingly, signeR works well
even when the model is mildly misspecified, due to the tendency of the signeR
to return a conservative estimate of the number of signatures, as discussed in
Drummond et al. (2023). SigProfiler is even more conservative than signeR,
frequently underestimating $K^{0}$ slightly in these simulations. In contrast,
SignatureAnalyzer and PoissonCUSP strongly overestimate $K^{0}$ when there is
overdispersion.

Figure 2(B) displays the precision and sensitivity, averaged over the 20
replicates, for each model in each setting. Here, we define _precision_ as the
proportion of estimated signatures that have a cosine similarity $\geq 0.9$
with at least one of the ground truth signatures, and _sensitivity_ is the
proportion of ground truth signatures for which there is an estimated
signature with cosine similarity $\geq 0.9$. Since this is a simulation study,
the ground truth signatures $\rho_{k}^{0}$ and $r_{k}^{0}$ are known. The
cutoff value of $0.9$ was chosen following Islam et al. (2022), but we also
vary it as described below. The contour lines in Figure 2(B) represent the
$F_{1}$ score, defined as
$2\times\mathrm{precision}\times\mathrm{sensitivity}/(\mathrm{precision}+\mathrm{sensitivity})$.

Points in the top-right corner of the plot in Figure 2(B) indicate better
performance. When the Poisson model is correct (no overdispersion), all of the
methods perform well, except that SigProfiler has lower sensitivity for
smaller sample sizes. CompNMF+cosmic consistently performs the best, with
precision and sensitivity close to one in all of these settings. CompNMF and
signeR also perform well, but with somewhat reduced precision and sensitivity
when $K^{0}$ is larger and when there is overdispersion. SigProfiler has good
precision but much lower sensitivity, particularly in the presence of
overdispersion. Meanwhile, SignatureAnalyzer and PoissonCUSP struggle in the
overdispersed scenarios, exhibiting severely degraded precision as well as low
sensitivity.

To see the range of performance exhibited by each method as the cosine
similarity cutoff varies, Figure 2(C) shows the average $F_{1}$ score versus
the cutoff value. As before, we see that CompNMF+cosmic is a clear standout,
generally exhibiting the best performance overall, and particularly excelling
at cutoff values between 0.9 and 0.98. CompNMF, signeR, and SigProfiler are
roughly comparable by this metric, whereas SignatureAnalyzer and especially
PoissonCUSP suffer from significantly lower $F_{1}$ scores across all cutoff
values.

Figure 3: Average computation time (in minutes) for each method. Times shown
are the average across all combinations of overdispersion $\tau$ and number of
signatures $K^{0}$ used in the simulation study.

Figure 3 shows the total computation time required by each method, as a
function of the number of samples $J$. All computations were performed on an
AMD Ryzen 3900-based dedicated server with 128GB of memory on Ubuntu 20.04, R
version 4.3.1 linked to Intel MKL 2019.5-075. Calculations were split across
20 cores via the foreach package, allocating one dataset per core for each
combination of $\tau$, $J$, and $K^{0}$, and running each method sequentially.
Not surprisingly, SignatureAnalyzer is the fastest method by an order of
magnitude, always taking under one minute to complete. PoissonCUSP is the
second fastest method, because the number of signatures in the model varies
adaptively within the Gibbs sampler, preventing unnecessary computations on
inactive factors. CompNMF is slightly slower than PoissonCUSP since no
adaptation in performed, but it is faster than signeR. SigProfiler was the
slowest among all methods by a wide margin.

Thus, in terms of computation time, both CompNMF and CompNMF+cosmic are
competitive with the next best-performing existing method, signeR.
Furthermore, signeR is written in C++, whereas CompNMF and CompNMF+cosmic are
written in R code, so they could be significantly faster if implemented in C++
or similar. Additionally, they could be sped up by updating only the non-
compressed factors when running the Gibbs sampler, following the adaptive
approach of PoissonCUSP.

The dominant performance of CompNMF+cosmic demonstrates the benefits of using
our informative prior. Notably, since CompNMF+cosmic includes all
$K_{\mathrm{pre}}=67$ cosmic signatures—whereas only $K_{\mathrm{pre}}^{0}=4$
cosmic signatures are used to generate the simulated data—this implies that
$63$ of these signatures are correctly compressed out of the model,
illustrating the effectiveness of our compressive hyperprior. signeR performs
impressively well in these simulations, tending to fall between CompNMF and
CompNMF+cosmic in nearly all metrics. However, as mentioned above, the CompNMF
methods would be significantly faster than signeR if implemented in the same
language.

The poor performance of SignatureAnalyzer and PoissonCUSP in the misspecified
settings appears to be due to overfitting. Indeed, the overestimation of
$K^{0}$ seen in Figure 2(A), along with the fact that SignatureAnalyzer and
PoissonCUSP yield the lowest root-mean-square error for the count matrices
(see Section S5), is a clear sign of overfitting. In the case of PoissonCUSP,
this is not surprising since the flexible nonparametric nature of the model
allows it to closely fit the data, although unfortunately in this case it is
fitting noise rather than signal.

##  6 Application: 21 breast cancer data

In this section, we apply our compressive NMF method to the benchmark 21
breast cancer dataset considered by Nik-Zainal et al. (2012) in their landmark
paper originating the study of mutational signatures in cancer. In particular,
we aim to assess the effect of our compressive hyperprior for small sample
sizes and evaluate whether our informative prior provides an advantage in
practice.

The dataset is based on whole-genome sequencing of $J=21$ patients, and
consists of mutation counts for each of the $I=96$ channels for each patient.
We obtained the data from the signeR package (Drummond et al., 2023). The
total mutation count is fairly homogeneous across patients, except for patient
PD4120a, for whom a large number of mutations was detected ($70,690$) compared
to the others ($18,871$ on average). This is attributable to the fact that
PD4120a was sequenced more deeply at nearly 200x coverage compared to around
30x for the others.

We evaluate the performance of our compressive NMF method (with and without
the informative cosmic-based prior) compared to the same alternatives as in
Section 5. Here, for both CompNMF and CompNMF+cosmic, we randomly initialize
by sampling from the prior, and run the Gibbs sampler for $12{,}000$
iterations, discarding the first $10{,}000$ as burn-in. This is repeated four
times, and we select the run yielding the highest average log-posterior.
Details of the settings used in all of the other methods are reported in
Section S6.

###  6.1 Application results

All of the methods fit the count matrix roughly equally well, with the
exception of SigProfiler, which exhibited significantly higher error.
Specifically, the root-mean-squared error (RMSE) between the mutation count
matrix $X$ and $\hat{R}\hat{\Theta}$ was 9.51 for CompNMF, 9.57 for
CompNMF+cosmic, 9.81 for signeR, 10.07 for PoissonCUSP, 10.08 for
SignatureAnalyzer, and 37.08 for SigProfiler.

Figure 4: (A) Cosine similarity between the signatures inferred by each method
and the best matching signature in cosmic. Vertical error bars indicate 90%
posterior credible intervals. (B) Signatures inferred by CompNMF+cosmic. (C)
Signatures inferred by CompNMF. The vertical grey lines in panels (B) and (C)
indicate univariate 90% posterior credible intervals for each mutational
channel. The numbers on the top-right corners of each panel denote the cosine
similarity with the best matching cosmic signature.

However, major differences were observed in the sets of estimated signatures.
For each method, Figure 4(A) shows the cosine similarity between each
estimated signature and the best matching cosmic signature in terms of cosine
similarity; see Section S6 for all estimated signatures. Overall,
CompNMF+cosmic recovers the most signatures with the highest cosine
similarities, which is not surprising since the cosmic signatures are used to
construct the informative prior used in CompNMF+cosmic. Nonetheless, this
compellingly demonstrates that the informative prior enables more signal to be
extracted from the data. For instance, CompNMF+cosmic finds evidence of the
presence of SBS8 and SBS40a, which are not recovered by any other method
except for SigProfiler, which finds a signature with moderate cosine
similarity to SBS40a. SBS8 appears in other breast cancers and may be
associated with homologous recombination deficiency. SBS40a is a rather flat
signature that has unknown etiology but appears in every cancer type
(Alexandrov et al., 2020).

CompNMF recovers all but two of the signatures found by CompNMF+cosmic,
although with somewhat lower cosine similarities. PoissonCUSP yields similar
results to CompNMF, but includes what appears to be a spurious match to SBS9,
since it has low cosine similarity and SBS9 is not found by any other method.
Next, signeR is also similar to CompNMF, but does not recover SBS34, which
does appear to be truly present since CompNMF+cosmic detects it in many
samples. SignatureAnalyzer recovers four of the same signatures as CompNMF,
with comparable cosine similarities; however, it misses SBS34 and produces two
signatures that are both matched to SBS2, one of which is an amalgamation of
SBS2 and SBS13. Lastly, SigProfiler recovers only three signatures. The strong
performance of the CompNMF methods demonstrates the utility of our compressive
hyperprior.

Figures 4(B) and (C) show the estimated signature vectors $r_{k}$ for our two
methods (CompNMF+cosmic and CompNMF); corresponding figures for all other
methods are in Section S6. The informative prior enables CompNMF+cosmic to
obtain clean signature estimates that are unambiguously matched to known
cosmic signatures, while still allowing for dataset-specific departures from
the cosmic signatures. In contrast, due to the small size of this dataset, all
of the other methods produce “merged” signatures that are combinations of two
or more cosmic signatures. For instance, CompNMF+cosmic perfectly
distinguishes SBS2 and SBS13, whereas these are merged into a single signature
by all other methods; compare SBS2 and SBS13 in Figure 4(B) to CompNMF SigC in
Figure 4(C) and the other methods’ results in Section S6. Both SBS2 and SBS13
are associated with APOBEC activity and appear frequently in breast cancers.
Similarly, CompNMF+cosmic perfectly separates SBS1 and SBS40a, whereas these
two are combined by all the other methods; see SigA in Figure 4(C) and see
Section S6.

Figure 5: Average posterior loadings (%) by patient. Left: CompNMF+cosmic.
Right: CompNMF.

The vertical bars in Figure 4(A,B,C) indicate model-based uncertainty in the
form of 90% credible intervals; there are no bars for SigProfiler and
SignatureAnalyzer since they only provide point estimates. This uncertainty
quantification is particularly useful for identifying signatures that may be
spurious. For instance, for both CompNMF+cosmic and CompNMF, the estimated
signature matched to SBS98 has high uncertainty in the cosine similarity and
in the signature vector itself, indicating that one should be skeptical about
the validity of the estimated signature and whether it should be matched to
SBS98.

Figure 5 shows the loadings $\theta_{kj}$ estimated by CompNMF+cosmic and
CompNMF; see Section S6 for all other methods. The CompNMF+cosmic loadings
provide insight into the interpretation of the other methods’ results. For
instance, in Figure 5, the loadings for SBS40a appear to have been split up
and added to the loadings for SBS1 and SBS3 to make the CompNMF loadings for
SigA and SigD, respectively. Indeed, looking at these signatures, it appears
that SBS40a and SBS1 (Figure 4(B)) were agglomerated to make SigA (Figure
4(C)), and that SBS40a and SBS3 were agglomerated to make SigD. Similarly,
CompNMF SigC appears to be a combination of SBS2 and SBS13. This is especially
apparent in patient PD4120a, for whom almost all mutations are attributed to
SigC by CompNMF, whereas CompNMF+cosmic splits the loading into equal
contributions from SBS2 and SBS13.

##  7 Discussion

This article introduces a novel Bayesian NMF method that obtains state-of-the-
art performance for mutational signatures analysis, with lower computational
burden. In particular, our compressive hyperprior provides a simple but
effective technique for determining the subset of active factors. This enables
the use of rich informative priors based on the cosmic database of known
signatures, significantly boosting the method’s precision and sensitivity for
recovering true signatures. The informative prior also disambiguates the
allocation of loadings to signatures, resulting in more accurate estimation of
the contribution of each signature to each sample, as well as clarifying the
results of other methods. Furthermore, the method provides posterior
uncertainty quantification, which helps distinguish real from spurious
signatures and can be used for downstream analyses.

There are several interesting directions for future work. First, one can
envision several extensions of the model. Following Grabski et al. (2023), it
would be interesting to jointly model multiple studies or multiple cancer
types using a hierarchical model with study-specific or cancer type-specific
parameters. Another useful extension of the model would be to include sample-
specific covariates, which could be helpful in improving targeted therapies
(Aguirre et al., 2018).

Additionally, the scope of applicability of the compressive hyperprior
technique is potentially broader than Poisson NMF models, and might prove
useful in other latent factorization models such as Gaussian factor models
(Bhattacharya and Dunson, 2011; Legramanti et al., 2020) or in user-item
recommendation (Gopalan et al., 2014), especially when prior information on
the factors is available.

Finally, while our method can handle mild misspecification in the form of
small overdispersion, it is fundamentally based on the assumption that the
counts are Poisson distributed – like all of the leading methods considered in
our empirical results (signeR, SigProfiler, SignatureAnalyzer, PoissonCUSP).
Consequently, larger departures from the assumed Poisson NMF model can be
expected to negatively impact the performance of all of these methods. Of
course, overdispersion can simply be handled by modeling the data as negative
binomial rather than Poisson (Lyu et al., 2020). However, there are many other
plausible ways in which the model could be misspecified, and any parametric
elaboration of the model will inevitably be misspecified in some way. Thus, an
important area for future work is providing improved robustness to
misspecification for mutational signatures analysis as well as non-negative
matrix factorization more generally.

## Acknowledgments

We would like to thank Scott Carter for suggesting this research direction. We
also thank Jonathan Huggins, Catherine Xue, and Giovanni Parmigiani for many
helpful discussions. J.W.M. was supported by the Department of Defense grant
W81XWH-18-1-0357 and the National Institutes of Health grant 5R01CA240299.

## Code availability

Code is publicly accessible via GitHub at
https://github.com/alessandrozito/CompressiveNMF.

## Supplementary material

The supplementary material contains proofs of the theoretical results, the
rationale for the Poisson NMF model likelihood, simulation details, and
additional empirical results.

## References

  * Abramowitz and Stegun (1972) Abramowitz, M. and I. A. Stegun (1972).  Handbook of Mathematical Functions with formulas, graphs, and mathematical tables.  U.S. Dept. of Commerce, National Bureau of Standards.
  * Aguirre et al. (2018) Aguirre, A. J., J. A. Nowak, N. D. Camarda, et al. (2018).  Real-time genomic characterization of advanced pancreatic cancer to enable precision medicine.  Cancer Discovery 8(9), 1096–1111.
  * Alexandrov et al. (2020) Alexandrov, L. B., J. Kim, N. J. Haradhvala, et al. (2020).  The repertoire of mutational signatures in human cancer.  Nature 578(7793), 94–101.
  * Alexandrov et al. (2013) Alexandrov, L. B., S. Nik-Zainal, D. C. Wedge, et al. (2013).  Signatures of mutational processes in human cancer.  Nature 500(7463), 415–421.
  * Armero and Bayarri (1997) Armero, C. and M. Bayarri (1997).  A Bayesian analysis of a queueing system with unlimited service.  Journal of Statistical Planning and Inference 58(2), 241–261.
  * Ayed and Caron (2021) Ayed, F. and F. Caron (2021).  Nonnegative Bayesian nonparametric factor models with completely random measures.  Statistics and Computing 31(5), 63.
  * Bhattacharya and Dunson (2011) Bhattacharya, A. and D. B. Dunson (2011).  Sparse Bayesian infinite factor models.  Biometrika 98(2), 291–306.
  * Bishop (2006) Bishop, C. M. (2006).  Pattern Recognition and Machine Learning.  Springer.
  * Carvalho et al. (2008) Carvalho, C. M., J. Chang, J. E. Lucas, J. R. Nevins, Q. Wang, and M. West (2008).  High-dimensional sparse factor modeling: Applications in gene expression genomics.  Journal of the American Statistical Association 103(484), 1438–1456.
  * Cemgil (2009) Cemgil, A. T. (2009).  Bayesian inference for nonnegative matrix factorisation models.  Computational Intelligence and Neuroscience 2009, 785152.
  * Drummond et al. (2023) Drummond, R. D., A. Defelicibus, M. Meyenberg, R. Valieris, E. Dias-Neto, R. A. Rosales, and I. T. da Silva (2023).  Relating mutational signature exposures to clinical data in cancers via signeR 2.0.  BMC Bioinformatics 24(1), 439.
  * Dunson and Herring (2005) Dunson, D. B. and A. H. Herring (2005).  Bayesian latent variable models for mixed discrete outcomes.  Biostatistics 6(1), 11–25.
  * Ferrari and Dunson (2021) Ferrari, F. and D. B. Dunson (2021).  Bayesian factor analysis for inference on interactions.  Journal of the American Statistical Association 116(535), 1521–1532.
  * Fischer et al. (2013) Fischer, A., C. J. Illingworth, P. J. Campbell, and V. Mustonen (2013).  EMu: probabilistic inference of mutational processes and their localization in the cancer genome.  Genome Biology 14(4), R39.
  * Frühwirth-Schnatter (2023) Frühwirth-Schnatter, S. (2023).  Generalized cumulative shrinkage process priors with applications to sparse Bayesian factor analysis.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 381(2247), 20220148.
  * Gopalan et al. (2014) Gopalan, P., F. J. Ruiz, R. Ranganath, and D. Blei (2014).  Bayesian nonparametric Poisson factorization for recommendation systems.  In S. Kaski and J. Corander (Eds.), Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, Volume 33 of Proceedings of Machine Learning Research, pp. 275–283.
  * Grabski et al. (2023) Grabski, I., L. Trippa, and G. Parmigiani (2023).  Bayesian multi-study non-negative matrix factorization for mutational signatures.  bioRxiv:10.1101/2023.03.28.534619, 1–12.
  * Greenman et al. (2007) Greenman, C., P. Stephens, R. Smith, et al. (2007).  Patterns of somatic mutation in human cancer genomes.  Nature 446(7132), 153–158.
  * Gulhan et al. (2019) Gulhan, D. C., J. J.-K. Lee, G. E. M. Melloni, I. Cortés-Ciriano, and P. J. Park (2019).  Detecting the mutational signature of homologous recombination deficiency in clinical samples.  Nature Genetics 51(5), 912–919.
  * Hoffman et al. (2010) Hoffman, M. D., D. M. Blei, and P. R. Cook (2010).  Bayesian nonparametric matrix factorization for recorded music.  In Proceedings of the 27th International Conference on International Conference on Machine Learning, pp. 439–446.
  * Islam et al. (2022) Islam, S. A., M. Díaz-Gay, Y. Wu, et al. (2022).  Uncovering novel mutational signatures by de novo extraction with SigProfilerExtractor.  Cell Genomics 2(11), 100179.
  * Kim et al. (2016) Kim, J., K. W. Mouw, P. Polak, et al. (2016).  Somatic ERCC2 mutations are associated with a distinct genomic signature in urothelial tumors.  Nature Genetics 48(6), 600–606.
  * Koh et al. (2021) Koh, G., A. Degasperi, X. Zou, S. Momen, and S. Nik-Zainal (2021).  Mutational signatures: emerging concepts, caveats and clinical applications.  Nature Reviews Cancer 21(10), 619–637.
  * Lal et al. (2021) Lal, A., K. Liu, R. Tibshirani, A. Sidow, and D. Ramazzotti (2021).  De novo mutational signature discovery in tumor genomes using SparseSignatures.  PLOS Computational Biology 17(6), 1–24.
  * Lee and Seung (2000) Lee, D. and H. S. Seung (2000).  Algorithms for non-negative matrix factorization.  In T. Leen, T. Dietterich, and V. Tresp (Eds.), Advances in Neural Information Processing Systems, Volume 13. MIT Press.
  * Legramanti et al. (2020) Legramanti, S., D. Durante, and D. B. Dunson (2020).  Bayesian cumulative shrinkage for infinite factorizations.  Biometrika 107(3), 745–752.
  * Liu et al. (2019) Liu, Y., W. Dong, W. Song, and L. Zhang (2019).  Bayesian nonnegative matrix factorization with a truncated spike-and-slab prior.  In 2019 IEEE International Conference on Multimedia and Expo (ICME), pp. 1450–1455.
  * Lyu et al. (2020) Lyu, X., J. Garret, G. Rätsch, and K.-V. Lehmann (2020).  Mutational signature learning with supervised negative binomial non-negative matrix factorization.  Bioinformatics 36(Issue Supplement_1), i154–i160.
  * Neal (1996) Neal, R. M. (1996).  Bayesian Learning for Neural Networks.  Number 118. Springer. Lecture Notes in Statistics.
  * Nebgen et al. (2021) Nebgen, B. T., R. Vangara, M. A. Hombrados-Herrera, S. Kuksova, and B. S. Alexandrov (2021).  A neural network for determination of latent dimensionality in non-negative matrix factorization.  Machine Learning: Science and Technology 2(2), 025012.
  * Nik-Zainal et al. (2012) Nik-Zainal, S., L. B. Alexandrov, D. C. Wedge, P. Van Loo, C. D. Greenman, K. Raine, D. Jones, J. Hinton, J. Marshall, L. A. Stebbings, A. Menzies, S. Martin, K. Leung, L. . Chen, and M. R. Stratton (2012).  Mutational processes molding the genomes of 21 breast cancers.  Cell 149(5), 979–993.
  * Nik-Zainal et al. (2016) Nik-Zainal, S., H. Davies, J. Staaf, et al. (2016).  Landscape of somatic mutations in 560 breast cancer whole-genome sequences.  Nature 534(7605), 47–54.
  * Ročková and George (2016) Ročková, V. and E. I. George (2016).  Fast Bayesian factor analysis via automatic rotations to sparsity.  Journal of the American Statistical Association 111(516), 1608–1622.
  * Rosales et al. (2016) Rosales, R. A., R. D. Drummond, R. Valieris, E. Dias-Neto, and I. T. da Silva (2016).  signeR: an empirical Bayesian approach to mutational signature discovery.  Bioinformatics 33(1), 8–16.
  * Rousseau and Mengersen (2011) Rousseau, J. and K. Mengersen (2011).  Asymptotic behaviour of the posterior distribution in overfitted mixture models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73(5), 689–710.
  * Rubanova et al. (2020) Rubanova, Y., R. Shi, C. F. Harrigan, R. Li, J. Wintersinger, N. Sahin, A. G. Deshwar, et al. (2020).  Reconstructing evolutionary trajectories of mutation signature activities in cancer using TrackSig.  Nature Communications 11(1), 731.
  * Tan and Févotte (2013) Tan, V. Y. and C. Févotte (2013).  Automatic relevance determination in nonnegative matrix factorization with the beta-divergence.  IEEE Transactions on Pattern Analysis and Machine Intelligence 35(7), 1592–1605.
  * Wang et al. (2020) Wang, S., Z. Tao, T. Wu, and X.-S. Liu (2020).  Sigflow: an automated and comprehensive pipeline for cancer genome mutational signature analysis.  Bioinformatics 37(11), 1590–1592.
  * Wang and Zhang (2013) Wang, Y.-X. and Y.-J. Zhang (2013).  Nonnegative matrix factorization: A comprehensive review.  IEEE Transactions on Knowledge and Data Engineering 25(6), 1336–1353.
  * Zhou (2018) Zhou, M. (2018).  Nonparametric Bayesian negative binomial factor analysis.  Bayesian Analysis 13(4), 1065–1093.
  * Zhou and Carin (2015) Zhou, M. and L. Carin (2015).  Negative binomial process count and mixture modeling.  IEEE Transactions on Pattern Analysis and Machine Intelligence 37(2), 307–320.

Supplementary Material for
“Compressive Bayesian non-negative matrix factorization for mutational
signatures analysis”

Section S1 contains the proofs of the theoretical results. Section S2 provides
additional results related to the inverse Kummer distribution and the marginal
distribution of the latent counts. Section S3 gives a detailed description of
the sampler for the NMF with informative priors. Section S4 provides
additional background on mutational signatures, including a first-principles
derivation of the model. Section S5 contains additional details of the
simulations in Section 5, including an adaptation of the CUSP model of
Legramanti_2020 to the setting of Poisson non-negative matrix factorization.
Finally, Section S6 provides additional results on the application.

##  S1 Proofs

#### Proof of Theorem 1

###### Proof.

The proof proceeds by directly marginalizing $r_{k}$ and $\theta_{k}$ from the
joint distribution. For notational simplicity, we first handle the general
case of $\mu_{k}\sim\mathrm{InvGamma}(a_{0},b_{0})$, and then plug in the
values of $a_{0}$ and $b_{0}$ for the compressive hyperprior. Under the model
in Equation 7, the joint density of $Y=(Y_{ijk})$, $\Theta=(\theta_{kj})$,
$R=(r_{ik})$, and $\mu=(\mu_{k})$ is

$\displaystyle\pi(Y,\Theta,R,\mu)\propto$
$\displaystyle\bigg{\\{}\prod_{i,j,k}e^{-r_{ik}\theta_{kj}}\frac{(r_{ik}\theta_{kj})^{Y_{ijk}}}{Y_{ijk}!}\bigg{\\}}\times$
(Latent Poisson counts) (S1)
$\displaystyle\bigg{\\{}\prod_{i,k}r_{ik}^{\alpha-1}\bigg{\\}}\times$
(Dirichlet prior on $r_{k}$)
$\displaystyle\bigg{\\{}\prod_{j,k}\Big{(}\frac{a}{\mu_{k}}\Big{)}^{a}\theta_{kj}^{a-1}e^{-a\theta_{kj}/\mu_{k}}\bigg{\\}}\times$
(Gamma prior on $\theta_{k}$)
$\displaystyle\bigg{\\{}\prod_{k}\mu_{k}^{-a_{0}-1}e^{-b_{0}/\mu_{k}}\bigg{\\}},$
(InvGamma prior on $\mu_{k}$)

dropping constants of proportionality. Since $\sum_{i}r_{ik}=1$ for all $k$,
we have $\prod_{i}e^{-r_{ik}\theta_{kj}}=e^{-\theta_{kj}}$. Thus, by Equation
S1, we have

$\displaystyle\pi(Y,\Theta,R,\mu)=f(Y,R)\,\bigg{\\{}\prod_{j,k}\theta_{kj}^{\sum_{i}Y_{ijk}+a-1}e^{-\theta_{kj}-a\theta_{kj}/\mu_{k}}\bigg{\\}}\bigg{\\{}\prod_{k}\mu_{k}^{-Ja-
a_{0}-1}e^{-b_{0}/\mu_{k}}\bigg{\\}}$ (S2)

where $f(Y,R)$ is a function that does not depend on $\Theta$ or $\mu$. Hence,

$\displaystyle\pi(\mu\mid Y)$
$\displaystyle\underset{\mu}{\propto}\pi(Y,\mu)=\int\int\pi(Y,\Theta,R,\mu)\,d\Theta\,dR$
$\displaystyle=\bigg{\\{}\int
f(Y,R)\,dR\bigg{\\}}\bigg{\\{}\prod_{k}\mu_{k}^{-Ja-
a_{0}-1}e^{-b_{0}/\mu_{k}}\bigg{\\}}\bigg{\\{}\prod_{j,k}\int\theta_{kj}^{\sum_{i}Y_{ijk}+a-1}e^{-\theta_{kj}-a\theta_{kj}/\mu_{k}}d\theta_{kj}\bigg{\\}}$
$\displaystyle\underset{\mu}{\propto}\bigg{\\{}\prod_{k}\mu_{k}^{-Ja-
a_{0}-1}e^{-b_{0}/\mu_{k}}\bigg{\\}}\bigg{\\{}\prod_{j,k}\frac{\Gamma(a+\sum_{i}Y_{ijk})}{(1+a/\mu_{k})^{a+\sum_{i}Y_{ijk}}}\bigg{\\}}.$

Since this factors over $k$ into products that depend on
$\mu_{1},\ldots,\mu_{K}$, respectively, it follows that

$\displaystyle\pi(\mu_{k}\mid Y)$ $\displaystyle\propto\mu_{k}^{-Ja-
a_{0}-1}e^{-b_{0}/\mu_{k}}\Big{(}1+\frac{a}{\mu_{k}}\Big{)}^{-Ja-\sum_{i,j}Y_{ijk}}$
$\displaystyle\propto\mu_{k}^{-(a_{0}-\sum_{i,j}Y_{ijk})-1}\Big{(}\frac{\mu_{k}}{a}+1\Big{)}^{-Ja-\sum_{i,j}Y_{ijk}}e^{-b_{0}/\mu_{k}}$
$\displaystyle\propto\mathrm{InvKummer}\Big{(}\mu_{k}\;\Big{|}\;a_{0}+Ja,\,b_{0},\,Ja+\sum_{i,j}Y_{ijk},\,a\Big{)}$

by Definition 1. The proof is completed by letting $a_{0}=aJ+1$ and
$b_{0}=\varepsilon aJ$. ∎

#### Proof of Theorem 2

To establish Theorem 2, we prove the following concentration result for
inverse Kummer distributions. Theorem 2 then follows immediately from Theorem
1 and Theorem 3.

######  Theorem 3.

Let $\varepsilon>0$, $a>0$, and $y_{n}\geq 0$ such that $y_{n}\to y$ for some
$y\in[0,\infty)$ as $n\to\infty$. If
$\mu_{n}\sim\mathrm{InvKummer}(2an+1,\varepsilon an,ny_{n}+an,a)$ then for any
$d>0$,

$\mathds{P}(|\mu_{n}-\mu_{*}|\leq d)\xrightarrow[n\to\infty]{}1$

where
$\mu_{*}=2a\varepsilon/(\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon}-(y-a+\varepsilon))$.

###### Proof.

We first prove concentration of $T_{n}=a/\mu_{n}$, and use this to show the
concentration of $\mu_{n}$. Applying this change of variables to Equation 10
and simplifying, we find that the density of $T_{n}$ is

$f_{n}(t)\propto t^{2an}(1+t)^{-ny_{n}-an}e^{-\varepsilon
nt}=\exp(-ng_{n}(t)),$

where $g_{n}(t)$ and its derivative are

$\displaystyle g_{n}(t)$
$\displaystyle=-2a\log(t)+(y_{n}+a)\log(1+t)+\varepsilon t$ (S3)
$\displaystyle g_{n}^{\prime}(t)$
$\displaystyle=-\frac{2a}{t}+\frac{y_{n}+a}{1+t}+\varepsilon.$ (S4)

By Lemma 1, $g_{n}^{\prime}(t)$ has a unique zero in $(0,\infty)$ at

$t_{n}^{*}=\frac{\sqrt{(y_{n}-a+\varepsilon)^{2}+8a\varepsilon}-(y_{n}-a+\varepsilon)}{2\varepsilon},$
(S5)

$g_{n}^{\prime}(t)$ is strictly monotone increasing on $(0,t_{n}^{*})$, and
for all $t_{0}>0$, we have

$\inf_{t\geq
t_{0}}g_{n}^{\prime}(t)\geq\min\\{\varepsilon,g_{n}^{\prime}(t_{0})\\}.$ (S6)

Further, $g_{n}^{\prime}(t)<0$ for all $t\in(0,t_{n}^{*})$, and
$g_{n}^{\prime}(t)>0$ for all $t>t_{n}^{*}$. Define
$t^{*}=\lim_{n\to\infty}t^{*}_{n}$, noting that by continuity and the fact
that $y_{n}\to y$ by assumption,

$t^{*}=\frac{\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon}-(y-a+\varepsilon)}{2\varepsilon}.$
(S7)

Let $\delta\in(0,t_{*})$ be arbitrary. Let $B_{r}(t^{*})=\\{t>0:|t-t^{*}|\leq
r\\}$ denote the ball of radius $r$ around $t^{*}$ in $(0,\infty)$, and let
$B_{r}(t^{*})^{c}=\\{t>0:|t-t^{*}|>r\\}$ denote its complement in
$(0,\infty)$. Then

$\mathds{P}(|T_{n}-t^{*}|>\delta)=\frac{\int_{B_{\delta}(t^{*})^{c}}e^{-ng_{n}(t)}\mathrm{d}t}{\int_{0}^{\infty}e^{-ng_{n}(t)}\mathrm{d}t}\leq\frac{\int_{B_{\delta}(t^{*})^{c}}e^{-ng_{n}(t)}\mathrm{d}t}{\int_{B_{\delta}(t^{*})}e^{-ng_{n}(t)}\mathrm{d}t}.$
(S8)

We seek a suitable upper bound for the numerator and a lower bound for the
denominator to show that Equation S8 converges to 0 as $n\to\infty$. Define
$\alpha_{n}=\min\\{g_{n}(t^{*}+\delta),g_{n}(t^{*}-\delta)\\}$ and
$\beta_{n}=\min\\{\varepsilon,g_{n}^{\prime}(t^{*}+\delta),-g_{n}^{\prime}(t^{*}-\delta)\\}$.
To bound the numerator, we claim that for all $n$ sufficiently large,

$\displaystyle g_{n}(t)$
$\displaystyle\geq\alpha_{n}+\beta_{n}\big{(}t-(t^{*}+\delta)\big{)}\text{ for
$t>t^{*}+\delta$,}$ (S9) $\displaystyle g_{n}(t)$
$\displaystyle\geq\alpha_{n}-\beta_{n}\big{(}t-(t^{*}-\delta)\big{)}\text{ for
$0<t<t^{*}-\delta$.}$ (S10)

To see this, choose $n$ large enough that $|t_{n}^{*}-t^{*}|<\delta$. Then by
Equation S6 with $t_{0}=t^{*}+\delta$, we have
$g_{n}^{\prime}(t)\geq\min\\{\varepsilon,g_{n}^{\prime}(t^{*}+\delta)\\}\geq\beta_{n}$
for all $t\geq t^{*}+\delta$. Hence, for $t\geq t^{*}+\delta$,

$g_{n}(t)=g_{n}(t^{*}+\delta)+\int_{t^{*}+\delta}^{t}g_{n}^{\prime}(\tau)\mathrm{d}\tau\geq\alpha_{n}+\int_{t^{*}+\delta}^{t}\beta_{n}\mathrm{d}\tau=\alpha_{n}+\beta_{n}\big{(}t-(t^{*}+\delta)\big{)},$

proving Equation S9. Meanwhile, since $g_{n}^{\prime}(t)$ is monotone
increasing on $(0,t_{n}^{*})$ and $t^{*}-\delta\leq t_{n}^{*}$, we have
$g_{n}^{\prime}(t)\leq g_{n}^{\prime}(t^{*}-\delta)$ for $0<t<t^{*}-\delta$.
Thus, for $0<t<t^{*}-\delta$,

$g_{n}(t)=g_{n}(t^{*}-\delta)-\int_{t}^{t^{*}-\delta}g_{n}^{\prime}(\tau)\mathrm{d}\tau\geq\alpha_{n}+\int_{t}^{t^{*}-\delta}\beta_{n}\mathrm{d}\tau=\alpha_{n}-\beta_{n}\big{(}t-(t^{*}-\delta)\big{)},$

proving Equation S10. Hence, for the numerator in Equation S8, for all $n$
sufficiently large,

$\displaystyle\int_{B_{\delta}(t^{*})^{c}}e^{-ng_{n}(t)}\mathrm{d}t$
$\displaystyle=\int_{0}^{t^{*}-\delta}e^{-ng_{n}(t)}\mathrm{d}t+\int_{t^{*}+\delta}^{\infty}e^{-ng_{n}(t)}\mathrm{d}t$
$\displaystyle\leq\int_{0}^{t^{*}-\delta}e^{-n\big{(}\alpha_{n}-\beta_{n}(t-(t^{*}-\delta))\big{)}}\mathrm{d}t+\int_{t^{*}+\delta}^{\infty}e^{-n\big{(}\alpha_{n}+\beta_{n}(t-(t^{*}+\delta))\big{)}}\mathrm{d}t$
$\displaystyle\leq
2e^{-n\alpha_{n}}\int_{0}^{\infty}e^{-n\beta_{n}t}\mathrm{d}t$
$\displaystyle\leq\frac{2e^{-n\alpha_{n}}}{n\beta_{n}}.$ (S11)

Next, we turn to finding a lower bound for the denominator of Equation S8.
Since $y_{n}\to y$, we have that $g_{n}(t)\to g(t)$ and $g_{n}^{\prime}(t)\to
g^{\prime}(t)$ pointwise for all $t>0$, where $g(t)$ is defined by replacing
$y_{n}$ with $y$ in Equation S3. Define $\alpha=\lim_{n\to\infty}\alpha_{n}$
and $\beta=\lim_{n\to\infty}\beta_{n}$. By Lemma 1 applied to $g^{\prime}(t)$,
we have $g^{\prime}(t)<0$ for all $t\in(0,t^{*})$, and $g^{\prime}(t)>0$ for
all $t>t^{*}$. Thus, $\beta>0$ and $\alpha>g(t^{*})$ since $t^{*}$ is the
unique minimizer of $g(t)$. Pick any $\gamma$ such that
$g(t^{*})<\gamma<\alpha$. By the continuity of $g(t)$, we can choose
$r\in(0,\delta)$ such that $\sup_{t\in B_{r}(t^{*})}g(t)<\gamma$. Note that
$g_{n}$ converges uniformly to $g$ on $B_{\delta}(t^{*})$, since for all $t\in
B_{\delta}(t^{*})$,

$|g_{n}(t)-g(t)|=|y_{n}-y|\log(1+t)\leq|y_{n}-y|\log(1+t^{*}+\delta).$

Therefore, for all $n$ sufficiently large, $\sup_{t\in
B_{r}(t^{*})}g_{n}(t)<\gamma$, and hence,

$\int_{B_{\delta}(t^{*})}e^{-ng_{n}(t)}\mathrm{d}t\geq\int_{B_{r}(t^{*})}e^{-ng_{n}(t)}\mathrm{d}t\geq\int_{B_{r}(t^{*})}e^{-n\gamma}\mathrm{d}t=2re^{-n\gamma}.$
(S12)

Applying the bounds in Equations S11 and S12 to Equation S8 yields

$\mathds{P}(|T_{n}-t^{*}|>\delta)\leq\frac{\int_{B_{\delta}(t^{*})^{c}}e^{-ng_{n}(t)}\mathrm{d}t}{\int_{B_{\delta}(t^{*})}e^{-ng_{n}(t)}\mathrm{d}t}\leq\frac{2e^{-n\alpha_{n}}/(n\beta_{n})}{2re^{-n\gamma}}=\frac{e^{-n(\alpha_{n}-\gamma)}}{n\beta_{n}r}\xrightarrow[n\to\infty]{}0$

since $\alpha_{n}\to\alpha>\gamma$ and $\beta_{n}\to\beta>0$. Since
$\delta\in(0,t^{*})$ was arbitrary, this implies
$\mathds{P}(|T_{n}-t^{*}|\leq\delta)\to 1$ for all $\delta>0$.

Now, define $\mu_{*}=a/t^{*}$ and let $d\in(0,\mu_{*})$. Then
$|\mu_{n}-\mu_{*}|\leq d$ if and only if $\mu_{*}-d\leq\mu_{n}\leq\mu_{*}+d$,
or equivalently, $a/(\mu_{*}+d)\leq a/\mu_{n}\leq a/(\mu_{*}-d)$. Choose
$\delta>0$ small enough that

$\frac{a}{\mu_{*}+d}<\frac{a}{\mu_{*}}-\delta=t^{*}-\delta<t^{*}+\delta=\frac{a}{\mu_{*}}+\delta<\frac{a}{\mu_{*}-d}.$

Then, recalling that $T_{n}=a/\mu_{n}$, we have
$\mathds{P}(|\mu_{n}-\mu_{*}|\leq d)\geq\mathds{P}(|T_{n}-t^{*}|\leq\delta)\to
1$ as $n\to\infty$. ∎

######  Lemma 1.

Let $a>0$, $y\geq 0$, and $\varepsilon>0$, and define
$h(t)=-2a/t+(y+a)/(1+t)+\varepsilon$ for $t>0$. Then

  1. 1.

$h(t)$ has a unique zero on $(0,\infty)$ at
$t^{*}=\big{(}\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon}-(y-a+\varepsilon)\big{)}/(2\varepsilon)$,

  2. 2.

$h(t)<0$ for $0<t<t^{*}$ and $h(t)>0$ for $t>t^{*}$,

  3. 3.

$h(t)$ is strictly monotone increasing on $(0,t^{*})$, and

  4. 4.

for all $t_{0}>0$, we have $\inf_{t\geq
t_{0}}h(t)\geq\min\\{\varepsilon,h(t_{0})\\}$.

###### Proof.

To find the zeros, we set $h(t)=0$ and solve. For $t>0$, we have $h(t)=0$ if
and only if

$0=t(t+1)h(t)=-2a(t+1)+(y+a)t+\varepsilon t(t+1)=\varepsilon
t^{2}+(y-a+\varepsilon)t-2a.$

By the quadratic formula, this equation has zeros at
$t^{*}=(-b\pm\sqrt{b^{2}+8a\varepsilon})/(2\varepsilon)$ where
$b=y-a+\varepsilon$. Since $|-b|=\sqrt{b^{2}}<\sqrt{b^{2}+8a\varepsilon}$,
there is exactly one positive solution $t^{*}$ at

$t^{*}=\frac{\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon}-(y-a+\varepsilon)}{2\varepsilon}.$

Therefore, $h(t)$ has a unique zero on $(0,\infty)$ at $t^{*}$, proving part
1. Since $h(t)\to-\infty$ as $t\to 0$, and $h(t)\to\varepsilon$ as
$t\to\infty$, the intermediate value theorem implies that $h(t)<0$ for
$0<t<t^{*}$ and $h(t)>0$ for $t>t^{*}$. This proves part 2.

For part 3, to show that $h(t)$ is strictly increasing on $(0,t^{*})$, let
$0<t_{1}<t_{2}<t^{*}$. Define
$\tilde{h}(t)=-2a/t+(y+a)/(1+t)+\tilde{\varepsilon}$ where
$\tilde{\varepsilon}=\varepsilon-h(t_{2})$, noting that
$\tilde{\varepsilon}>0$ since $h(t_{2})<0$ by part 2. Since
$\tilde{h}(t_{2})=0$ by the definition of $\tilde{h}$, by applying part 1 to
$\tilde{h}$ we know $t_{2}$ is the unique zero of $\tilde{h}$. Therefore, by
applying part 2 to $\tilde{h}(t)$, we know $\tilde{h}(t_{1})<0$ since
$t_{1}<t_{2}$. Since $\tilde{h}(t)=h(t)-h(t_{2})$, this implies
$h(t_{1})<h(t_{2})$ as desired.

To show part 4, we use proof by contradiction. Let $t_{0}>0$, and suppose
there exists $t_{1}>t_{0}$ such that
$h(t_{1})<\min\\{\varepsilon,h(t_{0})\\}$. Define
$\tilde{h}(t)=-2a/t+(y+a)/(1+t)+\tilde{\varepsilon}$ where
$\tilde{\varepsilon}=\varepsilon-h(t_{1})$, noting that
$\tilde{\varepsilon}>0$ since $h(t_{1})<\varepsilon$ by assumption. Further,
$\tilde{h}(t_{1})=0$ by the definition of $\tilde{h}$, and
$\tilde{h}(t_{0})>\tilde{h}(t_{1})$ since
$\tilde{h}(t_{0})-\tilde{h}(t_{1})=h(t_{0})-h(t_{1})>0$ by assumption. Thus,
by the intermediate value theorem, there exists $t^{\prime}\in(0,t_{0})$ such
that $\tilde{h}(t^{\prime})=\tilde{h}(t_{1})=0$, because
$\tilde{h}(t)\to-\infty$ as $t\to 0$. But then $\tilde{h}(t)$ would have two
zeros in $(0,\infty)$, which is a contradiction. Therefore,
$h(t)\geq\min\\{\varepsilon,h(t_{0})\\}$ for all $t>t_{0}$, which proves part
4. ∎

#### Proof of Corollary 1

###### Proof.

Suppose $\bar{Y}_{k}\to 0$. Then in Theorem 2, $y=0$ and therefore the
concentration point is
$\mu^{*}=2a\varepsilon/(\sqrt{(\varepsilon-a)^{2}+8a\varepsilon}-(\varepsilon-a))$.
It holds that $\mu^{*}<\varepsilon$, since
$(\varepsilon-a)^{2}+8a\varepsilon>(\varepsilon-a)^{2}+4a\varepsilon=(\varepsilon+a)^{2}$,
and hence, $\sqrt{(\varepsilon-a)^{2}+8a\varepsilon}-(\varepsilon-a)>2a$. Let
$d=C\varepsilon-\mu^{*}$, noting that $d>0$ since $C>1$ and
$\mu^{*}<\varepsilon$. Then by Theorem 2,

$\mathds{P}(\mu_{k}>C\varepsilon\mid Y)\leq\mathds{P}(|\mu_{k}-\mu^{*}|>d\mid
Y)=1-\mathds{P}(|\mu_{k}-\mu^{*}|\leq d\mid
Y)\xrightarrow[J\to\infty]{}0.\qed$

##  S2 Additional results

###  S2.1 Further properties of the Inverse Kummer

We study the relationship between the hyperparameter $\gamma\in\mathds{R}$ of
the inverse Kummer distribution (Definition 1), and the first moment of the
distribution. In particular, the following proposition implies that the mean
of an inverse Kummer with $\gamma>0$ is larger than the mean of the
corresponding inverse gamma distribution, when $\lambda>2$.

######  Proposition 1.

Let $\mu\sim\mathrm{InvKummer}(\lambda,\beta,\gamma,\delta)$. If $\lambda>2$,
then $\mathds{E}(\mu)$ is monotone increasing as a function of $\gamma$, for
$\gamma\in(0,\infty)$.

###### Proof.

Fix $\epsilon\in(0,1)$. Let
$\mu_{\gamma}\sim\mathrm{InvKummer}(\lambda,\beta,\gamma,\delta)$ and
$\mu_{\gamma+\epsilon}\sim\mathrm{InvKummer}(\lambda,\beta,\gamma+\epsilon,\delta)$,
and define $f(t)=t^{\lambda-1}(1+t)^{-\gamma}e^{-\beta t/\delta}$. Our aim is
to show that $\mathds{E}(\mu_{\gamma+\epsilon})\geq\mathds{E}(\mu_{\gamma})$.
By Equation 11, this is true if and only if

$\frac{\int_{0}^{\infty}\frac{1}{t(1+t)^{\epsilon}}f(t)\mathrm{d}{t}}{\int_{0}^{\infty}\frac{1}{(1+t)^{\epsilon}}f(t)\mathrm{d}{t}}\geq\frac{\int_{0}^{\infty}\frac{1}{t}f(t)\mathrm{d}{t}}{\int_{0}^{\infty}f(t)\mathrm{d}{t}}.$
(S13)

Let $X$ be the continuous random variable on $(0,\infty)$ with probability
density function $p(x)=f(x)/\int_{0}^{\infty}f(t)\mathrm{d}{t}$. Then, after
multiplying and dividing the left-hand side by
$\int_{0}^{\infty}f(t)\mathrm{d}{t}$, Equation S13 can be written in terms of
expectations as

$\frac{\mathds{E}\Big{(}\frac{1}{X(1+X)^{\epsilon}}\Big{)}}{\mathds{E}\Big{(}\frac{1}{(1+X)^{\epsilon}}\Big{)}}\geq\mathds{E}\Big{(}\frac{1}{X}\Big{)}.$
(S14)

or equivalently,

$\mathrm{Cov}\Big{(}\frac{1}{X},\frac{1}{(1+X)^{\epsilon}}\Big{)}\geq 0.$
(S15)

Define $g(x)=-1/x$ and $h(x)=-1/(1+x)^{\epsilon}$ for $x\in(0,\infty)$.
Observe that, by transformation of random variables,
$1/X\sim\mathrm{InvKummer}(\lambda,\beta,\gamma,\delta)$. Thus, by Equation
11, since $\lambda>2$ by assumption,

$\displaystyle\mathds{E}|g(X)|^{2}$
$\displaystyle=\mathds{E}\Big{(}\frac{1}{X^{2}}\Big{)}<\infty,$
$\displaystyle\mathds{E}|h(X)|^{2}$
$\displaystyle=\mathds{E}\bigg{(}\frac{1}{(1+X)^{2\epsilon}}\bigg{)}\leq\mathds{E}\Big{(}\frac{1}{X^{2\epsilon}}\Big{)}<\infty.$

Therefore, since $g(x)$ and $h(x)$ are monotone increasing and have finite
second moments, we have $\mathrm{Cov}(g(X),h(X))\geq 0$ by
schmidt2003covariance, which is equivalent to Equation S15. This completes the
proof. ∎

###  S2.2 Characterizing the concentration point of the inverse Kummer

We characterize the point at which the inverse Kummer concentrates, $\mu_{*}$,
in Theorem 2. Define

$\mu_{*}(\varepsilon,y,a)=\frac{2a\varepsilon}{\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon}-(y-a+\varepsilon)}$
(S16)

for $\varepsilon>0$, $y\geq 0$, and $a>0$.

######  Proposition 2.

For all $\varepsilon>0$ and $a>0$, $\mu_{*}(\varepsilon,y,a)$ is monotone
increasing as a function of $y$.

###### Proof.

Fix $\varepsilon>0$ and define

$g(y)=\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon},$

so that $\mu_{*}(\varepsilon,y,a)=2a\varepsilon/(g(y)-(y-a+\varepsilon))$.
Then

$\frac{\partial\mu_{*}}{\partial
y}=\frac{-2a\varepsilon\big{(}g^{\prime}(y)-1\big{)}}{\big{(}g(y)-(y-a+\varepsilon)\big{)}^{2}}.$

Differentiating $g$, we find that

$g^{\prime}(y)=\frac{y-a+\varepsilon}{\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon}}<1.$

Therefore, $\partial\mu_{*}/\partial y>0$, showing that $\mu_{*}$ is monotone
increasing as a function of $y$. ∎

Next, we derive Equation 12 using a first-order Taylor approximation. Fix
$y\geq 0$ and $a>0$, and define

$h(\varepsilon)=\sqrt{(y-a+\varepsilon)^{2}+8a\varepsilon}$ (S17)

for $\varepsilon>0$. Differentiating and simplifying, we find that

$h^{\prime}(\varepsilon)=\frac{y+\varepsilon+3a}{h(\varepsilon)}.$ (S18)

Thus, $h(0)=|y-a|$ and $h^{\prime}(0)=(y+3a)/|y-a|$. Hence, a first-order
Taylor approximation to $h$ at $\varepsilon=0$ yields

$h(\varepsilon)\approx
h(0)+h^{\prime}(0)\varepsilon=|y-a|+\frac{y+3a}{|y-a|}\varepsilon$ (S19)

when $\varepsilon$ is small relative to $|y-a|$. Plugging this into the
definition of $\mu_{*}$ in Equation S16, we obtain

$\mu_{*}(\varepsilon,y,a)\approx\frac{2a\varepsilon}{|y-a|+\frac{y+3a}{|y-a|}\varepsilon-(y-a+\varepsilon)}.$
(S20)

When $y>a$, we have $|y-a|=y-a$, so in this case Equation S20 becomes

$\mu_{*}(\varepsilon,y,a)\approx\frac{y-a}{2}.$

Meanwhile, when $0\leq y<1$, we have $|y-a|=a-y$, so in this case

$\mu_{*}(\varepsilon,y,a)\approx\frac{\varepsilon
a(a-y)}{(a-y)^{2}+(a+y)\varepsilon}$

by collecting and rearranging terms. Therefore, the first-order Taylor
approximation is

$\mu_{*}(\varepsilon,y,a)\approx\begin{cases}\displaystyle\frac{y-a}{2}&\text{
if $y>a$}\\\ \displaystyle\frac{\varepsilon
a(a-y)}{(a-y)^{2}+(a+y)\varepsilon}&\text{ if $0\leq y<a$}\end{cases}$ (S21)

as claimed in Equation 12.

###  S2.3 Relationship between the relevance weights and the latent counts

We further characterize the relationship between the latent counts $Y_{ijk}$
and relevance weights $\mu_{k}$. In particular, we derive the distribution of
$Y_{ijk}\mid\mu_{k}$, integrating out $r_{k}$ and $\theta_{k}$ in Equation 7.
We show that, appealingly, this distribution has a closed-form expression in
terms of hypergeometric functions. First, the distribution of
$Y_{ijk}\mid\mu_{k},r_{k}$, integrating out $\theta_{kj}$ in Equation 7, is
easily seen to be
$Y_{ijk}\mid\mu_{k},r_{k}\sim\mathrm{NegBin}\big{(}a,a/(a+r_{ik}\mu_{k})\big{)}$,
where the negative binomial is parametrized such that the mean and variance
are $\mathds{E}(Y_{ijk}\mid\mu_{k},r_{k})=\mu_{k}r_{ik}$ and
$\mathrm{Var}(Y_{ijk}\mid\mu_{k},r_{k})=\mu_{k}r_{ik}(a+\mu_{k}r_{ik})$. When
both $r_{k}$ and $\theta_{k}$ are integrated out, we obtain the following
result, where
$\prescript{}{2}{F}_{1}(a,b,;c,z)=\sum_{n=0}^{\infty}\frac{(a)_{n}(b)_{n}}{(c)_{n}}\frac{z^{n}}{n!}$
denotes the Gauss-hypergeometric function and $(a)_{n}=\Gamma(a+n)/\Gamma(a)$
is the ascending factorial abramowitz_stegun_1972.

######  Proposition 3.

The probability mass function of $Y_{ijk}\mid\mu_{k}$ under the model in
Equation 7 is

$\displaystyle\mathds{P}(Y_{ijk}$
$\displaystyle=y\mid\mu_{k})=\left(\frac{\mu_{k}}{a}\right)^{y}\frac{(a)_{y}(\alpha)_{y}}{y!(\alpha
I)_{y}}\prescript{}{2}{F}_{1}\Big{(}y+a,\,y+\alpha,\,y+\alpha
I,\,-\frac{\mu_{k}}{a}\Big{)},$

for $y\in\\{0,1,2,\ldots\\}$. Furthermore, the mean of this distribution is
$\mathds{E}(Y_{ijk}\mid\mu_{k})=\mu_{k}/I$.

###### Proof.

Recall that
$Y_{ijk}\mid\mu_{k},r_{k},\theta_{k}\sim\mathrm{Poisson}(r_{ik}\theta_{kj})$
and $\theta_{kj}\mid\mu_{k}\sim\mathrm{Gamma}(a,a/\mu_{k})$. Also,
$r_{ik}\sim\mathrm{Beta}(\alpha,\,\alpha I-\alpha)$ by marginalization
property of the Dirichlet distribution. Hence,

$\displaystyle\pi(Y_{ijk}=y\mid\mu_{k})=\int_{0}^{1}\bigg{\\{}\int_{0}^{\infty}\mathds{P}(Y_{ijk}=y\mid\mu_{k},r_{k},\theta_{k})p(\theta_{kj}\mid\mu_{k})d\theta_{kj}\bigg{\\}}p(r_{ik})dr_{ik}$
$\displaystyle=\int_{0}^{1}\bigg{\\{}\int_{0}^{\infty}(r_{ik}\theta_{kj})^{y}\frac{e^{-r_{ik}\theta_{kj}}}{y!}\frac{(a/\mu_{k})^{a}}{\Gamma(a)}\theta_{kj}^{a-1}e^{-(a/\mu_{k})\theta_{kj}}d\theta_{kj}\bigg{\\}}\frac{\Gamma(\alpha
I)}{\Gamma(\alpha)\Gamma(\alpha I-\alpha)}r_{ik}^{\alpha-1}(1-r_{ik})^{\alpha
I-\alpha-1}dr_{ik}$
$\displaystyle=\frac{(a/\mu_{k})^{a}}{y!\,\Gamma(a)}\frac{\Gamma(\alpha
I)}{\Gamma(\alpha)\Gamma(\alpha
I-\alpha)}\int_{0}^{1}\bigg{\\{}\int_{0}^{\infty}\theta_{kj}^{y+a-1}e^{-(r_{ik}+a/\mu_{k})\theta_{kj}}d\theta_{kj}\bigg{\\}}r_{ik}^{y+\alpha-1}(1-r_{ik})^{\alpha
I-\alpha-1}dr_{ik}$
$\displaystyle=\frac{(a/\mu_{k})^{a}}{y!\,\Gamma(a)}\frac{\Gamma(\alpha
I)}{\Gamma(\alpha)\Gamma(\alpha
I-\alpha)}\int_{0}^{1}\frac{\Gamma(y+a)}{(r_{ik}+a/\mu_{k})^{y+a}}r_{ik}^{y+\alpha-1}(1-r_{ik})^{\alpha
I-\alpha-1}dr_{ik}$
$\displaystyle=\frac{(a)_{y}(\mu_{k}/a)^{y}}{y!}\frac{\Gamma(\alpha
I)}{\Gamma(\alpha)\Gamma(\alpha
I-\alpha)}\int_{0}^{1}(t\mu_{k}/a+1)^{-(y+a)}t^{y+\alpha-1}(1-t)^{(y+\alpha
I)-(y+\alpha)-1}dt$
$\displaystyle=\Big{(}\frac{\mu_{k}}{a}\Big{)}^{y}\frac{(a)_{y}(\alpha)_{y}}{y!(\alpha
I)_{y}}\prescript{}{2}{F}_{1}\big{(}y+a,\,y+\alpha,\,y+\alpha
I,\,-\mu_{k}/a\big{)},$

where $(x)_{n}=\Gamma(x+n)/\Gamma(x)$ is the ascending factorial, and

$\prescript{}{2}{F}_{1}(a,b,c,z)=\frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_{0}^{1}t^{b-1}(1-t)^{c-b-1}(1-tz)^{-a}\mathrm{d}t,$

with $c>b>0$ and $a,z\in\mathds{R}$ is an alternative representation of the
Gauss-hypergeometric function; see abramowitz_stegun_1972. ∎

This result further explains the role of each $\mu_{k}$ in the mutational
process: it directly controls the contribution of signature $k$ in determining
the number of mutations $X_{ij}$ in channel $i$ for patient $j$, for given
values of hyperparameters $a$ and $\alpha$.

##  S3 Gibbs sampler for NMF with informative priors

In this section, we present the general Gibbs sampler for the model in
Equation 8 with priors as in Equation 9. Each step follows from simple semi-
conjugate prior updates, so we omit the derivations. Note that we set the
value of $\beta_{k}$ depending on the level of sparsity of each cosmic
signature. Specifically, for a given $s_{k}$, we calculate $\beta_{k}$ by
drawing 1000 samples
$\rho_{k}\sim\mathrm{Dirichlet}(\beta_{k}s_{1k},\ldots,\beta_{k}s_{Ik})$ for a
range of plausible $\beta_{k}$ values (from 10 to 5000, evenly spaced on a log
scale), and we select the value for which the median cosine similarity between
$s_{k}$ and the sampled $\rho_{k}$ vectors is closest to 0.975. This ensures
that all signatures have approximately equal variance under the prior. All
values are available at https://github.com/alessandrozito/CompressiveNMF.

Inference in the CompNMF+cosmic model is performed by iterating the following
steps.

  1. 1.

For $i=1,\ldots,I$ and $j=1,\ldots,J$, update the latent mutation counts by
drawing

$(Y_{ij}\mid-)\sim\mathrm{Multinomial}\big{(}X_{ij},(\tilde{q}_{ij1},\ldots,\tilde{q}_{ijK_{\mathrm{pre}}},q_{ij1},\ldots,q_{ijK_{\mathrm{new}}})\big{)}$

where $\tilde{q}_{ijk}=\rho_{ik}\omega_{kj}/Q_{ij}$ and
$q_{ijk}=r_{ik}\theta_{kj}/Q_{ij}$, with
$Q_{ij}=\sum_{k=1}^{K_{\mathrm{pre}}}\rho_{ik}\omega_{kj}+\sum_{k=1}^{K_{\mathrm{new}}}r_{ik}\theta_{kj}$.

  2. 2.

For $k=1,\ldots,K_{\mathrm{pre}}$, update the cosmic signatures by drawing

$(\rho_{k}\mid-)\sim\mathrm{Dirichlet}\bigg{(}\beta_{k}s_{1k}+\sum_{j=1}^{J}Y_{1jk},\ldots,\beta_{k}s_{Ik}+\sum_{j=1}^{J}Y_{Ijk}\bigg{)}.$

  3. 3.

For $k=K_{\mathrm{new}}+1,\ldots,K_{\mathrm{pre}}+K_{\mathrm{new}}$, update
the _de novo_ signatures by drawing

$(r_{k}\mid-)\sim\mathrm{Dirichlet}\bigg{(}\alpha+\sum_{j=1}^{J}Y_{1jk},\,\ldots,\,\alpha+\sum_{j=1}^{J}Y_{Ijk}\bigg{)}.$

  4. 4.

For $k=1,\ldots,K_{\mathrm{pre}}$ and $j=1,\ldots,J$, update the loadings
associated to the cosmic signatures by drawing

$(\omega_{kj}\mid-)\sim\mathrm{Gamma}\bigg{(}b+\sum_{i=1}^{I}Y_{ijk},\;\frac{b}{\tau_{k}}+1\bigg{)}.$

  5. 5.

For $k=K_{\mathrm{new}}+1,\ldots,K_{\mathrm{pre}}+K_{\mathrm{new}}$ and
$j=1,\ldots,J$, update the loadings associated to the _de novo_ signatures by
drawing

$(\theta_{kj}\mid-)\sim\mathrm{Gamma}\bigg{(}a+\sum_{i=1}^{I}Y_{ijk},\;\frac{a}{\mu_{k}}+1\bigg{)}.$

  6. 6.

For $k=1,\ldots,K_{\mathrm{pre}}$, update the relevance weights associated to
the cosmic signatures by drawing

$(\tau_{k}\mid-)\sim\mathrm{InvGamma}\bigg{(}2bJ+1,\;\varepsilon
bJ+b\sum_{j=1}^{J}\omega_{kj}\bigg{)}.$

  7. 7.

For $k=K_{\mathrm{new}}+1,\ldots,K_{\mathrm{pre}}+K_{\mathrm{new}}$, update
the relevance weights associated to the _de novo_ signatures by drawing

$(\mu_{k}\mid-)\sim\mathrm{InvGamma}\bigg{(}2aJ+1,\;\varepsilon
aJ+a\sum_{j=1}^{J}\theta_{kj}\bigg{)}.$

One important behavior we noticed is that, occasionally, the sampler either
(i) morphs a novel signature into a cosmic one even if that cosmic signature
has been specified in the prior, or (ii) morphs an existing cosmic signature
into another cosmic one. This is due to the multi-modal nature of the NMF
model and can be influenced by initialization. While carefully eliciting each
$\beta_{k}$ as above does help, such incoherence sometimes can hold for
relatively flat cosmic signatures, such as SBS3, SBS5, or SBS40a,b. To solve
the issue, we apply a label-switching step at 2/3 of the burn-in phase, where
the signatures that have not been compressed out of the model are re-matched
to the cosmic signatures via the Hungarian algorithm. This does not invalidate
the MCMC algorithm since it is only performed in the burn-in phase.

##  S4 Rationale for the model likelihood

###  S4.1 Types of base pair substitutions

In DNA, there are four bases: cytosine (C), thymine (T), adenine (A), and
guanine (G). Considering both strands of the double helix, cytosine always
pairs with guanine, and thymine always pairs with adenine. Thus, if we
distinguish one of the two strands of a given DNA molecule, there are four
possible base pairs at each point: C-G, G-C, T-A, and A-T.

When considering base pair substitutions at a given point, the convention is
to distinguish the strand containing the pyrimidine (C or T) before the
substitution has been made. Recall that cytosine (C) and thymine (T) are
_pyrimidines_ , whereas adenine (A) and guanine (G) are _purines_. With this
convention, there are six possible types of substitutions at any given point:

| before | after | abbreviation
---|---|---|---
1 | C-G | A-T | C>A
2 | C-G | G-C | C>G
3 | C-G | T-A | C>T
4 | T-A | A-T | T>A
5 | T-A | C-G | T>C
6 | T-A | G-C | T>G

Sometimes, these are abbreviated denoting only the pre-substitution pyrimidine
and what it changes to, as seen above.

These six classes can be further divided by considering the trinucleotide
context, that is, the bases directly adjacent to the base undergoing
substitution. The convention is to label the context in terms of the bases (C,
T, A, or G) on the 5’ and 3’ sides on the strand containing the pre-
substitution pyrimidine. For instance, in a substitution C>A, the C may be
flanked by a T on the 5’ side and a G on the 3’ side:

before | after | abbreviation
---|---|---
TCG | TAG | T[C>A]G
5’ 3’ |  5’ 3’ |

There are $4\times 4=16$ different contexts for each of the original six
substitution types. Therefore, there are $16\times 6=96$ single-base
substitution types when the trinucleotide context is taken into account.At
each position in the genome, one of the two strands contains a pyrimidine C or
T, flanked by bases on the 5’ and 3’ sides, say, X and Y, respectively: so the
trinucleotide context is either XCY or XTY. Thus, each position in the genome
can be in one of $2\times 4\times 4=32$ possible states. Since there are three
possible single-base substitutions at every position, we arrive at a total of
$32\times 3=96$ _mutational channels_.

###  S4.2 Continuous-time Markov process for substitutions

Focusing on one position $\ell$ in the genome, let us assume mutations at
$\ell$ occur as a time-homogeneous continuous-time Markov process, holding the
neighboring bases fixed. More precisely, when the current state is $a$, it
remains $a$ for an $\mathrm{Exponential}(|\Lambda_{aa}|)$ amount of time and
then transitions to $b\neq a$ with probability $\Lambda_{ab}/|\Lambda_{aa}|$,
where $\Lambda$ is a $32\times 32$ matrix such that (i) $\Lambda_{ab}\geq 0$
for $a\neq b$, and (ii) $\sum_{b}\Lambda_{ab}=0$. This is equivalent to saying
that transitions from $a$ to $b$ occur with rate $\Lambda_{ab}$; thus,
$\Lambda$ is called the transition rate matrix. See lawler2018introduction for
background.

Let $S^{t}_{\ell}$ denote the state at locus $\ell$ at time $t$, and let
$S^{0}_{\ell}$ be the state at $\ell$ for the normal (germline) genome of the
individual under consideration. Let $P^{t}_{ab}=\mathds{P}(S^{t}_{\ell}=b\mid
S^{0}_{\ell}=a)$ be the probability that the state is $b$ at time $t$ given
that the state is $a$ at time 0. From the theory of continuous-time Markov
processes, we have that

$P^{t}=\exp(t\Lambda)=\sum_{k=0}^{\infty}\frac{(t\Lambda)^{k}}{k!}$

where $\exp(\cdot)$ denotes the matrix exponential. Since the mutation rates
$\Lambda_{ab}$ are very small, it is reasonable to use a first-order Taylor
approximation, $P^{t}\approx I+t\Lambda$.

###  S4.3 Substitution counts are approximately Poisson distributed

Let $a_{i}$ and $b_{i}$ denote the starting and ending states, respectively,
for each of the substitution types $i=1,\ldots,96$. Let
$\lambda_{i}=\Lambda_{a_{i}b_{i}}$, and define
$X^{t}_{i}=\\#\\{\ell:S^{0}_{\ell}=a_{i},S^{t}_{\ell}=b_{i}\\}$, that is,
$X^{t}_{i}$ is the number of positions in the genome that undergo substitution
$i$, starting at state $a_{i}$ at time $0$ and ending at state $b_{i}$ at time
$t$.

Now, consider all of the positions $\ell$ that are in state $a$ at time 0, and
to simplify the math, let us assume that (a) no two of these positions are
adjacent, and (b) that substitions occur independently across positions. Of
the 32 states, only four of them can be reached from $a$: the state can remain
at $a$, or one of three substitutions can occur. Suppose these three
substitutions are $i=1,2,3$, so that the starting states are
$a_{1}=a_{2}=a_{3}=a$ and the ending states are $b_{1},b_{2},b_{3}$,
respectively. Let $s^{0}=(s_{\ell}^{0}:\ell=1,\ldots,L)$ be a fixed vector of
starting states for all positions $\ell$, and let
$n=\\#\\{\ell:s^{0}_{\ell}=a\\}$ be the number of positions starting in state
$a$ at time $0$. Let $\lambda_{0}=\lambda_{1}+\lambda_{2}+\lambda_{3}$ be the
sum of the rates for substitution types $1,2,3$. By the definition of
$P^{t}_{ab}$ and the assumption of independence across positions, letting
$X_{0}^{t}=X_{1}^{t}+X_{2}^{t}+X_{3}^{t}$, the vector
$(X_{1}^{t},X_{2}^{t},X_{3}^{t},n-X_{0}^{t})$ follows a multinomial
distribution. Specifically, for non-negative integers $x_{1},x_{2},x_{3}$ such
that $x_{0}:=x_{1}+x_{2}+x_{3}\leq n$, we have

$\displaystyle\mathds{P}(X^{t}_{1:3}=x_{1:3}\mid S^{0}=s^{0})$
$\displaystyle=\frac{n!}{(n-x_{0})!x_{1}!x_{2}!x_{3}!}(P^{t}_{aa})^{n-x_{0}}\prod_{i=1}^{3}(P^{t}_{a_{i}b_{i}})^{x_{i}}$
$\displaystyle\approx\frac{n!}{(n-x_{0})!x_{1}!x_{2}!x_{3}!}(1-t\lambda_{0})^{n-x_{0}}(t\lambda_{1})^{x_{1}}(t\lambda_{2})^{x_{2}}(t\lambda_{3})^{x_{3}}$
(S22)

by the first-order Taylor approximation, $P^{t}\approx I+t\Lambda$. Since the
genome is large and mutation rates are small, it is natural to assume that $n$
is large and $t\lambda_{0}=O(1/n)$. Thus, letting $c=nt\lambda_{0}$ we have
$(1-t\lambda_{0})^{n}=(1-c/n)^{n}\approx e^{-c}=\exp(-nt\lambda_{0})$ and
$(1-t\lambda_{0})^{-x_{0}}=(1-c/n)^{-x_{0}}\approx 1$ when $x_{0}\ll n$, which
is the case with high probability. Plugging these approximations into Section
S4.3 yields

$\displaystyle\approx\frac{n!}{(n-x_{0})!x_{1}!x_{2}!x_{3}!}\exp(-nt\lambda_{0})(t\lambda_{1})^{x_{1}}(t\lambda_{2})^{x_{2}}(t\lambda_{3})^{x_{3}}$
$\displaystyle=\frac{n!\,n^{-x_{0}}}{(n-x_{0})!}\prod_{i=1}^{3}\exp(-nt\lambda_{i})\frac{(nt\lambda_{i})^{x_{i}}}{x_{i}!}$

since $x_{0}=x_{1}+x_{2}+x_{3}$ and
$\lambda_{0}=\lambda_{1}+\lambda_{2}+\lambda_{3}$ by definition. By Stirling’s
approximation,

$\frac{n!\,n^{-x_{0}}}{(n-x_{0})!}\sim\frac{\sqrt{2\pi
n}\,(n/e)^{n}n^{-x_{0}}}{\sqrt{2\pi(n-x_{0})}\,((n-x_{0})/e)^{n-x_{0}}}=\sqrt{\frac{n}{n-x_{0}}}\frac{e^{-x_{0}}n^{n}}{(n-x_{0})^{n}}\frac{(n-x_{0})^{x_{0}}}{n^{x_{0}}}\longrightarrow
1$

as $n\to\infty$ with $x_{0}$ fixed, since $(1-x_{0}/n)^{n}\to e^{-x_{0}}$.
Hence, we have

$\displaystyle\mathds{P}(X^{t}_{1:3}=x_{1:3}\mid
S^{0}=s^{0})\approx\prod_{i=1}^{3}\mathrm{Poisson}(x_{i}\mid nt\lambda_{i})$

when $n$ is large, $x_{0}\ll n$, and $t\lambda_{0}=O(1/n)$.

For each of the 32 distinct possible starting states $a$, the same
approximation applies to the set of positions starting in state $a$. Modeling
these 32 sets of positions independently, we have

$\displaystyle\mathds{P}(X^{t}_{1:96}=x_{1:96}\mid
S^{0}=s^{0})\approx\prod_{i=1}^{96}\mathrm{Poisson}(x_{i}\mid
n_{i}t\lambda_{i})$

where $n_{i}=\\#\\{\ell:s^{0}_{\ell}=a_{i}\\}$. In other words, the counts of
the $96$ substitution types are approximately distributed as independent
Poisson random variables with rates $n_{i}t\lambda_{i}$.

The preceding derivation ignores the fact that a substitution at one position
changes the context of the two adjacent positions. However, since it is rare
for single-base substitutions to occur at two adjacent positions, the effect
of ignoring this should be negligible.

###  S4.4 Multiple mutational processes

Suppose $x_{ij}$ is the number of mutations of substitution type $i$ for
subject $j$, for $i=1,\ldots,I$ and $j=1,\ldots,J$, where $I=96$. The
derivation above justifies modeling these mutation counts as

$X_{ij}\sim\mathrm{Poisson}(n_{ij}\lambda_{ij}t_{j})$

independently, where $t_{j}$ is the age or exposure time of subject $j$,
$\lambda_{ij}$ is the mutation rate for substitution type $i$ in subject $j$,
and $n_{ij}$ is the number of positions that are in state $a_{i}$ in the
normal genome of subject $j$, out of all positions that were measured. The
positions measured may be a subset of the genome due to whole-exome/targeted
sequencing or low sequencing depth, for example.

From birth, each subject is exposed to many mutational processes, such as
environmental exposures, replication errors, defective DNA repair mechanisms,
and so on. Each mutational process causes each substitution type to occur at a
given rate, and the profile of rates across the 96 substitution types can be
expected to vary depending on the mutational process. Since rates are additive
in a continuous-time Markov process, it is natural to model the subject-
specific mutation rates $\lambda_{ij}$ as linear combinations of these
mutational process rate profiles, with non-negative weights depending on the
exposure of the subject to each process. Further, assuming the opportunity
counts $n_{ij}$ are constant (or nearly constant) across all subjects $j$, one
can absorb $n_{ij}$ into $\lambda_{ij}$, which changes the interpretation of
$\lambda_{ij}$ by reparametrizing it. This leads to using a representation of
the form

$n_{ij}\lambda_{ij}t_{j}=\sum_{k=1}^{K}r_{ik}\theta_{kj}$

where the weight $\theta_{kj}\geq 0$ is the exposure of subject $j$ to process
$k$, and $(r_{1k},\ldots,r_{Ik})$ is the mutation rate profile for mutational
process $k$, which is referred to as its mutational signature. Thus, we arrive
at the Poisson non-negative matrix factorization model in Equation 1,

$X_{ij}\sim\mathrm{Poisson}\Big{(}\sum_{k=1}^{K}r_{ik}\theta_{kj}\Big{)}.$

A statistical issue with this representation is that there is a non-
identifiability between the $r_{ik}$’s and $\theta_{kj}$’s, since arbitrary
multiplicative constants $c_{k}$ can be moved between them. We deal with this
by normalizing the mutational signatures to sum to $1$, that is, by enforcing
the constraint $\sum_{i=1}^{I}r_{ik}=1$ for all $k$.

##  S5 Simulation details and additional results

###  S5.1 Details of PoissonCUSP

The cusp model described in Legramanti_2020 is a spike-and-slab shrinkage
prior that enables automatic selection of the number of latent factors in
Gaussian factorization models. We adapt it to the Poisson factorization model
as follows. We begin by specifying the following prior structure for the
signatures and the loadings:

$(r_{1k},\ldots,r_{Ik})\sim\mathrm{Dirichlet}(\alpha,\ldots,\alpha),\quad\theta_{kj}=\vartheta_{kj}\mu_{k},\quad\vartheta_{kj}\sim\mathrm{Gamma}(a,a).$

Then, we let

$\mu_{k}\sim(1-\pi_{k})\mathrm{Gamma}(a_{0},b_{0})+\pi_{k}\delta_{\mu_{\infty}},$

where $\delta_{x}$ denotes the point mass at $x$ and $\pi_{k}$ is the prior
probability of sampling the spike $\mu_{\infty}=0.01$, which is modeled as

$\pi_{k}=\sum_{\ell=1}^{k}\phi_{\ell},\quad\phi_{\ell}=v_{\ell}\prod_{m=1}^{\ell-1}(1-v_{m}),\quad
v_{\ell}\sim\mathrm{Beta}(1,\alpha).$

Hence, $\pi_{k}$ increases as $k$ increases. As a spike-and-slab prior, this
enables automatic selection of the number of signatures. Each iteration of the
Gibbs sampler for the basic PoissonCUSP model consists of the following steps:

  1. 1.

For $i=1,\ldots,I$ and $j=1,\ldots,J$, sample auxiliary variables
$Y_{ij}=(Y_{ij1},\ldots,Y_{ijK})$ according to
$Y_{ij}\sim\mathrm{Multinomial}\big{(}X_{ij},(q_{ij1},\ldots,q_{ijK})\big{)},$
where
$q_{ijk}=r_{ik}\theta_{kj}/\sum_{\kappa=1}^{K}r_{i\kappa}\theta_{\kappa,j}$.

  2. 2.

Sample the individual loadings $\vartheta_{kj}$ from
$\vartheta_{kj}\sim\mathrm{Gamma}\big{(}a+\sum_{i=1}^{I}Y_{ijk},\,a+\mu_{k}\big{)}.$

  3. 3.

Sample the signatures $r_{k}=(r_{1k},\ldots,r_{ik})$ from
$r_{k}\sim\mathrm{Dirichlet}\big{(}\alpha+\sum_{j=1}^{J}Y_{1jk},\ldots,\alpha+\sum_{j=1}^{J}Y_{ijk}\big{)}.$

  4. 4.

For each $k=1,\ldots,K$, sample the categorical auxiliary variables $Z_{k}$ as
follows

$\mathds{P}(Z_{k}=\ell\mid-)=\begin{cases}\phi_{\ell}\,\mu_{\infty}^{\sum_{i,j}Y_{ijk}}\exp(-\mu_{\infty}\sum_{j}\vartheta_{kj})&\text{
if }1\leq\ell\leq k\\\\[8.0pt]
\displaystyle\phi_{\ell}\,\frac{b_{0}^{a_{0}}}{\Gamma(a_{0})}\frac{\Gamma(a_{0}+\sum_{ij}Y_{ijk})}{(b_{0}+\sum_{j}\vartheta_{kj})^{a_{0}+\sum_{ij}Y_{ijk}}}&\text{
if }k<\ell\leq K\end{cases}$

where $\phi_{\ell}=\mathds{P}(Z_{k}=\ell)$ is the prior probability of the
auxiliary variables. These probabilities are derived by integrating out the
parameters $\mu_{k}$, relying on the conjugacy of the model.

  5. 5.

For $\ell=1,\ldots,K-1$, sample the sticks from

$v_{\ell}\sim\mathrm{Beta}\bigg{(}1+\sum_{k=1}^{K}\mathds{1}(Z_{k}=\ell),\;\alpha_{\pi}+\sum_{k=1}^{K}\mathds{1}(Z_{k}>\ell)\bigg{)},$

and fix $v_{K}=1$.

  6. 6.

Calculate $\phi_{1},\ldots,\phi_{K}$ via the stick-breaking construction,
namely $\phi_{\ell}=v_{\ell}\prod_{m=1}^{\ell-1}(1-v_{m})$.

  7. 7.

For $k=1,\ldots,K$, if $Z_{k}\leq k$ then set $\mu_{k}=\delta_{\mu_{\infty}}$,
otherwise sample

$\mu_{k}\sim\mathrm{Gamma}\bigg{(}a_{0}+\sum_{i=1}^{I}\sum_{j=1}^{J}Y_{ijk},\;b_{0}+\sum_{j=1}^{J}\vartheta_{kj}\bigg{)}.$

Notice that here we assume a multiplicative structure for $\theta_{kj}$, while
our compressive hyperprior approach assumes a hierarchical one. In principle,
one could specify $\theta_{kj}\sim\mathrm{Gamma}(a,a/\mu_{k})$ and
$\mu_{k}\sim(1-\pi_{k})\mathrm{InvGamma}(a_{0},b_{0})+\pi_{k}\delta_{\mu_{\infty}}$
to mimic our model. However, in practice, we found that this approach did not
yield the desired shrinkage effect, since the Gibbs sampler preferred to
sample from the prior slab rather than allocating the signatures to the spike.
This happened both when $a_{0}=1$ and $b_{0}=1$ were fixed and when they were
chosen according to our compressive model as $a_{0}=aJ+1$ and
$b_{0}=a\mu_{\infty}J$ with $\mu_{\infty}=0.01$. This is likely due to a
mixing issue and to the strong multimodal nature of the resulting posterior.
Instead, we found that the sampler described above worked better.

Inference for the number of signatures $K^{*}$ in PoissonCUSP can be performed
using the same adaptive Metropolis sampler as in Algorithm 2 in
Legramanti_2020. This automatically tunes the number of columns in the
factorization in a random manner, by eliminating the columns that in a given
iteration fall within the spike (that is, the signatures for which $Z_{k}\leq
k$) and potentially adding novel ones. Refer to Legramanti_2020 for a
description. We implement their algorithm with their choice of tuning
hyperparameters.

###  S5.2 Additional simulation results

We now present additional results for the simulations in Section 5. In
particular, we compare the models by assessing the root mean squared error
(RMSE) for (i) the observed counts $X$, (ii) the true mean matrix
$\Lambda^{0}=(\lambda^{0}_{ij})$, (iii) the true matrix of signatures $R^{0}$,
and (iv) the true loadings matrix $\Theta^{0}$.

First, we consider the RMSE for $X$ and $\Lambda^{0}$, calculated as
$\textsc{rmse}(X,\hat{\Lambda})=\big{(}\sum_{ij}(X_{ij}-\hat{\lambda}_{ij})^{2}/IJ\big{)}^{1/2}$
and
$\textsc{rmse}(\Lambda^{0},\hat{\Lambda})=\big{(}\sum_{ij}(\lambda^{0}_{ij}-\hat{\lambda}_{ij})^{2}/IJ\big{)}^{1/2}$,
respectively, where
$\hat{\lambda}_{ij}=\sum_{k=1}^{K^{*}}\hat{r}_{ik}\hat{\theta}_{kj}$ and
$K^{*}$ is the estimated number of signatures for each method. Figures S1 and
S2 display the comparison across all models and replicate data sets for
overdispersion $\tau=0$ and $\tau=0.15$, respectively. In the correctly
specified case ($\tau=0$), no method performs overwhelmingly better than all
the others. This is expected since all models correctly estimate the true
number of signatures in this case; see Figure 2(A). One exception is
SigProfiler, which performs noticeably worse than the others, particularly in
terms of the RMSE for $\Lambda^{0}$. The performance of SigProfiler does
improve somewhat as the sample size increases, suggesting that the algorithm
might struggle in small dimensions or on relatively low mutation counts.

The situation changes, however, in the overdispersed (negative binomial) case
with $\tau=0.15$, displayed in Figure S2. Here, SignatureAnalyzer and
PoissonCUSP obtain the lowest $\textsc{rmse}(X,\hat{\Lambda})$ across all
values of $J$, however, their corresponding
$\textsc{rmse}(\Lambda^{0},\hat{\Lambda})$ shows the opposite trend. This is a
clear indication of overfitting, which is reinforced by the fact that both
models overestimate the number of signatures; see Figure 2(A).

Figure S1: RMSE between (top) $\hat{\Lambda}$ and the count matrix $X$, and
(bottom) $\hat{\Lambda}$ and the true mean matrix $\Lambda^{0}$, for each
method across 20 replicates, when data are generated with overdispersion
$\tau=0$, i.e., the Poisson model. Figure S2: RMSE between (top)
$\hat{\Lambda}$ and the count matrix $X$, and (bottom) $\hat{\Lambda}$ and the
true mean matrix $\Lambda^{0}$, for each method across 20 replicates, when
data are generated with overdispersion be $\tau=0.15$, i.e., negative
binomial.  Figure S3: RMSE between (top) the estimated loadings $\hat{\Theta}$
and true loadings $\Theta^{0}$, and (bottom) the estimated signatures
$\hat{R}$ and true signatures $R^{0}$, over 20 replicate datasets, when data
are generated with overdispersion $\tau=0$, i.e., the Poisson model.  Figure
S4: RMSE between (top) the estimated loadings $\hat{\Theta}$ and true loadings
$\Theta^{0}$, and (bottom) the estimated signatures $\hat{R}$ and true
signatures $R^{0}$, over 20 replicate datasets, when data are generated with
overdispersion $\tau=0.15$, i.e., negative binomial.

To enable a direct comparison of performance in terms of estimating the true
signature and loadings matrices $R^{0}$ and $\Theta^{0}$, we first perform a
matching step to maximize the total pairwise cosine similarity between the
estimated $\hat{R}$ and true $R^{0}$ using the Hungarian algorithm. If the
true number of signatures and the estimated number differ, we also pad the
smaller matrix with zeros for the non-matched signatures, which penalizes
incorrect estimation of $K$. The rows of $\hat{\Theta}$ and $\Theta^{0}$ are
also permuted accordingly, and are also padded with zeros to make the
dimensions the same. We then calculate
$\textsc{rmse}(R^{0},\hat{R})=\big{(}\sum_{ik}(r^{0}_{ik}-\hat{r}_{ik})^{2}/IK\big{)}^{1/2}$
and
$\textsc{rmse}(\Theta^{0},\hat{\Theta})=\big{(}\sum_{kj}(\theta^{0}_{kj}-\hat{\theta}_{kj})^{2}/KJ\big{)}^{1/2}$,
where $\hat{r}_{ik}$ and $\hat{\theta}_{kj}$ are the estimated signatures and
loadings, obtained by averaging posterior samples in the case of the Bayesian
models. The results for $\tau=0$ and $\tau=0.15$ are shown in Figures S3 and
S4, respectively. The best performance is attained by CompNMF+cosmic. While
this model uses information from the true cosmic signatures in the prior, it
correctly filters out the signatures that are not needed. Moreover, the
improved estimation of the cosmic signatures improves the estimation of the
associated loadings. The performance of all the other models, with the
exception of SigProfiler, is virtually identical when $\tau=0$. Finally, when
$\tau=0.15$, PoissonCUSP and SignatureAnalyzer perform poorly due to
overfitting, as before.

##  S6 Additional results from the 21 breast cancer data

This section presents the complete sets of signatures inferred by each method
on the 21 breast cancer data from Section 6. The models we run are the same as
in Section 5, with only the following changes:

  1. (i)

CompNMF: $K=15$ signatures, $\varepsilon=0.01$.

  2. (ii)

CompNMF+cosmic: $K_{\mathrm{new}}=10$, _de novo_ and $K_{\mathrm{pre}}=67$
cosmic v3.4 signatures, $\varepsilon=0.01$.

  3. (iii)

PoissonCUSP: starting at $K=15$,

  4. (iv)

signeR: we set estimate_hyper = TRUE as in Rosales_2016, and let $K$ range
from $2$ to $15$. Unlike Rosales_2016, however, we do not include the
opportunity counts in the model.

  5. (v)

SignatureAnalyzer: $K=15$, selection method set to L1W.L2H (same as
Alexandrov_2020).

  6. (vi)

SigProfiler: $K$ ranging from $2$ to $15$.

Both compressive methods and PoissonCUSP are run for 12,000 iterations,
discarding the first 10,000 as burn-in.

We obtain the following results. SigProfiler estimates only three signatures
(matched to SBS2, SBS3, and SBS40a), and the cosine similarity is high only
for SBS3; see Figure S5. This is likely a consequence of the small sample
size, which can lead to the merging of two or more signatures due to
insufficient signal to distinguish them. Interestingly, SignatureAnalyzer
estimates five signatures (matched to SBS1, SBS2, SBS2, SBS3, and SBS13), but
there is duplication since one of them appears to be a combination of SBS2 and
SBS13, and ends up being matched to SBS2; see Figure S6.

We find that signeR infers five signatures (matched to SBS1, SBS2, SBS3,
SBS13, and SBS96), similar to the results in the signeR paper Rosales_2016 but
with slight differences since we do not account for the opportunity matrix in
order to provide a consistent comparison between methods; see Figure S7. In
the signeR results, SBS1 (which is sparse) appears to have been merged with a
flatter signature and thus is retrieved with lower cosine similarity. The
method also infers a signature that is matched to SBS96, but with low cosine
similarity.

CompNMF and PoissonCUSP estimate a larger number of signatures (six and seven,
respectively); see Figures S8 and S9. In both cases, SBS1, SBS2, SBS3, and
SBS13 are inferred with cosine similarities comparable to the other methods.
They also both infer SBS34 and SBS98, but these have larger credible
intervals, indicating greater uncertainty. SBS98 also has particularly low
cosine similarity in both cases, suggesting that it may be spurious.
PoissonCUSP also estimates a signature matched to SBS9, but again with
relatively high uncertainty and low cosine similarity. SBS9 has been found in
other breast cancer types Alexandrov_2020, but its current hypothesized
etiology (polymerase eta somatic hypermutation in lymphoid cells) has not been
validated.

Finally, our CompNMF+cosmic model estimates eight signatures, all with cosine
similarity near $1$, except for SBS98; see Figure S10. The estimated signature
matched to SBS98 has very high uncertainty and low cosine similarity,
suggesting it is probably a spurious match. As in the simulations, the
informative prior appears to provide significantly improved sensitivity to
detect the presence of signatures, while still allowing for departures from
the cosmic signatures.

The loadings estimated by PoissonCUSP, SignatureAnalyzer, SigProfiler, and
signeR are displayed in Figure S11 as percentages.

Figure S5: Mutational signatures inferred by SigProfiler on the 21 breast
cancer dataset. Figure S6: Mutational signatures inferred by SignatureAnalyzer
on the 21 breast cancer dataset. Figure S7: Mutational signatures inferred by
signeR on the 21 breast cancer dataset. Figure S8: Mutational signatures
inferred by PoissonCUSP on the 21 breast cancer dataset. Figure S9: Mutational
signatures inferred by CompNMF on the 21 breast cancer dataset. Figure S10:
Mutational signatures inferred by CompNMF+cosmic on the 21 breast cancer
dataset. Figure S11: Loading of each signature (as a percentage of the total
loading) for each patient, for PoissonCUSP, signeR, SignatureAnalyzer, and
SigProfiler on the 21 breast cancer dataset.

chicago references