WO2006087240A1

WO2006087240A1 - Method, system and software for carrying out biological interpretations of microarray experiments

Info

Publication number: WO2006087240A1
Application number: PCT/EP2006/002052
Authority: WO
Inventors: Audrey Kauffmann; Philippe Dessen; Vladimir Lazar
Original assignee: Institut Gustave Roussy (IGR)
Current assignee: Institut Gustave Roussy (IGR)
Priority date: 2005-02-15
Filing date: 2006-02-14
Publication date: 2006-08-24
Anticipated expiration: 2007-08-15
Also published as: EP1691311A1

Abstract

The invention relates to DNA microarray data analysis. A p-value of each gene of a cluster, such as a pathway given by a database, is calculated on the basis of fluorescence signals issued from a microarray chip, corresponding to a probability that the gene is reported as being differentially expressed between several cell types, under a certain condition. Clusters from pre-existing databases containing a given number of genes of the microarray are formed, and, for each cluster, significant genes are selected with respect to their p-values. According to the invention, genes of the microarray are randomly drawn in order to evaluate a Z-score indicating a significance of each cluster.

Description

Method, system and software for carrying out biological interpretations of microarrav experiments

This invention relates to searching biological interpretations from microarray experiments. Specifically the present invention relates to statistical methods for functional analysis and identification of differentially expressed genes of microarray data.

The present invention also relates to statistical methods for determining the molecular functions of genes and/or the metabolic pathways and/or the interaction networks and/or biological processes involving the gene products and/or the sub-cellular localization of said gene products and/or the chromosomal localization of said genes.

Biological systems are controlled by hundreds of genes working in concert. Microarray technology allows to monitor simultaneously thousands genes and identify a subset of differentially expressed gene in a particular cell type under a given biological or pathological condition. Changes of gene expression levels illustrate modification in gene functions.

The basic concept of microarray analysis is simple. Nucleic acids from a cell type are harvest and labeled to generate targets. Said nucleic acids are hybridized to probe DNA sequences (cDNA or oligonucleotides) corresponding to specific genes that have been affixed onto a solid matrix in a known configuration. Hybridization, based on base-pairing between probes and targets, provides a quantitative measurement of the abundance of a particular gene of the target nucleic acids. This information is captured digitally and subject to analysis. Comparison of hybridization pattern enables the identification of nucleic acids that differs in abundance in two or more target samples. Simulations are powerful tools for understanding dynamic complexity and envisioning a variety of outcomes through experimentation.

