Discovery and Scoring of Protein Interaction Subnetworks Discriminative of Late Stage Human Colon Cancer*S

We used a systems biology approach to identify and score protein interaction subnetworks whose activity patterns are discriminative of late stage human colorectal cancer (CRC) versus control in colonic tissue. We conducted two gel-based proteomics experiments to identify significantly changing proteins between normal and late stage tumor tissues obtained from an adequately sized cohort of human patients. A total of 67 proteins identified by these experiments was used to seed a search for protein-protein interaction subnetworks. A scoring scheme based on mutual information, calculated using gene expression data as a proxy for subnetwork activity, was developed to score the targets in the subnetworks. Based on this scoring, the subnetwork was pruned to identify the specific protein combinations that were significantly discriminative of late stage cancer versus control. These combinations could not be discovered using only proteomics data or by merely clustering the gene expression data. We then analyzed the resultant pruned subnetwork for biological relevance to human CRC. A number of the proteins in these smaller subnetworks have been associated with the progression (CSNK2A2, PLK1, and IGFBP3) or metastatic potential (PDGFRB) of CRC. Others have been recently identified as potential markers of CRC (IFITM1), and the role of others is largely unknown in this disease (CCT3, CCT5, CCT7, and GNA12). The functional interactions represented by these signatures provide new experimental hypotheses that merit follow-on validation for biological significance in this disease. Overall the method outlines a quantitative approach for integrating proteomics data, gene expression data, and the wealth of accumulated legacy experimental data to discover significant protein subnetworks specific to disease.

We used a systems biology approach to identify and score protein interaction subnetworks whose activity patterns are discriminative of late stage human colorectal cancer (CRC) versus control in colonic tissue. We conducted two gel-based proteomics experiments to identify significantly changing proteins between normal and late stage tumor tissues obtained from an adequately sized cohort of human patients. A total of 67 proteins identified by these experiments was used to seed a search for proteinprotein interaction subnetworks. A scoring scheme based on mutual information, calculated using gene expression data as a proxy for subnetwork activity, was developed to score the targets in the subnetworks. Based on this scoring, the subnetwork was pruned to identify the specific protein combinations that were significantly discriminative of late stage cancer versus control. These combinations could not be discovered using only proteomics data or by merely clustering the gene expression data. We then analyzed the resultant pruned subnetwork for biological relevance to human CRC. A number of the proteins in these smaller subnetworks have been associated with the progression (CSNK2A2, PLK1, and IGFBP3) or metastatic potential (PDGFRB) of CRC. Others have been recently identified as potential markers of CRC (IFITM1), and the role of others is largely unknown in this disease (CCT3, CCT5, CCT7, and GNA12). The functional interactions represented by these signatures provide new experimental hypotheses that merit follow-on validation for biological significance in this disease. Overall the method outlines a quantitative approach for integrating proteomics data, gene expression data, and the wealth of accumulated legacy experimental data to discover significant protein subnetworks specific to disease.

