An automated approach for determining the number of components in non-negative matrix factorization with application to mutational signature learning

NMF receives as input a matrix and a number of factors K, and outputs such that

In its most general formulation, NMF solves the following minimization problem:

for some cost function . A useful class of distance functions can be obtained by using Bregman divergence, i.e for a univariate convex smooth function φ

is a cost function. Some useful φ functions are given in table 2.1.

Table 2.1. Examples for widely used NMF variants.

Description	φ(x)
Frobenius norm		1
Kullback–Leibler divergence (KL)	xlogx
Itakura-Saito divergence (IS)	–logx

In practice, sometimes a weight matrix is given to us, e.g. reflecting entries with available data or the confidence in each cell, leading to the following extension:

In order to optimize NMF, which is known to be NP hard [13], we start from a random guess for W, H and apply some iterative update rules. For the initialization we use the random scheme used in scikit-learn; for the updates we use a scalar block coordinate descent (sBCD) [14], which we found to be the simplest and most efficient algorithm for a general cost function . We adapt both the initialization (lines 3–5) and sBCD (lines 6, 9–16) to work with a weight matrix as in (2). To check for convergence, we use a minimal relative improvement threshold ε and a maximum iteration number T (lines 7–8, 17) parameters. This is summarized in algorithm 1, which was adapted from [14]. In practice, to avoid bad local minimum, we initialize the model many times, train each for a few hundred iterations and continue training the best performing model.

In order to score different values of K, we hide some of the original matrix entries and test how well we can recover them. To this end, we set the weight matrix O to be binary, acting as a mask on the matrix values. To score the imputation of W, H we first normalize W, H, V, then average the L¹ norm over all missing values. The entire scheme is summarized in algorithm 2. In lines 2–4 we move the weights from the rows of H to the columns of W without changing the product WH. To reduce the variance caused by high/low values to impute, we normalize the rows of V and W to sum to one (lines 5–6, 7–8).

Algorithm 1 NMF(V, O, K)
1: Input:
2: Parameters: φ convex and univariate function, max iterations
3:
4:
5:
6: E = V − WH
7:
8: While do
9:
10: for do
11:
12: for do
13:
14: for do
15:
16:
17:
18: return W, H

Algorithm 2 Imputation Error (V, O, W, H)
1: Input:
2:
3:
4: W = WDH
5:
6:
7:
8:
9:
10: return

Algorithm 3 CV2K ()
1: Input:
2: Parameters: p = fraction to hide, T = number of runs for each K
3: for do
4: for do
5:
6: W, H = NMF(V, O, k)
7: S[k, t] = ImputationError(V, O, W, H)
8:
9: do
10:
11:
12: While
13: return K*

Our CV2K scheme is described in algorithm 3 and consists of two parts:

• (a)  
  (lines 3–8) Evaluate all applicable values for K using imputation with a fraction p of missing values by repeatedly running Algorithm 2 on random O matrices (line 5). Note we require O to have at least one non-zero entry in every row and column. We let K* be the value of K attaining the best median error over all runs. An example application of this criterion is shown in figure 1.
• (b)  
  (lines 9–12) Use Wilcoxon rank-sum test to compare the obtained value K* to smaller values, and reduce the value if they are not significantly different. This parsimony consideration is used to resolve cases in which the imputation error reaches a plateau, as exemplified in figure 2. Table 2 summarizes the effect of this part on the performance of our method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As we show below, the results of CV2K are quite robust to p across a wide range, yet we hypothesize that fractions between 1/N and 1/M will work best. This requires further investigation in more size-varying data sets. In order to efficiently identify K*, we use a binary search algorithm, starting from an arbitrary guess of 20. Once we find K*, we scan a small window (±3) around it to test for a better value. If there is a close K with a better median, we adopt it and repeat, otherwise we use Wilcoxon rank sum test to check the possibility of decreasing K and repeat.

