System and Method for Classification of Patients

A system (100) for classifying a biological test sample, including a database (112) populated with reference expression data. The reference expression data includes expression levels of a plurality of molecules (polynucleotides or polypeptides), including a set of marker molecules, in a plurality of reference samples. Each reference sample has a pre-assigned value for each of one or more clinically significant variables. The system includes at least one processor (110) and at least one storage medium containing program instructions for execution by said processor (110). The program instructions cause the processor to accept (122) input expression data including a test vector of expression levels of the marker molecules in the biological test sample; and pass the input expression data to one or more analysis programs (130a, 130b, 35). The analysis programs include at least one statistical classification program (135) for assigning a value of at least one of said clinically significant variables to the test sample.

1-43. (canceled)

44. A method of determining the risk of breast cancer recurrence in a breast cancer patient, comprising the steps of:

isolating a nucleic acid sample from an isolated biological sample from the patient; and

testing and measuring expression levels of a set of nucleic acid marker molecules in the nucleic acid sample, wherein the nucleic acid marker molecules comprise the nucleic acids listed in Table 5, thereby determining the risk of breast cancer recurrence in the subject.

45. The method of claim 44, wherein the method further comprises testing additional clinical covariates for the subject selected from patient age, grade, nodes, tumour size or ER status.

46. The method of claim 45, further comprising calculating a prognostic index according to Formula 1:

PI=∑i=1200wixi=0.139601(grade)+0.64644(ER)+0.938702(nodes)+0.010679(size(mm))+0.023595(age)+0.243639.

47. The method of claim 44, wherein the expression data are generated using a platform selected from the group consisting of cDNA microarrays, oligonucleotide microarrays, microRNA (miRNA) arrays, and high-throughput quantitative polymerase chain reaction (qPCR).

48. The method of claim 44, wherein the expression data is at least partly assessed according to the distribution across reference samples of one or more statistics derived from reference data, wherein the statistics are selected from the group consisting of background intensity, percentage of molecules above detection threshold, ratio of 3′ expression level to 5′ expression level, slope of RNA degradation curve, normalization factor, and log(base 10) ratio of mean intensity to mean background intensity.

49. The method of claim 44 further comprising normalizing the distribution of the expression data to be comparable with the distribution of reference expression data.

50. A system for performing the method of claim 44, comprising:

at least one processor; and

at least one storage medium containing program instructions for execution by said processor, said program instructions causing said processor to execute steps comprising:

accepting input data in the form of expression levels of the set of the nucleic acid marker molecules in the biological sample; and

assigning a clinical annotation or pre-assigned value of one or more clinically significant variables to the biological sample,

wherein the clinical annotation is selected from the group consisting of anatomical system, tissue of origin, tumour subtype and risk of breast cancer recurrence; and wherein the one or more clinically significant variables is selected from the group consisting of disease state, disease prognosis, and treatment response.

51. A method of classifying a biological sample isolated from a cancer patient, comprising the step of:

isolating a nucleic acid sample from the isolated biological sample from the patient;

testing and measuring expression levels of a set of nucleic acid marker molecules in the nucleic acid sample to assign a clinical annotation to the biological sample, wherein the nucleic acid marker molecules are any combination of 100 or more genes in Table 4;

comparing the expression levels of the set of nucleic acid marker molecules in the biological sample to expression levels of the set of nucleic acid marker molecules in a set of reference samples, each member of the set of reference samples having a pre-assigned value for each of one or more clinically significant variables selected from the group consisting of disease state, disease prognosis, and treatment response;

the comparing of expression levels in the biological sample to the expression levels in the reference samples comprising one or more analysis methods, the analysis methods comprising at least one statistical classification method trained to distinguish among said pre-assigned values on the basis of that part of the reference data corresponding to expression levels of the nucleic acid marker molecules; and

classifying the biological sample for at least one of the clinically significant variables according to the pre-assigned values of reference samples using the statistical classification program.

52. The method of claim 51, wherein the at least one statistical classification algorithm is selected from the group consisting of k-nearest neighbors (kNN), linear discriminant analysis, principal components analysis, nearest centroid classification, and support vector machines.

53. The method of claim 51, wherein the one or more clinically significant variables are organized according to a hierarchy, wherein the levels of the hierarchy are selected from the group consisting of anatomical system, tissue type, and tumour subtype.

54. The method of claim 51, wherein the cancer is breast cancer.

55. The method of claim 51, wherein the expression data are generated using a platform selected from the group consisting of cDNA microarrays, oligonucleotide microarrays, microRNA (miRNA) arrays, and high-throughput quantitative polymerase chain reaction (qPCR).

56. The method of claim 51, wherein the expression data is at least partly assessed according to the distribution across reference samples of one or more statistics derived from reference data, wherein the statistics are selected from the group consisting of background intensity, percentage of molecules above detection threshold, ratio of 3′ expression level to 5′ expression level, slope of RNA degradation curve, normalisation factor, and log(base 10) ratio of mean intensity to mean background intensity.

57. The method of claim 51, further comprising normalizing the distribution of the expression data to be comparable with the distribution of reference expression data.

58. A system for performing the method of claim 51, comprising:

at least one processor; and

at least one storage medium containing program instructions for execution by said processor, said program instructions causing said processor to execute steps comprising:

accepting input data in the form of expression levels of the set of the nucleic acid marker molecules in the biological sample; and

assigning a clinical annotation or pre-assigned value of one or more clinically significant variables to the biological sample,

wherein the clinical annotation is selected from the group consisting of anatomical system, tissue of origin, tumour subtype and risk of breast cancer recurrence; and wherein the one or more clinically significant variables is selected from the group consisting of disease state, disease prognosis, and treatment response.