Molecular & Cellular Proteomics 8:827-845, 2009.
A fundamental presumption of the -omics revolution is that high dimensional data sets resulting, for example, from pro-teomics and genomics experiments should be integrated with functional annotations to give a more complete account of the cellular changes underlying the etiology of human disease. Nevertheless the accumulation of specific gene annotations and experimental protein or gene expression data is presently outpacing data integration. Network modeling of protein-protein interactions provides a context for such data integration (1)(2)(3). These modeling approaches can build networks using databases created from literature curation, inference by homology, high throughput data, or a combination of these (4). Network generation, analysis, and modeling are clearly fundamental to a new generation of systems biology approaches that promise an improved understanding of the causes of human disease as well as providing novel biomarkers of its progression.
We undertook a systems biology approach to identify protein "signatures" that were significantly discriminative of late stage human colorectal cancer (CRC) 1 versus control. CRC continues to be the second leading cause of cancer death in adult Americans (5). Although a great deal of research is focused on the early detection of CRC, comparably less attention has been paid to understanding the pathophysiology of a late stage (Duke's D) phenotype. As the prognosis of a late stage diagnosis is significantly poorer (Ͻ10% long term survivability (5)) than one following early detection (Ͼ90% (5)), identifying significant network-level changes in a late stage cohort holds the possibility of more clearly elucidating the mechanisms of tumorigenesis specific to this phenotype.
Proteomics studies of CRC using tumor and adjacent normal tissue obtained by biopsy from human patients have been conducted (6 -8). Even more studies have profiled protein expression changes in colon cancer cell lines (9 -13). However, these tissue-based studies have used either a sample cohort of mixed pathologic stage or a cohort size that was smaller than optimal. As colon cancer is a disease that progresses over a number of years and is marked by distinct pathologic stages of increasing severity (Duke's A-D), it is reasonable to expect changes in the proteome that are associated with particular stages of disease. Hence a cohort of homogenous pathology may improve the detection of stagespecific changes. Further we expected that the most dramatic changes in protein expression, in terms of quantity and magnitude, would be detectable between control and tumor using a statistically robust late stage cohort. Fig. 1 outlines our overall experimental design and emphasizes an integrated -omics approach to understanding the pathophysiology of stage D colon cancer. We began by quantitative proteomics profiling of a late stage CRC cohort and used the differentially expressed proteins to seed a search for protein interaction subnetworks possibly involved in this disease ( Fig. 1, steps 1-11). Our list of differentially expressed proteins was imported into a bioinformatics data mining tool that permits a search of a database comprised of tens of thousands of manually curated protein-protein interactions (14). This provided us with a list of subnetworks that we were able to rank by significance and reduce to a subsequently manageable number on which to focus our analysis ( Fig. 1, steps [12][13][14]. We seeded our initial search for subnetworks with proteomics data (versus transcriptomics data) because changes in both protein expression and isoform abundance provide the most direct functional readout of the cell. As such, we expected our seed proteins would represent "fence posts" within the subnetworks, each with one or more functional roles. These subnetworks represent expansions of the cancer proteome in the sense that the algorithm we used (see "Experimental Procedures") builds an extended interactome comprised of many targets around a smaller set of seed proteins. This extended interactome, although it provides clues to the regulatory connections that drive the observed abundance changes, is merely qualitative and inferential. Thus, a potential criticism of this type of approach is that a small set of proteins that is potentially causative with respect to a disease state has simply been expanded to a larger set whose members may or may not be important to the phenotype. If the myriad of protein interaction networks (and there are many tools to build interactomes from a set of seed targets) are to be useful to researchers for informing new hypotheses, new methods are needed to quantitatively evaluate the significance of the targets within the subnetwork that is generated. To address this, we adapted a method described by Ideker and co-workers (15) to score our subnetworks, and then we systematically searched within each one using the metric of mutual information (MI) to identify statistically significant combinations of proteins that were highly discriminative of stage D cancer versus control. Gene expression data were an ideal basis for scoring our subnetworks because of their complete coverage, i.e. we were able to assign an mRNA expression value to every target.
Overall our guiding principle is that protein is the immediate effector of phenotype. Therefore, profiling changes in the proteome is likely to provide the most direct evidence for cellular changes causing, or resulting from, a disease. However, proteomics experiments typically have incomplete coverage of the proteome. In particular, gel-based expression experiments are most likely to detect high abundance proteins. These high abundance targets may include subnetwork nodes that participate in larger subnetworks of protein interactions and may also be regulated at the level of transcription. If so, patterns of mRNA expression can be useful for discriminating between disease and non-diseased states within these "discovered" subnetworks. Of course, mRNA expression data have the characteristic of whole genome coverage, and these data can enable queries for subnetworks of interest that are "saturated." Here we present an integrated approach to cancer biology, one that shows how proteomics data, genomics data, and a vast database of legacy experimental data can be integrated with MI scoring schemes to reveal protein signatures significantly discriminative of disease. The signatures are useful for focusing follow-on experiments to verify their functional role in a disease phenotype. In addition, our approach is very general for use with existing public data sets as well as newly generated data and can be applied in the context of multiple types of protein interaction networks.

EXPERIMENTAL PROCEDURES
Sample Preparation pI 3-10 Experiment-Tissue samples were procured from a human tissue repository at the Case Comprehensive Cancer Center (supplemental Data S2). In addition to a tumor biopsy during surgical resection, a normal biopsy adjacent to the patient's tumor was also taken, typically Ͼ10 cm from tumor. Validation of the tissue as normal or tumor (including stage of the tumor) was performed by a pathologist. Tissues were immediately frozen and stored at Ϫ80°C. A 50-mg sample provided an adequate mass of protein for the 2D DIGE experiment. The tissue was weighed and placed in lysis buffer (4% CHAPS, 7 M urea, 2 M thiourea, 30 mM Tris) on ice, and the cells were disrupted by a three-cycle sonication protocol in a 4°C cold room. A protease inhibitor mixture (Sigma-Aldrich, catalog number P8340) and a wide spectrum phosphatase inhibitor (Roche Applied Science) were added to the buffer at the manufacturer's suggested concentration to inhibit protein degradation and dephosphorylation, respectively. The homogenate was centrifuged at 12,000 rpm for 10 min, and the protein fraction was withdrawn by pipette. Protein concentration was quantified by a colorimetric assay, similar to the Bradford assay, using the 2D-Quant kit (GE Healthcare). Aliquots were stored at Ϫ80°C. pI 4 -7 Experiment-Protein fractions from the prior experiment were thawed and cleaned with the 2D-Cleanup kit (GE Healthcare), and the concentration was redetermined as before. Aliquots were re-stored at Ϫ80°C.

2D Gel Electrophoresis
We used the 2D DIGE system available from GE Healthcare (formerly Amersham Biosciences) described by Marouga et al. (39). This system provides two distinct advantages over conventional 2D PAGE. First, it allows for up to three distinct samples to be labeled by spectrally resolvable fluorophores (CyDyes Cy2, Cy3, and Cy5) and multiplexed in a single gel. Second, by using one of these CyDyes (typically Cy2) to label a pooled sample, constituted by a proportional amount of every sample in the experiment, the Cy2 dimension is useful as an internal standard. This internal standard is crucial in the image analysis phase to a confident assessment of real biological variation from gel to gel as distinct from changes arising from variance in protein loading. For the purpose of detection by image analysis, 50 g of protein is sufficient for labeling by each of the CyDyes. Additionally gels intended to be used for spot excision were loaded with an additional 350 g of an unlabeled, pooled sample sufficient for tryptic digestion and detection of the peptides by LC-MS 2 .
First Dimension-Each minimal CyDye was reconstituted in fresh N,N-dimethylformamide, and a 400-pmol quantity was used to label 50 g of protein at pH 8 -9. Cy2 was used to label the pooled internal standard as described above. Cy3 and Cy5 were used to label the normal and tumor samples, and we alternately swapped the dyes on subsequent sample pairs to alleviate dye-specific effects that could bias image analysis. Labeling proceeded for 30 min in the dark and was quenched with 10 mM lysine. Samples were then mixed with an equal volume of 2ϫ sample buffer (8 mM urea, 4% CHAPS, 2% DTT, 2% Pharmylyte, pH 3-10 or 4 -7, non-linear), placed on ice for 10 min, then loaded onto non-linear pH 3-10 (or 4 -7) Immobiline DryStrips (GE Healthcare), placed in a strip holder, and focused with an IPGphor system using a step gradient protocol ranging from 30 to 8000 V for approximately 27 h. The strips were then stored at Ϫ80°C, ready for the second dimension. Additionally for the first experiment (pI 3-10), two pooled, unlabeled 350-g samples were prepared and focused separately to be subsequently separated in the second dimension on separate gels intended for spot excision. By contrast, for the second experiment (pI 4 -7), the unlabeled, pooled sample was mixed with the labeled samples and run on the same gel. This is possible because the Deep Purple gel stain (GE Healthcare) we used to stain the unlabeled sample is spectrally resolvable from the CyDye fluorophores. This reduces the number of gels required for the experiment.
Second Dimension-For separation by molecular weight we used the Ettan DALT Twelve apparatus. The DryStrips were rehydrated in 10 l of re-equilibration buffer (8 M urea, 100 mM Tris-HCl, pH 6.8, 30% glycerol, 1% SDS, 45 mg/ml iodoacetamide (to reduce streaking)) for 10 min, laid across the top of a homogeneous 12.5% polyacrylamide gel "sandwiched" between two glass plates submerged in running buffer (40% bisacrylamide, 1.5 Tris, 10% SDS, 10% ammonium persulfate, 10% TEMED), and then covered with a 0.5% agarose solution. Separation proceeded at 15°C at 0.5 watts/gel and then 1.0 watt/gel for 15 h. Separation was stopped when the bromphenol dye front reached the bottom of the gel.
Gel Fixation-Gels to be used for spot excision were previously secured to one glass plate in the "sandwich" with silane (Bind-Silane, GE Healthcare). After the experiment was stopped, these gels were fixed in 50% methanol and 7.5% acetic acid and subsequently stained with Deep Purple Total Protein Stain (GE Healthcare).

Image Analysis
Gels were scanned using a Typhoon 9400 variable mode imager (GE Healthcare). During this phase each CyDye fluorophore is independently excited by laser light specific to its particular excitation spectrum. Emission sensitivity, i.e. photo multiplier tube, was adjusted until the most intense spot on the gel approached saturation. This tuning was performed at a 1000-m resolution, and once the photo multiplier tube value was optimal, a final high resolution scan was performed at 100 m. The gel(s) intended to be used for spot excision was poststained with Deep Purple (GE Healthcare), imaged using 532/560 nm wavelength light, and then stored in the dark at 4°C. In general, by following the recommended settings for the Typhoon imager, our experience indicates that dye-specific biases in spot intensity are eliminated or reduced below significance. After imaging the three fluorophores for each gel, the images were imported to the DeCyder image analysis software (GE Healthcare) for spot detection, spot matching (intragel), and determination of statistically significant biological variation (intergel) based on the measurement of relative abundance change after background subtraction and normalization to the internal standard. Typically about 90% of the spots on a gel will fall within 2 standard deviations of the mean and not show a significant -fold change (the null hypothesis), although this is highly sample-dependent.
Image analysis is a time-consuming part of this experiment. A statistically significant spot is one whose mean -fold change is greater than or equal to Ϯ50% (depending on statistical power) and paired t test is less than or equal to 0.05. Each spot that passes significance must then be manually checked to ensure that it is likely a protein spot and not a gel artifact. Satisfying these criteria, a pick list is generated and exported to the software controlling the Ettan robotic spot picker (GE Healthcare). Spots were excised with a 3-mm core from the poststained gel and loaded to a 96-well plate for digestion.

In-gel Digestion
Excised gel plugs were washed four times for 10 min with 50 l of both 25 mM ammonium bicarbonate (ABC) and 50% ACN, removing the liquid between each wash. 10 mM DTT freshly prepared in 30 l of 25 mM ABC was added to each gel plug. The samples were then incubated for 45 min at 56°C. Following incubation, gel plugs were cooled at room temperature for 20 min. The DTT was removed, and 30 l of 55 mM iodoacetamide was added. The samples were incubated in the dark for 45 min at room temperature. The iodoacetamide was removed, and samples were washed four times with 50 l of both 25 mM ABC and 50% ACN. Gel plugs were covered with 10 l of a 100-ng trypsin solution, incubated at room temperature for 10 min to allow absorption of trypsin, then covered with 15 l of 25 mM ABC, and placed in a 37°C water bath overnight for digestion. The reaction was quenched the following day with 7 l of 1% (final concentration) formic acid. Extraction of the peptides from the gel plugs was completed by adding 30 l of 50% ACN, 5% formic acid, vortexing for 30 min, spinning the samples, and finally sonication for 5 min.

Mass Spectrometry and Database Software
Most samples were analyzed by tandem mass spectrometry using an LTQ mass spectrometer (Thermo Electron Corp., Bremen, Germany) equipped with an Ettan multidimensional LC system (GE Healthcare). Six samples were run on a Finnigan LTQ FT hybrid mass spectrometer (Thermo Electron Corp.) operated in positive ion mode. 2.5 l of tryptic peptides were desalted on a C 18 pre column (PepMap 100, 300 ϫ 5-m particle size, 100 Å, Dionex) and then separated on a reverse-phase column (C 18 , 75 m ϫ 150 nm, 3 m, Dionex) using mobile phases A (0.1% formic acid) and B (84% acetonitrile, 0.1% formic acid) with a linear gradient of 2%/min, beginning with 100% A. Peptides were subsequently infused at a flow rate of 300 nl/min via a Pico Tip emitter (New Objective, Inc., Woburn, MA) at a voltage of 1.8 kV. Mass spectra were recorded in the ion trap, and MS 2 spectra were acquired for the five most intense ions in the LTQ using a collision energy of 35 eV and an isolation width of 2.5 Da.
Bioworks version 3.2 (Thermo Electron Corp.) using the SEQUEST software was used to search against an indexed human database with a peptide mass tolerance of 2.5 Da and a fragment tolerance of 1.0 Da. Search parameters included partial methionine oxidation, complete carbamidomethylation of cysteine, and two missed cleavage sites. Statistically significant peptides were those satisfying p Ͻ 0.001 and cross-correlation (Xcorr) values of 1.9, 2.5, and 3.0 for 1ϩ, 2ϩ, and 3ϩ charged ions, respectively. The protein probability cutoff was p Ͻ 0.001, and each "hit" necessarily required at least three peptides for consideration with rare exceptions for low molecular weight proteins. Surviving that filter, each protein call was manually "rationalized" to the gel; that is the theoretical pI and molecular weight were compared with the observed values on the gel image. Cleavage products and post-translational modifications were considered in this step.

Statistical Power Analysis
The power of a statistical test (1 Ϫ ␤) is a measure of the probability of correctly rejecting the null hypothesis, H 0 , if it is false. Low power studies consequently have a higher rate of false negatives. Formally ␤ is functionally related to sample size (n), the standard deviation of the distribution (), the difference in the means being tested ( Ϫ 0 ), and the area under the standard curve at a given significance level (␣) (Z 100 ).
Prior to our second experiment we estimated the average spot variance () by considering all spots on all 12 gels under the assumption that the source of variance was primarily biological and that experimental variance was relatively minimal. Because the samples had been prepared, labeled, and separated at once under near identical conditions we thought this assumption reasonable. Next, at a fixed level of significance (␣ ϭ 0.05) we calculated the relationship between power and fold change ( 0 Ϫ 1 ) at three different sample levels (n ϭ 3, 6, and 12). This provided an estimate of the minimum number of paired samples required to measure a particular minimum fold change with a power of 0.8 (see supplemental Data S1).

Gene Expression Data
mRNA expression was measured by cDNA microarray on 171 human colon tissue samples of various stages of colon cancer (normal ϭ 16, stage B ϭ 41, stage C ϭ 25, stage D ϭ 50, metastatic ϭ 39) using the Affymetrix Human EXON 1.0 ST chip. Expression values were generated with the Expression Console program from Affymetrix (Affymetrix, Santa Clara, CA) using the probe logarithmic intensity error (PLIER) algorithm to minimize the effect of outliers. Expression values for all 171 samples for select genes in our networks, plus the decoy database of 1000 genes, were generously provided to us by the Case Comprehensive Cancer Center. The decoy genes were randomly chosen from Ͼ17,800 probe sets with core evidence. The distribution of the decoy was evaluated to ensure representation across all 23 (1-22 plus X) chromosomes (supplemental Data S3) and was verified not to overlap any genes in our four networks. The data set is available upon request.

Protein Interaction Network Database and Subnetwork Build Algorithm
We used MetaCore from GeneGo Inc. (version 4.6 build 12332) to search for protein-protein interaction networks. MetaCore uses a protein interaction database comprising tens of thousands of protein interactions that have been manually curated based on a thorough reading of evidence reported in the literature. MetaCore covers 2400 journals and does not use natural language processing algorithms. In essence, the database represents a vast wealth of legacy experimental evidence that can be quickly mined in a number of ways for proteins and interactions relevant to a particular disease. These data can be usefully represented by directed graphs ("networks") that illustrate not only which proteins interact with each other but the functional nature of the interaction between them (binding, cleavage, phosphorylation, etc.). We will use the term "network" to refer to the entire database of protein interactions, and "subnetwork" will mean any network smaller than the whole network. There are a number of algorithms available in MetaCore with which to build subnetworks around a set of differentially expressed targets ("seed"). We chose an algorithm that would extend a subnetwork around our seed while minimizing the number of outgoing and/or incoming connections needed to enclose the seed in a "cloud" of interactions topologically constrained by the shortest path. As this subnetwork is likely to be very large, it is subsequently divided into smaller subnetworks by maximizing the saturation of the seed targets in each while obeying our input constraint of subnetwork size (n ϭ 50). We further constrained the search by species (human) and tissue type (colon); all other prefilter options retained their default values. The end result was a list of subnetworks (13) ranked by p value and zScore. The reported p value is calculated assuming a hypergeometric distribution, and it represents the probability of a particular mapping arising by chance given the numbers of genes in the set of total networkable genes (i.e. genes or network objects that have at least one annotated functional interaction), all genes on maps/subnetworks/processes, genes on a particular map/subnetwork/process, and genes in our experiment. The zScore is a statistical measure of the concentration of the seed targets in the subnetwork. zScore ϭ where N equals the total number of nodes in MetaCore database, R equals the number of the objects of the subnetwork corresponding to the genes in the import list, n equals the total number of nodes in each subnetwork generated from the import list, and r equals the number of nodes with data in each subnetwork generated from the import list. A white paper providing additional details of network construction algorithms is available upon request.

Network Scoring and Significance Tests
A flow chart of the scoring scheme is outlined in Fig. 5. We obtained global gene expression data as measured by cDNA microarray (for detail see "Gene Expression Data") for every gene product in each of the four subnetworks chosen for analysis. Importantly although the search criteria constrained each subnetwork to 50 proteins or less, some proteins in the subnetworks were complexes involving multiple gene products, subunits, or isoforms. We chose to include all these gene products when we exported the subnetwork list of genes from MetaCore. Consequently and depending on the particular subnetwork, the number of genes may exceed 50. Additionally the microarray data set also included five experiments that had used microdissected epithelial cells from colonic crypts. We considered these to be controls for the normal tissue samples because of the homogeneous cell type, whereas the normal tissue samples had detectable levels of stromal markers (e.g. vimentin). Hence genes with an average expression value less than 40 (below the detection limit of quantitative PCR) across the crypt samples were considered unexpressed in the epithelium layer and removed from consideration during scoring.
Mutual Information-MI is a concept from information theory used to measure the dependence of two random, discrete variables, say X and Y, based on their joint and marginal distributions. A high MI score (0 Յ MI Յ 1) indicates that X and Y are non-randomly associated to each other, whereas in the limit an MI value of 0 indicates that the two variables are statistically independent. To apply the concept to our problem, we retrieved the mRNA expression for every gene product in each of the four networks and used these data to populate two distributions of network activity values, one for a set of normal samples (X) and one for a set of stage D samples (Y). With these two distributions we were able to calculate an MI score between normal and stage D for each network and also use it as an optimization metric to search within the network for combinations of proteins (signatures) that would maximize this score. Intuitively a high MI score would indicate that the corresponding proteins non-randomly associate between normal and stage D, inferring their functional importance in the network in late stage cancer and suggesting new experiments for elucidating mechanisms of tumorigenesis.
The raw mRNA expression values in each subnetwork were first normalized by subtracting the population mean and dividing that difference by the population standard deviation. Next a subnetwork activity score was determined for each sample (column) in each network by summing the corresponding normalized mRNA expression values (rows). The activity score across phenotypes is a random, continuous variable that we discretized using a binning procedure. The number of bins (Fig. 5, step 3) was determined by log 2 (number of samples) ϩ 1, which is Sturge's rule, in the same manner as described by Ideker and co-workers (15). The range of the bins was determined from the range of the normalized expression values plus or minus a small adjustment to ensure all values fell within a bin. Finally the marginal and joint distributions were determined for the two phenotypes being compared (normal and stage D), and the mutual information value was computed (Fig. 5, steps 4 and 5). Note that the example in Fig. 5 indicates that network activity values were calculated over every patient sample (column), which as a matter of course we did do, but we only calculated the MI between normal and stage D consistent with the proteomics comparison. The other values, however, were useful for testing hypothesis 2 (H2; see Fig. 6).
Significance Testing-To test the hypothesis that the genes in a given subnetwork were not significantly discriminative of phenotype compared with a random selection of n genes (where n equals the number of genes in the subnetwork being measured), we randomly selected 1,000,000 combinations of n genes from the decoy database to create the null distribution and then evaluated the actual MI score of the subnetwork on the cumulative distribution function (CDF). To test the second hypothesis, which is that the genes in our network do not associate with a particular phenotype, we permuted the phenotypes (columns) of the relevant array 100,000 times and evaluated significance in the same way. The evaluation on the CDF is expressed as a percent value. A value of, say, 95% indicates there is a 5% chance (p ϭ 0.05) of observing a higher MI value, assuming the null hypothesis is true. The programs required to import the expression data, organize it for analysis, visualize it, and perform the scoring and optimization search as well as the hypothesis testing were written using Matlab and are available on request.

Label-free Mass Spectrometry: Protein -Fold Change Determination
Sample Preparation-50 g of total protein derived from colonic tissue lysate was precipitated with acetone (Ϫ20°C) at Ϫ80°C for 20 min followed by centrifugation at 12,000 rpm for 10 min at 4°C. The samples were then dried, rebuffered in 20 l of 0.2% ProteasMax surfactant (Promega Corp., Madison, WI), and gently shaken for 30 min. The buffer volume was then increased to 93.5 l with 50 mM ammonium bicarbonate and incubated with 1 l of 0.5 M DTT for 20 min at 56°C followed by incubation with 2.7 l of 0.55 M iodoacetamide for 15 min in the dark. 1 l of 1.0% ProteasMax was added to the buffer followed by 1.8 l of trypsin (Promega Corp.) that had been dissolved in 50 mM ammonium bicarbonate to a final concentration of 1 g/l. Digestion was carried out for 3 h at 37°C. Peptides were concentrated on a 100-l C 18 UtraMicroTip column (Net Group, Inc.), eluted in 20 l of 0.1% formic acid in 60% acetonitrile, then diluted with UltraPure water to a final concentration of 500 ng/l, and stored at Ϫ80°C.
LC-MS/MS-Analyses were performed on an LC Packings/Dionex Ultimate 3000 HPLC-Orbitrap XL (Finnigan, San Jose, CA) system. The HPLC system is equipped with two independent ternary gradient pumps suitable for high throughput dual column parallel HPLC mode applications. A standard injection volume of 10 l was used for all the samples, giving a total of 1 g of digest on the column. The data collection method incorporated a 30,000 resolution Orbitrap full scan in the FT mode using profile mode data collection followed by data-dependent mode MS 2 acquisition of the top five precursors from each full scan (centroid mode). CID mode fragmentation was chosen for generating the MS 2 spectra in the linear ion trap with a standardized value of 30% normalized collision energy being chosen for the fragmentation of the peptides. The LC method included a slow, 95-min acetonitrile ramp from an initial 4.8% until a final composition of 50.2% was achieved. Further high organic elution was performed (5 min) to complete the elution of peptides from the analytical column followed by equilibration of the column for succeeding analyses.
Analysis-Raw files were searched by Mascot against the IPI_Human database (ipi.HUMAN.v3.28.fasta). For each protein of interest one tryptic peptide with tandem MS evidence in at least one of six replicate samples was used to measure the relative expression change between normal and tumor. Peptide abundance was determined by the area under the elution curve that was extracted from the total chromatogram using a mass window adequate to capture all isotopes of the observed monoisotopic mass (Յ10 ppm). The curves were smoothed by an 11-point Gaussian filter and base line-subtracted using the Xcalibur software (Thermo Electron Corp.). -Fold change between normal and tumor was calculated as the ratio of the integrated curve areas. Three replicate runs of a single sample pair (patient 507, normal/tumor) were used to estimate technical variance. The coefficient of variation for -fold change ranged from 6 to 39%. The -fold change for IGFBP3 was determined by densitometry from Western blot analysis.

2D DIGE Discovery Proteomics
The features of our cohort and the overall design of our proteomics experiments are shown in Fig. 1. Twenty-four tissue samples (12 normal and 12 tumor, each pair from the same patient) were prepared using standard procedures (see "Experimental Procedures"). The tissue samples had been vetted by a pathologist to establish tumor grade. The experimental design involved alternately labeling tumor and control samples with Cy3 and Cy5 dyes, whereas Cy2 was used to label a 50-g pooled fraction that served as an internal standard for each gel. The usefulness of this standard cannot be overemphasized as it assists in providing a confident assessment of real biological variation by controlling for variance in protein loading (16). Each tripartite sample was first separated by IEF over a broad pH range (3)(4)(5)(6)(7)(8)(9)(10) and then by molecular weight using 12.5% homogenous SDS-polyacrylamide gels. All the samples were labeled and separated simultaneously under identical conditions to minimize experimental variation (see "Experimental Procedures"). Each gel yields three images (Fig. 1, step 4), and along with the poststained gel used for spot excision, a total of 37 images were imported to the DeCyder software (GE Healthcare) for differential image analysis (i.e. spot matching) followed by statistical analysis of biological variation. Fig. 2 is a Cy5 image representative of a typical analytical gel (patient 5144) indicating the significant spots matched. In total for this experiment 58 spots were identified as significantly (p Յ 0.05, mean fold change Ͼ50%) changing between normal and tumor. For the majority of these spots DeCyder was able to detect a match on greater than 30 of the images. In no case did that number fall below 20 or fail to be matched on the poststained gel used for spot excision. The spots were robotically excised from gel, digested by trypsin overnight, and submitted to reverse-phase LC-MS 2 followed by database search. Twenty-three spots were confidently (p Յ 0.001, peptide and protein) identified by database search (Table  I and Fig. 2, annotated). Thirteen proteins were up-regulated in cancer, and seven were down-regulated.
The IEF range chosen for this experiment resulted in spot overlay in a number of regions of interest. We anticipated that we could improve separation and focus more spots by using an IEF range of pI 4 -7. Additionally using a measure of spot variance from the first experiment, we found we only needed six sample pairs to capture -fold changes greater than Ϯ30% while maintaining statistical power of 0.8 (see "Statistical Power Analysis" under "Experimental Procedures"). Accordingly we performed a second experiment using a subset of six sample pairs from the first experiment (numbers 145, 321, 362, 468, 480, and 602). The protein fractions were thawed and cleaned with a kit (2D-Cleanup kit, GE Healthcare) to remove impurities known to interfere with proper separation. The protein concentration was redetermined as before by colorimetric assay. We used even more stringent criteria in the image analysis phase as compared with experiment 1. Spots not only had to satisfy the same statistical criteria, but additionally, for a spot to be considered for picking, it had to have been matched by DeCyder on every one of the 19 gel images. Using these criteria we identified 150 significantly (p Յ 0.05, mean -fold change Ͼ30%) changing spots (Fig. 3). Activating the false discovery rate filter in DeCyder (based on the method of Benjamini and Hochberg (17)), over 40 of these spots retained their significance (q Յ 0.05).The indicated spots were excised from the gel, digested by trypsin, and submitted to reverse-phase LC-MS 2 followed by database search. Of 150 spots, we confidently identified 67 proteins. Thirty-five spots were up-regulated in cancer, and 32 were down-regulated (Table II), indicative of a lack of bias toward upor down-regulated proteins.  (1); unlabeled samples were run on separate gels for poststaining or mixed in the analytical gels (see "Experimental Procedures"). Samples were separated by isoelectric focusing (2) and then by molecular weight (3), and each fluorophore was imaged independently (4). Using the DeCyder software, spots were matched on an intragel basis with differential image analysis (DIA) (5) and on an intergel basis with biological variation analysis (BVA) to assess biological variation (6). Significant spots (⌬fold dependent on statistical power) were selected for robotic excision (7) and digested by trypsin (8), and the peptides were separated by reverse-phase (RP) chromatography and detected by tandem mass spectrometry (9). MS 2 spectra were searched using SEQUEST (10), and identified proteins were imported to MetaCore to search for relevant networks (11). Significant protein signatures are scored by mutual information using gene expression profiles (12); signatures are extended out one hop to infer functional relevance (13). Resultant subnetworks were analyzed for biological relevance to CRC and new hypothesis generation (14).

Highly Significant Late Stage CRC Signatures
As stated above, our guiding hypothesis was that the differentially expressed proteins ("targets") found in our experiments represented nodes upstream and/or downstream of other nodes in one or more functional subnetworks that may be dysregulated in stage D colon cancer. Accordingly we used the set of unique targets from both experiments (n ϭ 67) to seed a search for functional subnetworks. A detailed description of the protein interaction database we searched and the subnetwork construction algorithm we used is provided under "Experimental Procedures." The search returned 13 subnetworks, each of which contained a variable number of between one and 12 of the seed targets. We limited our attention to four subnetworks judged most significant by a combination of p value and zScore. One of these subnetworks is annotated in Fig. 4 with a breakdown of the most significant gene ontological processes, their percent representation in the subnetwork, the p value, zScore, and subnetwork size, i.e. the total number of gene products it contains. The remaining three are provided in supplemental Data S4.
To test our hypothesis that our 67 proteomic targets were significant for late stage CRC, we implemented a quantitative method to score the subnetworks. Fig. 5 outlines our approach. A more detailed explanation of the scoring procedure is provided under "Experimental Procedures." Briefly we obtained mRNA expression data for every gene product in each of the four subnetworks from a set of unpublished microarray experiments (Affymetrix) performed on a large cohort of clinical tissue samples of varying CRC stage (Fig. 5, step 1). With these data we computed an activity value for each subnetwork over each normal experiment (n ϭ 16) and each stage D experiment (n ϭ 50) (Fig. 5, step 2). The activity values represent a continuous variable and for the purpose of computing MI need to be discretized by a binning procedure (Fig. 5,  step 3). With discrete values we computed the relevant distributions (Fig. 5, step 4) and then the MI between normal and stage D for each subnetwork. The scores are shown in Table  III (last column).
There are two null hypotheses relevant to evaluate the significance of the MI score. The first, which we will call H1, states that genes in a subnetwork are not discriminative of  Table I. the disease phenotype compared with a random set of genes. For example, if the subnetwork contained 10 genes, then any 10 genes taken at random would produce an MI score at least as good as the real subnetwork of 10 genes under the null hypothesis. The second, call it H2, is subtly different; it states that the expression levels of the genes in the subnetwork do not associate with a particular phenotype. For example, if the network contains 10 genes that produce a high MI score between normal and stage D, then under the null hypothesis scores at least as high will be found for random permutations of phenotypes, i.e. by disrupting the real association between patient and gene expression.
From the microarray we obtained mRNA expression data for 1000 random genes ("decoys") and ensured that these genes had no overlap with any of the genes in the four subnetworks. A null distribution was estimated for testing H1 by evaluating an MI score for 1,000,000 combinations of n random genes selected from the decoy data set where n equals the size of the particular subnetwork being assessed. The null distribution for H2 was estimated from 100,000 permutations of the phenotypes (columns) in the two-dimensional array representing a subnetwork. Significance was then determined by evaluating the MI score on the CDF of the respective null distribution. 1 Ϫ CDF indicates the probability of finding a higher MI score. Probability values of 1% or less were considered to be significant.
By this measure neither the null hypothesis for H1 nor that for H2 could be rejected for any of the four subnetworks, i.e. when all the gene products in the subnetwork were used to compute its activity value (Table III, last column). We reasoned that this result could be attributed to how the subnetworks were built and scored. Although all four subnetworks were discovered by proteomics profiling and were judged statistically significant, the individual interactions in a subnetwork are nevertheless based in large part on a diverse set of experiments performed in vitro and in vivo in a variety of different tissues. Consequently although many of the subnetworks are indicated to be active in colonic tissue, they are not necessarily important to a metastatic cancer phenotype. Also we observed that the subnetworks, in terms of the nodes and edges they contain, were sensitive to the parameters used to search for them, even including the version of the database software. This results in a certain degree of arbitrariness in the overall topology of the subnetwork. Given these caveats the statistically insignificant MI scores were not very surprising.
However, we further hypothesized that the activity of specific protein combinations within the subnetwork(s) would be highly discriminative of disease. Thus, we performed an exhaustive combinatorial search over each subnetwork for combinations that would maximize the MI score between control and stage D. We did this for up to six combinations (readily accomplished on a conventional desktop computer). For each subnetwork except one (subnetwork 2), the MI score steadily increased with combination number (Table III). MI scores for these specific combinations (signatures) were much higher than those for any of the subnetworks taken as a whole, and more importantly, in each case the MI score was highly significant with respect to both null hypotheses, H1 and H2 (Fig. 6, compare with MI in column labeled "signature 6" in Table III). Notably the top scoring signatures included proteins for which we had independently found direct proteomics evidence (e.g. CCT2, HSP90AB1, SERPINA1, and CapG) as did certain other signatures of fewer combinations. Fig. 7 highlights the gene products (gray) from each subnetwork participating in signature 6. To extend the potential functional importance of these signatures, we added to each directed graph those proteins that were one hop away from the signature proteins, i.e. those from the corresponding parent subnetwork immediately up-or downstream. We briefly discuss the proteins and functional interactions of these expanded sub-networks in more detail in the following section. Additionally similar to the conclusion of Ideker and co-workers (15), we found that the signature genes did not cluster in a dendrogram computed using traditional distance metrics, e.g. Euclidean or Spearman (data not shown), and would likely have been overlooked by conventional gene classification techniques.
As mentioned above, with the exception of parent network 2, MI scores increased or were constant for successively larger combinations of proteins. Further we found that combinations of less than six proteins (signatures 1-5) were also significant when tested against the appropriate H1 and H2 null distributions (data not shown). Indeed the relevant distributions for smaller combinations of proteins had characteristics (e.g. mean and variance) similar to those for signature 6. If the proteins appearing in successively larger combinations were completely different this would suggest that our method was not sensitive to small variations in the underlying network activity patterns. However, this is not FIG. 3. Representative gel from experiment 2. Polygons indicate spots significantly changing between normal and cancer as determined by DeCyder. Magenta polygons indicate that the spot passed a multiple comparison filter test (false discovery rate) in DeCyder. Labeled spots were identified by mass spectrometry and appeared in one or more of four significant MetaCore networks. See Table II. what we observed. As indicated in Table III, for each network all signatures (1-6) frequently show the repeated contribution of either the same protein or proteins involved in the same complex of proteins, suggesting that the method is sensitive to the specific interactions of functionally similar subnetwork proteins. Even for the signatures derived from parent subnetwork 2, the contribution of proteins capable of phospholipase activity (PLA2 and PRDX6) appears in five of six signatures.

Biological Relevance to CRC of Signature Proteins in Extended Subnetworks
It merits emphasis that each of the most significant extended subnetworks (Fig. 7) contained targets for which we had direct proteomics evidence (CCT2, HSP90AB1, SER-PINA1, and CapG), indicating that gene products significant by their contribution to MI maintain their significance at the level of the proteome in late stage CRC. The most significant targets (highlighted gray) are generally classified according to those with a known role in CRC or a role in other human cancers and those with no known role in cancer. The fact that we found significant genes with a known role in CRC can be understood as a positive control of our analytical method.
Genes with a Role in CRC-IGFBP3 (also known as IBP3) is an insulin-like growth factor-binding protein that was recently identified to cause apoptosis in a tumor necrosis factor-related apoptosis-inducing ligand (TRAIL)-mediated fashion in relevant CRC cells (18). Additionally a large association study found that paired polymorphisms in IGFBP3 and its substrate predicted a significant increase in risk for CRC (19). The integrin family of proteins has been well studied in CRC. They are generally responsible for cell-cell and cell-matrix adhesion. Loss of expression of certain subunits in this family has been associated with increased neoplastic transformation in colonic epithelium, and the specific loss of ␤1 (ITGB1) chains was associated with benign to malignant transformations (20). Notably integrins are active as heterodimers, and although only ITGB1 contributed to the significant MI score, the subnetwork indicates that the dimer ITGB1/ITGA4 is of particular interest. IFITM1 (IFI17) is a member of the family of interferoninducible transmembrane proteins. A recent study (21) proposed it as a possible marker for human colorectal tumors. The subnetwork also revealed it to be regulated at the level of transcription by the PBAF complex, certain members of which we also found significant. SMARCA4 (also known as BRG-1) plays a key role in the chromatin-remodeling complex in mammals. One study (22) showed in vivo evidence that BRG-1 interacts with ␤-catenin and induces the transcription of T-cell factor (TCF) target genes. Mutant forms of BRG-1 lacking ATPase activity disrupted this induction. TCF target gene activation is the final consequence of the WNT signaling pathway, mutations in which are well known to cause tumorigenesis in the human colon. The role of platelet-derived growth factor receptor (PDGFR) has been well studied in CRC along with other receptors capable of tyrosine kinase activity and downstream signaling. One recent study (23) found a significant association between the stromal expression of the B subunit (PDGFRB) and the metastatic potential of CRC tumors. In addition to being overexpressed in a number of human cancers, casein kinase II (CSNK2A2) in CRC suppresses apoptosis by desensitizing cells to TRAIL in a caspase-dependent manner but independent of NF-␤ (24). It has also been shown to promote survival of colon cancer cells by increasing the expression of survivin via the canonical transcription pathway hyperactive in CRC (TCF/ LEF (lymphoid enhancer binding factor)) (25). PLA2G12A is a member of the family of secreted phospholipases, many of which display distinct patterns of expression in adenocarcinomas (26). Genes with a Role in Other Human Cancers-CapG, a gelsolin-like capping protein, has been identified as a possible tumor suppressor gene (27), although our proteomics screen revealed it to be up-regulated in cancer in agreement with the mRNA expression. A closer look at this study revealed that the authors had measured a near complete loss of the CapG protein in a variety of primary human cancer tissues but not colon tissue. Our evidence that the CapG message and protein are up-regulated in CRC indicates that it may have oncogenic activity in the colon. Further we actually identified CapG at two closely spaced but distinct spots on the gel, suggesting that post-translational modification may be important to its activity. The human gene PLK1 (or PLK), a serine/threonine protein kinase, was characterized many years ago (28), and its expression was found to strongly correlate with the mitotic activity of a variety of tumor cell lines, including those derived from human colon. Notably the study found that PLK1 was not expressed in a variety of the normal human tissues with the exception of normal colon tissue. More recently, PLK1 was found to be overexpressed in primary CRC tumors (29), identified as a prognostic factor for CRC (30), and when knocked down or inhibited in human adenocarcinoma cells (RKO) lead to dramatic mitotic arrest (31), thus showing promise as a possible drug target. Lastly driver mutations in PLK1 are not unknown as was recently revealed by a large screen for somatic mutations on over 500 protein kinases covering a large cohort of human cancers (32). RPS2 is a gene encoding a ribosomal protein in the small 40 S subunit. RPS2 was found by a proteomics screen to be a novel kinase substrate differentially phosphorylated in breast cancer development (33). the remaining significant subnetwork targets in human cancer (HSP90AB1, HIST1H2AB, TUBA4A, TUBB3, GNA12,  TRAP1, DYNLT3, CCT3, CCT5, CCT7, and POLR2D). Guided by the evaluation of interactions on the subnetwork along with select proteomics evidence, these targets may merit follow-on experiments to discover their role, if any, in late stage CRC. For example, the chaperone containing t-complex proteins (CCTs) play a role in protein folding in eukaryotes and are widely expressed in the cytosol. One study did find a significant elevation of the CCT transcript in human colon carcinomas and validated the change in protein expression by immunohistochemistry, whereas our proteomics screen revealed it to be down-regulated in the cancer tissue. Notably that study had not vetted the samples for tumor stage, highlighting the importance of stagespecific studies. Additionally CCT3 and CCT5 were identified by microarray analysis to be significantly differentially expressed in the epithelium of other human cancer tissues (esophageal, breast, ovarian, and lung), but their functional role in cancer, CRC in particular, is unknown.  (1). Normalized values were then summed to produce an activity value for each sample (2); activity scores are continuous and need to be assigned into discrete bins for an MI calculation (3). A joint distribution matrix is calculated between two sets of samples, e.g. normal (N) and stage D (4). MI is calculated as shown where p(x) is the marginal distribution of normal activity, p(y) is the marginal distribution of stage D activity, and p(x,y) is the joint distribution of x and y (5).

Significant Targets in the Developmental Process Subnetwork Are Coordinately Regulated
We used a label-free mass spectrometry approach (see "Experimental Procedures") to verify the relative expression of four of the significant targets in the developmental process subnetwork (Fig. 7, panel 1) in a new cohort of clinical tissue samples. The differential expression of IGFBP3 was determined by Western blot. The relative expression change between normal and stage D at the level of mRNA for each of these targets was computed from the microarray and used for comparison. Most all of the targets were up-regulated in cancer in all patients at the level of mRNA and protein (Fig. 8).
Overall these data indicate coordinated regulation at the level of mRNA and protein but also highlight the relatively large variation of expression of both mRNA and protein across patients. An interpretation consistent with this observation is that subtle changes in the transcription of one or more targets may have a synergistic effect on the activity of the other targets in maintaining the phenotype, something the measure of mutual information is well suited to capture and that is consistent with our guiding hypothesis.

The Advantages of an Integrated -Omics Approach to Cancer
Colon cancer has a strong genetic basis due to the accumulation of somatic mutations in oncogenes and tumor suppressor genes. However, it is also widely accepted that because of the resiliency of mammalian cells single gene mutations are usually insufficient to cause this disease (34). Although a great deal of work has been done identifying genes involved in colon cancer as well as the canonical pathways to which they resolve (35)(36)(37), comparatively little work has been done to evaluate the functional protein interactions derived directly from proteomics data. It is in fact not known how genomics, transcriptomics, or proteomics perspectives may differently inform our understanding of colon cancer onset or progression. Classification of disease phenotype using candidate gene, candidate RNA, or candidate protein target approaches has been the bedrock of modern -omics research. However, in some cases these single gene/protein models of disease have been disappointing in follow-on studies (3). Alternative approaches, using network and subnetwork classifiers, are currently under examination. In this study, we searched for protein subnetworks by leveraging a database built on a very large number of legacy experiments using proteomics data as a seed to discover subnetworks discriminative of late stage CRC. It was thought that this approach would quickly lead us to significant protein-protein interaction subnetworks that would reveal the functional cause, or consequence, of stage-specific phenotype(s). We then developed a novel approach for searching within these subnetworks for particular signatures that are significant discriminators of the disease phenotype using a scoring process based on gene expression data. Computational and biostatistical methods to unify proteomics and gene expression data are a powerful way of identifying novel interaction signatures involved in late stage cancer. Although gene lists combined with rigorous statistical analysis can identify significant genes whose expression profiles cluster, the resultant gene lists alone provide no functional information of the post-transcriptional mechanism(s) of dysregulation. One criticism of our approach is that gel-based proteomics experiments, which provided the seed proteins for our search, typically identify highly abundant proteins as differentially expressed. Many of these are either so-called "housekeeping" genes with a role in metabolism or in any case may often be considered unimportant to a disease phenotype such as cancer as they may lack transcription factors or receptors as protein classes. However, our integrated approach was able to locate these seed proteins within regulatory subnetworks of great interest. Each of the four subnetworks we scored included between eight and 12 proteomic targets that were directly identified. This underscores the usefulness of a network-based approach that identifies specific and significant functional interactions possibly relevant to the pathophysiology of cancer. It also revealed the large diversity of subnetwork interactions in which these high expressers are evidently involved. However, as the scoring shows, the entire subnetwork(s) is not statistically significant for classifying phenotype. Only the end product of our quantitative approach, which identified FIG. 6. Estimated null distributions (probability density function) for hypotheses H1 and H2. For the H1 null distribution, an array of pseudo, sixgene signatures was computed from 1,000,000 combinations of genes randomly selected six at a time from the decoy data set. The activity value for each pseudosignature was computed between normal and stage D followed by the MI scores, which were then used to populate the distribution. For the H2 null distributions, as each signature (best six) comprises different genes, a separate H2 null distribution was computed for each. First we computed an array of 100,000 (100K) random permutations of phenotypes (n ϭ 171). Then using the six genes for each signature (Table III), we computed an activity value for 16 pseudonormals and 50 pseudostage D samples (refer to array in Fig. 5). We then computed 100,000 MI scores to populate the H2 distributions. H1 and H2 were modeled by the generalized extreme value distribution in Matlab.
FIG. 7. Expanded subnetworks from corresponding signatures. Signature 6 proteins were expanded by one hop inside the corresponding parent subnetwork(s) (Fig. 4 and supplemental Data S4) to infer functional relevance. Signature 6 proteins are highlighted gray; overlapping gray circles indicate multiple members or subunits of a complex participating in the signature. See also Table III, column labeled Signature 6. Horizontal lines demark cellular compartments. Panels 1-4 are pruned versions of the parent subnetworks, Fig. 4 and supplemental S4 (subnetworks 2-4), respectively. root nodes with functional interactions significant for phenotype, presents a focused set of testable hypotheses suitable for validation by perturbation experiments. Additionally we acknowledge that different protein interaction databases are likely to return different subnetworks given the same target seed. However, this is most likely because present day databases represent an undersampling of the human interactome, coverage of which has been recently estimated at less than 1% (38), and not because of any inherent arbitrariness attributable to our approach. Indeed as coverage of the human interactome continues to improve, interaction databases are likely to converge with respect to subnetwork selectivity.
Using the measure of mutual information to score the networks had an advantage over other classification methods in that there is no requirement that the underlying data be normally distributed. This made the method particularly well suited to examining gene expression data that, for many of the genes in our networks, exhibited non-normal distributions of expression for particular stages of cancer. Pairing this approach with exhaustive combinatorial search, versus a greedy search, reduces the possibility that the signatures represent a local rather than global maximum. A complete exhaustive search of the expression landscape for even larger combinations of gene products (Ͼ6) is limited only by computer power. Finally some of the proteins identified in our signatures did not, independently, have a significant change in mRNA expression at least not enough to be considered significant by simple gene expression profiling. But it is certainly conceivable that the cumulative effect of small changes in network activity (mRNA expression) may lead to significant changes in the proteome, and this was a guiding hypothesis of our study.
As high throughput methods continue to produce more genomics and proteomics data, it will become increasingly important to find new ways to integrate these data and to provide precise, quantitatively significant classifications of human disease stage. These classifiers will likely be critical to the assessment of individual phenotype important for development of personalized medicine.  507, 534, and 540). mRNA values are the difference between the means measured using the normal samples (n ϭ 16) and stage D samples (n ϭ 50) obtained from the microarray. Protein -fold change was determined by label-free mass spectrometry except for IGFBP3, which was determined by Western blot. As most targets were up-regulated in the tumor, carbonic anhydrase I (CA1) is included as a loading control to indicate that the observed up-regulation of most targets was not merely due to a greater amount of tumor digest on column versus normal.