To eliminate the fraction hyperparameter, we propose a variant (CV2K-x) which is described in Algorithm 4 (showing only the difference from Algorithm 3). The idea is to combine cross validation on the rows of V to learn the H matrix (signatures) and cross validation on the columns to learn the W matrix (exposures; for the latter we use SciPy's non-negative least squares [15]). Finally, reconstruction errors are summed up from each column and row folds to produce a score for a given rank. When applying this method with number of folds ranging from 5 to 30, we noticed that the results were the same or deviated by at most 1 across this range. Henceforth, we used 10-fold cross validation (F = 10).

Algorithm 4 CV2K-x ()—replaces Algorithm 3: lines 2–7
2: Parameter: F = number of sample folds
3: Randomly split V samples to F groups
4: for do
5: for do
6:
7: Vtest  = V[f]
8:
9:
10:
11: for do
12:
13:
14:
15:
16: Htest  = H[m]
17:

Two previous methods, Bi-CV [6] and scattered holdout [9], have explored the use of cross validation for learning the number of components in NMF. To demonstrate the utility of our new method, we also tested it against these two previous methods. We found that the Bi-CV method is not stable, producing results with no single local minimum. For example, on a sample simulated data set ((SA) ii]) in which the true number of components is 39, the Bi-CV method gives multiple local minima that are far from the true value (for K = 15, 64, 76, 52, 29, in this order). To test the scattered holdout approach we used Alex Williams' implementation (http://alexhwilliams.info/itsneuronalblog/2018/02/26/crossval/). We adjusted the implementation (abbreviated scatCV below) to use median values (rather than mean) to score each rank, thus decreasing the problem of outliers.

Lin et al [10] also suggests using imputation to choose K. Their approach, implemented in the NNLM package, is based on masking 30% of the input matrix entries in cross validation and picking K to minimize the reconstruction error. See https://rdrr.io/cran/NNLM/f/inst/doc/Fast-And-Versatile-NMF.pdf for additional details.

2.6.1. Simulated data

2.6.2. Real data

We analyzed breast cancer (BRCA), chronic lymphocytic leukemia (CLLE), and malignant lymphoma (MALY) mutation datasets from whole-genome sequences (WGS) from the International Cancer Genome Consortium (ICGC). We downloaded mutations from the ICGC data portal (https://docs.icgc.org/). For BRCA, we used release 22 to match Nik-Zainal et al [17], while for the other datasets we used the most recent release (release 27). These data sets had 560, 100 and 100 samples, respectively, and more than 1000 mutations for most samples.

In addition, we downloaded mutations from whole-exome-sequencing (WXS) data of 411 ovarian cancer (OV) patients from The Cancer Genome Atlas [18]. There are about 100 mutations per sample in this collection. We extracted pancreas (PANC) and prostate (PROST) mutations from WGS PCAWG data set linked from COSMIC. These data sets have 326 and 286 samples, respectively, with more then 1000 mutations for most samples.

Zoom In Zoom Out Reset image size

For evaluation, we use three figures of merit: (i) Average percent deviation is the average over , emphasizing errors in data sets with a small number of signatures. (ii) Average distance is the average over , emphasizing large deviations from the true value especially in data sets with a high number of signatures. (iii) Overall wins is the number of datasets where the predicted rank is closest to the true one. In addition, we use a signed version of the deviation, which is the average over , to examine whether we over- or under-estimate the true rank.

Our first set of results is focused on simulated data from [16] of different characteristics, accounting for different numbers of signatures and different similarity relationships between the signatures. We examine our performance on simulated data over a range of fractions p of hidden entries, as summarized in table 4, observing that the results are relatively robust across a wide range of values for this parameter and that—as a general rule—our method tends to under-estimate K, with smaller fraction values leading to closer estimations compared to larger values. We further observe that the use of the Wilcoxon rank sum step improves the performance of the algorithm, as summarized in table 2.

We compare our performance against two widely used tools for mutational signatures extraction, SigProfiler (SP) and SignatureAnalyzer (SA). We report the results of the authors of these tools on the simulated data, also published in [16]. As SigProfiler and SignatureAnalyzer miss results on few data sets (4 and 2, respectively), to aid the comparison we fill all missing values with the optimal value (marked with an asterisk in table 3). We further compare to a general model selection approach, the Akaike information criterion (AIC), as well as to the scatCV and NNLM methods reviewed above. The results are summarized in table 3 and depicted in figure 4 (left panel). They show the superiority of our approach (regardless of the choice of the hidden fraction p to use) over the existing ones.

Zoom In Zoom Out Reset image size

Next, we evaluated the different methods on six real mutation data sets described above. As the true number of components or signatures in these applications is not known, we used the number of signatures suggested in COSMIC [19] as a surrogate. For SA we report in these tests the most frequent choice across 10 runs. SP was not used in these tests due to its high running time (using the original matlab code, we could not execute the new python code on our servers) and lack of automation. The results are summarized in table 5 and depicted in figure 4 (right panel), showing that our method yields numbers of signatures that are in higher agreement with the current knowledge on the corresponding data sets as represented by COSMIC.