METHOD FOR ROBUST COMPARISON OF DATA

This invention pertains in general to the field of bioinformatics. More particularly the invention relates to a computer-implemented method for analyzing multivariate data from genetic studies.

Many modern experimental techniques make it possible to study a large number of biological variables in a single experiment. For example, with the microarray technique researchers can monitor the activity of tens of thousands of genes simultaneously. Other measurement techniques, such as next generation sequencing, provide the possibility to study many aspects of genomics and proteomics and provide data sets with even higher number of variables. Clearly, to make sense of such a vast data set, it is necessary to be able to rank the variables in order of importance for a particular application. Having such a ranking, the researcher can select a number of top-ranked genes that are strongly linked to the studied condition and that should therefore be studied in more depth. Methods for generating variable rankings from experimental results are abundant. For example, the genes can be ranked by their relative increase (or decrease) in activity among breast cancer patients compared to healthy persons, in order to get a deeper understanding of risk factors and causes of breast cancer. Obviously, to make it possible for the researcher to draw generalizable and biologically relevant conclusions the ranking of the genes should not depend significantly on which particular cancer patients and healthy controls that are included in the study. In practice, however, it has been shown in several studies that there is a strong dependence on the selected sample (Fortunel et al, 2003; Michiels et al, 2005; Fan et al, 2006; Abraham et al, 2010). For example, Abraham et al (2010) studied data from breast cancer patients classified as high risk patients or low risk patients based on the time to occurrence of distant metastasis (data set from Ivshina et al, 2006). They concluded that the ranked list of genes that were differentially expressed between high and low risk patients was extremely sensitive to the choice of patients. A single gene could easily be ranked anywhere between position 1 and 1,500 (among 22,215 genes) if different subsets of the patients were used as the basis for the gene ranking. Similarly, Soneson and Fontes (2012) noted, using a data set comparing lung cancer patients with good and poor outcome, respectively (Bhattacharjee et al, 2001) as well as a data set comparing breast cancer patients with different mutation status (Hedenfalk et al, 2001) that the gene list obtained with a variety of different ranking methods depended considerably on the selection of included patients. Clearly, this ambiguity makes correct determination of the top-ranked genes from a single sample almost impossible, and biological interpretation of the results extremely difficult. It has been suggested that part of the reason for the apparent instability of gene rankings is that there are many genes that have similar functions in the cell, and that only one gene from a specific pathway needs to be altered in any given patient. As a consequence, the set of top-ranked genes will depend on the constitution of the particular sample.

The goal of a biological study is almost never to just rank the genes in order of importance, but rather to compare the results to those obtained from other studies, or to databases containing current knowledge. Apparently, to obtain reliable results from such comparisons it is important to have a stable and reliable ranking of the genes from the experiment. Furthermore, to be able to make such comparisons, it is necessary to have a unified way of representing the results from different sources. For example, while the results of a microarray experiment are generally represented by the ranked list of genes, the contents of a database can be represented e.g. as small, unordered collections of genes, where the genes in each collection share similar functions in the cell. By unifying the representation, we can compare any two gene lists, whether ranked or unranked and independently of the number of included genes, within the same framework.

Comparing the results from a microarray experiments to those from other studies, or to databases, helps the researcher to interpret the results in a biological context. Of course, stability and reproducibility are as important as the gene ranking. By studying the biological meaning, i.e. the indication with regards to a certain biological state, of the results instead of only the actual ranking, the results are slightly more stable and the dependence on the selection of patients is somewhat reduced. However, using current methods the stabilization, the effect is still not sufficient to allow biologically relevant conclusions to be drawn.

There is thus a need for an improved method for analyzing multivariate data from genetic studies and particularly for a method enabling a determination of relationship between samples, which relationship is indicative of a biological state, and which enables improved analysis of the biological state.

An objective of the present invention is to provide relatively robust comparison of data, such as gene expression data, by providing a computer-implemented method and a computer-readable medium, that analyze data according to the appended patent claims.

A general solution according to embodiments of the invention is to use information extracted from the experimental data to obtain a more stable list, which is represented in a form allowing straightforward comparisons to other lists of predetermined data, both in ordered and unordered form.

According to a first aspect, a computer-implemented method for identifying a set having a relationship between a plurality i of bioinformatics data samples of each of a plurality of measurement variables g, based on a ranking method for sorting said data, said relationship being indicative of a biological state. Said computer-implemented method comprising sorting the bioinformatics data into a list l according to the ranking provided by the ranking method; computing an exchangeability score of subsets of variables of the ordered bioinformatics data, resulting in an exchangeability matrix V_l, carrying information about the exchangeability between the measurement variables g; expanding the ordered lists l of the data, based on the exchangeability matrix V_l; and comparing the expanded experimental list vector l_lwith another list vector comprising a plurality of measurement variables g and thereby determining the relationship between the samples, thereby reducing the plurality i into an identified set of data samples being indicative of the biological state, which enables improved analysis of the biological state. From the expanded list vector it is straightforward to extract biologically relevant information that is stabilized compared to the biological information obtained from the original experiment.

The method according to an embodiment may be summarized as follows.

A computer-implemented method is provided for identifying a set having a relationship between a plurality i of bioinformatics data samples of each of a plurality of measurement variables based on a ranking method for sorting said data, said relationship being indicative of a biological state. The computer-implemented method comprises sorting the bioinformatics data into a list l according to the ranking provided by the ranking method. The computer-implemented method further comprises computing an exchangeability score of each pair (or subsets) of variables. The exchangeability score may be defined as a total exchangeability variation score (PS), mean exchangeability score (ES(mean)), maximal exchangeability score (ES(max)), one-sided mean exchangeability score (oES(mean)), one-sided maximal exchangeability score (oES(max)), or normalized variants thereof. These different types of exchangeability scores are defined as follows:

Consider a probability triple (Ω, , P) and let X₁, . . . , X_mdenote random variables on n, taking values in some space . Given X₁, . . . , X_mwe define the multivariate random variable

X1×…×Xm:Ω→×…×m

by X₁× . . . ×X_m(ω)=(X₁(ω), . . . , X_m(ω)). To each random variable X₁× . . . ×X_mthere is an associated measure Pr_{X1× . . . ×Xm}defined by

Pr_{X1× . . . ×Xm}(A)=P({ω∈Ω;X₁× . . . ×X_m(ω)∈A})

