*K-means and cluster models for cancer signatures

We present *K-means clustering algorithm and source code by expanding statistical clustering methods applied in https://ssrn.com/abstract=2802753 to quantitative finance. *K-means is statistically deterministic without specifying initial centers, etc. We apply *K-means to extracting cancer signatures from genome data without using nonnegative matrix factorization (NMF). *K-means’ computational cost is a fraction of NMF’s. Using 1389 published samples for 14 cancer types, we find that 3 cancers (liver cancer, lung cancer and renal cell carcinoma) stand out and do not have cluster-like structures. Two clusters have especially high within-cluster correlations with 11 other cancers indicating common underlying structures. Our approach opens a novel avenue for studying such structures. *K-means is universal and can be applied in other fields. We discuss some potential applications in quantitative finance.


Introduction and summary
Every time we can learn something new about cancer, the motivation goes without saying. Cancer is different. Unlike other diseases, it is not caused by "mechanical" breakdowns, biochemical imbalances, etc. Instead, cancer occurs at the DNA level via somatic alterations in the genome structure. A common type of somatic mutations found in cancer is due to single nucleotide variations (SNVs) or alterations to single bases in the genome, which accumulate through the lifespan of the cancer via imperfect DNA replication during cell division or spontaneous cytosine deamination [1,2], or due to exposures to chemical insults or ultraviolet radiation [3,4], etc. These mutational processes leave a footprint in the cancer genome characterized by distinctive alteration patterns or mutational signatures.
If we can identify all underlying signatures, this could greatly facilitate progress in understanding the origins of cancer and its development. Therapeutically, if there are common underlying structures across different cancer types, then a therapeutic for one cancer type might be applicable to other cancers, which would be a great news. 2 However, it all boils down to the question of usefulness, i.e., is there a small enough number of cancer signatures underlying all (100+) known cancer types, or is this number too large to be meaningful or useful? Indeed, there are only 96 SNVs, 3 so we cannot have more than 96 signatures. 4 Even if the number of true underlying signatures is, say, of order 50, it is unclear whether they would be useful, especially within practical applications. On the other hand, if there are only a dozen or so underlying signatures, then we could hope for an order of magnitude simplification.
To identify mutational signatures, one analyzes SNV patterns in a cohort of DNA sequenced whole cancer genomes. The data is organized into a matrix G is , where the rows correspond to the N = 96 mutation categories, the columns correspond to d samples, and each element is a nonnegative occurrence count of a given mutation category in a given sample. Currently, the commonly accepted method for extracting cancer signatures from G is [5] is via nonnegative matrix factorization (NMF) [6,7]. Under NMF the matrix G is approximated via G ≈ W H, where W iA is an N × K matrix, H As is a K × d matrix, and both W and H are nonnegative. The appeal of NMF is its biologic interpretation whereby the K columns of the matrix W are interpreted as the weights with which the K cancer signatures contribute into the N = 96 mutation categories, and the columns of the matrix H are interpreted as the exposures to the K signatures in each sample. The price to pay for this is that NMF, which is an iterative procedure, is computationally costly and depending on the number of samples d it can take days or even weeks to run it. Furthermore, it does not automatically fix the number of signatures K, which must be either guessed or obtained via trial and error, thereby further adding to the computational cost. 5 Some of the aforesaid issues were recently addressed in [8], to wit: (i) by aggregating samples by cancer types, we can greatly improve stability and reduce the number of signatures; 6 (ii) by identifying and factoring out the somatic mutational noise, or the "overall" mode (this is the "de-noising" procedure of [8]), we can further greatly improve stability and, as a bonus, reduce computational cost; and (iii) the number of signatures can be fixed borrowing the methods from statistical risk models [9] in quantitative finance, by computing the effective rank (or eRank) [10] for the correlation matrix Ψ ij calculated across cancer types or samples (see below). All this yields substantial improvements [8].
In this paper we push this program to yet another level. The basic idea here is quite simple (but, as it turns out, nontrivial to implementsee below). We wish to apply clustering techniques to the problem of extracting cancer signatures. In fact, we argue in Section 2 that NMF is, to a degree, "clustering in disguise". This is for two main reasons. The prosaic reason is that NMF, being a nondeterministic algorithm, requires averaging over many local optima it produces. However, each run generally produces a weights matrix W iA with columns (i.e., signatures) not aligned with those in other runs. Aligning or matching the signatures across different runs (before averaging over them) is typically achieved via nondeterministic clustering such as k-means. So, not only is clustering utilized at some layer, the result, even after averaging, generally is both noisy 7 and nondeterministic! I.e., if this computationally costly procedure (which includes averaging) is run again and again on the same data, generally it will yield different looking cancer signatures every time!
The second, not-so-prosaic reason is that, while NMF generically does not produce exactly null weights, it does produce low weights, such that they are within error bars. For all practical purposes we might as well set such weights to zero. NMF requires nonnegative weights. However, we could as reasonably require that the weights should be, say, outside error bars (e.g., above one standard deviationthis would render the algorithm highly recursive and potentially unstable or computationally too costly) or above some minimum threshold (which would still further complicated as-is complicated NMF), or else the noncompliant weights are set to zero. As we increase this minimum threshold, the matrix W iA will start to have more and more zeros. It may not exactly have a binary cluster-like structure, but it may at least have some substructures that are cluster-like. It then begs the question: are there cluster-like (sub)structures present in W iA or, generally, in cancer signatures?
To answer this question, we can apply clustering methods directly to the matrix G is , or, more, precisely, to its de-noised version ′ G is (see below) [8]. The naïve, brute-force approach where one would simply cluster G is or ′ G is does not work for a variety of reasons, some being more nontrivial or subtle than others. Thus, e.g., as discussed in [8], the counts G is have skewed, long-tailed distributions and one should work with log-counts, or, more precisely, their de-noised versions. This applies to clustering as well. Further, following a discussion in [11] in the context of quantitative trading, it would be suboptimal to cluster denoised log-counts. Instead, it pays to cluster their normalized variants (see Section 2 hereof). However, taking care of such subtleties does not alleviate one big problem: nondeterminism! 8 If we run a vanilla nondeterministic algorithm such as k-means on the data however massaged with whatever bells and whistles, we will get random-looking disparate results every time we run k-means with no stability in sight. We need to address nondeterminism! Our solution to the problem is what we term *K-means. The idea behind *K-means, which essentially achieves determinism statistically, is simple. Suppose we have an N × d matrix X is , i.e., we have N dvectors X i . If we run k-means with the input number of clusters K but initially unspecified centers, every run will generally produce a new local optimum. *K-means reduces and in fact essentially eliminates this indeterminism via two levels. At level 1 it takes clusterings obtained via M independent runs or samplings. Each sampling produces a binary N × K matrix Ω iA , whose element equals 1 if X i belongs to the cluster labeled by A, and 0 otherwise. The aggregation algorithm and the source code therefor are given in [11]. This aggregationfor the same reasons as in NMF (see above)involves aligning clusters across the M runs, which is achieved via k-means, and so the result is nondeterministic. However, by aggregating a large number M of samplings, the degree of nondeterminism is greatly reduced. The "catch" is that sometimes this aggregation yields a clustering with K′ < K clusters, but this does not pose an issue. Thus, at level 2, we take a large number P of such aggregations (each based on M samplings). The occurrence counts of aggregated clusterings are not uniform but typically have a (sharply) peaked distribution around a few (or manageable) number of aggregated clusterings. So this way we can pinpoint the "ultimate" clustering, which is simply the aggregated clustering with the highest occurrence count. This is the gist of *K-means and it works well for genome data.
So, we apply *K-mean to the same genome data as in [8] consisting of 1389 (published) samples across 14 cancer types (see below). Our target number of clusters is 7, which was obtained in [8] using the eRank based algorithm (see above). We aggregated 1000 samplings into clusterings, and we constructed 150,000 such aggregated clusterings (i.e., we ran 150 million k-means instances). We indeed found the "ultimate" clustering with 7 clusters. Once the clustering is fixed, it turns out that within-cluster weights can be computed via linear regressions (with some bells and whistles) and the weights are automatically positive. That is, we do not need NMF at all! Once we have clusters and weights, we can study reconstruction accuracy and withincluster correlations between the underlying data and the fitted data that the cluster model produces.
We find that clustering works well for 10 out the 14 cancer types we study. The cancer types for which clustering does not appear to work all that well are Liver Cancer, Lung Cancer, and Renal Cell Carcinoma. Also, above 80% within-cluster correlations arise for 5 out of 7 clusters. Furthermore, remarkably, one cluster has high within-cluster correlations for 9 cancer types, and another cluster for 6 cancer types. These appear to be the leading clusters. Together they have high withincluster correlations in 11 out of 14 cancer types. So what does all this mean?
Additional insight is provided by looking at the within-cluster correlations between signatures Sig1 through Sig7 extracted in [8] and our clusters. High within-cluster correlations arise for Sig1, Sig2, Sig4 and Sig7, which are precisely the signatures with "peaks" (or "spikes" -"tall mountain landscapes"), whereas Sig3, Sig5 and Sig6 do not have such "peaks" ("flat" or "rolling hills landscapes"); see Figs. 14 through 20 of [8]. The latter 3 signatures simply do not have cluster-like structures. Looking at Fig. 21  In Section 2 we (i) discuss why applying clustering algorithms to extracting cancer signatures makes sense, (ii) argue that NMF, to a degree, is "clustering in disguise", and (iii) give the machinery for building cluster models via *K-means, including various details such as what to cluster, how to fix the number of clusters, etc. In Section 3 we discuss (i) cancer genome data we use, (ii) our application of *K-means to it, and (iii) the interpretation of our empirical results. Section 4 contains some concluding remarks, including a discussion of potential applications of *K-means in quantitative finance, where we outline some concrete problems where *K-means can be useful. Appendix A contains R source code for *K-means and cluster models.

