Biclustering and Co-Clustering

Biclustering did not begin with a cloud of points. It began with two rectangles. One was a yeast (Saccharomyces cerevisiae) cell-cycle matrix, 2,884 genes by 17 conditions. The other was a human B-cell lymphoma matrix, 4,026 genes by 96 conditions. Cheng and Church brought biclustering to expression data on those two matrices in 2000.

In a matrix like that, a gene is not simply similar or dissimilar to another gene. Two genes can move together across a handful of the 17 conditions and go their own way across the rest. A method that must score similarity across all 17 columns at once averages that structure away.

Biclustering searches for local row–column relationships instead. The unit of structure is a submatrix: a set of rows named together with the set of columns on which those rows cohere. It is not a partition that must explain every feature equally. The same reasoning transfers wherever data arrives as a rectangle — documents by terms, users by items, tumour samples by genes.

One object can belong to different patterns depending on the variables considered. In every case below the rows and the relevant columns were found together, on a named matrix.

The difference is whether every feature contributes to the same notion of similarity. The size of that difference has been measured rather than asserted. In 2003 Dhillon and colleagues put information-theoretic co-clustering against ordinary one-dimensional document clustering. On the CLASSIC3 collection the two methods are close, and the authors say so exactly: “Observe that co-clustering extracted the original clusters almost correctly resulting in a micro-averaged-precision of 0.9835 while 1D-clustering led to a micro-averaged-precision of 0.9432.”

Their Table 5 is where the gap opens. On the 20-Newsgroups subsets, co-clustering against 1D-clustering: Binary 0.98 vs 0.64, Binary subject 0.96 vs 0.67, Multi5 0.87 vs 0.34, Multi5 subject 0.89 vs 0.37, Multi10 0.56 vs 0.17, Multi10 subject 0.54 vs 0.19. The harder the collection, the more of the answer sits in the coupling between documents and words, and the less a single global feature space recovers.

Those numbers became a reference baseline. Thirteen years later another group reprinted them in their own hierarchical co-clustering paper, as the “ICC [14]” column: 0.87 on Multi5, 0.89 on Multi5s, 0.56 on Multi10, 0.54 on Multi10s. The same table is also the honest counterweight. Row-only agglomerative clustering on Multi10 reaches 0.24 with single-link, but 0.60 with complete-link — above the co-clustering figure. Joint structure is not a guarantee. It is a hypothesis about where the structure lives, and it wins by margins that depend on the collection.

The three families differ in what they commit to. Ordinary row clustering groups rows using all selected columns in one feature space: one membership per row, a global feature weighting to justify, local subspace patterns it can miss, and the practical advantage that a new row is easy to assign. Biclustering finds row subsets coherent on selected column subsets: it captures local interactions, can produce overlapping memberships, complicates interpretation, and needs coverage and redundancy controls. Matrix factorization represents rows and columns through continuous latent factors: graded participation, an ability to reconstruct missing entries, factors that may be harder to name, and identifiability that has to be treated with caution. FABIA, factor analysis for bicluster acquisition, is the factorization member of the twelve algorithms Eren and colleagues compared.

Evaluate both the internal block and the collection of blocks. Five steps, three of which have published instruments rather than good intentions.

1. Define matrix semantics — document what rows, columns, zeros, missing entries and exposure represent. Cheng and Church had to make this decision in public. The human B-cell lymphoma matrix carried 4,026 genes and 96 conditions with 47,639 missing values, 12.3% of all matrix elements. Their choice was deliberate rather than defensive: “Missing data in the matrices were replaced with random numbers. The expectation was that these random values would not form recognizable patterns and thus would be the leading candidates to get removed in node deletion.” The lesson is not that random filling is right. It is that 12.3% of that matrix was a modelling decision, written down, with a stated expectation of what the algorithm would do with it.

2. Normalize appropriately — control row scale, column scale, sparsity, and batch effects before fitting. The choice is part of the method, not preparation for it. scikit-learn's SpectralBiclustering exposes three normalizations: independent row and column normalization, bistochastization, and log. The block structure you recover depends on which you pick.

3. Inspect block coherence — review values, residuals, exemplars, and domain meaning inside each bicluster. A δ-bicluster is defined by a residual score, so the score and the biology have to be read together.

4. Measure collection quality — quantify coverage, overlap, redundancy, and unassigned evidence. Cheng and Church did this for their first 100 biclusters and reported it as a number: “In the yeast data experiment, the 100 biclusters covered 2,801, or 97.12% of the genes, 100% of the conditions, and 81.47% of the cells in the matrix.” On the human B-cell lymphoma matrix the same 100 blocks covered 3,687 genes, 91.58%, and 100% of the conditions — but only 36.81% of the cells. Same method, same number of blocks, less than half the cell coverage. A count of biclusters says nothing about how much of a matrix has been explained. Only the coverage figure does.