The present invention relates to classification of patients on the basis of expression of multiple biological markers. It is particularly suited to expression data from microarrays and other high-throughput platforms, although it will be appreciated that the invention may have wider applicability.

It has long been recognised that diagnosis and treatment of disease on the basis of epidemiologic studies may not be ideal, especially when the disease is a complex one having multiple causative factors and many subtypes with possibly wildly varying outcomes for the patient. This has recently led to an increased emphasis on so-called “personalised medicine”, whereby specific characteristics of the individual are taken into account when providing care.

An important development in the move towards personalised care has been the ability to identify molecular markers which are associated with a particular disease state or which are predictive of the individual's response to a particular treatment.

For example, in relation to breast cancer, the estrogen receptor (ER) or HER2/neu (ErbB-2) status of a tumour can be used in determining a patient's suitability for therapies that target these molecules in the tumor cells. These molecular markers are examples of “companion diagnostics” which are used in conjunction with traditional tests such as histological status in order to guide treatment regimes.

In cancer cases where a tumour has metastasized, it is important to determine the tissue of origin of the tumour. The current diagnostic standard in such cases includes imaging, serum tests and immunohistochemistry (IHC) using one or more of a panel of known antibodies of different tumour specificity (Pavlidis et al, Eur J Cancer 39, p 1990 (2003); Burton et al, JAMA 280, p 1245 (1998); Varadhachary et al, Cancer 100, p 1776 (2004)). For approximately 3-5% of all cases, known as Cancer of Unknown Primary (CUP), these conventional approaches do not reach a definitive diagnosis, although some may eventually be solved with further, more extensive investigations (Horlings et al, J Clin Oncol 26, p 4435 (2008); Raab et al, Cancer 104, p 2205 (2005)). The range of tests able to be performed can depend not only on an individual patient's ability to tolerate potentially invasive, costly and time consuming diagnostic procedures, but also on the diagnostic tools at the clinician's disposal, which may vary between hospitals and countries.

To date, most diagnostic protocols are primarily reliant on microscopy, single gene or protein biomarkers (IHC) and imaging techniques such as MRI and PET Scan. Unfortunately, these techniques all have limitations and may not on their own provide adequate information to diagnose widely metastasized tumours, poorly differentiated malignancies, rare subtypes or unusual presentations of common cancers.

It has been hypothesized that the information gained from gene expression profiling can be used as a companion diagnostic to the above protocols, helping to confirm or refine the predicted primary origin in a focused and efficient manner.

Since the advent of various robotic and high throughput genomic technologies, including RT-PCR and microarray, several groups (van Laar et al, Int J Cancer 125, p 1390 (2009); Rosenfeld et al, Nature Biotechnology 26, p 462 (2008); Tothill et al, Cancer Res 65, p 4031 (2005); Bloom et al, Am J Pathol 164, p 9 (2004); Monzon et al, J Clin Oncol 27, p 2503 (2009); Ramaswamy et al, PNAS 98, 15149 (2001)) have investigated the use of gene expression data to predict the primary origin of a metastatic tumor. Prediction accuracies in the literature range from 78% to 89%.

A number of gene expression based, commercial diagnostic services have arisen since the sequencing of the human genome, offering a range of personalized diagnostic and prognostic assays. These services represent a significant advance in patient access to personalized medicine. However the requirement of shipping fresh or preserved human tissue to an interstate or international reference laboratory has the potential to expose sensitive biological molecules to adverse weather conditions and logistical delays. In some parts of the world it may also be prohibitively expensive to ship human tissue to a reference laboratory in a timely fashion, thus limiting access to this new technology.

Most current commercially available gene-expression based cancer tests use a proprietary “diagnostic” microarray or PCR-based assay (van Laar et al; Rosenfeld et al; Dumur et al, J Mol Diagn 10, p 67 (2008)). Such arrays allow assaying of a small set of genes chosen for a particular purpose and are custom manufactured for that purpose. Because of the limited set of genes that are quantified by these existing assays, the data generated generally cannot be used for multiple diagnostic or prognostic analyses if a different set of genes is required. Furthermore, whatever data is generated, it is generally not accessible to the clinician requesting the test should it be desired to conduct further investigations or compile a custom database of gene expression data for research purposes.

In view of the above deficiencies, it is desirable to provide a more flexible and efficient method and system for diagnosis and prognosis of a patient based on expression of multiple biological markers.

Accordingly, in a first aspect, the present invention provides a system for classifying a biological test sample, including:

a database populated with reference expression data, the reference expression data including expression levels of a plurality of molecules (polynucleotides or polypeptides) in a plurality of reference samples, the molecules including a set of marker molecules, each reference sample having a pre-assigned value for each of one or more clinically significant variables;

at least one processor; and

at least one storage medium containing program instructions for execution by said processor, said program instructions causing said processor to execute the steps of:

accepting input expression data, the input expression data including a test vector of expression levels of the marker molecules in the biological test sample; and

passing the input expression data to one or more analysis programs, the analysis programs including at least one statistical classification program which has been trained to distinguish among said pre-assigned values on the basis of that part of the reference data corresponding to expression levels of the marker molecules; and

assigning one of said pre-assigned values to the test sample for at least one of said clinically significant variables using the statistical classification program.

By providing a reference data set with known clinical annotation in a single database in combination with the ability to accept input data from a user of the system, it is possible to have a centralised repository of disease classification which can be used to conduct different diagnostic or prognostic analyses (using different classification programs) with different sets of marker molecules. The system thus provides flexibility in that different tests may be conducted using the same reference data and input data without needing to re-assay the biological test sample.