for ⊂× . . . ×.

Throughout this specification exchangeability of pairs of discrete random variables with values in ={1, . . . , M} are mainly considered. In this special case, since X₁×X₂is the reflection of X₂×X₁with respect to the line {(x,y)∈²; y=x}, a total exchangeability variation score (PS), mean exchangeability score (ES(mean)), maximal exchangeability score (ES(max)) is defined as follows. First the corresponding distances are defined according to

P_X1×X2^Var=|Pr_X1×X2−Pr_X1×X1|, where |μ| denotes the total variation of a (real-valued) measure μ

and

EDX1×X2max=HDρ(suppX1×X2,suppX2×X1)ρ((1,1),(M,M)),

where HD denotes the Hausdorff distance between the two sets defined by

HDρ(A,B):=max(supa∈Adistρ({a},B),supb∈Bdistρ({b},A)).andEDX1×X2mean=X1×X2[distρ(X1×X2,suppX2×X1)]ρ((1,1),(M,M))

The one-sided mean exchangeability distance for a pair of discrete random variables (X₁, X₂) is defined by

oEDX1×X2mean:=X1×X2[distρ(X1×X2,suppX2×X1⋂R(X1×X2))]ρ((1,2),(M-1,M))

if supp X₂×X₁∩R(X₁×X₂)(ω)≠ø for all ω∈Ω, and oED_X1×X2^mean=1 otherwise.

It is also possible to define a one-sided variant of the maximal exchangeability distance (oED_Xi×Xj^max) in an analogous manner. Due to the normalization factors introduced in the estimates of weak exchangeability, all measures introduced above attain only values in [0, 1].

This allows us to define similarity measures (exchangeability scores) for pairs of genes as follows:

PS_Xi×Xj^Var=1−P_Xi×Xj^Var,ES_Xi×Xj^mean=1−ED_Xi×Xj^mean,ES_Xi×Xj^max=1−ED_Xi×Xj^max,

oES_Xi×Xj^mean=1−oED_Xi×Xj^mean,oES_Xi×Xj^max=1−oED_Xi×Xj^max.

Finally, normalized values of the exchangeability scores are defined by comparing them to the corresponding values for pairs of random variables with some pre-specified distribution representing a null hypothesis of no association. In this specification, the main focus is on weak exchangeability of pairs of discrete random variables, in which case it is natural to compare to pairs of random variables (Y₁, Y₂) uniformly distributed on a set S⊂× with cardinality equal to that of suppX₁×X₂. We show only the normalization for oES_Xi×Xj^mean, as the other scores can be normalized analogously.

The normalized one-sided mean exchangeability score for a pair of discrete random variables (X₁,X₂) is defined by

noESX1×X2mean=(oESX1×X2mean-oESY1×Y2mean1-oESY1×Y2mean)+

where Y₁×Y₂is a uniformly distributed random variable on a set S⊂× where |S|=|suppX₁×X₂|, and (a)₊=max(a,0).

This finally resulting in an M×M exchangeability matrix V_l, wherein each entry (V_l)_ijis the normalized exchangeability score of measurement variables g_iand g_j, and carrying information about the exchangeability between each of the plurality of measurement variables g.

In an embodiment of the computer-implemented method according to the first aspect, the third step of expanding comprises the steps of computing a diagonal position matrix A_l, by defining diagonal elements as:

(A_i)_ii=u(r_i)

where r_iis the ranking statistic of variable i used in sorting the data into an ordered list l, and u:→ is a monotone function; computing a diagonal global weight matrix W_l, by weighting the data; and computing an expanded experimental list vector l_laccording to a function:

f(A_l,V_l,W_l)=l_l.

In an embodiment of the computer-implemented method according to the first aspect, the computing of an expanded experimental list vector l_lis done according to the formula:

l_l=((l_l)₁, . . . ,(l_l)_M)

by letting (l_l)_i=h((G_l)_i), where G_l=A_lV_lW_l; and (G_l)_iis the i:th column of G_land h:^M→.

In an embodiment of the computer-implemented method according to the first aspect, the subsets of variables are pairs of variables.

The pairs of variables may be pairs of variables with an exchangeability score exceeding a predefined threshold value.

In an embodiment of the computer-implemented method according to the first aspect, the ranking method is signal-to-noise ratio (SNR), fold change of average expression value, t-test, ANOVA, and non-parametric tests.

In an embodiment of the computer-implemented method according to the first aspect, the second step of computing an exchangeability score of each pair of variables is total exchangeability variation score, PS, mean exchangeability score, ES(mean), maximal exchangeability score, ES(max), one-sided mean exchangeability score, oES(mean), one-sided maximal exchangeability score, oES(max), or normalized variants thereof.

In an embodiment of the computer-implemented method according to the first aspect, computing the global weight matrix is according to the formula:

(W)=log(N+1Ni+1)

where N is the number of lists used as reference data and N, is the number of reference lists comprising g_i.

In an embodiment of the computer-implemented method according to the first aspect, h is h₁(x)=Σ_i=1^Mx_i; h₂(x)=min_1≦i≦M|x_i|; or h₃(x)=∥x∥ for some norm ∥•∥ on ^M.

In an embodiment of the computer-implemented method according to the first aspect, the fourth step of comparing the expanded experimental list vector l_l, called l_l1, with the other list vector comprising a plurality of measurement variables g, called l_l2, is done by computing a similarity coefficient s(l₁, l₂) according to the formula

s(1,2)=l1·l2l1l2;

and a dissimilarity coefficient d(l₁, l₂) according to the formula

d(l₁,l₂)=1−s(l₁,l₂).

In an embodiment of the computer-implemented method according to the first aspect, the measurement variables are of the same kind.

In another embodiment of the computer-implemented method according to the first aspect, the measurement variables are of different kinds.

The bioinformatics data may be any kind of data representing a plurality i of bioinformatics data samples of each of a plurality of measurement variables g. Preferably, the bioinformatics data is bioinformatics expression data, such as DNA microarray data, RNA-seq data or real-time PCR data; microRNA data, such as RNA microarray data, RNA-seq data or real-time PCR data; DNA methylation data, such as DNA methylation microarray data; or protein expression data, such as protein microarray data, antibody array data, 2-D gel data or LC-MS data. As is appreciate by a person skilled in the art, the context of the data is what makes the data carry biological information, hence bioinformatics data. This is the subject of the ranking method, which ultimately enables indication of a biological state.