In particular, when microarray experiments are to be analyzed, the genes are filtered on a criteria based, for example, on p-value, fold change, log(ratio) and/or intensity. Then, various statistical tests can be performed to refine the selection so that the found genes can fit correctly. Finally, common properties between these genes defined as having an interest can be found. This last step often consists in clustering genes by molecular functions or biological processes according to the Gene Ontology consortium (Ashburner et al., Nature Genetics, 2000; 25: 25-29), for instance, FatiGO (Al-Shahrour et al., Bioinformatics, 2004; 20: 578-580), GOminer (Zeeberg et al., Genome Biology, 2003; 4:R28), MAPPFinder (Doniger et al., Genome Biology, 2003; 4:R7) or by pathways from KEGG (Ogata et al., Nucleic Acids Research, 1999; 27: 29-34.) with DAVID (Dennis et al., Genome Biology, 2003; 4:P3) or from BioCarta (http://www.biocarta.com/genes/index.asp).

While analyzing the approaches of the prior art, it has been realized that an "a priori" that can either be right or wrong was imposed during the first two steps. For instance, genes with a low fold change were simply not taken into account, thus ruling out potentially important regulation pathways.

Statistical analysis of microarray data is probably the most difficult problem associated with the use of this technology. The aim is to apply statistical approaches to determine gene expression change significance under a given biological or pathological condition, thus enabling the extraction of significant biological information.

It is therefore an object of the present invention to provide statistical methods for identifying significant changes of gene expression, to cluster genes, and to determine whether the differences among individual genes and/or clusters of genes are significant in identifying differentially expressed genes in different cell types, under different biological or pathological conditions.

The present invention aims a method as claimed in claim 1. Embodiments of the method are aimed in dependent claims. The invention aims also a system for performing the method. The invention aims also a computer program product for implementing all or some of the steps of the method.

Other advantages and characteristics of the invention appear on reading the following detailed description of an implementation given below by way of example with reference to the accompanying figures, in which:

- Figure 1 A represents a microarray fluorescence readout;

- Figure 1 B shows an outline of the algorithm of the present invention;

- Figure 2a shows an expression data set of a microarray experiment wherein for each sample of interest, each gene is associated with a measurement of fluorescence intensity generated from the microarray fluorescence readout as shown in Figure 1A;

- Figure 2b shows a table wherein each gene is assigned to a p-value;

- Figure 2c shows a table wherein genes are clustered into pathways according to information retrieved from a pre-existing database; - Figure 2d shows a table wherein each pathway of Figure 2c is assigned to p-values (listed in the table of figure 2b) associated to its constituting genes;

- Figure 2e shows a table wherein each pathway is assigned to a percentage of significant genes; - Figure 2e' shows a similar table than that of Figure 2e wherein each pathway is assigned to a percentage of significant genes when a Fisher's test is carried out;

- Figure 2f shows a table wherein each pathway is assigned to randomized p-values; - Figure 2g shows a table wherein each pathway is assigned to a percentage of significant genes found randomly;

- Figure 2h shows a table wherein each pathway is characterized by a Z- score value; - Figure 2h' shows a similar table than that of Figure 2h wherein each pathway is characterized by a p-value determined according to a Fisher's test;

- Figure 3 shows an example of representative results table consisting of a rank list of clusters generated by the implementation of the present invention;

- Figure 4a shows an expression data set of a microarray experiment as shown in Figure 2a;

- Figure 4b shows a table wherein each pair of genes is assigned to a correlation estimator; - Figure 4c shows a table wherein known and unknown genes are grouped with respect to the correlation estimators;

- Figure 4d shows a table wherein known genes are clustered into pathways which are characterized by Z-score values;

- Figure 4e shows a table wherein unknown genes are assigned to clusters of known genes;

- Figure 5 shows a system for performing the invention;

- Figure 6 shows a representative flowchart of an embodiment of a main algorithm of the software according to the invention;

- Figure 7 shows a flowchart algorithm for calculating a Z-score value (z- scored assigned to a cluster PJ; and

- Figure 8 shows a flowchart algorithm for implementing the flowchart of Figure 7 to clusters Pj of a pre-existing database and for ranking the clusters according to their Z-score value.

Pattern of gene expression in a sample of interest, for example a cell, is characteristic of its current conditions. The conditions may be biological conditions, pathological conditions, or diagnostical or prognostical conditions. The types of cells may also be cultured cells and primary cells isolated from an organism. The effect of a compound or composition on a type of cells may also be assessed using the method of the present invention.

Signals issued from a microarray chip may correspond to measurements of fluorescent intensities. Accordingly, for each gene, there will be as many individual signals as cell types under a given condition.

Figure 1A represents a microarray fluorescence readout consisting of columns of genes (HUGO or Unigene symbols as column name) and rows of samples of interest under different conditions. Complementary DNA (cDNA) were prepared from the samples of interest under condition 0 (the last 8 samples) and condition 1 (the first 4 samples) and were labeled with red and green fluorescence dyes, respectively. Each square corresponds to the fluorescence intensity of a gene in a specific sample. Due to a lack of coloring, the light grey and dark grey squares reflect samples that have been labeled with green and red fluorescent dyes, respectively. This fluorescence is based on the hybridization of the labeled cDNA from said specific sample to the control cDNA that has been affixed onto the microarray chip. Dual-channel laser excites the corresponding dye, which fluoresces in proportion to the degree of hybridization that has occurred. Fluorescence intensity is proportional to the level of gene expression in each sample.

Microarray data may be interpreted within the context of gene function and functional relationship between genes. Usually, microarray data are related to existing biological databases.

The database that can be used in the context of the invention is selected from the group consisting of, for example, the Gene Ontology consortium, the

KEGG database or the Biocarta database. Each database gathers genes/gene products that are functionally related into clusters. A gene product might have one or more molecular functions, and may be used in one or more biological process; it might be associated with one or more cellular components. Accordingly, the Gene Ontology consortium that describes how genes products behave in cellular context gathers genes products into clusters according to three principles, i.e. molecular functions, biological processes and cellular components. It is also possible to gather genes from the databases according to their chromosomal localization.

Therefore, it is proposed to seek biological interpretation from microarray experiments directly by searching functions of interest in a whole data set. Instead of being the final step of the analysis, the step of searching functions is considered in the present invention as a starting point to focus on few groups of genes in the rest of the study.

In the context of the present invention, searching biological interpretation from microarray experiments is a freely available tool which allows the identification of pathways and ontologies that contain a relevant number of genes that are differentially expressed between two tested conditions by DNA microarray. The present invention can be applied, for example, to human or mouse studies and makes use of the KEGG pathways, the BioCarta networks and/or the Gene Ontology clusters. The present invention relies on the assessment of p-values associated to t-tests and the calculation of Z-scores after randomization of the genes.

A significance of the results obtained when implementing the invention can be fixed by user. A significance threshold can be chosen preferably between 95% and 99.9 %. Accordingly, the genes that are characterized by a p-value between 5.10^"2 and 10~³ will be deemed as being significantly differentially expressed. Similarly, a significant Z-score value is superior to 2.5 and more preferably superior to 3, thus indicating that the cluster of selected genes associated to the Z-score contains genes that are significantly differentially expressed.

In figure 1B, pre-existing databases that may be used are the KEGG Pathways, BioCarta Pathways or Gene Ontology. Each database may be organized into groups of gene clusters according to their common involvement in interactive pathways, biological process, molecular function or cellular component. For example, the BioCarta pathways database consists of individual genes, such as G1 , G3, G5...Gx or G6, G7...Gy, that are classified into gene clusters according to their common involvement in a pathway, such as pathways x and y, respectively.

The data Set may consist of the log intensities or log ratios that are based on the microarray fluorescence intensity readout shown in Figure 1A. These log intensities or ratios are determined for each gene in order to generate its p- value using a student t-test, or more generally an ANOVA analysis. The p- value corresponds to the probability that a gene, is reported as differentially expressed between conditions 0 and 1.

With each pathway of the database, the program of the present invention searches the data set for the genes of the pathway and stores their associated p-value within said pathway. The pathways for which no gene was identified in the data set are suppressed from the rest of the study. Then, for each pathway, the program counts the number of genes found in the data set and the number of genes that have a p-value lower than a predetermined threshold and expressed it as a percentage value (Ps). To test the significance of the percentage of differentially expressed genes by pathways, the program randomly selects X genes from the data set and calculates the percentage of genes with a p-value lower than the threshold compared with X (Pr). This operation is repeated N times (the number N of iterations is preferably determined by the user) and a Z-score is computed for each pathway as followed: _v

_Z - ^Ps -Pr

, under the null hypothesis HO: Z ~N (0,1 ) P wherein P_s is the percentage of significant genes found in the data set for a given annotation; P_r is the percentage of significant genes found randomly, Pr the mean of the N P_r and σ_Pr the square of the variance of the N P_r.

Figures 2a to 2h refer to a first embodiment of the method according to the invention. The method generates a results table such as the one shown in figure 3. On one hand, an expression data set of microarray experiments is collected (figure 2a). Said data set preferably consists of fluorescence intensities that are obtained from the microarray fluorescence read out (figure 1A) and may be expressed either in log(ratio), ratio, log(intensity), fold-change or any numerical expression related to the fluorescence. Missing values or unusable experimental values, in the submitted data set, can be represented by a blank value. It is possible to choose how many "blank values" per row can be authorized and keep the corresponding gene within the analysis. For a given gene, if the number of "blank values" is over a threshold, this given gene will not be taken into account; if it is lower, a t-test or ANOVA will be performed on the non missing values. For each gene that is kept in the analysis, the method according to the present invention performs a statistical test (either ANOVA or t-test in the case of two conditions) that compares the mean of expression level between the studied conditions. ANOVA analysis tests whether the null hypothesis that the difference between the means of expression level in each condition is zero. If the difference between the means is among those that would occur rarely by chance when the null hypothesis is true, the null hypothesis is rejected and the results is considered as statistically significant. Using said test, each gene is assigned to a p-value in order to ascertain whether said gene is significantly differentially expressed between the tested conditions (figure 2b).

On the other hand, files containing information from database such as KEGG, Biocarta and Gene Ontology are collected (figure 2c). These files are automatically updated. For each cluster of genes (of an ontology term or pathway), the program according to the present invention searches the genes that are present in both said cluster and the submitted data set (figure 2d). Then, the program counts for each cluster the number of these genes that have an associated p-value lower than a predetermined threshold (fixed-by- user), thus identifying the genes of the cluster that are differentially expressed between the tested conditions. When the program is applied to each cluster of a given database, the clusters that have a significant number of differentially expressed genes are identified (figure 2e). Then, for X genes in a given cluster X, a step of the method according to the invention consists of picking out randomly X genes in the data set (as shown in the example of figure 2b) and, counting how many of these genes are differentially expressed (figure 2f). This step is repeated N times, preferably a thousand times (N=1000). Then, clusters that have a significant number of genes differentially expressed at random are known (figure 2g). Finally, for each given cluster, a Z-score (figure 2h) is calculated by comparing the observed values (figure 2e) to the randomized values (figure 2g).

Figure 3 shows a representative results table consisting of a rank list of clusters from pre-existing database. In the example of the present results table, genes were clustered according to their cellular functions or activities (pathway names listed in the box) within the Gene Ontology database with only two genes of each listed (right box). Each pathway cluster has a calculated Z-score (Figure 1B) that is used to rank them in decreasing order of value. Only the pathway clusters with Z-scores that are greater than a predetermined threshold, having preferably a value of 3 in the given example, are deemed to be significant. The pathway cluster with the highest Z-score is the most significant, which indicates that it has the most differential expression between conditions 0 and 1.

Figures 4a to 4e refer to a second embodiment of the method according to the present invention, which allows the determination of a function of an unknown gene as described hereafter. From the data set (figure 4a), which is preferably the same than the data set of figure 2a, the software calculates correlation estimators between each possible pair of genes (figure 4b) using a correlation test. The correlation estimator may either be a distance or a correlation coefficient, or any estimator of the similarity in expression level between two genes. Thanks to these correlation estimators, groups of genes are clustered according to the similarity in their expression level (figure 4c). The method according to this embodiment may also combine different types of correlation estimators (figure 4b), such as distances and correlation coefficients to group genes according to different criteria (figure 4c). In any cases, associations are made between the results (figure 4d) and the groups obtained (as shown in figure 4c) in order to assign unknown genes (figure 4e) to clusters of known genes (the pathways shown in the example of figure 4e) that share similar level of expression (in figure 4c). Accordingly, the method according to this embodiment assigns putative functions, processes, interactions to unknown genes or pathways (as shown in figure 4e).

A system for performing the invention is shown in Figure 5. Said system includes a processing unit PU composed of:

- an interface INT for receiving signals issued from the microarray data set,

- a processor MP, - a memory for storing results and the signals measurements obtained from the microarray,

- a working memory RAM (random access memory).

More specifically, the memory MEM is arranged for storing a computer program product in order to perform the steps according to the method of the invention. The software may also be loaded from a CD-ROM player and/or from the internet. In addition, the system may also include a communication interface COM, such as a modem or an equivalent device in order to import information from pre-existing databases.

Figure 6 shows a representative flowchart of an embodiment of the algorithm of the software according to the invention. In step 60, the software activates a connection to a pre-existing database DB to download a number of ni genes that are involved in a cluster Pj. The cluster Pj may be a pathway, a function, a process or a cellular component. In step 61, a p-value is determined for each of the nj genes involved in the cluster Pj as described in Figure 1B. In step 62, the software assigns a variable n_j that is equal to the number of nj genes involved in the cluster Pj. In step 63, the value ofia floating variable j is set at j=1 and a recurrent treatment is implemented for the floating variable j. The test 64 compares the p-value (p-val)_j of gene G_j of the cluster P₁ with a predetermined threshold THR. For example, the value of the threshold may be equal to 10^'3. If the p-value (p-val)_j is greater than said threshold THR (arrow y), the floating variable nj as defined in step 62 is decremented by one unit. If the p-value (p-val)j is significant, i.e. less than the threshold THR, step 66 verifies whether the floating variable j reaches the value of the number of genes n-, of the cluster Pj. If test 66 is positive (arrow y), the percentage (Per)j of genes that have a significant p-value in the cluster Pj is calculated in step 67. Otherwise (if variable j has not reached the number ni : arrow n), the value of the floating variable j is incremented and step 64 is repeated. The flowchart algorithm shown in Figure 6 consists of a main routine PROC that is intended to be applied to both the genes involved in a cluster from a pre-existing database and the randomly drawn genes. Therefore the Z-score value (z-score)j assigned to each cluster Pj can be calculated by using a software having a flowchart algorithm as represented in Figure 7. Step 70 of Figure 7 consists of applying the main routine PROC to the cluster Pj from a pre-existing database with the exact genes given by the database REF[Pj]. In step 71 , a percentage Preή is obtained by running the main routine PROC as previously described with reference to Figure 6. On the other hand, a floating variable k is set in step 72. The floating variable k corresponds to an index of random draws amongst N=IOOO draws in the given example. In step 73, the number of nj genes that are involved in the cluster Pj is obtained from the database. In step 74, nj genes are randomly drawn. In step 75, the p-value of each randomly drawn gene is used to perform the main routine PROC on the nj randomly drawn genes (step 76). Then, the percentage of genes amongst the randomly drawn genes that have a significant p-value is determined as previously described with reference to Figure 6. Test 78 verifies whether the number of draws reaches a value of N=IOOO. If the floating variable k equals to N (arrow y), a Z-score value is calculated and assigned to each cluster Pj in step 79. Otherwise (arrow n), the value of the floating variable k is incremented. Each Z-score value is calculated according to the formula:

/ x Prefj -X \z - score)_! = , ^σPR wherein:

representing the mean of the N (Per)i_k values over the N draws (step 81 ), and σpR is related to the variance of the draws (step 82). Referring to Figure 8, the flowchart algorithm that is represented in Figure 7 is applied to all the clusters Pj obtained from the pre-existing database. In step 90, a Z-score value is calculated and assigned to each cluster Pj. In step 91 , the clusters Pj are ranked according to their respective Z-score value. The clusters that have a Z-score value greater than a predetermined threshold are selected as significant. In step 92, a results table such as the one represented in Figure 3 is finally obtained.

The results obtained by implementing the present invention have been validated by laboratory experiments, in particular on human and mouse samples. The tests show that the bigger the submitted microarray data sets are, the more relevant the results obtained by the present invention are.

Certainly, the present invention should not be limited to the above-described embodiment, but may also include others embodiments.

Thus, as an alternative embodiment of the above-described Z-score test, it is also possible to carry out a Fisher's test. In the Fisher's test, each percentage of significantly differentially expressed genes that are identified per category or pathway is compared with the percentage of significantly differentially expressed genes of the complete data set. The corresponding p-value is determined using the following formula

_ (a + b)\(c + d)\(a + c)l(b + d)\ ^{P ~} (a + b + c + d)\a\b\c\d\

wherein a is the number of significant genes in the data set, b is the number of non significant genes in the data set ; c is the number of significant genes in a given category or pathway and d is the number of non significant genes in said given category or pathway. In such embodiment, Figures 2e' and 2h' replace the above-described Figures 2e and 2h of the Z-score test's embodiment. In addition, the steps corresponding to Figures 2f and 2g of the Z-score test's embodiment can be deleted if the Fisher's test is carried out.

Claims

1. Method implemented by computer means for analyzing DNA microarray data given by a DNA microarray set that generates: - for a plurality of genes, on one hand,

- and for a plurality of cell types under different biological and/or pathological condition, on the other hand, signals corresponding to an expression level of a gene in a cell type under a predetermined condition, said process comprising: a) for each gene, on the basis of said signals, determining a statistical parameter that corresponds to a probability that said gene is reported as being differentially expressed between said plurality of cell types and under said predetermined condition, b) assigning to each gene of the microarray at least one cluster of a preexisting database, said cluster containing a given number of genes of the microarray, and c) for each cluster, selecting the genes based on a significance of their statistical parameter and determining a number of the selected genes, wherein the method further comprises the steps of: d) randomly draw N times a given number of genes as defined in step b) from the microarray and associate to each of said random genes their statistical parameter for determining, for each draw, a number of random genes that have a significant statistical parameter, and evaluate a mean value, over the N draws, of said numbers of random genes that have a significant statistical parameter, e) compare said mean value with said number of selected genes as determined in step c), a result of said comparison indicating a significance of said cluster.