5. Test stability — refit across samples and preprocessing choices, aligning similar blocks before comparison. The alignment is the hard part, and it has a standard instrument: the consensus score of Hochreiter and colleagues. scikit-learn computes it by taking pairwise Jaccard indices, solving the one-to-one assignment with scipy.optimize.linear_sum_assignment (a modified Jonker-Volgenant algorithm), and dividing by the size of the larger set. It is the measure on which their own method was judged. Their abstract reports that “On 100 simulated datasets with known true, artificially implanted biclusters, FABIA clearly outperformed all 11 competitors”, while on three microarray datasets with known subclusters “it was two times the best and once the second best method”.

Joint methods differ in whether rows, columns, or local submatrices may overlap. Each option below is a published model with shipped code behind it.

Co-clustering partitions rows and columns into a global checkerboard of blocks. That word is a specification, not a picture. Kluger and colleagues set it out in Genome Research in 2003, and scikit-learn states the assumption together with its arithmetic: “The SpectralBiclustering algorithm assumes that the input data matrix has a hidden checkerboard structure. The rows and columns of a matrix with this structure may be partitioned so that the entries of any bicluster in the Cartesian product of row clusters and column clusters are approximately constant. For instance, if there are two row partitions and three column partitions, each row will belong to three biclusters, and each column will belong to two biclusters.” A second, independent implementation describes the same model in the same terms. The R package biclust says its BCSpectral method “Performs Spectral Biclustering as described in Kluger et al., 2003. Spectral biclustering supposes that normalized microarray data matrices have a checkerboard structure that can be discovered by the use of svd decomposition in eigenvectors, applied to genes (rows) and conditions (columns).”

Non-overlapping biclustering finds row–column blocks without assigning one element to several blocks. scikit-learn's SpectralCoclustering is the strict case, citing Dhillon 2001 and a bipartite-graph normalized cut. Each row and each column ends in exactly one bicluster.

Overlapping biclustering allows a row or column to participate in multiple local patterns. The Cheng and Church δ-bicluster, found one block at a time by node deletion, is of this kind. So is the order-preserving submatrix.

The factor model view represents the matrix through latent row and column factors rather than discrete blocks. FABIA sits here. Eren and colleagues placed all of these side by side in Briefings in Bioinformatics in 2013: “Twelve algorithms were chosen for comparison in this article” — Cheng and Church, OPSM, xMOTIFs, QUBIC, BiMax, ISA, COALESCE, Plaid, BBC, CPB, FABIA and Spectral. The last of those “uses singular value decomposition to find a checkerboard pattern in the data in which each bicluster is up- or downregulated”. Their conclusion is a design principle for the whole map: “Algorithms that model the entire dataset seem more resilient to noise than algorithms that seek individual biclusters”.

The richer question is “which rows behave similarly with respect to which columns?” That formulation can reveal interactions hidden by global feature weighting. The extra flexibility raises the burden on whoever has to read the result. Every reported block should justify its coherence, distinctiveness, coverage, and intended use.

That flexibility was priced in 2000, in the paper that introduced it. Cheng and Church's Theorem 1 shows that finding the largest square δ-bicluster is NP-hard, by reduction from BALANCED COMPLETE BIPARTITE SUBGRAPH: “Theorem 1. The problem of finding the largest square δ-bicluster (|I| = |J|) is NP-hard”. The paper is candid about the cost. NP-hardness, they note, “implies that there may be sizable biclusters with good scores evading the search by any efficient algorithm”.

The price is not an artifact of their particular definition. Ben-Dor and colleagues changed the definition entirely and met the same wall. Their 2003 abstract states the model and the bill in one sentence: “Specifically, we look for order-preserving submatrices (OPSMs), in which the expression levels of all genes induce the same linear ordering of the experiments (we show that the OPSM search problem is NP-hard in the worst case).” Their Section 3, “OPSM IS NP-COMPLETE”, proves “Theorem 3.1. OPSM is NP-complete” by reduction from the balanced complete bipartite subgraph problem — the same reduction family Cheng and Church used. They ran the model anyway, on a breast tumour matrix of 3,226 genes and 22 samples. Seven years on, the difficulty had not softened. Trapp and Prokopyev opened their 2010 paper in INFORMS Journal on Computing by restating it — “This problem is known to be NP-hard” — and attacked it with integer programming, which is what one does when no efficient algorithm is available.

The comparison came six years after Cheng and Church. Prelić and colleagues put five biclustering algorithms against one deliberately simple binary reference model, Bimax, and against a hierarchical baseline, on synthetic and real data for Saccharomyces cerevisiae and Arabidopsis thaliana; the work appeared in Bioinformatics in 2006. Biclustering “in general has advantages over a conventional hierarchical clustering approach”. Yet “there are considerable performance differences between the tested methods”, which is the part a report has to state. And the third finding is the one that keeps a project honest: “already the simple reference model delivers relevant patterns within all considered settings”.