According to a second aspect, a computer program product is provided, comprising computer program code means for executing the method according the first aspect, when said computer program code means are run by an electronic device having computer capabilities.

According to a third aspect, a computer readable medium is provided, having stored thereon a computer program product comprising computer program code means for executing the method according to the first aspect, when said computer program code means are run by an electronic device having computer capabilities.

According to a fourth aspect, a computer is provided, configured to perform the method according the first aspect.

According to a fifth aspect, an apparatus for identifying a set having a relationship between a plurality i of bioinformatics data samples of each of a plurality of measurement variables g, based on a ranking method for sorting said data, said relationship being indicative of a biological state is provided. Said apparatus comprising a memory for storing data and instructions and a controller, wherein said controller is configured to sort the bioinformatics data into a list l according to the ranking provided by the ranking method; compute an exchangeability score of subsets of variables of the ordered bioinformatics data, resulting in an exchangeability matrix V_l, carrying information about the exchangeability between the measurement variables g; expand the ordered lists l of the data, based on the exchangeability matrix V_l; and compare the expanded experimental list vector l_lwith another list vector comprising a plurality of measurement variables g and thereby determining the relationship between the samples and thereby reducing the plurality i into an identified set of data samples being indicative of the biological state, which enables improved analysis of the biological state.

According to a sixth aspect, an apparatus for identifying a set having a relationship between a plurality i of bioinformatics data samples of each of a plurality of measurement variables g, based on a ranking method for sorting said data, said relationship being indicative of a biological state is provided. Said apparatus comprising a memory for storing data and instructions; a controller for executing said instructions; a sorter for sorting the bioinformatics data into a list l according to the ranking provided by the ranking method; a determiner for determining an exchangeability score of subsets of variables of the ordered bioinformatics data, resulting in an exchangeability matrix V_l, carrying information about the exchangeability between the measurement variables g; an expander for expanding the ordered lists l of the data, based on the exchangeability matrix V_l; and a comparator for comparing the expanded experimental list vector l_lwith another list vector comprising a plurality of measurement variables g and thereby determining the relationship between the samples and thereby reducing the plurality i into an identified set of data samples being indicative of the biological state for use in a diagnostic step for making a diagnosis on the biological state, which enables improved analysis of the biological state.

Embodiments of the present invention have the advantage over the prior art that it enables relatively robust conclusions to be drawn from the data, since they provide stability with respect to re-sampling of data. A technical advantage of this is ability to determine relationship between samples, which is indicative of a biological state, and which enables improved analysis of the biological state.

Several embodiments of the present invention will be described in more detail below with reference to the accompanying drawings in order for those skilled in the art to be able to carry out the invention. The invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art. The embodiments do not limit the invention, but the scope of the invention is only limited by the appended patent claims. Furthermore, the terminology used in the detailed description of the particular embodiments illustrated in the accompanying drawings is not intended to be limiting of the invention.

Exchangeability of Random Variables

Consider a probability triple (Ω, , P) and let X₁, . . . , X_mdenote random variables on n, taking values in some space . Given X₁, . . . , X_mwe define the multivariate random variable

X1×…×Xm:Ω->×…×m

by X₁× . . . ×X_m(ω)=(X₁(ω), . . . , X_m(ω)). To each random variable X₁× . . . ×X_mthere is an associated measure Pr_{X1× . . . ×Xm}defined by

Pr_{X1× . . . ×Xm}(A)=P({ω∈Ω;X₁× . . . ×X_m(ω)∈A})

for ⊂× . . . ×.

The concept of exchangeability of random variables was made popular by de Finetti in the 1930's. Conventionally, a finite sequence (X₁, . . . , X_m) of random variables is called exchangeable if their joint distribution is invariant under permutation of X₁, . . . , X_m, i.e. if

Pr_{X1× . . . ×Xm}=Pr_{Xπ(1)× . . . ×Xπ(m)}

for each π∈S_m(the group of permutations of {1, . . . , m}). This means that from a statistical point of view the order of the variables in the product is completely irrelevant. From this definition, it is clear that any sequence of independent and identically distributed (i.i.d.) random variables is exchangeable, but the reverse implication is false. The definition of exchangeability given above is rather strong, and we introduce a much weaker notion of exchangeability as follows:

The finite sequence of random variables (X₁, . . . , X_m) is weakly exchangeable if the null sets of Pr_{X1× . . . ×Xm are invariant under permutations, i.e.}

Pr_{Xπ(1)× . . . ×Xπ(m)}<<Pr_{Xτ(1)× . . . ×Xτ(m)}

for all π, τ∈S_m. Here, μ<<ν denotes that the positive measure μ is absolutely continuous with respect to the positive measure ν.

It is clear that a finite sequence of random variables (X₁, . . . , X_m) that is exchangeable is weakly exchangeable, but that the opposite implication is false in general.

Measure of Exchangeability

In this section we will discuss some ways to quantify the degree of exchangeability for a sequence of random variables.

Given a finite sequence of random variables (X₁, . . . , X_m), the total exchangeability variation is given by

PX1×…×XmVar:=1Sm-1∑π∈SmPrXπ(1)×…×Xπ(m)-1Sm∑τ∈SmPrXτ(1)×…×Xτ(m).

Here, |μ| denotes the total variation of the (real-valued) measure μ.

The total exchangeability variation is a positive measure and we note that P_{X1× . . . ×Xm}^Var(Ω)=0 if and only if the sequence (X₁, . . . , X_m) is exchangeable.

We now turn to a discrete probability space (Ω, , P), where Ω is a finite set, =2^Ω is the σ-algebra consisting of all events and P:→[0,1] is a probability measure.

We let X₁, . . . , X_mbe random variables on Ω taking values in :={1, . . . , M}. The support of the random variable X₁× . . . ×X_mis defined by

suppX₁× . . . ×X_m:={(q_h, . . . ,q_m)∈× . . . ×;Pr_{X1× . . . ×Xm}({(q₁, . . . ,q_m)})>0}

For finite sets of discrete random variables, Definition 1 implies that (X₁, . . . , X_m) is weakly exchangeable if

supp(X₁× . . . ×X_m)=supp(X_π(1)× . . . ×X_π(m))

for all Therefore, to quantify the degree of weak exchangeability for a set of discrete random variables we will compare the support of the joint distributions.