2. The method according to claim 1 , wherein the number N is approximately one thousand.

3. The method according to claim 1 or claim 2, wherein said comparison of step e) comprises a Z-score calculation according to a formula:

_Cr , wherein: ^σp_r

- P_s is a percentage of significant genes selected in the microarray for a given cluster i;

- P^r is a mean value, over N draws, of N percentages P_r of randomly drawn significant genes, and

- σ_Pr is the square of a variance of said N percentages P_r.

4. The method according to any of the precedent claims, wherein the selection of step c) is based on a comparison of the statistical parameter of each gene and a threshold value.

5. The method according to any of the precedent claims, wherein said statistical parameter is a p-value in an analysis of variance such as ANOVA.

6. The method according to any of the precedent claims, wherein said statistical parameter is a p-value in a mean comparison analysis such as a Student t-test.

7. The method according to claim 4, taken in combination with any of claims 5 and 6, wherein a selection of a gene in step c) is based on whether its p-value is inferior to a threshold value.

8. The method according to claim 7, wherein said threshold value is preferably 1(T³.

9. The method according to any of the precedent claims, wherein: - said pre-existing database comprises a plurality of clusters, each cluster comprising a given number of genes of the microarray,

- steps d) and e) are performed for each cluster, the method further comprising a step of: f) ranking the clusters of said plurality of clusters with respect to their comparison results.