Cluster models
The chief objective of this paper is to introduce a novel approach to identifying cancer signatures using clustering methods. In fact, as we discuss below in detail, our approach is more than just clustering. Indeed, it is evident from the get-go that blindly using nondeterministic clustering algorithms, 9 which typically produce (unmanageably) large numbers of local optima, would introduce great variability into the resultant cancer signatures. 10 On the other hand, deterministic algorithms such as agglomerative hierarchical clustering 11 typically are (substantially) slower and require essentially "guessing" the initial clustering, 12 which in practical applications 13 can often turn out to be suboptimal. So, both to motivate and explain our new approach employing clustering methods, we firstso to speak -"break down" the NMF approach and argue that it is in fact a clustering method in disguise!

"Breaking down" NMF
The current "lore"the commonly accepted method for extracting K cancer signatures from the occurrence counts matrix G is (see above) [5] is via nonnegative matrix factorization (NMF) [6,7]. Under NMF the matrix G is approximated via G ≈ W H, where W iA is an N × K matrix of weights, H As is a K × d matrix of exposures, and both W and H are nonnegative. However, not only is the number of signatures K not fixed via NMF (and must be either guessed or obtained via trial and error), NMF too is a nondeterministic algorithm and typically produces a large number of local optima. So, in practice one has no choice but to execute a large number N S of NMF runswhich we refer to as samplingsand then somehow extract cancer signatures from these samplings. Absent a guess for what K should be, one executes N S samplings for a range of values of K (say, K min ≤ K ≤ K max , where K min and K max are basically guessed based on some reasonable intuitive considerations), for each K extracts cancer signatures (see below), and then picks K and the corresponding signatures with the best overall fit into the underlying matrix G. For a given K, different samplings generally produce different weights matrices W. So, to extract a single matrix W for each value of K one averages over the samplings. However, before averaging, one must match the K cancer signatures across different samplingsindeed, in a given sampling X the columns in the matrix W iA are not necessarily aligned with the columns in the matrix W iA in a different sampling Y. To align the columns in the matrices W across the N S samplings, once often uses a clustering algorithm such as k-means. However, since k-means is nondeterministic, such alignment of the W columns is not guaranteed toand in fact does notproduce a unique answer. Here one can try to run multiple samplings of k-means for this alignment and aggregate them, albeit such aggregation itself would require another level of alignment (with its own nondeterministic clustering such as k-means). 14 And one can do this ad infinitum. In practice, one must break the chain at some level of alignment, either ad hoc (essentially by heuristically observing sufficient stability and "convergence") or via using a deterministic algorithm (see footnote 14 ). Either way, invariably all this introduces (overtly or covertly) systematic and statistical errors into the resultant cancer signatures and often it is unclear if they are meaningful without invoking some kind empirical biologic "experience" or "intuition" (often based on already well-known effects of, e.g., exposure to various well-understood carcinogens such as tobacco, ultraviolet radiation, aflatoxin, etc.). At the end of the day it all boils down to how usefulor predictivethe resultant method of extracting cancer signatures is, including signature stability. With NMF, the answer is not at all evident…