Let ρ:(× . . . ×)×(× . . . )→ be a metric and define the distance between two sets A, B⊂(× . . . ×) by

dist_ρ(A,B):=min_a∈A,b∈Bρ(a,b).

Furthermore, define the Hausdorff distance between the two sets by

HD_ρ(A,B):=max(sup_a∈Adist_ρ({a},B),sup_b∈Bdist_ρ({b},A)).

Given a finite sequence of discrete random variables (X₁, . . . , X_m), the maximal exchangeability distance is given by

EDX1×…×Xmmax=∑π∈Sm∑τ∈SmHDρ(suppXπ(1)×…×Xπ(m),suppXτ(1)×…×Xτ(m))ρ(1m,Mm)Sm(Sm-1)

and the mean exchangeability distance is given by

EDX1×…×Xmmean=∑π∈Sm∑τ∈SmXπ(1)×…×Xπ(m)[distρ(Xπ(1)×…×Xπ(m),suppXτ(1)×…×Xτ(m))]ρ(1m,Mm)Sm(Sm-1)

Here, l_m=(l, . . . , l) and M_m=(M, . . . , M).

Clearly, ED_{X1× . . . ×Xm}^max=ED_{Xπ(1)× . . . ×Xπ(m)}^maxfor any π∈S_m, and ED_X1^{× . . . ×X}_m^max=0 if and only if (X₁, . . . , X_m) is weakly exchangeable, and the same is true for ED_X1^{× . . . ×X}_m^max.

For the purpose of this application, we will mainly consider exchangeability of pairs of discrete random variables with values in ={1, . . . , M}. In this special case, since X₁×X₂is the reflection of X₂×X₁with respect to the line {(x,y)∈²; y=x}, we get

PX1×X2Var=PrX1×X2-PrX2×X1EDX1×X2max=HDρ(suppX1×X2,suppX2×X1)ρ((1,1),(M,M))EDX1×X2mean=X1×X2[distρ(X1×X2,suppX2×X1)]ρ((1,1),(M,M))

The Exchangeability Plot

To illustrate the degree of weak exchangeability of a pair of discrete random variables visually we propose the exchangeability plot. The exchangeability plot for the random variables X₁and X₂is obtained by depicting both suppX₁×X₂and suppX₂×X₁in the same figure. A pair of random variables (X₁,X₂) is weakly exchangeable if the two sets overlap completely. FIG. 2 shows exchangeability plots for two pairs of random variables. The pair in FIG. 2A is weakly exchangeable, while the pair in FIG. 2B is not.

The Exchangeability of Genes

In an embodiment, wherein the bioinformatics data is gene expression data we assume that we are given a universal set of M genes, denoted ={g₁, . . . , g_m}. Data in the form of a population (e.g. cancer patients and healthy control subjects) and a variable ranking method (e.g. a t-test contrasting the two groups in the population) is the starting point, called experiment, of the method. The sample space Ω consists of all possible rankings of the M genes, and the random variables X₁, . . . , X_M:Ω→{1, . . . , M} represent the ranking positions of the genes in . A finite set of genes g_i1, . . . , g_im is said to be (weakly) exchangeable if the corresponding sequence of random variables (Xi₁, . . . , Xi_m) is (weakly) exchangeable. Intuitively, a set of genes is exchangeable if their positions in the variable ranking can be interchanged without changing the biological interpretation of the ranking.

Of course, we do not know the measures Pr_Xi for the variables in practice, so these have to be estimated. In this description we use subsampling to generate a collection of B data sets for which we compute gene rankings. From the B gene rankings, we construct a position vector S_ifor each gene g_iby collecting all positions of the gene in the B rankings. The elements of a position vector S_iare then samples of the random variable X_i. Combining two position vectors S_iand S_jgives samples of the variables X_i×X_jand X_j×X_iwhich can be used to obtain estimates of P_Xi×Xj^Var, ED_Xi×Xj^maxand ED_Xi×Xj^mean.

We now introduce another measure of exchangeability which is especially adapted to the study of ranked gene lists. The rationale behind this measure is that we may not want two genes to obtain a small exchangeability distance if they always appear in the same order in the rankings. For example, a gene which is always ranked first is not highly exchangeable with a gene that is always ranked second, since the first gene is clearly more strongly related to the response than the second. The new measure penalizes such situations in the computation of the exchangeability distance. For a pair of random variables (X₁, X₂) we define a new set-valued random variable on Ω by

R(X₁×X₂)(ω):={(x,y)∈×;sign(x−y)=(X₁(ω)−X₂(ω))}.

The one-sided mean exchangeability distance for a pair of discrete random variables (X₁, X₂) is defined by

oEDX1×X2mean=X1×X2[distρ(X1×X2,suppX2×X1⋂R(X1×X2))]ρ((1,2),(M-1,M))

if suppX₂×X₁∩R(X₁×X₂)(ω)≠ø for all ω∈Ω, and oED_X1×X2=1 otherwise.

It is also possible to define a one-sided variant of the maximal exchangeability distance (oED_Xi×Xj^max) in an analogous manner. We note that due to the normalization factors introduced in the estimates of weak exchangeability, all measures introduced above attain only values in [0, 1]. This allows us to define similarity measures (exchangeability scores) for pairs of genes as follows:

PS_Xi×Xj^Var=1−P_Xi×Xj^Var,ES_Xi×Xj^mean=1−ED_Xi×Xj^mean,ES_Xi×Xj^max=1−ED_Xi×Xj^max,

oES_Xi×Xj^mean=1−oED_Xi×Xj^mean,oES_Xi×Xj^max=1−oED_Xi×Xj^max.

Finally, we define normalized values of the exchangeability scores by comparing them to the corresponding values for pairs of random variables with some pre-specified distribution representing a null hypothesis of no association. In this paper, the main focus is on weak exchangeability of pairs of discrete random variables, in which case it is natural to compare to pairs of random variables (Y₁, Y₂) uniformly distributed on a set S⊂× with cardinality equal to that of suppX₁×X₂. We show only the normalization for oES_Xi×Xj^mean, the other scores can be normalized analogously.

The normalized one-sided mean exchangeability score for a pair of discrete random variables (X₁,X₂) is defined by

noESX1×X2mean=(oESX1×X2mean-oESY1×Y2mean1-oESY1×Y2mean)+