10. The method according to claim 9, wherein a selection of the clusters that have a comparison result which is higher than a predetermined threshold, as the clusters the most significant, is performed at step f).

11. The method according to claim 10, wherein said predetermined threshold is fixed by user. i

12. The method according to any of claims 10 and 11 , taken in combination with claim 3, wherein the ranking of step f) is performed on Z-score values associated to clusters of said plurality of clusters, respectively, and wherein said most significant clusters have a Z-score value, as a comparison result, which is higher than said predetermined threshold.

13. The method of claim 12, wherein said Z-score threshold value is preferably fixed at 3.

14. The method according to any of the precedent claims, wherein:

- in step a), a correlation test is further performed on each pair of genes of the microarray and a correlation estimator is calculated for each pair of genes, - groups of genes having a high correlation estimator between them are clustered, and

- the content of a cluster given by the pre-existing database in step b) is enhanced with at least one unknown gene which belongs to a same group than a known gene already indexed in said given cluster of said database.

15. The method according to claim 14, wherein said correlation estimator represents a distance between the expression levels of two genes of a pair.

16. The method according to claim 15, wherein said groups are clustered with respect to a similarity criteria between expression levels of two genes.

17. The method according to any of the precedent claims, wherein said signals are derived from florescence intensity measurements.

18. The method according to any of the precedent claims, wherein said clusters are related to metabolic pathways and/or molecular functions and/or interaction networks and/or biological processes and/or cellular components, and/or chromosomal localization.

19. A computer system for performing the method according to any of the precedent claims, comprising:

- an interface module for collecting the signals issued from a DNA microarray set, - a communication module for obtaining said clusters from said pre-existing database,

- a unit for processing said signals, and

- a memory medium for storing a computer program intended to be run by said processing unit for performing all or some of the steps of the method according to any of the precedent claims.

20. A computer program product intended to be stored in a memory of a processor unit or on a removable memory medium suitable for co-operating with a reader of said processor unit, or downloadable from a distant site, wherein the program product comprises instructions for implementing all or some of the steps of the method according to any of claims 1 to 18.