In this section we describe in detail our technique for finding clusters like the example in Figure . A gene expression matrix is an N
matrix of real-valued measurements X
. The rows are genes, the columns are samples, and xij
is the measured (log) expression, relative to a baseline. Typically there are missing entries in X.
We use a technique described in [17
], an iterative algorithm based on the singular value decomposition, for imputing missing expression values;our analysis here assumes that X
has no missing values.
Let Sk be the indices of a cluster of k genes, and
be the collection of p column averages of the expression values for this cluster. Then for each cluster size k, gene shaving seeks a cluster Sk having the highest variance of the column averages:
The important question of how to choose the cluster size k is addressed in the next section.
Our procedure generates a sequence of nested clusters Sk, in a top-down manner, starting with k = N, the total number of genes, and decreasing down to k = 1 gene. At each stage the largest principal component of the current cluster of genes is computed. This eigen gene is the normalized linear combination of genes with largest variance across the samples. We then compute the inner product (essentially the correlation) of each gene with the eigen gene, and discard ('shave off') a fraction of the genes having lowest (absolute) inner product. The process is repeated on the reduced cluster of genes. The shaving algorithm is depicted in Figure and given in detail in Box .
Schematic of the gene shaving process.
The Shaving algorithm.
There are a number of duplicate genes in the dataset. In some cases the sequence for a given gene appears on the microarray more than once, either by design or by accident. In other cases, more than one different EST (expressed sequence tag) is present for the same gene. Gene shaving can be affected by duplicate genes, since they are highly correlated with each other. We therefore averaged expression profiles for the duplicate genes, leaving 3624 unique gene profiles.
The sequence of operations 1-5 in Box gives the first cluster of rows - the first ribbon in Figure . Step 6 orthogonalizes the data to encourage discovery of a different (uncorrelated) second cluster. Note that although we work with the orthogonalized matrix in the shaving process for the second and subsequent clusters, the derived clusters and their averages involve the original genes.
The shaving process requires repeated computation of the largest principal component of a large set of variables. If naively implemented, this requires the construction of a N × N sample covariance matrix ∑ of the genes, and the computation of its largest eigenvector. We can avoid the computational burden by working in the dual space, where the matrices have dimension p × p. Furthermore, as we require only the largest eigenvector, the computations can be reduced even further by using the power method, using the eigenvector of the previous cluster as a starting value.
The three resulting clusters are shown in Figure and again in Figure . Figure shows the pairwise scatterplots of each of the three column averages ('super genes') from the clusters. The absolute correlations range from 0.27 to 0.68. The full gene names for the members of the first three clusters are given in Table .
The three gene clusters from unsupervised shaving
It is useful to evaluate how much of the dimensionality of the gene expression variation is captured by the clusters derived from gene shaving. We can approximate the expression profile for each gene in the complete dataset as a linear combination of the three super genes from each cluster (which are simple averages of the genes in each cluster). The percent variance explained by the first j = 1,2, ... 10 super genes is shown in Figure .
Figure 7 Percent of gene variance explained by first j gene shaving column averages (j = 1,2,... 0) (solid curve), and by first j principal components (broken curve). For the shaving results, the total number of genes in the first j clusters is also indicated. (more ...)
Thus the three gene clusters, involving a total of 33 genes, explain about 20% of the variation. The percent variance explained by the first j principal components (broken curve) exceeds that from gene shaving. Each principal component gives a non-zero weight to all 3624 genes, however.
The gap estimate of cluster size
Each shaving sequence produces a nested set of gene clusters Sk, starting with the entire expression matrix and then proceeding down to some minimum cluster size (typically 1). If we applied this procedure to null data, in which the rows and columns were independent of each other, we could still find some interesting-looking patterns in the resulting blocks. Hence, we need to calibrate this process so that we can differentiate real patterns from spurious ones. Even in the case of real structure, it is unlikely that a distinct set of genes is correct for a cluster, and the rest not. More likely there is a graduation of membership eligibility, and we have to decide where to draw the line. Here we describe a procedure based on randomization that helps us select a reasonable cluster size.
Our method requires a quality measure for a cluster. We favor both high-variance clusters, and high coherence between members of the cluster. As the generation of the cluster sequence was driven strongly by the former, we focus on the latter in selecting a good cluster. By analogy with the analysis of variance for grouped data, we define the following measures of variance for a cluster Sk:
The between variance is the variance of the (signed) mean gene. The within variance measures the variability of each gene about the cluster average, also averaged over samples. As this can be small if the overall variance is small, a more pertinent measure is the between-to-within variance ratio VB/VW, or alternatively, the percent variance explained
A large value of R2 implies a tight cluster of coherent genes. This is the quality measure we use to select a cluster from the shaving sequence Sk.
Let Sk index the clusters of a given shaving sequence (with k being the number of genes). Let Dk be the R2 measure for the kth member of sequence. We wish to know whether Dk is larger than we would expect by chance, if the rows and columns of the data were independent.
be a permuted data matrix, obtained by permuting the elements within each row of X.
We form B
such matrices, indexed by b
= 1,2,... B
. Let Dk*b
be the R2
measure for cluster Sk*b
. Denote by k*
the average of Dk*b
. The Gap
function is defined by
We then select as the optimal number of genes that value of k producing the largest gap:
The idea is that at the value
the observed variance is the most ahead of expected. Multiple clusters are produced for the randomized data just like for the original data, and the gap test is used repeatedly to select the cluster size at each stage.
For the DLCL data, the maximum for the first cluster occurs at eight genes. Figure shows the percent-variance curves, Dk, for both the original and randomized tumor data as a function of size, and the gap curves used to select the specific cluster sizes in Figure .
Figure 8 (a) Variance plots for real and randomized data. The percent variance explained by each cluster, both for the original data, and for an average over three randomized versions. (b) Gap estimates of cluster size. The gap curve, which highlights the difference (more ...)
Predicting patient survival
One important motivation for developing gene shaving was the wish to identify distinct sets of genes whose variation in expression could be related to a biological property of the samples. In the present example, finding genes whose expression correlates with patient survival is an obvious challenge. Group factors g1
were created by splitting each gene cluster in Figure into two groups of 24 patients. We used each of these groupings as a factor in Cox's proportional hazards model for predicting overall survival [18
]. Of the group factors only g2
was significant, at the 0.05 level (p
], a cluster of 380 genes was chosen on the basis of their large variation over samples, and their 'germinal center B-like' or 'activated B-like' expression profiles. Using these 380 genes, a hierarchical clustering produced two groups of patients which were (just) statistically different in survival. Close inspection shows that 18 of the 23 genes in the second cluster above also fall into this cluster of 380 genes. Hence, gene shaving can find clinically and biologically relevant subdivisions in gene expression data in an unsupervised fashion.
It may be fortuitous that one of these groupings correlates with survival, as the clusters were not chosen with survival in mind. We next describe a modification of gene shaving that explicitly looks for clusters that are related to patient survival.
The methods discussed so far have not used information about the columns to 'supervise' the shaving of rows. Here we generalize gene shaving to incorporate full or partial supervision.
As in Equation (1), we consider a cluster of genes Sk
having column average vector
. Let y
, ... yp
) be a set of auxiliary measurements available for the samples. For example each yj
, might be a survival time for the patient corresponding to sample j
or a class label for each sample, such as a diagnosis category. Supervised shaving maximizes a weighted combination of column variance and an information measure J
for fixed 0 ≤ α ≤ 1. The value α = 1 gives full supervision; values between o and 1 provide partial supervision.
Choice of the measure J
) depends on the nature of the auxiliary information y
. If the y
codes class labels, J
) can be taken as the sum of squared differences between the category averages
. For censored survival times y
, think of
as a covariate in a Cox (proportional hazards) model. If the score vector from this model is g,
we set J
) = ggT
, a p
matrix. Computationally we first scale the genes so that the within-risk set variance is 1.
When fully supervised, the shaving procedure reduces to simply ranking the genes from largest to smallest Cox model score test. Thus there is no role for principal components or top-down shaving in this case. However, to encourage coherence within the gene clusters, it can be better to use a partially supervised criterion, which does use principal components and top-down shaving. This is demonstrated in the example below. One can also include other prognostic factors in the model, and direct shaving to find a gene cluster whose column average is a strong predictor in the model. This can be done with other models, for example a linear regression model for a quantitative measurement. Operationally, all of these choices for J
are quadratic functions of the column averages
, and gene shaving can be carried out via principal components of a suitably modified variance matrix.
We applied supervised shaving to the set of 3624 genes from the DLCL samples. Figure examines the effect of different values of the supervision weight α, showing the average (absolute) gene correlation and Cox model p value for each choice. From this plot we chose the value α = 0.10, which gives a good mix of high gene correlation and low p value. Partially supervised gene shaving produced a cluster with 234 genes, shown in Figure and in Table .
Average (absolute) gene correlation and Cox model p value, for clusters of size 200 from supervised shaving and for different values of α. The value of Qa = 0.1 seems best, and is used in the gene shaving procedure.
Cluster of 234 genes from supervised shaving.
Cluster from supervised shaving applied to full set of 3624 genes
Some of the genes are close together in the hierarchical clustering order (indicated by the first number in Table ), many are not. Some genes have a negative sign, and others have no sign. We will call these 'negative' and 'positive' genes respectively. The negative genes have their expression values flipped before being averaged with other gene expression profiles. Figure shows the gap curve, suggesting a cluster size of 35. However, further analysis below suggests the better cluster size of 234.
Figure 11 (a) Gap curve for supervised shaving. (b) Survival curves in the two groups defined by the low or high expression of the 234 genes. Group I has high expression of positive genes, and low expression of negative genes; group 2 has low expression of positive (more ...)
The cluster of 234 genes contains many of the strongest individual genes for predicting survival. For example, 130 of the strongest 200 genes are in the cluster. Some weaker genes are, however, also in the cluster, the weakest being the 1332nd strongest gene among the full list of 3624. Figure shows the survival curves stratified by the low and high expression of this gene cluster, using the median of the cutoff point. The two resulting groups are shown in Figure .
The two groups of samples that showed highest and lowest expression of the gene cluster associated with survival.
Using this grouping as a predictor in the Cox model for survival gave a highly significant p value (0.00042). However, this p value must be scrutinized. Figure , shows the Cox model p value as a function of the cluster size, for both partially and fully supervised shaving. We will call these the 'training set p values'. As the gene shaving process was supervised by the survival times, the training set p values will be biased downward, and it is important to adjust them appropriately. We permuted the survival times, re-ran the shaving process and computed the corresponding p values. This was repeated 100 times, and for each cluster size we computed the proportion of times the permutation p values were less than or equal to the training set p values. These proportions are shown in Figure ,, and are unbiased estimates of the true p values. Fully supervised shaving is too aggressive, and the upward adjustment of the p values is large. As a result the p value is around 0.05 for the full set of genes, but none of the smaller clusters is significant at the 0.05 level. For partially supervised shaving, many of the p values are below 0.05, and from this we chose the cluster size of 234 near the minimum.
Figure 13 Supervised gene shaving from full gene set. (a,c) Partially supervised with α = 0.10; (b,d) fully supervised (α = 1). (a,b) Training set p values; (c,d) permutation p values for the cluster average as a function of cluster size. The chosen (more ...)
Using the full set of genes, flipping each to have positive correlation with survival, averaging their expression values and finally cutting at the median, gave a grouping nearly the same as Groups 1 and 2 in Figure . The only change was a swap between DLCL-oo14 and DLCL-oo18, and these two samples are right at the median cutpoint between the two groups in Figure .
The international prognostic index (IPI) A score was also available for these patients. Components of the IPI include age, level of the enzyme lactate dehydrogenase (LDH) and the number of extranodal sites. As in [14
], we dichotomized IPI scores into low (0, 1 or 2) and high (3, 4 or 5). The resulting grouping seems to be about as predictive as the IPI index, and is quite independent from it, as Table indicates.
Cross-tabulation of gene shaving groups with IPI index
When added to a Cox model containing IPI, this grouping had a training set p value of 0.0006. Figure shows the Kaplan-Meier survival curves for each group, stratified by low and high IPI.
Figure 14 Kaplan-Meier survival curves in the two groups defined by the cluster of 234 genes shown in Figure , stratified by IPI. Group I has high expression of positive genes and low expression of negative genes in Figure , (more ...)
], two patient groups were defined from a hierarchical clustering tree grown from a 38o-gene cluster. As a predictor, the grouping was just significant in the low IPI group only, at the 0.05 level. Table gives a cross-tabulation of that grouping with the one used in this paper in Figure .
Table 4 A comparison of the patient groups obtained by gene shaving with those obtained previously 
Thus 25/36 = 69% of the patients are classified the same way by both groupings. The patient grouping of Alizadeh et al.
] was based on a cluster of 380 genes, chosen for their large variation over the samples. Our cluster of 234 genes has 38 genes in common with this cluster of 380, and they are indicated by an asterisk in Table . Five of the 234 genes also appear in the unsupervised clusters found earlier, in the second of the three clusters.