where Y₁×Y₂is a uniformly distributed random variable on a set S⊂× where |S|=|suppX₁×X₂|, and (a)₊=max(a,0).

We note that the measures of exchangeability depend on the number of genes in the ranking (M). For two genes having the exchangeability plot shown in FIG. 2A we get noES_X1×X2^mean=1.0 irrespective of the number of genes since in this case, the distance between any value of X₁×X₂and suppX₂×X₁is zero. For the exchangeability plot in FIG. 2B we obtain noES_X1×X2^mean=0.17 if M=15 and noES_X1×X2^mean=0.99 if M=1,000.

General Idea

In this section, we present a general framework for list representation and comparison, with a focus on an embodiment of the present invention applicable to analysis of biologic or bioinformatics data, and in particular analysis of gene expression data.

However, as appreciated by a person skilled in the art, the bioinformatics data may be any kind of data representing a plurality i of bioinformatics data samples of each of a plurality of measurement variables g. Preferably, the bioinformatics data is bioinformatics expression data, such as DNA microarray data, RNA-seq data or real-time PCR data; microRNA data, such as RNA microarray data, RNA-seq data or real-time PCR data; DNA methylation data, such as DNA methylation microarray data; or protein expression data, such as protein microarray data, antibody array data, 2-D gel data or LC-MS data. As is appreciate by a person skilled in the art, the context of the data is what makes the data carry biological information, hence bioinformatics data. This is the subject of the ranking method, which ultimately enables indication of a biological state.

The lists are represented as vectors in ^M, where the entry in position i gives the contribution of gene g_i. The vector representation allows us to compare both ranked and unranked gene lists within the same framework, using one of the many similarity or dissimilarity measures available to compare vectors in ^M. This is an advantage compared to existing methods for list comparison, which are specifically designed to compare certain types of lists. Our framework also provides a way to determine which genes are most important for explaining the similarity between two lists. Assume for example that some measure based on the scalar product in ^Mis used to measure the similarity between two vectors x and y. Then the value of x_iy_iis a measure of the influence of the i′th variable on the similarity between the two lists. Finally, ordering the genes by their weights in the vector provides a new ranking of the genes.

In an embodiment, as above, we have a universal set of M genes, ={g₁, . . . , g_M}, where the genes are indexed in a fixed (but otherwise arbitrary) fashion. The universal set can be e.g. all genes on a microarray chip. An ordered (or unordered) gene list is then an ordered (or unordered) subset of the universal set.

Let l⊂ denote a list. If l is ordered and if gene g_iis contained in l, we denote its position by π_l(i). If g_i∉l, we define π_l(i)=0. For an unordered list, we let π_l(i):=χ_l(g_i), where χ_lis the characteristic function of the set l. Given a list l we define a corresponding list matrix G_las the product of three basic M×M matrices;

G_l:=A_lV_lW_l

The three basic matrices are designed to represent different characteristics of l. We call A_lthe position matrix, V_lthe exchangeability matrix and W_lthe global weight matrix. From the list matrix we form a list vector l_l:=((l_l)₁, . . . , (l_l)_M), by letting (l_l)_i:=h((G_l)_i) where (G_l)_iis the i:th column of G_land h:^M→ is a summarization function, e.g. a norm. The list vector will be used as the vector representation of the list. Once all lists of interest are represented by vectors in ^M, we can define the similarity between them e.g. as the cosine of the angle between the corresponding list vectors, i.e.

s(1,2)=l1·l2l1l2

where · denotes the inner product in ^M, and we can obtain a dissimilarity coefficient as

d(l₁,l₂)=1−s(l₁,l₂)

Choosing A_l, V_l, W_l, h and the (dis)similarity coefficient on ^Msuitably, most methods currently available for list comparison fit into this general framework.

The Position Matrix

The position matrix A_lis defined as a diagonal matrix that contains information about the type of list (ordered or unordered) and the positions of the genes within the list. We define the diagonal elements (the position values) by (A_l)_ii=u (r_i) where r_iis the ranking statistic of gene g_i, and u:→ is a monotone function.

For all genes not in the list, we define (A_l)_ii=0. In this paper we mainly use rank-based functions, i.e. letting (A_l)_ii=ũ(π_l(i)) where ũ:→ is a decreasing function of π_l(i) if π_l(i)>0, and ũ(0)=0.

This means that the diagonal elements corresponding to genes in the top of the list l are high, while the genes further down in the list obtain lower values. Furthermore, all genes not in the list have position value zero. We note that in some cases, other choices of position values may be better suited for unordered lists, where it may be desirable to give the genes different weights, e.g. depending on some external criterion, even though there is no specified ordering.

The Exchangeability Matrix

The exchangeability matrix V_lcarries information about the exchangeability between the genes in in the specific experiment giving the list l. In an embodiment, we define the entry (V_l)ij to be the estimated normalized one-sided mean exchangeability score of g_iand g_j(i.e. _Xi×Xj^mean), so the diagonal elements are always 1. If V_lis diagonal, i.e. V_l=IM, then the only non-zero elements in the list vector are those corresponding to genes that are actually contained in l and consequently only the genes that are present in the list affect its vector representation. However, if V_lis not diagonal, there is a possibility that the vector representation of the list is extended, i.e. that it contains non-zero entries for genes which are not themselves present in the list, but are exchangeable with some gene in the list. The (absolute) weight of a gene in the list can also be increased if it is highly exchangeable with a gene with a higher (absolute) position value, since the high exchangeability indicates that the genes could as well have switched positions without changing the interpretation of the list.

We note that the general framework for list representation supports any matrix of gene similarities in the place of V_l. Thus, in an embodiment, V_lcould be defined from some kind of expert biological knowledge, e.g. concerning which genes are related to the same biological function. This could be used for example when comparing lists from different experiments to each other. In another embodiment, the positive part of the correlation matrix is used in place of V_l, since a high correlation between the expression levels of two genes is often considered to indicate similar biological functions of the genes.

The Global Weight Matrix

The global weight matrix W_lis a diagonal matrix that permits weighting the influence of the genes differently, depending on some informativeness or reliability estimate. For example, we may wish to downweight the influence of a gene that has a high probability to be present in an arbitrarily chosen list, since this gene is unspecific and may not give much relevant information about the similarity between a pair of lists.