Clustering in disguise?
So, in practice, under the hood, NMF already uses clustering methods. However, it goes deeper than that. While NMF generically does not produce vanishing weights for a given signature, some weights are (much) smaller than others. E.g., often one has several "peaks" with high concentration of weights, with the rest of the mutation categories having relatively low weights. In fact, many weights can even be within the (statistical plus systematic) error bars. 15 Such weights can for all practical purposes be set to zero. In fact, we can take this further and ask whether proliferation of low weights adds any explanatory power. One way to address this is to run NMF with an additional constraint that the weights (obtained via averagingsee above) should be higher than either (i) some multiple of the corresponding error bars 16 or (ii) some preset fixed minimum weight. This certainly sounds reasonable, so why is this not done in practice? A prosaic answer appears to be that this would complicate the already nontrivial NMF algorithm even further, require additional coding and computation resources, etc. However, arguendo, let us assume that we require, say, that the weights be higher than a preset fixed minimum weight w min or else the weights are set to zero. As we increase w min , the so-modified NMF would produce more and more zeros. This does not mean that the resulting matrix W iA would have a binary cluster structure, i.e., that = W w δ i G i A iA ( ), , where δ AB is a Kronecker delta and G : {1, …, N} ↦ {1, …, K} is a map from N = 96 mutation categories to K clusters. Put another way, this does not mean that in the resulting matrix W iA for a given i (i.e., mutation category) we would have a nonzero element for one and only one value of A (i.e., signature). However, as we gradually increase w min , generally the matrix W iA is expected to look more and more like having a binary cluster structure, albeit with some "overlapping" signatures (i.e., such that in a given pair of signatures there are nonzero weights for one or more mutations). We can achieve a binary structure via a number of ways. Thus, a rudimentary algorithm would be to take the matrix W iA (equally successfully before or after achieving some zeros in it via nonzero w min ) and for a given value of i set all weights W iA to zero except in the signature A for which W iA = max(W iA |A = 1, …, K). Note that this might result in some empty signatures (clusters), i.e., signatures with W iA = 0 for all values of i. This can be dealt with by (i) ether simply dropping such signatures altogether and having fewer K′ < K signatures (binary clusters) at the end, or (ii) augmenting the algorithm to avoid empty clusters, which can be done in a number of ways we will not delve into here. The bottom line is that NMF essentially can be made into a clustering algorithm by reasonably modifying it, including via getting rid of ubiquitous and not-too-informative low weights. However, the downside would be an even more contrived algorithm, so this is not what we are suggesting here. Instead, we are observing that clustering is already intertwined in NMF and the question is whether we can simplify things by employing clustering methods directly.

Making clustering work
Happily, the answer is yes. Not only can we have much simpler and apparently more stable clustering algorithms, but they are also computationally much less costly than NMF. As mentioned above, the biggest issue with using popular nondeterministic clustering algorithms such as k-means 17 is that they produce a large number of local optima. For definiteness in the remainder of this paper we will focus on kmeans, albeit the methods described herein are general and can be applied to other such algorithms. Fortunately, this very issue has already been addressed in [11] in the context of constructing statistical industry classifications (i.e., clustering models for stocks) for quantitative trading, so here we simply borrow therefrom and further expand and adapt that approach to cancer signatures.

K-means
A popular clustering algorithm is k-means [12][13][14][15][16][17][18]. The basic idea behind k-means is to partition N observations into K clusters such that each observation belongs to the cluster with the nearest mean. Each of the N observations is actually a d-vector, so we have an N × d matrix X is , i = 1, …, N, s = 1, …, d. Let C a be the K clusters, C a = {i|i ∈ C a }, a = 1, …, K. Then k-means attempts to minimize 18 are the cluster centers (i.e., cross-sectional means), 19 and n a = |C a | is the number of elements in the cluster C a . In (1) the measure of "closeness" is chosen to be the Euclidean distance between points in R d , albeit other measures are possible.
One "drawback" of k-means is that it is not a deterministic algorithm. Generically, there are copious local minima of g in (1) and the algorithm only guarantees that it will converge to a local minimum, not the global one. Being an iterative algorithm, unless the initial centers are preset, k-means starts with a random set of the centers Y as at the initial iteration and converges to a different local minimum in each run. There is no magic bullet here: in practical applications, typically, trying to "guess" the initial centers is not any easier than "guessing" where, e.g., the global minimum is. So, what is one to do? One possibility is to simply live with the fact that every run produces a different answer. In fact, this is acceptable in many applications. However, in the context of extracting cancer signatures this would result in an exercise in futility. We need a way to eliminate or greatly reduce indeterminism.

Aggregating clusterings
The idea is simple. What if we aggregate different clusterings from multiple runswhich we refer to as samplingsinto one? The question is how. Suppose we have M runs (M ≫ 1). Each run produces a clustering with K clusters. Let ). Here we are assuming that somehow we know how to properly order (i.e., align) the K clusters from each run. This is a nontrivial assumption, which we will come back to momentarily. However, assuming, for a second, that we know how to do this, we can aggregate the binary matrices Ω r ia into a single matrix ia . Now, this matrix does not look like a binary clustering matrix. Instead, it is a matrix of occurrence counts, i.e., it counts how many times a given mutation was assigned to a given cluster in the process of M samplings. What we need to construct is a map G such that one and only one mutation belongs to each of the K clusters. The simplest criterion is to map a given mutation to the cluster in which ∼ Ω ia is maximal, i.e., where said mutation occurs most frequently. A caveat is that there may be more than one such clusters. A simple criterion to resolve such an ambiguity is to assign said mutation to the cluster with most cumulative occurrences (i.e., we assign said mutation to the cluster with the largest ia ). Further, in the unlikely event that there is still an ambiguity, we can try to do more complicated things, or we can simply assign such a mutation to the cluster with the lowest value of the index atypically, there is so much noise in the system that dwelling on such minutiae simply does not pay off.
However, we still need to tie up a loose end, to wit, our assumption that the clusters from different runs were somehow all aligned. In practice each run produces K clusters, but (i) they are not the same clusters and there is no foolproof way of mapping them, especially when we have a large number of runs; and (ii) even if the clusters were the same or similar, they would not be ordered, i.e., the clusters from one run generally would be in a different order than the clusters from another run.
So, we need a way to "match" clusters from different samplings. Again, there is no magic bullet here either. We can do a lot of complicated and contrived things with not much to show for it at the end. A simple pragmatic solution is to use k-means to align the clusters from different runs. Each run labeled by r = 1, …, M, among other things, produces a set of cluster centers Y r as . We can "bootstrap" them by row into a (KM) takes values 17 Which are preferred over deterministic ones for the reasons discussed above. 18 Below we will discuss what X is should be for cancer signatures. 19 Throughout this paper "cross-sectional" refers to "over the index i". 20 Note that here the superscript r in Ω r ia , G r (i) and N a r (see below) is an index, not a power.
Z. Kakushadze, W. Yu Biomolecular Detection and Quantification 13 (2017) 7-31 = … a 1, , (KM) ͠ . We can now cluster ∼ Y as ͠ into K clusters via k-means. This will map each value of a͠ to {1, …, K} thereby mapping the K clusters from each of the M runs to {1, …, K}. So, this way we can align all clusters. The "catch" is that there is no guarantee that each of the K clusters from each of the M runs will be uniquely mapped to one value in {1, …, K}, i.e., we may have some empty clusters at the end of the day. However, this is fine, we can simply drop such empty clusters and aggregate (via the above procedure) the smaller number of K′ < K clusters. I.e., at the end we will end up with a clustering with K′ clusters, which might be fewer than the target number of clusters K. This is not necessarily a bad thing. The dropped clusters might have been redundant in the first place. Another evident "catch" is that even the number of resulting clusters K′ is not deterministic. If we run this algorithm multiple times, we will get varying values of K′. Malicious circle?