Preferably, one of said analysis programs includes instructions for assessing the quality of the input expression data. The quality of the input expression data may be at least partly assessed according to the distribution across reference samples of one or more statistics derived from the reference data, the statistics including (for example) background intensity, percentage of molecules above detection threshold, ratio of 3′ expression level to 5′ expression level, slope of RNA degradation curve, normalisation factor, and log(base 10) ratio of mean intensity to mean background intensity.

Providing a quality control module as one of the analysis programs allows the clinician or other user to check that the data, as a whole, fall within acceptable ranges so that low-quality data are not passed to the classifier or classifiers. Use of low-quality data could lead to a diagnosis which is inconsistent with other tests which may have been conducted, such as imaging or immunohistochemistry.

One of the classification programs may be a prediction of patient gender. This serves as a further quality check since, for example, if a female patient is predicted as male (on the basis of comparison to the reference data which are stratified according to gender), the fidelity of the data is cast into doubt.

In one embodiment, one of said analysis programs includes instructions for normalising the distribution of the input expression data to be comparable with the distribution of the reference expression data. This can help to increase the likelihood that differences between the input and reference data are due to real biological differences, and not due to mere statistical artifacts or to differences in the laboratory protocols used in generating the two data sets.

In a particularly preferred embodiment, each analysis program is executed on a different one of said processors. This can vastly improve the speed of the analysis.

In another aspect, the present invention provides a method for classifying a biological test sample, including the steps of:

choosing a set of marker molecules;

providing a database populated with reference expression data, the reference expression data including expression levels of a plurality of molecules in a plurality of reference samples, the plurality of molecules including at least the marker molecules, each reference sample having a pre-assigned value for each of one or more clinically significant variables;

accepting input expression data, the input expression data including a test vector of expression levels of the marker molecules in the biological test sample; and

assigning one of said pre-assigned values to the test sample for at least one of said clinically significant variables by passing the test vector to a statistical classification program;

wherein the statistical classification program has been trained to distinguish among said pre-assigned values on the basis of that part of the reference data corresponding to expression levels of the marker molecules.

The database may be in communication with a server computer which is interconnected to at least one client computer by a data network, said server computer being configured to accept the input expression data from the client computer.

Hosting the database on a server and allowing remote upload can improve the speed and efficiency of diagnosis. The clinician, having conducted a biopsy and assayed the sample (either themselves, or via a service laboratory located on site or nearby) to obtain a data file containing the expression levels of the marker molecules, can then simply upload the data file to the server for analysis and receive the test results within a short space of time, possibly within seconds. The server may reside on an internal network to which the clinician has access, or may be located on a wide area network, for example in the form of a Web server. The latter is particularly advantageous as it allows hosting and maintenance of a server accessing a large database of samples in one location, while a clinician located anywhere in the world and having access to relatively modest local resources can upload a data file to obtain a diagnosis based on a comprehensive set of annotated samples, such an analysis otherwise being inaccessible to the clinician.

The or each clinically significant variable may be selected from the group including disease state, disease prognosis, and treatment response. For example, the disease may be cancer, and the clinically significant variables may be organised according to a hierarchy, the levels of which may be selected from the group consisting of anatomical system, tissue type and tumour subtype. In that case, the classification program may include a multi-level classifier which classifies the test sample according to anatomical system, then tissue type, then tumour subtype. This provides a multi-marker, multi-level classification which is analogous to, but independent of, traditional approaches to diagnosis of tumour origin.

The marker molecules may include any combination of 100 or more of the polynucleotides listed in Table 4. We have found that sets of 100 or more of these molecules can provide a classification accuracy of greater than 85% for anatomical system and greater than 75% for tissue type.

In another embodiment, the disease is breast cancer, in which case the clinically significant variable may be risk of recurrence of the disease. The marker molecules in this embodiment may include the polynucleotides listed in Table 5. This is a prognostic, rather than diagnostic, application of the invention.

The invention is further applicable to other contexts in which predictive analysis is desired. For example, if a reference data set including expression levels for cancer patients having undergone one or more of various drug treatments is available, and the patients are annotated according to response to treatment, it would be possible to build and train a classifier to predict response of a patient who had not yet undergone the treatment, based on the expression levels of marker molecules in that patient.

In a particularly preferred embodiment, the reference expression data may be generated using a platform selected from the group including cDNA microarrays, oligonucleotide microarrays, protein microarrays, microRNA (miRNA) arrays, and high-throughput quantitative polymerase chain reaction (qPCR).

Oligonucleotide microarrays are particularly preferred for use in the present invention. If this type of microarray is used, each molecule being assayed is a polynucleotide, which may either be represented by a single probe on the microarray or by multiple probes, each probe having a different nucleotide sequence corresponding to part of the polynucleotide. If multiple probes are present, one of said analysis programs might include instructions for summarising the expression levels of the multiple probes into a single expression level for the polynucleotide.

Oligonucleotide microarrays such as those manufactured by Affymetrix, Inc and marketed under the trademark GeneChip currently represent the vast majority of microarrays in use for gene (and other nucleotide) expression studies. As such, they represent a standardised platform which particularly lends itself to collation of large databases of expression data, for example from cancer patients, in order to provide a basis for diagnostic or prognostic applications such as those provided by the present invention.

Preferably, the input expression data are generated using the same platform as the reference expression data. If the input expression data are generated using a different platform, then the identifiers of the molecules in the input data are matched to the identifiers of the molecules in the reference data prior to performing classification, for example on the basis of sequence similarity, or by any other suitable means such as on the basis of GenBank accession number, Refseq or Unigene ID.

Preferably, the statistical classification program includes an algorithm selected from the group including k-nearest neighbours (kNN), linear discriminant analysis, principal components analysis, nearest centroid classification and support vector machines.

In a further aspect of the present invention, there is provided a method of classifying a biological test sample from a cancer patient, including the step of:

comparing expression levels of a set of marker molecules in the test sample to expression levels of said set of marker molecules in a set of reference samples, each member of the set of reference samples having a known clinical annotation, to assign a clinical annotation to the test sample,

wherein the clinical annotation is selected from the group including anatomical system, tissue of origin, tumour subtype and risk of breast cancer recurrence.

In a yet further aspect, the present invention provides a system for classifying a biological test sample from a cancer patient, including:

a database populated with reference data, the reference data including expression levels of a set of marker molecules in a set of reference samples, each member of the set of reference samples having a known clinical annotation;

at least one processor; and

at least one storage medium containing program instructions for execution by said processor, said program instructions causing said processor to execute steps including:

accepting input data in the form of expression levels of the set of marker molecules in the test sample; and

assigning a clinical annotation to the test sample on the basis of the similarity of the input data to the expression levels of the set of marker molecules in the reference data;

wherein the clinical annotation is selected from the group including anatomical system, tissue of origin, tumour subtype and risk of breast cancer recurrence.

The marker molecules may include any combination of 100 or more of the polynucleotides listed in Table 4, or may include the polynucleotides listed in Table 5.

In the following discussion, embodiments of the invention will be described mostly by reference to examples employing Affymetrix GeneChips. However, it will be understood by the skilled person that the methods and systems described herein may be readily adapted for use with other types of oligonucleotide microarray, or other measurement platforms.

The terms “gene”, “probe set” and “molecule” are used interchangeably for the purposed of the preferred embodiments described herein, but are not to be taken as limiting on the scope of the invention.

Referring to FIGS. 1 and 2, there is shown in schematic form a system 100 and method 200 for classifying a biological test sample. The sample is acquired 220 by a clinician and then treated 230 to extract, fluorescently label and hybridise RNA to microarray 115 according to standard protocols prescribed by the manufacturer of the microarray. Following hybridisation, the surface of the microarray is scanned at high resolution to detect fluorescence from regions of the surface corresponding to different RNA species. In the case of Affymetrix arrays, each scanned “feature” region contains hundreds of thousands of identical oligonucleotides (25 mers), which hybridise to any complementary fluorescently labelled molecules present in the test sample. The fluorescence intensity detected from each feature region is thus correlated with the abundance (expression level) of the complementary sequence in the test sample.

The scanning step results in the production of a raw data file (a CEL file), which contains the intensity values (and other information) for each probe (feature region) on the array. Each probe is one of the 25 mers described above and forms part of one of a multiplicity of “probe sets”. Each probe set contains multiple probes, usually 11 or more for a gene expression microarray. A probe set usually represents a gene or part of a gene. Occasionally, a gene will be represented by more than one probe set.

Once the CEL file is obtained, the user may upload it (step 120 or 240) to server 110.

In the preferred embodiments, the system is implemented using a network including at least one server computer 110, for example a Web server, and at least one client computer. Software running on the Web server can be used to accept the input data file (CEL file) containing the multiple molecule abundance measurements (probe signals) for a particular patient from the client computer over a network connection. This information is stored in the system user's dedicated directory on a file server, with upload filenames, date/time and other details stored in a relational database 112 to allow for later retrieval.

The Web server 110 subsequently allows the user to select individual CEL files for analysis by a list of available diagnostic and prognostic methods, the list being able to be configured to add new methods as they are implemented. Results from the specific analysis requested, in the format of text, numbers and images, are also stored in the relational database 112 and delivered to the user via the Web server 110. All data generated by a particular user is linked to a unique identifier and can be retrieved by the user by logging into to the Web server 110 using a username and password combination.

When an analysis is requested by the user, at step 122, the raw data from the CEL file are passed to a processor, which executes a program 130a contained on a storage medium, which is in communication with the processor.

Accepting Clinical Data Input

In conjunction with the file that contains the multiple molecule abundance measurements (probe signals) for a particular patient, the user can also be asked to input other information about the patient. This information can be used for predictive, prognostic, diagnostic or other data analytical purposes, independently or in association with the molecular data. These variables can include patient age, gender, tumor grade, estrogen receptor status, Her-2 status, or other clinicopathological assessments. An electronic form can be used to collect this information, which the user can submit to a secure relational database.

Algorithms that combine ‘traditional’ clinical variables or patient demographic data and molecular data can result in more statistically significant results than algorithms that use only one or the other. The ability to collect and analyse all three types of data is a particularly advantageous aspect of at least some embodiments of the invention.

Program 130a is a low-level analysis module, which carries out steps of background correction, normalisation and probe set summarisation (grouped as step 250 in FIG. 2).

Background adjustment is desirable because the probe signals (fluorescence intensities) include signal from non-biological sources, such as optical and electronic noise, and non-specific binding to sequences which are not exactly complementary to the sequence of the probe. A number of background adjustment methods are known in the art. For example, Affymetrix arrays contain so-called ‘MM’ (mismatch) probes which are located adjacent to ‘PM’ (perfect match) probes on the array. The sequence of the MM probe is identical to that of the PM probe, except for the 13^thbase in its sequence, and accordingly the MM probes are designed to measure non-specific binding. A number of known methods use functions of PM-MM or log₂(PM)-log₂(MM) to derive a background-adjusted probe signal, for example the Ideal Mismatch (IM) method used by the Affymetrix MAS 5.0 software (Affymetrix, “Statistical Algorithms Description Document” (2002), Santa Clara, Calif., incorporated herein in its entirety by reference). Other methods ignore MM, for example the model-based adjustment of Irizarry et al (Biostatistics 4, p 249 (2003)), or use sequence-based models of non-specific binding to calculate an adjusted probe signal (Wu et al, JASA 99, p 909 (2004)).