The specific embodiment following description focuses on an embodiment of the present invention applicable to analysis of biologic data, and in particular analysis of gene expression data. However, it will be appreciated that the invention is not limited to this application but may be applied to many other data including for example meteorology and astronomy.

In an embodiment according to FIG. 1, bioinformatics data comprising a plurality of samples was obtained from a microarray study of lung adenocarcinoma patients, as provided in Bhattacharjee, A., Richards, W.G., Staunton, J., Li, C., Monti, S., Vasa, P., Ladd, C., Beheshti, J., Bueno, R., Gillette, M., Loda, M., Weber, G., Mark, E. J., Lander, E. S., Wong, W., Johnson, B. E., Golub, T. R., Sugarbaker, D. J., and Meyerson, M. (2001), Classification of human lung carcinomas by mRNA expression profiling reveals distinct adenocarcinoma subclasses. A ranking method, signal-to-noise ratio, was also provided.

The abovementioned bioinformatics data set consists of gene expression profiles from 31 patients with “good outcome” and 31 patients with “poor outcome”. The data set was downloaded from http://www.broadinstitute.org/gsea/datasets.jsp (Lung_Boston_collapsed.gct), on Jan. 11, 2011. In the downloaded file, original Affymetrix probe set IDs have been replaced with gene symbols and all probe sets matching to the same gene have been collapsed into one variable.

The embodiment according to FIG. 1 is further visualized in FIG. 10. Here, the large group of persons on the left side indicates the population that we are interested in studying. For example, it can be patients with breast cancer and healthy controls, or lung cancer patients with good or poor prognosis, respectively. It can also be for example rats treated with different drugs, or a homogeneous collection of plants that we want to characterize. From the population, we draw a sample of patients, rats, etc. to include in the experimental study. Reference (A) represents the ranked list of variables obtained in step 110.

In FIG. 10, step 120 corresponds to the exchangeability computation step. The subsamples are subcollections of the whole sample. As will be appreciated by a person skilled in the art, the right way to construct them depends on the application, specifically it depends on what we want to stabilize the ranked list with respect to. The subsamples are used to create ranked variable lists (B_—1 to B_—20). From these ranked variable lists we compute the exchangeabilities among the variables. These are collected into the exchangeability matrix.

Step 130 represents the expansion step. From the original ranked list (A) we compute a position matrix for the variables. We are also given a weight matrix (possibly the identity matrix) that contains prior knowledge about the reliability of the variables. The position matrix, the exchangeability matrix, and the weight matrix are multiplied and the values in the resulting matrix are summarized into an expanded list vector (C).

In step 140 the expanded list vector (C) is compared to other list vectors (D_—1 to D_m). These list vectors can be obtained in the same way as (C) or in any other way. An important special case is when the list vectors are used to represent gene sets. In this case, the list vectors (D_—1) to (D_m) are binary vectors. By comparing the expanded list vector (C) to the list vectors (D_—1) to (D_m) new biological insights concerning the population can be generated. The goal of the process 110-140 is to obtain a stabilization effect so that these biological insights are not sensitive to the specific choice of the sample.

In the data set, 62 patient observations (31 with “good” outcome and 31 with “poor” outcome) thus gave 62 samples in form of observation vectors, comprising 5,217 measurement variables g, i.e. studied genes.

Signal-to-noise ratio (SNR) for each variable in the data was computed as

SNR=mpoor-mgoodσpoor+σgood

where m_poorand m_goodis the mean value of the expression values in the “poor outcome” and “good outcome” group, respectively and σ_poorand σ_goodis the standard deviation of the expression values in the “poor outcome” and “good outcome” group, respectively. The signal-to-noise ratios provide information about the biological importance of the variables with respect to the current task.

The variables were then sorted 110 by decreasing value of the computed statistic. In this way, each variable i obtains a position in the ranking. The variable with the highest (positive) SNR value obtained position 1, and the variable with the lowest (most negative) SNR value obtained position 5217. The experimental list is now the set of variables sorted in the indicated order, i.e. an ordered list l. The goal is to analyze the obtained list to enable robust biological conclusions. In FIG. 10, the experimental ranked list is indicated by (A). The top-ranked genes in this ordered list are DBP, ALDH2, PBP, CYP4B1, CAT, EEF1A1, KIAA0256, IL12RB1, ATP5A1 and RPLP1. Table 1 shows the top-ranked genes in this ordered list. The ranking of the variables by their SNRs is unstable and sensitive to the initial choice of patients to include in the study. Thus, in order to provide biological insights of a general character, it needs to be stabilized.

Next, exchangeabilities between each pair of variables were computed 120 and collected in a 5217×5217 matrix. The exchangeabilities were computed in two steps.

Firstly, position vectors were computed by generating 20 slightly modified data sets by sub sampling the original data set, each time keeping ⅔ of the patients (randomly chosen) from each outcome group. For each modified data set, a signal-to-noise ratio (SNR) is computed for each variable as shown above, and used to rank the variables. In FIG. 10, the ranked lists are indicated by B_—1, . . . , B_—20. Table 2 shows the top-ranked genes in three of the ranked lists obtained in this way. Note that there are considerable discrepancies between the top-ranked genes in these lists, further indicating that the variable ranking is highly dependent on the choice of included samples and needs stabilization.

The positions that a variable obtains in the 20 rankings are collected in a position vector for the variable. The position vector for variable i is denoted by S_i(1≦i≦5217).

Secondly, for each pair of variables, the exchangeability score was estimated according to the following. For each pair of variables the joint position vectors were formed as:

S_i×j{((S_i)_k,(S_j)_k)}_k=1²⁰

and

S_j×i={((S_j)_k,(S_i)_k)}_k=1²⁰

These vectors describe the observed positions for the ordered pairs of variables. Then, the one-sided mean exchangeability distance between the two variables was estimated by

o.EDXi×Xjmean=∑k=120120distρ({((Si)k,(Sj)k)},{((Sj)v,(Si)v)}v=120⋂Rij(k))ρ((1,2),(5216,5217))

In this equation, ρ is the standard Euclidean metric in ²and, for two sets A and B in ², dist_ρ(A, B)=min_a∈A,b∈Bρ(a, b) were defined.

Finally,

R_ij(k)={(x, y)∈×; sign(x−y)=sign ((S_i)_k−(S_j)_k)} denotes the points in {1, . . . 5217}×{1, . . . 5217} that are on the same side of the diagonal as the k′th point. In case {((S_j)_ν,(S_i)_ν)}_ν=1^B∩R_ij(k)=0 for some k, the one-sided mean exchangeability distance between the two variables was defined to be 1.