Fixing the "ultimate" clustering
Not really! There is one other trick up our sleeves we can use to fix the "ultimate" clustering thereby rendering our approach essentially deterministic. The idea above is to aggregate a large enough number M of samplings. Each aggregation produces a clustering with some K′ ≤ K clusters, and this K′ varies from aggregation to aggregation. However, what if we take a large number P of aggregations (each based on M samplings)? Typically there will be a relatively large number of different clusterings we get this way. However, assuming some degree of stability in the data, this number is much smaller than the number of a priori different local minima we would obtain by running the vanilla kmeans algorithm. What is even better, the occurrence counts of aggregated clusterings are not uniform but typically have a (sharply) peaked distribution around a few (or manageable) number of aggregated clusterings. In fact, as we will see below, in our empirical genome data we are able to pinpoint the "ultimate" clustering! So, to recap, what we have done here is this. There are myriad clusterings we can get via vanilla k-means with little to no guidance as to which one to pick. 21 We have reduced this proliferation by aggregating a large number of such clusterings into our aggregated clusterings. We then further zoom onto a few or even a unique clustering we consider to be the likely "ultimate" clustering by examining the occurrence counts of such aggregated clusterings, which turns out to have a (sharply) peaked distribution. Since vanilla k-means is a relatively fast-converging algorithm, each aggregation is not computationally taxing and running a large number of aggregations is nowhere as time consuming as running a similar number (or even a fraction thereof) of NMF computations (see below).

What to cluster?
So, now that we know how to make clustering work, we need to decide what to cluster, i.e., what to take as our matrix X is in (1). The naïve choice X is = G is is suboptimal for multiple reasons (as discussed in [8]).
First, the elements of the matrix G is are populated by nonnegative occurrence counts. Nonnegative quantities with large numbers of samples tend to have skewed distributions with long tails at higher values. I.e., such distributions are not normal but (in many cases) roughly log-normal. One simple way to deal with this is to identify X is with a (natural) logarithm of G is (instead of G is itself). A minor hiccup here is that some elements of G is can be 0. We can do a lot of complicated and even convoluted things to deal with this issue. Here, as in [8], we will follow a pragmatic approach and do something simple insteadthere is so much noise in the data that doing convoluted things simply does not pay off. So, as the first cut, we can take This takes care of the G is = 0 cases; for G is ≫ 1 we have R is ≈ ln(G is ), as desired. Second, the detailed empirical analysis of [8] uncovered what is termed therein the "overall" mode 22 unequivocally present in the occurrence count data. This "overall" mode is interpreted as somatic mutational noise unrelated to (and in fact obscuring) the true underlying cancer signatures and must therefore be factored out somehow.
Here is a simple way to understand the "overall" mode. Let the correlation matrix Ψ ij = Cor(X is , X js ), where Cor(·, ·) is serial correlation. 23 I.e., ii are variances, and the serial covariance matrix 24 ij between different mutation categories is nonzero and is in fact high for most cancer types we study. This is the aforementioned somatic mutational noise that must be factored out. If we aggregate samples by cancer types (see below) and compute the correlation matrix Ψ ij for the so-aggregated data (across the n = 14 cancer types we studysee below), 25 the average correlation ρ is over whopping 96%. Another way of thinking about this is that the occurrence counts in different samples (or cancer types, if we aggregate samples by cancer types) are not normalized uniformly across all samples (cancer types). Therefore, running NMF, a clustering or any other signature-extraction algorithm on the vanilla matrix G is (or its "log" X is defined in (3)) would amount to mixing apples with oranges thereby obscuring the true underlying cancer signatures.
Following [8], factoring out the "overall" mode (or "de-noising" the matrix G is ) therefore most simply amount to cross-sectional (i.e., across the 96 mutation categories) demeaning of the matrix X is . I.e., instead of X is we use ′ X is , which is obtained from X is by demeaning its columns: 26 We should note that using ′ X is instead of X is in (1) does not affect clustering. Indeed, g in (1) is invariant under the transformations of the form X is → X is + Δ s , where Δ s is an arbitrary d-vector, as thereunder we also have Y as → Y as + Δ s , so X is − Y as is unchanged. In fact, this is good: this means that de-noising does not introduce any additional errors into clustering itself. However, the actual weights in the matrix W iA are affected by de-noising. We discuss the algorithm for fixing W iA below. However, we need one more ingredient before we get to determining the weights, and with this additional ingredient de-noising does affect clustering.

Normalizing log-counts
As was discussed in [11], clustering X is (or equivalently ′ X is ) would be suboptimal. 27 The issue is this. Let ′ σ i be serial standard deviations, references therein) corresponding to the overall movement of the broad market, which affects all stocks (to varying degrees)cash inflow (outflow) into (from) the market tends to push stock prices higher (lower). This is the market risk factor, and to mitigate it one can, e.g., hold a dollar-neutral portfolio of stocks (i.e., the same dollar holdings for long and short positions). 23 Throughout this paper "serial" refers to "over the index s". 24 The overall normalization of C ij , i.e., d − 1 (unbiased estimate) vs. d (maximum likelihood estimate) in the denominator in the definition of C ij in (4), is immaterial for our purposes here. 25 So, in this case d = n = 14 in (4). 26 For the reasons discussed above, we should demean X is , not G is . 27 More precisely, the discussion of [11] is in the financial context, to wit, quantitative trading, which has its own nuances (see below). However, some of that discussion is quite is , where, as above, Cov(·, ·) is serial covariance. Here we assume that samples are aggregated by cancer types, so s = 1, …, d with d = n = 14. Now, ′ σ i are not cross-sectionally uniform and vary substantially across mutation categories. The density of ′ σ i is depicted in Fig. 1  X is , this variability in ′ σ i will not be accounted for.
A simple solution is to cluster normalized demeaned log-counts This way we factor out the nonuniform (and skewed) standard deviation out of the log-counts. Note that now denoising does make a difference in clustering. Indeed, if we use (1) and (2), the quantity g (and also clusterings) will be different.

Fixing cluster number
Now that we know what to cluster (to wit, ′ ∼ X is ) and how to get to the "unique" clustering, we need to figure out how to fix the (target) number of clusters K, which is one of the inputs in our algorithm above. 29 In [8] it was argued that in the context of cancer signatures their number can be fixed by building a statistical factor model [9], i.e., the number of signatures is simply the number of statistical factors. 30 So, by the same token, here we identify the (target) number of clusters in our clustering algorithm with the number of statistical factors fixed via the method of [9].