Normalisation is generally required in order to remove systematic biases across arrays due to non-biological variation. Methods known in the art include scaling normalisation, in which the mean or median log probe signal is calculated for a set of arrays, and the probe signals on each array adjusted so that they all have the same mean or median; housekeeping gene normalisation, in which the probe or probe set signals for a standard set of genes (known to vary little in the biological system of interest) in the test sample are compared to the probe signals of that same set of genes in the reference samples, and adjusted accordingly; and quantile normalisation, in which the probe signals are adjusted so that they have the same empirical distribution in the test sample as in the reference samples (Bolstad et al, Bioinformatics 19, p 185 (2003)).

If the arrays contain multiple probes per probe set, then these can be summarised by program 130a in any one of a number of ways to obtain a probe set expression level, for example by calculating the Tukey biweight of the log (PM-IM) values for the probes in each probe set (Affymetrix, “Statistical Algorithms Description Document” (2002)).

Once the low-level analysis is completed, the background-corrected, normalised and, if necessary, summarised, data are passed (step 124) to program 130b, which is a quality control (QC) module. The execution of program 130b is depicted as step 260 in FIG. 2.

Quality data from an individual array can be used to infer the reliability and reproducibility of the entire molecular/genomic profile. One way to do this is to establish ranges for each quality metric that correspond to acceptable, warning and unacceptable levels. By analysing a large number of genomic profiles from reference samples comprising disparate tissue types and laboratory locations, a large body of quality data can be accumulated and stored in database 112.

The data for each of the quality metrics used herein approximately follow a log-normal distribution, as illustrated schematically in FIG. 4. Acceptable, warning and unacceptable ranges for each metric are thus calculated by determining the 25^thpercentile (Q1, indicated by 410), 75^thpercentile (Q3, indicated by 430) and corresponding interquartile range (IQR, indicated by 420) of the log-transformed values. Acceptable values are defined as those which lie between Q1−1.5*IQR and Q3+1.5*IQR.

Values in the ranges 405, 435 corresponding to ranges (Q1−1.5*IQR) to (Q1−3.0*IQR) or (Q3+1.5*IQR) to (Q3+3.0*IQR) are referred to as outliers, and are given a warning label. Values which lie to the left 403 or right 437, respectively, of those ranges are referred to as “far outliers” and are deemed to be unacceptably outside of the range of values used to develop and validate the gene expression test for which the test sample is being submitted.

The median, Q1/Q3 and IQR rather than mean and standard deviation are used to determine thresholds as the former are robust to outliers. This prevents the ranges from being overly influenced by a small number of samples that may not be representative of the true general distribution.

Table 1 is an example output from program 130b which describes each quality measurement (QC1 to QC8) and shows the value determined from the specific array being investigated. It also identifies the acceptable range and a variable classifier (Ok/Warning/Reject) column, which may change colour based on the contents of each cell. This allows the end user to rapidly determine if their input data is suitable for further analysis.

If a test sample passes the QC checks of program 130b it can then proceed (step 270) to predictive analysis as carried out by statistical classification program 135, which is used to assign a value of a clinically relevant variable to the sample. Such clinical parameters could include:

• • The primary tissue of origin for a biopsy of metastatic cancer;
  • The molecular similarity to patients who do or do not experience disease relapse with a defined time period after their initial treatment;
  • The molecular similarity to patients who respond poor or well to a particular type of therapeutic agent;
  • The status of clinicopathological markers used in disease diagnosis and patient management, including ER, PR, Her2, angiogenesis markers (VEGF, Notch), Ki67 etc.;
  • Possible chromosomal aberrations, including deletions and amplifications of part or whole of a chromosome;
  • The molecular similarity to patients who respond poor or well to a particular type of radiotherapy;
  • Other methods that may be developed by 3^rdparty developers and implemented in the system via an Application Programming Interface (API).

The predictive algorithms used in at least some embodiments of the present invention function by comparing the data from the test sample, to the series of reference samples for which the variable of interest is confidently known, usually having been determined by other more traditional means. The series of known reference samples can be used as individual entities, or grouped in some way to reduce noise and simplify the classification process.

Algorithms such as the K-nearest neighbor (KNN) algorithm use each reference sample of known type as separate entities. The selected genes/molecules (probe sets) are used to project the known samples into multi-dimensional gene/molecule space as shown in FIG. 5, in which the first three principal components for each sample are plotted. The number of dimensions is equal to the number of genes. The test sample is then inserted into this space and the nearest K reference samples are determined, using one of a range of distance metrics, for example the Euclidean or Mahalanobis distance between the points in the multi-dimensional space. Evaluating the classes of the nearest K reference samples to the test sample and determining the weighted or non-weighted majority class present can then be used to infer the class of the test sample.

The variation of classes present in the K nearest neighbors can also be used as a confidence score. For example, if 4 out of 5 of the nearest neighbor samples to a given test sample were of the same class (eg Ovarian cancer) the predicted class of the test sample would be Ovarian cancer, with a confidence score of 4/5=80%.

Other methods of prediction rely on creating a template or summarized version of the data generated from the reference samples of known class. One way this can be done is by taking the average of each selected gene across clinically distinct groups of samples (for example, those individuals treated with a particular drug who experience a positive response compared to those with the same disease/treatment who experience a negative or no response). Once this template has been determined, the class of a test sample can be inferred by calculating a similarity score to one or both templates. The similarity score can be a correlation coefficient.

Classifiers such as the nearest centroid classifier (NCC), linear discriminant analysis (LDA) or support vector machines operate on this basis (SVM). LDA and SVM carry out weighting of the genes/molecules when creating the classification template, which can reduce the impact of outlier measurements and spread the classification workload evenly over all genes/molecules selected, rather than relying on a subset to contribute to a majority of the total index score calculated. This can be the case when using a simple correlation coefficient as a predictive index.