From the one-sided mean exchangeability distance a one-sided exchangeability score was formed by

{dot over (o)}ES_X1×Xj^{mean 10}=1−oED_X1×Xj^mean.

The one-sided mean exchangeability score is normalized to give a normalized one-sided mean exchangeability score, by

n.oESXi×Xjmean=max(0,o.EDXi×Xjmean-o.EDYi×Yjmean1-o.EDYi×Yjmean).

The normalization constant ėES_Yi×Yj^meanis the estimated value of the one-sided mean exchangeability score for variable pairs uniformly distributed on a set S⊂{1, . . . 5217}×{1, . . . 5217} with 20 elements. The estimated value {dot over (n)}oES_Xi×Xj^meanis then put in position (i,j) in the exchangeability matrix V_l. A part of the exchangeability matrix according to an embodiment, containing the exchangeabilities among the first 10 variables is shown in FIG. 11. Note that the exchangeability between a variable and itself is always 1.

Next, the ordered list l of the data is expanded 130, thus combining the computed exchangeabilities with the information about the position of the variables in the experimental list to obtain a representation of the experimental list as a vector in ⁵²¹⁷. In an embodiment, this is done in several steps.

First, a diagonal position matrix A_lis computed 130a, where the entry in the i′th row and column is derived from the position of the i′th variable in the list of experimental data. For all variables with a positive SNR as computed above (i.e. the variables in the top of the list),

Aii=350(π(i)-1)2+350

is defined where π(i) is the position of variable i in the list. For variables with a negative SNR,

Aii=-350(5217-π(i))2+350

is defined. FIG. 12 shows the part of the position matrix corresponding to the first 10 variables according to an embodiment. Note that some variables have positive position values, indicating that they have positive SNRs, and others have negative position values.

Second, a global weight matrix W_lis computed 130b, which permits inclusion of a weight for each variable, which can be used to put more emphasis on e.g. genes that are known to be reliably measured. Since there is no such information in this example, the global weight matrix is taken to be the identity matrix.

Third, an expanded experimental list vector is computed 130c. The expanded list vector l_lis taken as the vector representation of the experimental list in ⁵²¹⁷, according to the formula:

l_l=((l_l)₁, . . . ,(l_l)_M)

• • by letting

(l_l)_i=h((G_l)_i)

where G_l=A_lV_lW_l, and where (G_l)_iis the i′th column of G_land h:→ is defined by letting the i′th element of the list vector be the entry with the largest magnitude from column i in the list matrix. The list vector is called expanded since the value in each position is constructed using information from all the variables in the data set. FIG. 12 shows the first 10 elements in the expanded list vector according to an embodiment. These elements indicate the stabilized importance of the first 10 genes in the considered application. Table 3 shows the 50 genes corresponding to the highest values in the expanded list matrix. This ranking can be interpreted as a stabilized variant of the original gene ranking, cf. (A) in FIG. 10. Compared to Table 1, some genes (e.g. IL12RB1, RPLP1 and PPP3CB) are ranked higher in the stabilized ranking, indicating that in this more general ranking, these genes obtain increased importance for the discrimination between good and poor outcome lung cancer patients. Conversely, some genes (e.g. CAT, EEF1A1 and ATP5A1) are ranked higher in the original ranking than in the stabilized ranking. This suggests that their top-ranking positions in the original ranking could be an effect of the specific sample and not as suitable to generalize to the population.

Finally, the expanded experimental list vector l_lis compared 140 to another list vector comprising a plurality i of measurement variables g, i.e. the gene sets.

To compare the experimental list vector l_lto the gene sets from MSigDB, vector representations of the gene sets are constructed as well. Since the gene sets are not generated directly from an experiment as the experimental list, there is not enough information to construct the exchangeability matrix as for the experimental list. Moreover, the gene sets are not ranked, i.e. the elements in a gene set do not have a specific order. Each gene set is therefore represented by a vector with 5217 elements, where the elements corresponding to the genes in the gene set are set to 1 and the others to 0. Thus, the list vector for a gene set is similar to the list vector in FIG. 12, but the elements are either 0 or 1, indicating inclusion or exclusion of the corresponding gene in the gene set.

Now, both the experimental list and the gene sets are represented by vectors of length 5217. A dissimilarity score is computed between two vectors by taking 1 minus the cosine of the angle between them, which is easily computed directly from the vectors. The dissimilarity score lies in the interval [0,2], where 0 is obtained only for identical list vectors. The output from this step is a list of dissimilarity scores (one for each gene set). The gene sets with the smallest dissimilarity score are most closely related to the genes in the top of the list (i.e. to the genes over expressed in the group with poor outcome relative to the group with good outcome), while the gene sets with the largest dissimilarity score are most closely related to the genes in the bottom of the list (i.e. to the genes over expressed in the group with good outcome relative to the group with poor outcome). In this way, conclusions can be drawn concerning which gene sets (and thereby which cellular functions) are most closely related to the two groups of observations.

These conclusions may be important for different therapeutical and/or diagnostical reasons, or for research purposes. The gene sets that are closest to the experimental list indicate e.g. biological pathways that are associated with the outcome of the patients. Thus, the results can give rise to new biological insights concerning the population, since a determination of a specific relationship between samples is enabled, which relationship is indicative of a biological state, and which enables improved analysis of the biological state.

We will now estimate the robustness of the biological conclusions drawn in the abovementioned embodiment as described above. The effect of the stabilization will be estimated by repeating the same analysis without introducing the exchangeabilities, and both methods will be contrasted with what is obtained from Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005). The robustness is estimated in the following way:

First, we generate ten slightly modified data sets from the original data through bootstrapping, i.e. sampling 31 observations from each outcome group (“good” and “poor”) with replacement. Hence, some patient may be included more than once in the same modified data set.

Second, we compute two list vectors for each bootstrapped data set. The first list vector, the expanded list vector, is computed by sorting 110, computing 120 an exchangeability score and expanding 130, as described above. The second list vector, the non-expanded list vector, is obtained in the same way but using the identity matrix instead of the exchangeability matrix V_l. In this way, each entry in the list vector is affected only by the position of the corresponding variable.