Effective rank
So, following [9,8], we set 31 Here eRank(Z) is the effective rank [10] of a symmetric semi-positivedefinite (which suffices for our purposes here) matrix Z. It is defined as (8)(9) where λ (a) are the L positive eigenvalues of Z, and H has the meaning of the (Shannon a.k.a. spectral) entropy [34,35]. Let us emphasize that in (6) the matrix Ψ ij is computed based on the demeaned log-counts 32 ′ X is . The meaning of eRank(Ψ ij ) is that it is a measure of the effective dimensionality of the matrix Ψ ij , which is not necessarily the same as the number L of its positive eigenvalues, but often is lower. This is due to the fact that many d-vectors ′ X is can be serially highly correlated (which manifests itself by a large gap in the eigenvalues) thereby further reducing the effective dimensionality of the correlation matrix.

How to compute weights?
The one remaining thing to accomplish is to figure out how to compute the weights W iA . Happily, in the context of clustering we have significant simplifications compared with NMF and computing the weights becomes remarkably simple once we fix the clustering, i.e., the matrix Ω iA = δ G(i),A (or, equivalently, the map G : {i} ↦ {A}, i = 1, …, N, A = 1, …, K, where for the notational convenience we use K to denote the number of clusters in the "ultimate" clusteringsee above). Just as in NMF, we wish to approximate the matrix G is via a product of the weights matrix W iA and the exposure matrix H As , both of which must be nonnegative. More precisely, since we must remove the "overall" mode, i.e., de-noise the matrix G is , following [8], instead of G is we will approximate the re-exponentiated demeaned log-counts matrix ′ X is : We can include an overall normalization by taking s is is (recall that X s is the vector of column means of X issee Eq. (5)), etc., to make it look more like the original matrix G is ; however, this does not affect the extracted signatures. 33 Also, technically speaking, after re-exponentiating we should "subtract" the extra 1 we added in the definition (3) (assuming we include one of the aforesaid overall normalizations). However, the inherent noise in the data makes this a moot point.
So, we wish to approximate ′ G is via a product W H. However, with clustering we have = W w δ i G i A iA ( ), , i.e., we have a block (cluster) structure where for a given value of A all W iA are zero except for i ∈ J Fig. 1. Horizontal axis: serial standard deviation ′ σ i for N = 96 mutation categories (i = 1, …, N) of cross-sectionally demeaned log-counts ′ X is across n = 14 cancer types (for samples aggregated by cancer types, so s = 1, …, d, d = n). Vertical axis: density using R function density(). See Section 2.4.1 for details.
(footnote continued) general and can be adapted to a wide variety of applications. 28 Qu. = Quartile, SD = Standard Deviation, MAD = Mean Absolute Deviation. 29 A variety of methods for fixing the number of clusters have been discussed in other contexts, e.g., [22][23][24][25][26][27][28][29]. 30 In the financial context, these are known as statistical risk models [9]. For a discussion and literature on multifactor risk models, see, e.g., [30,31] and references therein. For prior works on fixing the number of statistical risk factors, see, e.g., [32,33]. 31 Here Round(·) can be replaced by floor(·) = ⌊·⌋. 32 Note that using normalized demeaned log-counts ′ ∼ X is gives the same Ψ ij . 33 This is because each column of W, being weights, is normalized to add up to 1.

Z. Kakushadze, W. Yu
Biomolecular Detection and Quantification 13 (2017)  (A) = {j|G(j) = A}, i.e., for the mutation categories labeled by i that belong to the cluster labeled by A. Therefore, our matrix factorization of G is into a product W H now simplifies into a set of K independent factorizations as follows: So, there is no need to run NMF anymore! Indeed, if we can somehow fix H As for a given cluster, then within this cluster we can determine the corresponding weights w i (i ∈ J(A)) via a serial linear regression: where ε is are the regression residuals. I.e., for each A ∈ {1, …, K}, we 2 . Furthermore, if H As are nonnegative, then the weights w i are automatically nonnegative as they are given by: Now, we wish these weights to be normalized: This can always be achieved by rescaling H As . Alternatively, we can pick H As without worrying about the normalization, compute w i via (13), rescale them so that they satisfy (14), and simultaneously accordingly rescale H As . Mission accomplished!

Fixing exposures
Well, almost… We still need to figure out how to fix the exposures H As . The simplest way to do this is to note that we can use the matrix Ω iA = δ G(i),A to swap the index i in ′ G is by the index A, i.e., we can take That is, up to the normalization constants η͠ A (which are fixed via (14)) we simply take cross-sectional means of ′ G is in each cluster. (Recall that n A = J(A).) The so-defined H As are automatically positive as all ′ G is are positive. Therefore, w i defined via (13) are also all positive. This is a good newsvanishing w i would amount to an incomplete weights matrix W iA (i.e., some mutations would belong to no cluster).
So, why does (15) make sense? Looking at (12), we can observe that, if the residuals ε is cross-sectionally, within each cluster labeled by A, are random, then we expect that ∑ i∈J(A) ε is ≈ 0. If we had an exact equality here, then we would have (15) with η A = 1 (i.e., = η n ͠ A A ) assuming the normalization (14). In practice, the residuals ε is are not exactly "random". First, the number n A of mutation categories in each cluster is not large. Second, as mentioned above, there is variability in serial standard deviations across mutation types. This leads us to consider variations.

A variation
Above we argued that it makes sense to cluster normalized demeaned log-counts ′ = ′ ′ ∼ X X σ / i is is due to the cross-sectional variability (and skewness) in the serial standard deviations ′ σ i . We may worry about similar effects in ′ G is when computing H As and w i as we did above. This can be mitigated by using normalized quantities is are serial variances. That is, we can define 35 where ν A = ∑ i∈J(A) 1/ω i . So, 1/ω i are the weights in the averages over the clusters.

Another variation
Here one may wonder, considering the skewed roughly log-normal distribution of G is and henceforth ′ G is , would it make sense to relate the exposures to within-cluster cross-sectional averages of demeaned logcounts ′ X is as opposed to those of ′ G is ? This is easily achieved. Thus, we can define (this ensures positivity of H As ): Exponentiating we get I.e., instead of an arithmetic average as in (15), here we have a geometric average. As above, here too we can introduce nontrivial weights. Note that the form of (17) is the same as (13), it is only H As that is affected by the weights. So, we can introduce the weights in the geometric means as follows: is . Thus, we have: So, the weights are the exponents ′ μ σ 1/ A i . Other variations are also possible.

Implementation
We are now ready to discuss an actual implementation of the above algorithm, much of the R code for which is already provided in [8,11]. The R source code is given in Appendix A hereof.