To make clinically useful predictions about a specimen of biological material that has been collected from an individual patient, a large database of reference data from patients with the same condition is desirable. The reference samples are preferably processed using similar, more preferably identical, laboratory processes and the reference data are ideally generated using the same type of measurement platform, for example, an oligonucleotide microarray, to avoid the need to match gene identifiers across different platforms.

The reference data can be generated from tissue specifically collected or obtained for the diagnostic test being created, or from publically available sources, such as the NCBI Gene Expression Omnibus (GEO: http://www.ncbi.nlm.nih.gov/geo/). Clinical details about each patient can be used to determine whether the finished database accurately reflects the targeted patient population, for example with regard to age/sex/ethnicity and other relevant parameters specific to the disease of interest.

Clinical annotations can be used for analysis of the same input data at different levels. For example, cancer can be classified using a hierarchy of annotations. These begin at the system level, and then progress to unique tissues and subtypes, which are defined on the basis of pathological or molecular characteristics. The NCI Thesaurus is a source of hierarchical cancer classification information (http://nciterms.nci.nih.gov/NCIBrowser/Dictionary.do).

All data generated or obtained can be stored in organized flat files or in relational database format, such as Microsoft Access or Microsoft SQL Server. In this format it can be readily accessed and processed by analytical algorithms trained to use all or part of the data to predict the status of a clinically relevant parameter for a given test sample.

Following execution of classification program 135, the clinical predictions are stored in relational database 112. An interface 111 from the server 110 to database 112 can be used to deliver online and offline results to the end user. Online results can be delivered in HTML or other dynamic file format, whereas portable document format (PDF) can be used for creating permanent files that can be downloaded from the interface 111 and stored indefinitely. Result information in the form of text, HTML or PDF can also be delivered to the user by electronic mail.

AJAX Web 2.0 technologies can be used to streamline the presentation of online results and general functionality of the Web site.

A single processor may be used to execute each of the programs 130a, 130b, 135 and any other analysis desired. However, it is advantageous to configure the system 100 such that each analysis module is managed by a separate processor. This allows parallel execution of different user requests to be performed simultaneously, with the results stored in a single centralized relational database 112 and structured file system.

In this embodiment, illustrated schematically in FIG. 3, each module is programmed to monitor 320 a specific network directory (“trigger directory”). When the system operator requests 305 an analysis, either by uploading a new data file or requesting an additional analysis on a previously uploaded data file, the Web server 110 creates a “trigger file” in the directory 325 being monitored by the processing application. This trigger file contains the operator's unique identifier and the unique name of the data file on which to carry out the analysis.

When the classification module 135 detects (step 330) one or more trigger files, the contents of the file are read and stored temporarily in memory. The processing application then performs its preconfigured analysis routine, using the data file corresponding to the information contained in the trigger file. The data file is retrieved from the user's data directory (residing on a storage medium in communication with the server or other network-accessible computer) and read into memory in order to perform the requested calculations and other functions. Once the analysis routine is complete, the trigger file is deleted and the module 135 returns to monitoring its trigger directory for the next trigger file.

Multiple versions of the same classification module 135 can run simultaneously on different processors, all configured to monitor the same trigger directory and write or save their output to the same relational database 112 and file storage system. Alternatively, different modules in addition to classification module 135 could be run on different processors at the same time using the same input data. For processes that take several minutes (eg initial chip processing and Quality Module 130a) this enables analysis requests 305 that are submitted, while an existing request is underway, to be commenced before the completion of the first.

It will be appreciated that many other types of analysis (diagnostic, predictive, prognostic or other) may be conducted within the framework of the system provided by the present invention. When a new analysis program is created, it can be added to the list of analysis modules selectable by a user for execution on one or more input data files.

Additional modules can be added to the system by creating additional ‘trigger’ directories, monitored by analysis scripts. These can, of course, be used in conjunction with existing modules, such as the quality module described above.

A molecular profile can be adapted for use with the system by providing

• • Details of the technology used to measure the status of the molecules necessary to perform the test (e.g. genes, proteins, antibodies);
  • A list of molecular identifiers (probe or probe set identifiers, or gene or protein databank accession numbers, for example) specific to the platform used to develop, and to be used for future application of, the test;
  • A reference set of data from patients with the target disease (or other clinical identification) derived from the same class of patients; and
  • A statistical equation which describes how data corresponding to the molecular identifiers and reference dataset are used to predict the status of a test sample.

A custom results interface can then be created and incorporated into the system, linked to the underlying databases and results delivery mechanisms described previously. Technology-specific quality control measurements can also be incorporated, if they are not adequately represented by those contained in the quality module described previously.

The expO data, NCBI GEO accession number GSE2109, generated by the International Genomics Consortium, was used as a reference data set to train a tumor origin classifier.

Downloaded CEL files corresponding to the reference samples were pre-processed with the algorithms from Affymetrix MAS 5.0 software and compiled into BRB ArrayTools format, with housekeeping gene normalization applied. Using the associated clinical information from GSE2109, samples were classified at 3 levels of clinical annotation; (1) anatomical system (n=13), (2) tissue (n=29) and (3) subtype (n=295), as shown in FIG. 6. For Level 1 and 2 annotations, a minimum class size of three was set. The mean class sizes for the three levels of sample annotation were: (1) 149, (2) 66 and (3) 6, correlating with number of neighbors used in the kNN algorithm (r²=0.99).

Predictive gene expression models were developed using BRB ArrayTools and translated to automated scripts in the R statistical language, incorporating functions from the Bioconductor project (Gentleman et al, Genome Biology 5, R80 (2004)). The Web service was constructed in the Microsoft ASP.net language (Redmond, USA; version 3.5) with supporting relational databases developed in Microsoft SQL Server 2008. Statistical analysis of internal cross validation and independent validation series results was performed using Minitab (Minitab Inc. State College Pa., version 15.1.3) and MedCalc (MedCalc Software, Mariakerke, Belgium).

Most cells in the human body express under most circumstances, at comparatively constant levels, a set of genes referred to as “housekeeping genes” for their role in maintaining structural integrity and core cellular processes such as energy metabolism. The Affymetrix U133 Plus 2.0 GeneChip (NCBI GEO accession number GPL 570) contains 100 probe sets that correspond to known housekeeping genes, which can be used for data normalization and quality control purposes. For normalization purposes, the 100 housekeeping genes present on a given array within the reference data set were compared to those of a specific normalization array. To select a normalization array for this test, BRB-ArrayTools was used to identify the “median” array from the entire reference data set. The algorithm used was as follows:

• • Let N be the number of arrays, and let i be an index of arrays running from 1 to N.
  • For each array i, compute the median log-intensity of the array (denoted M_i)
  • Select a median M from the [M₁, . . . , M_N] values. If N is even, then the median M is the lower of the two middle values.
  • Choose as the median array the one for which the median log-intensity M_iequals the overall median M.

Housekeeping gene normalization was applied to each array in the reference data set. The differences between the log₂expression levels for housekeeping genes in the array and log₂expression levels for housekeeping genes in the normalization array were computed. The median of these differences was then subtracted from the log₂expression levels of all 54,000 probe sets, resulting in a normalized whole genome gene expression profile.

To select probe sets for the prediction of tumor origin, ‘one-v-all’ comparisons (t-tests) were performed for each tissue type in the training set (n=29) to identify probe sets which were differentially expressed in each tissue type compared to the rest of the data set. The probe sets identified by this procedure provide a characteristic gene expression signature for tumours originating in each tissue type.

In each comparison, genes that had a p-value less than 0.01 for differential expression, and a minimum fold change of 1.5 in either direction (upregulated or downregulated) were identified as marker probe sets. The analysis was performed using BRB ArrayTools (National Institute of Health, US). The 29 sets of marker probe sets were combined into a single list of 2221 unique probe sets, which are listed in Table 4.

The normalized expression data corresponding to these marker probe sets was retrieved from the complete 1942 reference sample×54000 probe set reference data, and this subset was passed to a kNN algorithm at both Level 1 (Anatomical-system, 5NN (nearest neighbours) used) and Level 2 (Tissue, 3NN used) clinical annotation.

To evaluate whether a smaller set of probe sets would achieve lower misclassification rates, leave-one-out cross validation (LOOCV) of the level 1 and 2 classifiers was performed using multiples of 100 probe sets from 10 to 2220, after ranking in descending order of variance. For each cross-validation test, the percentage agreement between the true and predicted classes was recorded and this is shown in FIGS. 7(a) and 7(b). The maximum classification accuracy obtained was 90% for Level 1 and 82% for Level 2. Reducing the number of marker probe sets used did not significantly improve computation speed.

CEL files from 22 independent Affymetrix datasets (all Affymetrix U133 Plus 2.0) containing a total of 1,710 reference samples were downloaded from NCBI GEO and processed as previously described. These datasets represent a broad range of primary and metastatic cancer types, contributing institutes and geographic locations, as detailed in Table 2.

Of 1,461 primary tumor validation samples that passed all QC checks, the Level 1 and Level 2 classifiers predicted 92% and 82% correctly. Tumor subtype data were not available for most validation datasets; therefore percentage accuracy of this level (3) of the classifier was not calculated. The difference observed between Level 1 and Level 2 classifier accuracy is largely influenced by ovary/endometriod and colon/gastric misclassifications. As with all comparisons of novel diagnostic methods with clinically derived results, the percentage agreement is dependent on multiple factors, including the accuracy of the clinical annotation, integrity of the sample annotations and data files as well as the performance characteristics of the method itself.

General linear model analysis was performed on the proportion of correct level 1 and level 2 predictions, including tissue type (n=10) and geographic location (n=3) in a regression equation to determine if these variables were factors in overall result accuracy. For Level 1 predictions (anatomical system), no significant difference in result accuracy was observed for tissue type (P=0.13) or geographic location (P=0.86). For Level 2 predictions (tissue type), a marginally significant difference was observed with tissue type (P=0.049) but no significant difference associated with location (P=0.38). The significant difference associated with tissue type at Level 2 is most likely associated with the small sample size of some tumor types.

The total set of 2,775 U133 Plus 2.0 arrays used in the training and validation stages of this analysis was used to derive acceptable ranges, as discussed above, for 8 different QC parameters. The ranges are shown in Table 3.

The Bioconductor package ‘SimpleAffy’ (Wilson and Miller, Bioinformatics 21, p 3683 (2005)) was used to generate measurements of background intensity, percentage of probe sets detected and 3′/5′ ratios. RNA degradation slopes were computed using the ‘AffyRNAdeg’ function in the ‘affy’ package (Gautier et al, Bioinformatics 20, p 307 (2004)).

The quality module also includes two assessments of data normalization. These are the MAS5 scaling factor and the log(base 2) of the housekeeping gene set normalization factor (i.e. the median difference between the log expression levels of housekeeping genes in a given sample and those of the housekeeping genes in the reference data as a whole).

The final assessment that is made is the signal-to-noise ratio (SNR), which is the log 10 ratio of mean probe set intensity divided by the mean background intensity. This metric is designed to ensure there is a sufficiently large difference between probe and background hybridization, which will not occur if the RNA is heavily degraded or problems with procedures such as RNA labelling or chip washing have occurred.

As an additional data quality control measure, a predictive Diagonal Linear Discriminant Analysis (DLDA) algorithm can be used for gender classification. Genes that were differentially expressed with a p-value less than 0.001 and minimum fold change of 2 between the 1,453 female and 695 male patients (regardless of cancer type) were selected as those which could distinguish males from females. A test sample, which is submitted for the purpose of other classification analyses, can be passed to the DLDA algorithm, which predicts the gender of the test sample based on the expression levels of the gender-discriminating genes thus identified.

The trained DLDA classifier for patient gender consists of 183 probe sets. During 3×3 fold cross validation, the gender of 97% of the 2,148 samples was correctly predicted, with a sensitivity of 97% and specificity of 95% from this internal validation exercise.

Reflecting the nature of existing diagnostic workflows for metastatic tumors, a novel 3-tiered approach to predicting the origin of a metastatic tumor biopsy was developed. For each test sample analysed, 3 rounds of kNN classification were performed, using the 3 levels of annotation previously described, i.e. (1) anatomical system, (2) tissue and (3) histological subtype, with k=5, 3 and 1 respectively. The decreasing value of k with increasing specificity of tissue annotation was chosen based on the decreasing mean class size at each tier of the classifier, with which it is highly correlated (r^2=0.99).

A measurement of classifier confidence was generated for Level 1 (k=5) and Level 2 (k=3) results by determining the relative proportion of a test sample's 5 or 3 neighbors (respectively) that contribute to the winning class. The Level 3 prediction (k=1) identifies the specific individual tumor from the reference database that is closest to the test sample, in multi-dimensional gene expression space. As such, it is not possible to calculate a weighted confidence score for this level of classifier.

To determine the internal cross validation performance of the reference data and 3-tier algorithm, leave-one-out cross validation (LOOCV) was performed on the reference data set, using annotation levels 1 and 2. Results were tallied and overall percentage agreement and class-specific sensitivities and specificities were determined. The R/Bioconductor package “class” was used for kNN classification and predictive analyses.

Two training data sets from untreated breast cancer patients (GEO accession numbers GSE4922 and GSE6352), including a total of 425 samples hybridized to Affymetrix HG-U133A arrays (GEO accession number GPL96) were downloaded in CEL file format. Clinical data were available for age, grade, ER status, tumour size, lymph node involvement, and follow-up data for up to 15 years after diagnosis were also available. An independent validation data set, consisting of samples from 128 Tamoxifen-treated patients hybridized to Affymetrix HG-U133Plus2 arrays with age, grade, ER status, nodal involvement and tumour size data, was also obtained.

A semi-supervised method substantially in line with the method described by Bair and Tibshirani (PLoS Biology 2, p 511 (2004), incorporated herein in its entirety by reference) was used, with algorithm settings of k=2 (number of principal components for the “supergenes”), p-value threshold of 0.001 for significance of a probe set being univariately correlated with survival, 10-fold cross-validation, and age, grade, nodes, tumour size and ER status used as clinical covariates. The method identified 200 prognostic marker probe sets, shown in Table 5, and gave the following model for risk of recurrence (Formula 1):

PI=∑i=1200wixi=0.139601(grade)+0.64644(ER)+0.938702(nodes)+0.010679(size(mm))+0.023595(age)+0.243639

In Formula 1, w_iis the weight of the i^thprobe set, x_iis its log expression level, and PI is prognostic index.

FIGS. 8(a) and 8(b) show Kaplan Meier analysis of 10-fold cross validation predictions made for the 425-sample training set. Log rank tests were used to compare the survival characteristics of the two risk groups identified.

Evaluation of the cross-validation predictions made for the training set revealed a highly statistically significant difference in the survival characteristics of the high and low risk groups. Of the 425 patients, 297 (70%) were classified as high-risk and 128 (30%) as high risk. The p-value of the Kaplan Meier analysis log-rank test was P<0.0001 and the hazard ratio of the classifier was 3.75 (95% confidence interval 2.47 to 5.71).

In the training set, 85% of patients classified as low risk were disease-recurrence free at 5 years after treatment. In the high-risk group, 41% of patients experienced disease recurrence within this same time period.

FIGS. 9(a) and 9(b) show survival characteristics of the high and low risk groups for the independent validation data set. The groups identified in this cohort are more similar to each other up to 3 years after diagnosis. This is likely attributable to the use of Tamoxifen in these patients. After this time point survival characteristics are significantly different.

Kaplan Meier analysis and log-rank testing was performed on the independent validation set. The P-value associated with the log rank test was P=0.0007. A hazard ratio of 4.90 (95% confidence interval 1.96 to 12.28) was observed. These figures indicate that the classifier was able to stratify the patients into two groups with markedly different survival characteristics.

Overall those individuals in the high-risk group are 4.9 times more likely to experience disease recurrence than those in the low risk group in the 10 years after diagnosis. Three quarters of the independent validation patients are classified as low risk (n=97) and of these, 90% are recurrence-free after 5 years.

Additionally, multivariate Cox Proportional Hazards analysis was performed on the 128 sample independent validation set. Two models were built and tested, one including the clinical variables only, and the other including the clinical variables and classifier prediction variable (high/low risk). The significance level of the clinical-only model was P=0.0291, whilst for the clinical+classifier model it was P=0.0126. The classifier remained independently prognostic in the second model (P=0.048).

These results indicate that the classifier (comprised of 200 genes+5 clinical variables) is able to stratify patients into high and low risk groups for disease recurrence. Furthermore, the stratification of patients is more statistically significant than the use of clinical variables alone. The prognostic significance of the classifier has been evaluated in patients who do and do not receive Tamoxifen treatment following their initial diagnosis and surgical procedure.