Third, we compute the distances between the two list vectors and each of the gene sets by comparing 140 them for each of the ten bootstrapped data sets. The ten gene sets that have the shortest distance to the respective list vector and the ten gene sets that have the longest distance to the respective list vector are stored. For each pair of bootstrapped data sets, we compute the overlap between the collections of gene sets with shortest distances as well as the overlap between the collections of gene sets with longest distances to each of the two types of experimental list vectors (expanded and non-expanded).

For each of the ten bootstrapped data sets, we also apply GSEA using the R package provided by the authors at http://www.broadinstitute.org/gsea/downloads.jsp. For each data set we find the gene sets most related to each of the outcome groups and then we compute the overlap of these collections of gene sets for each pair of bootstrapped data sets. The GSEA algorithm uses a weighted Kolmogorov-Smirnov statistic to define an enrichment score (ES) for each gene set. These are normalized to eliminate the possible effect of the size of the gene sets, yielding normalized enrichment scores (NES) for the gene sets (denoted GSEANES).

The result is seen in FIGS. 3-5. For each gene set collection (related to top (A), i.e. poor outcome, and bottom (B), i.e. good outcome, of the experimental list, respectively), the figures show the size of the overlap between the extracted gene sets for any two of the bootstrapped data sets. FIG. 3 is the result according to an embodiment of the invention, FIG. 4 is the result without expanding 130 the list and FIG. 5 is the result obtained using GSEA. The numbers represent degree of overlap (1-10, where 10 is maximum overlap). As can be seen, the robustness of the gene set rankings obtaining by comparing the expanded list vectors in this respect is remarkably high compared to that of the other two methods meaning that sampling variations have a much smaller impact on the resulting biological conclusions.

The gene sets are ranked by the distance to the extended list vector. By studying the gene sets with the shortest (or longest, respectively) distance to the experimental list we find the biological processes most related to the experimental contrast. FIG. 6 gives the closest and farthest gene sets with the method according to an embodiment of the invention “The invention” as well as the genesets with the largest positive and negative normalized enrichment score from GSEA, “GSEA”. These are computed using the R implementation of the GSEA algorithm.

For all ranking metrics we also generate concordance plots for the gene set rankings obtained by ordering the gene sets by decreasing or increasing values of the respective ranking statistic, i.e. taking the direction of the association with the response into account. For each k the concordance plot shows the number of gene sets which are among the top-k ranked gene sets in all boot strapped data sets. The concordance plot where top gene sets are those positively associated with poor outcome is shown in FIG. 7, and the concordance plots for the rankings where the gene sets positively associated to good outcome are ranked in the top is shown in FIG. 8. In both plots, the number of shared gene sets is on the y axis and cutoff point is on the x axis.

In both FIGS. 7 and 8, it is seen that the gene set rankings obtaining by comparing the expanded list vectors, represented by a line (ExtendedSimilarity) are more stable than the other two methods, GSEA represented by dotted line (GSEANES) and unexpanded list represented by a dashed line (NotExtendedSimilarity).

Thus, the ranking metrics based on the similarity to the extended list vector consistently provide reproducible rankings.

Even if stability is a tractable property of a gene list from an experiment, it is not the only one. We can obtain a perfectly robust gene list by always assigning each gene the same position, but the resulting list is useless from the point of view of biological interpretation. In this part we extract gene lists in three different ways and compare the biological information contained in the top- and bottom-ranked genes in each list.

We compute three gene lists as follows. First, the non-expanded gene list which is the original list obtained by sorting the genes according to the SNR values. Then, the expanded gene list which is computed from the expanded list vector (obtained as described above) by sorting the genes according to their value in the expanded list vector. Finally, we compute an aggregated gene list by computing the median position in the position vector across 100 rankings from subsampled data sets for each gene and sorting the genes according to this, with the gene with lowest mean position in the top.

For each of the three gene lists, we extract a subset of the original data set consisting only of the top-k and bottom-k genes (for different choices of k). This small data set is fed into a centroid classifier (Schölkopf and Smola (2002): Learning with kernels, p. 4-5) and we estimate the classification ability of the subset by means of ten-fold cross-validation. Table 1 shows the estimated classification accuracy for top and bottom genes from the expanded gene list, the non-expanded gene list and the aggregated gene list, as well as the mean classification accuracy for 20 random gene lists. The best performing subset for each k is highlighted in bold.

Apparently, the increased stability does not come at the price of reduced biological meaning since the top and bottom genes from the expanded gene list perform best in the classification task.

Implementation

In an embodiment, a computer program product comprising computer program code means for executing the method according to some embodiments, when said computer program code means are run by an electronic device having computer capabilities.

In an embodiment according to FIG. 9, a computer readable medium 300 is disclosed, having stored thereon a computer program product comprising computer program code means for executing the method according some embodiments, when said computer program code means are run by an electronic device 200 having computer capabilities.

Computer program, software program, program product, or software, in the present context mean any expression, in any programming language, code or notation, of a set of instructions intended to cause a system having computer processing capabilities to perform a particular function directly after conversion to another language, code or notation.

The electronic device 200 having computer capabilities may be any device normally used for performing the involved tasks, e.g. a hardware, such as a processor with a memory. The processor may be any of variety of processors, such as Intel or AMD processors, CPUs, microprocessors, Programmable Intelligent Computer (PIC) microcontrollers, Digital Signal Processors (DSP), etc. However, the scope of the invention is not limited to these specific processors. The memory may be any memory capable of storing information, such as Random Access Memories (RAM) such as, Double Density RAM (DDR, DDR2), Single Density RAM (SDRAM), Static RAM (SRAM), Dynamic RAM (DRAM), Video RAM (VRAM), etc. The memory may also be a FLASH memory such as a USB, Compact Flash, SmartMedia, MMC memory, MemoryStick, SD Card, MiniSD, MicroSD, xD Card, TransFlash, and MicroDrive memory etc. However, the scope of the invention is not limited to these specific memories.

The invention may be implemented in any suitable form including hardware, software, firmware or any combination of these. However, preferably, the invention is implemented as computer software running on one or more data processors and/or digital signal processors. The elements and components of an embodiment of the invention may be physically, functionally and logically implemented in any suitable way. Indeed, the functionality may be implemented in a single unit, in a plurality of units or as part of other functional units. As such, the invention may be implemented in a single unit, or may be physically and functionally distributed between different units and processors.