Data summary
In our empirical analysis below we use the same genome data (from published samples only) as in [8]. This data is summarized in Table S1 (borrowed from [8]), which gives total counts, number of samples and the data sources, which are as follows: A1 = [36] [52]. Sample IDs with the corresponding publication sources are given in Appendix A of [8]. In our analysis below we aggregate samples by the 14 cancer types. The resulting data is in Tables S2 and S3. For tables and figures labeled S★ see Supplementary Materials (see Appendix C for a web link).

Structure of data
The underlying data consists of a matrixcall it G iswhose elements are occurrence counts of mutation types labeled by i = 1, …, N = 96 in samples labeled by s = 1, …, d. More precisely, we can work 34 The superscript T denotes matrix transposition. 35 I.e., here we assume that ε is /ω i are approximately random in (12). with one matrix G is which combines data from different cancer types; or, alternatively, we may choose to work with individual matrices [G(α)] is , where: α = 1, …, n labels n different cancer types; as before, i = 1, …, N = 96; and s = 1, …, d(α). Here d(α) is the number of samples for the cancer type labeled by α. The combined matrix G is is obtained simply by appending (i.e., bootstrapping) the matrices [G(α)] is together column-wise. In the case of the data we use here (see above), this "big matrix" turns out to have 1389 columns. Generally, individual matrices [G(α)] is and, thereby, the "big matrix", contain a lot of noise. For some cancer types we can have a relatively small number of samples. We can also have "sparsely populated" data, i.e., with many zeros for some mutation categories. As mentioned above, different samples are not necessarily uniformly normalized. Etc. The bottom line is that the data is noisy. Furthermore, intuitively it is clear that the larger the matrix we work with, statistically the more "signatures" (or clusters) we should expect to get with any reasonable algorithm. However, as mentioned above, a large number of signatures would be essentially useless and defy the whole purpose of extracting them in the first placewe have 96 mutation categories, so it is clear that the number of signatures cannot be more than 96! If we end up with, say, 50+ signatures, what new or useful does this tell us about the underlying cancers? The answer is likely nothing other than that most cancers have not much in common with each other, which would be a disappointing result from the perspective of therapeutic applications. To mitigate the aforementioned issues, at least to a certain extent, following [8], we can aggregate samples by cancer types. This way we get an N × n matrix, which we also refer to as G is , where the index s = 1, …, d now takes d = n values corresponding to the cancer types. In the data we use n = 14, the aggregated matrix G is is much less noisy than the "big matrix", and we are ready to apply the above machinery to it.

Genome data results
The 96 × 14 matrix G is given in Tables S2 and S3 is what we pass into the function bio.cl.sigs() in Appendix A as the input matrix x. We use: iter.max = 100 (this is the maximum number of iterations used in the built-in R function kmeans()we note that there was not a single instance in our 150 million runs of kmeans() where more iterations were required); 36 num.try = 1000 (this is the number of individual k-means samplings we aggregate every time); and num.runs = 150000 (which is the number of aggregated clusterings we use to determine the "ultimate"that is, the most frequently occurringclustering). So, we ran k-means 150 million times. More precisely, we ran 15 batches with num.runs = 10000 as a sanity check, to make sure that the final result based on 150,000 aggregated clusterings was consistent with the results based on smaller batches, i.e., that it was insample stable. 37 Based on Table S4, we identify Clustering-A as the "ultimate" clustering (cf. Clustering-B/C/D).
We give the weights for Clustering-A, Clustering-B, Clustering-C and Clustering-D using unnormalized and normalized regressions with exposures computed based on arithmetic averages (see Section 2.6) in Tables 1, 2, S5-S10, and Figs. 2 through Fig. 15 and S1 through S40. We give the weights for Clustering-A using unnormalized and normalized regressions with exposures computed based on geometric averages (see Section 2.6) in Tables 3, 4, and Figs. S41 through S54. The actual mutation categories in each cluster for a given clustering can be read off the aforesaid tables with the weights (the mutation categories with nonzero weights belong to a given cluster), or from the horizontal axis labels in the aforesaid figures. It is evident that Clustering-A, Clustering-B, Clustering-C and Clustering-D are essentially variations of each other (Clustering-D has only 6 clusters, while the other 3 have 7 clusters).

Reconstruction and correlations
So, based on genome data, we have constructed clusterings and weights. Do they work? I.e., do they reconstruct the input data well? It is evident from the get-go that the answer to this question may not be binary in the sense that for some cancer types we might have a nice clustering structure, while for others we may not. The aim of the following exercise is to sort this all out. Here come the correlations…

Within-cluster correlations
We have our de-noised 38 matrix ′ G is . We are approximating this matrix via the following factorized matrix: We can now compute an n × K matrix Θ sA of within-cluster cross-sectional correlations between ′ G is and G* is defined via (xCor(·, ·) stands for "cross-sectional correlation" to distinguish it from "serial correlation" Cor(·, ·) we use above) 39 We give this matrix for Clustering-A with weights using normalized regressions with exposures computed based on arithmetic means (see Section 2.6) in Table 5. Let us mention that, with exposures based on arithmetic means, weights using normalized regressions work a bit better than using unnormalized regressions. Using exposures based on geometric means changes the weights a bit, which in turn slightly affects the within-cluster correlations, but does not alter the qualitative picture.

Overall correlations
Another useful metric, which we use as a sanity check, is this. For each value of s (i.e., for each cancer type), we can run a linear crosssectional regression (without the intercept) of ′ G is over the matrix W iA . So, we have n = 14 of these regressions. Each regression produces multiple R 2 and adjusted R 2 , which we give in Table 5. Furthermore, we can compute the fitted valuesĜ * is based on these regressions, which are given by where (for each value of s) F As are the regression coefficients. We can now compute the overall cross-sectional correlations (i.e., the index i runs over all N = 96 mutation categories) These correlations are also given in Table 5 and measure the overall fit quality.

Interpretation
Looking at Table 5 a few things become immediately evident. Clustering works well for 10 out the 14 cancer types we study here. The cancer types for which clustering does not appear to work all that well are Breast Cancer (labeled by X4 in Table 5), Liver Cancer (X8), Lung 36 The R function kmeans() produces a warning if it does not converge within iter.max. 37 We ran these 15 batches consecutively, and each batch produced the same top-10 (by occurrence counts) clusterings as in Table S4; however, the actual occurrence counts are different across the batches with slight variability in the corresponding rankings. The results are pleasantly stable. 38 De-noising per se does not affect cross-sectional correlations. Adding extra 1 in (3) (recall that we obtain ′ G is by cross-sectionally demeaning X is and then re-exponentiating) has a negligible effect. So, in the correlations below we can use the original data matrix G is instead of ′ G is . 39 Due to the factorized structure (22), these correlations do not directly depend on H As .

Z. Kakushadze, W. Yu
Biomolecular Detection and Quantification 13 (2017)  Cancer (X9), and Renal Cell Carcinoma (X14). More precisely, for Breast Cancer we do have a high within-cluster correlation for Cl-5 (and also Cl-4), but the overall fit is not spectacular due to low within-cluster correlations in other clusters. Also, above 80% within-cluster correlations 40 arise for 5 clusters, to wit, Cl-1, Cl-3, Cl-4, Cl-5 and Cl-6, but not for Cl-2 or Cl-7. Furthermore, remarkably, Cl-1 has high within-cluster correlations for 9 cancer types, and Cl-5 for 6 cancer types. These appear to be the leading clusters. Together they have high within-cluster correlations in 11 cancer types. So what does all this mean? Additional insight is provided by looking at the within-cluster correlations between the 7 cancer signature extracted in [8] and the clusters we find here. Let iα W be the weights for the 7 cancer signatures from Tables 13 and 14 of [8]. We can compute the following withincluster correlations (α = 1, …, 7 labels the cancer signatures of [8], which we refer to as Sig1 through Sig7): These correlations are given in Table 6. High within-cluster correlations arise for Cl-1 (with Sig1 and Sig7), Cl-5 (with Sig2) and Cl-6 (with Sig4). And this makes perfect sense. Indeed, looking at Figs. 14 through 20 of [8], Sig1, Sig2, Sig4 and Sig7 are precisely the cancer signatures that have "peaks" (or "spikes" -"tall mountain landscapes"), whereas Sig3, Sig5 and Sig6 do not have such "peaks" ("flat" or "rolling hills landscapes"). No wonder such signatures do not have high withincluster correlationsthey simply do not have cluster-like structures. Looking at Fig. 21 in [8], it becomes evident why clustering does not work well for Liver Cancer (X8)it has a whopping 96% contribution from Sig5! Similarly, Renal Cell Carcinoma (X14) has a 70% contribution from Sig6. Lung Cancer (X9) is dominated by Sig3, hence no cluster-like structure. Finally, Breast Cancer (X4) is dominated by Sig2,  Table 2 with the next 48 mutation categories) for the 7 clusters in Clustering-A (see Table S4) based on unnormalized (columns 2-8) and normalized (columns 9-15) regressions (see Section 2.6 for details). Each cluster is defined as containing the mutations with nonzero weights. (The mutations are encoded as follows: XYZW = Y>W: XYZ. Thus, GCGA = C>A: GCG.) For instance, cluster Cl-2 contains 8 mutations GCGA, TCGA, ACGG, GCCG, GCGG, TCGG, GTCA, GTCG. In each cluster the weights are normalized to add up to 100% (up to 2 digits due to the aforesaid rounding). In Tables 1 through S10 "weights based on unnormalized regressions" are given by (13), (14) and (15), while "weights based on normalized regressions" are given by (17), (14) and (16), i.e., the exposures are calculated based on arithmetic averages (see Section 2.6 for details). which has a high within-cluster correlation with Cl-5, which is why Breast Cancer has a high within-cluster correlation with Cl-5 (but poor overall correlation in Table 5 Cancer (X13) apparently all have the Cl-5 cluster structure embedded in them substantially. Furthermore, note the overlap between these two lists, to wit, Esophageal Cancer(X6), Ovarian Cancer (X11), Pancreatic Cancer (X12) and Prostate Cancer (X13). We obtained this result purely statistically, with no biologic input, using our clustering algorithm and other statistical methods such as linear regression to obtain the actual weights. It is too early to know whether this insight will aid any therapeutic applications, but that is the hopesimilarities in the underlying genomic structures of different cancer types raise hope that therapeutics for one cancer type could perhaps be applicable to other cancer types. On the other hand, our findings above relating to Liver Cancer, Lung Cancer and Renal Cell Carcinoma (and possibly also Breast Cancer, albeit the latter does appear to have a not-so-insignificant overlap with Cl-5, which differentiates it from the aforesaid 3 cancer types) suggest that these cancer types apparently stand out.

Concluding remarks
Clustering ideas and techniques have been applied in cancer research in various incarnations and contexts aplentyfor a partial list of works at least to some extent related to our discussion here, see, e.g., [52,53,54,55,40,56,5,36, and references therein. As mentioned above, even in NMF clustering is used at some (perhaps not-so-evident) layer. What is new in our approachand hence new resultsis that: (i) following [8], we apply clustering to aggregated by cancer types and   de-noised data; ii) we use a tried-and-tested in quantitative finance bag of tricks from [11], which improves clustering; and (iii) last but not least, we apply our *K-means algorithm to cancer genome data. As mentioned above, *K-means, unlike vanilla k-means or its other commonly used variations, is essentially deterministic, and it achieves determinism statistically, not by "guessing" initial centers or as in  agglomerative hierarchical clustering, which basically "guesses" the initial (e.g., 2-cluster) clustering. Instead, via aggregating a large number of k-means clusterings and statistical examination of the occurrence counts of such aggregations, *K-means takes a mess of myriad vanilla k-means clusterings and systematically reduces randomness and indeterminism without ad hoc initial "guesswork".      As mentioned above, consistently with the results of [8] obtained via improved NMF techniques, Liver Cancer, Lung Cancer and Renal Cell Carcinoma do not appear to have clustering (sub)structures. This could be both good and bad news. It is a good news because we learned something interesting about these cancer typesand in two complementary ways. However, it could also be a bad news from the therapeutic standpoint. Since these cancer types appear to have little in common with others, it is likely that they would require specialized therapeutics. On the flipside, we should note that it would make sense to exclude these 3 cancer types when running clustering analysis. However, it would also make sense to include other cancer types by utilizing the International Cancer Genome Consortium data, which we leave for future studies. (For comparative reasons, here we used the same data as in [8], which was limited to data samples published as of the date thereof.) This paper is not intended to be an exhaustive empirical study but a proof of concept and an opening of a new avenue for extracting and studying cancer signatures beyond the tools that NMF provides.
And we do find that 11 out of the 14 cancer types we study here have clustering structures substantially embedded in them and clustering overall works well for at least 10 out of these 11 cancer types. 41 Now, looking at Fig. 14 of [8], we see that its "peaks" are located at ACGT, CCGT, GCGT and TCGT. The same "peaks" are present in our cluster Cl-1 (see Figs. 2 and 3). Hence the high within-cluster correlation between Cl-1 and Sig1. On the other hand, Sig1 of [8] is essentially the same as the mutational signature 1 of [40,36], which is due to spontaneous cytosine deamination. So, this is what our cluster Cl-1 describes. Next, looking at Fig. 15 of [8], we see that its "peaks" are located at TCAG, TCTG, TCAT and TCTT. The first two of these "peaks" TCAG and TCTG are present in our Cl-5 (see Figs. 10 and 11), the third "peak" TCAT is present in our Cl-1 (see Figs. 2 and 3), while the fourth "peak" TCTT is present in our Cl-4 (see Figs. 8 and 9), which is consistent with the high within-cluster correlations between Sig2 and Cl-4 Table 3 Weights (in the units of 1%, rounded to 2 digits) for the first 48 mutation categories for the 7 clusters in Clustering-A (see Table S4) based on unnormalized (columns 2-8) and normalized (columns 9-15) regressions with the exposures computed via geometric means (see Section 2.6 for details). Here "weights based on unnormalized regressions" are given by (13), (14) and (19), while "weights based on normalized regressions" are given by (17), (14) and (21). Other conventions are the same as in Table 1.   Table 5 The within-cluster cross-sectional correlations Θ sA (columns 2-8), the overall correlations Ξ s (column 11) based on the overall cross-sectional regressions, and multiple R 2 and adjusted R 2 of these regressions (columns 9 and 10). See Section 3.3 for details. Cancer types are labeled by X1 through X14 as in and Cl-5, albeit its within-cluster correlation with Cl-1 is poor. Note that Sig2 of [8] is essentially the same as the mutational signatures 2 + 13 of [40,36], which are due to APOBEC mediated cytosine deamination. In fact, it was reported as a single signature in [36], however, subsequently, it was split into 2 distinct signatures, which usually appear in the same samples. 42 Our clustering results indicate that grouping TCAG and TCTG into one signature makes sense as they belong to the same cluster Cl-5. However, grouping TCAT and TCTT together does not appear to make much sense. Looking at the figures for Clustering-A, Clustering-B, Clustering-C and Clustering-D, we see that the TCAT "peak" invariably appears together with the ACGT, CCGT, GCGT and TCGT "peaks" as in Cl-1 in Clustering-A, Cl-2 in Clustering-B, Cl-1 in Clustering-C, and Cl-1 in Clustering-D, but never with TCTT. So, our clustering approach tells us something new beyond the NMF "intuition". This may have an important implication for Breast Cancer, which, as mentioned above, is dominated by Sig2. Thus, based on our results in Table 5, we see that Breast Cancer has high within-cluster correlations with Cl-4 and Cl-5, but not with Cl-1. This may imply that clustering simply does not work well for Breast Cancer, which would appear to put it in the same "stand-alone" league as Liver Cancer, Lung Cancer and Renal Cell Carcinoma. In any event, clustering invariably suggests that the TCAT "peak" belongs in Cl-1 with the 4 "peaks" ACGT, CCGT, GCGT and TCGT related to spontaneous cytosine deamination, rather than those related to APOBEC mediated cytosine deamination. Now, let us check the remaining two signatures of [8] with "tall mountain landscapes" (see above), to wit, Sig4 and Sig7. Looking at Fig. 17 of [8], we see that its "peaks" are at CTTC, TTTC, CTTG and TTTG. The same peaks appear in our Cl-6 (see Figs. 12 and 13). Hence the high within-cluster correlation between Cl-6 and Sig4. Note that Sig4 is essentially the same as the mutational signature 17 of [40,36], and its underlying mutational process is unknown. Next, looking at Fig. 20 of [8], we see that its "peaks" for the C > G mutations are essentially the same as in Cl-1. Hence the high within-cluster correlation between Cl-7 and Sig1. So, there are no surprises with Sig1, Sig4 and Sig7. However, based on our clustering results, as we discuss above, with Sig2 we do findwhat we feel is a pleasantsurprise, that splitting it into two signatures (see above) might be inadequate and the TCAT "peak" might really belong with the Sig1 "peaks" (spontaneous v. APOBEC mediated cytosine deamination). This is exciting as it might be an indication of the limitations of NMF (or clustering…). 43 In Introduction we promised that we would discuss some potential applications of *K-means in quantitative finance, and so here it is. Let us mention that *K-means is universal, oblivious to the input data and applicable in a variety of fields. In quantitative finance *K-means a priori can be applied everywhere clustering methods are used with the added bonus of (statistical) determinism. 44 One evident example is statistical industry classifications discussed in [11], where one uses clustering methods to classify stocks. In fact, *K-means is an extension of the methods discussed in [11]. One thing to keep in mind is that in *K-means one sifts through a large number P of aggregations, which can get computationally costly when clustering 2000+ stocks into 100+ clusters. 45 Another potential application is in the context of combining alphas (trading signals)see, e.g., [79]. Yet another application is when we have a term structure, such as a portfolio of bonds (e.g., U.S. Treasuries or some other bonds) with varying maturities, or futures (e.g., Eurodollar futures) with varying deliveries. These cases resemble the genome data more in the sense that the number N of instruments is relatively small (typically even fewer than the number of mutation categories). Another example with a relatively small number of instruments would be a portfolio of various futures for different FX (foreign exchange) pairs (even with the uniform delivery), e.g., USD/ EUR, USD/HKD, EUR/AUD, etc., i.e., FX statistical arbitrage. One approach to optimizing risk in such portfolios is by employing clustering methods and a stable, essentially deterministic algorithm such as *Kmeans can be useful. Hopefully *K-means will prove a valuable tool in cancer research, quantitative finance as well as various other fields (e.g., image recognition).

Conflict of interest
Authors declare no conflict of interest.

Appendix A. R source code
In this appendix we give the R (R Package for Statistical Computing, http://www.r-project.org) source code for computing the clusterings and weights using the algorithms of Section 2. The code is straightforward and self-explanatory. 46 The main function is bio.cl.sigs(x, iter.max = 100, num.try = 1000, num.runs = 10000). Here: x is the N × d occurrence counts matrix G is (where N = 96 is the number of mutation categories, and d is the number of samples; or d = n, where n is the number of cancer types, when the samples are aggregated by cancer types); iter.max is the maximum number of iterations that are passed into the R built-in function kmeans(); num.try is the number M of aggregated clusterings (see Section 2.3.2); num.runs is the number of runs P used to determine the most frequently occurring clustering (the "ultimate" clustering) obtained via aggregation (see Section 2.3.3). The function bio.erank.pc() is defined in Appendix B of [8]. The function qrm.stat.ind.class() is defined in Appendix A of [11]. This function internally calls another function qrm.calc.norm.ret(), which we redefine here via the function bio.calc.norm.ret(). 47 The output is a list, whose elements are as follows: res$ind is an N × K binary matrix Ω iA = δ G(i),A (i = 1, …, N, A = 1, …, K, the map G : {1, …, N} ↦ {1, …, K}see Section 2), which defines the K clusters in the "ultimate" clustering; 48 res$w is an N-vector of weights obtained via unnormalized regressions using arithmetic means for computing Table 6 The within-cluster cross-sectional correlations Δ αA between the weights for 7 cancer signatures Sig1 through Sig7 of [8] and the weights (using normalized regressions with exposures based on arithmetic averages) for 7 clusters in Clustering A (see Section 3.3 for details). All quantities are in the units of 1% rounded to 2 digits. The values above 80% are given in bold font.  43 Or both… Alternativelyand that would be truly excitingperhaps there is a biologic explanation. In any event, it is too early to tellyet another possibility is that this is merely an