Abstract:
Techniques for processing gene expression data and predicting gene relationships are provided. More specifically, a method for processing gene expression ratios may include discretizing a first expression ratio for a first gene and a second expression ratio for a second gene for each of a plurality of experiments into one of three indications, calculating the probability that the combination of the discretized expression ratios for each of the plurality of experiments arises by random chance, and determining whether the first gene and the second gene are related.

Description:
CROSS-REFERENCES TO RELATED APPLICATIONS 
     This application claims priority to U.S. Provisional No. 60/338,745 filed Nov. 5, 2001, which is incorporated by reference herein. 
    
    
     STATEMENT AS TO RIGHTS TO INVENTIONS MADE UNDER FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT 
     Not Applicable 
     REFERENCE TO A “SEQUENCE LISTING,” A TABLE, OR A COMPUTER PROGRAM LISTING APPENDIX SUBMITTED ON A COMPACT DISK 
     Not Applicable 
     COPYRIGHT NOTICE 
     A portion of this application contains computer codes, which are owned by California Institute of Technology. All rights have been preserved under the copyright protection, California Institute of Technology ©2002. 
     BACKGROUND OF THE INVENTION 
     The present invention is related to techniques for processing gene expression data and predicting gene relationships. 
     Since ancient time, humans have been searching for early diagnosis and effective treatment of diseases. In the last several hundred years, medical techniques, including blood sample analysis and physical imaging analysis, greatly improved the ability for diagnosis and treatment. But no cure has yet been found for many deadly and debilitating diseases, including cancer. To further improve humans&#39; physical health, gene diagnosis and gene therapy have been proposed and experimented. The human genome sequencing has been largely completed and the focus of medical research has been shifted to unveiling functions of various genes and their relations. 
     Relationship of genes may be inferred from their expression data, and may be organized in various forms, including clustering, dendrogram, and relevance network. Under the clustering method, genes are grouped into clusters based on their proximity in a multi-dimensional expression space, as measured by, among others, Euclidean distance, linear correlation, and non-linear correlation. Consequently, all genes are organized into a hierarchical structure. In contrast, under the dendrogram method, genes are comprehensively compared with a metric, and then added to a binary tree in order of decreasing correlation such that pairs with the highest correlation are closest in the tree. Both the clustering method and the dendrogram method suffers from several drawbacks. First, the clustering method based on Euclidean distance cannot handle genes with missing data, because incomplete expression vectors for these genes cannot be accurately oriented in the multi-dimensional space. Second, the dendrogram method and the clustering method based on Euclidean distance cannot identify negatively correlated genes and therefore ignore some important biological relations. Third, many clustering methods do not allow genes to belong to multiple clusters, and thus cannot accurately describe genes that are under the control of two or more different regulatory factors. Similarly, the dendrogram method has the same problem. 
     Under the relevance network method, a probability function, such as that based on mutual information or combinatorics, is used to estimate the probability that genes are independent. If the probability is low, the genes are predicted as having significant relation. Consequently, they are placed as connected in distinct gene networks with varying number of elements. Unlike the dendrogram method and many clustering methods, the relevance network method incorporates only genes with significant relations. In addition, a probability function is non-metric, so that the probability that two genes are related does not need to satisfy the “triangle inequality.” Under the “triangle inequality,” the distance, such as Euclidean distance, between genes A and B cannot exceed the sum of distance between genes A and C and distance between genes B and C. This requirement limits the ability of the metric method to describe gene relationships. Genes A and B may be weakly related or unrelated, even though both genes are regulated in part by the same gene C. Hence distance between genes A and C and distance between genes B and C may be small, while distance between genes A and B is large, exceeding the sum of distance between genes A and C and distance between genes B and C. Such relationships among genes A, B, and C cannot be adequately described by the metric method, but a non-metric method such as the relevance network using a probability function can provide an adequate description. 
     Some relevance network methods utilize a 2×2 matrix to calculate the probability that genes are independent, as exemplified in Walker, M. G. et al., Prediction of Gene Function by Genome-Scale Expression Analysis: Prostate Cancer-Associated Genes,  Genome Research  9(12): 1198-1203, 1999, Walker, M. G. et al., Pharmaceutical Target Discovery Using Guilt-by-Association: Schizophrenia and Parkinson&#39;s Disease Genes,  Proceedings of the International Conference on Intelligent Systems for Molecular Biology  146: 282-286, 1999, and Walker, M. G., Drug Target Discovery by Gene Expression Analysis: Cell Cycle Genes,  Current Cancer Drug Targets  1(1), 2001. These methods use binary expression data representing presence or absence of genes in a particular cell sample, and analyze gene relations based on presence or absence of genes in a common set of cell samples. These methods do not study continuous expression data representing regulatory effect of genes, but these data contain important information on gene relations. 
     As an improvement, Liang proposed a relevance network method that discretizes continuous gene expression data into binary states, i.e., on or off states of regulatory effect, as described in Liang S., Reveal, A General Reverse Engineering Algorithm for Inference of Genetic Network Architectures,  Pacific Symposium on Biocomputing  3:18-29 (1998). This method uses the binary representation for genes, which does not fully represent richness of expression data. Expression data may show that genes are up regulated, down regulated, or unchanged. 
     In contrast, Butte et al. used a relevance network method that discretizes the continuous expression data representing gene regulatory effect into n sub-ranges, e.g., 10 sub-ranges, as described in Butte A. J. el al., Mutual Information Relevance Networks: Functional Genomic Clustering Using Pairwise Entropy Measurements, in Altman, R., Dunker, K., Hunter L., Lauderdale K., Klein T. eds.,  Pacific Symposium on Biocomputing , at 418-429, (2000), Hawaii, World Scientific. This method uses narrow sub-ranges; therefore it cannot effectively filter out measurement noises that are associated with gene expression data, and may yield inaccurate gene correlation. 
     BRIEF SUMMARY OF THE INVENTION 
     According to the present invention, techniques for processing gene expression data and predicting gene relationships are provided. Merely by way of example, a method for predicting gene relationships may include providing a control sample and a plurality of experimental samples, and performing one experiment on each of the plurality of experimental samples. In addition, the method measures a first gene control signal corresponding to a first gene from the control sample, and measures a first gene experimental signal corresponding to the first gene from each of the plurality of experimental samples. Further, the method calculates a first expression ratio corresponding to the first gene for each of the plurality of experimental samples and discretizes the first expression ratio for each of the plurality of experimental samples into one of three indications representing the first gene up-regulated, not changed, or down-regulated. Similarly, the method repeats the above measurements and analyses for a second gene to obtain a discretized second expression ratio for each of the plurality of experimental samples. Based on the discretized expression ratios for the first gene and the second gene, the method calculates the probability that the combination of the discretized expression ratios for each of the plurality of experimental samples arises by random chance. Finally, the method determines whether the first gene and the second gene are related. 
     In another embodiment, a method for processing gene expression ratios may include gathering gene expression ratio data from a database. The gene expression data comprises a first expression ratio for a first gene for each of a plurality of experiments and a second expression ratio for a second gene for each of the plurality of experiments. The method discretizes the first expression ratio for each of the plurality of experiments into one of three indications, the three indications representing the first gene up-regulated, not changed, or down-regulated, and similarly discretizes the second expression ratio for each of the plurality of experiments. Based on the discretized expression ratios, the method calculates the probability that the combination of the discretized expression ratios for each of the plurality of experiments arises by random chance. Further, the method determines whether the first gene and the second gene are related. 
     In yet another embodiment, a method for processing gene expression ratios may include gathering a first plurality of gene expression data for a first plurality of experiments with a first control sample and calculating a first probability in response to the first plurality of gene expression data that a pair genes are related. In addition, the method includes gathering a second plurality of gene expression data for a second plurality of experiments with a second control sample and calculating a second probability in response to the second plurality of gene expression data that said pair genes are related. Further, the method includes calculating a combined probability in response to the first probability and the second probability, and determining whether said pair of genes are related. 
     In yet another embodiment of the present invention, a method for processing gene expression ratios includes gathering gene expression data. The gene expression data include a first gene control signal and a first gene experimental signal for a first gene for each of a plurality of experiments, and a second gene control signal and a second gene experimental signal for a second gene for each of the plurality of experiments. In addition, the method includes discretizing the first expression ratio for each of the plurality of experiments into one of three indications, the three indications representing the first gene up-regulated, not changed, or down-regulated, and discretizing the second expression ratio for each of the plurality of experiments into one of three indications, the three indications representing the second gene up-regulated, not changed, or down-regulated. Moreover, the method includes calculating the probability that the combination of the discretized first expression ratio for each of the plurality of experiments and the discretized second expression ratio for each of the plurality of experiments arises by random chance. Further, the method includes determining whether the first gene and the second gene are related. 
     In yet another embodiment of the present invention, a computer-readable medium including instructions for processing gene expression ratios includes one or more instructions for gathering gene expression data. The gene expression data comprises a first expression ratio for a first gene for each of a plurality of experiments and a second expression ratio for a second gene for each of the plurality of experiments. Moreover, the computer-readable medium includes one or more instructions for discretizing the first expression ratio for each of the plurality of experiments into one of three indications, the three indications representing the first gene up-regulated, not changed, or down-regulated. Similarly, the computer-readable medium also includes one or more instructions for discretizing the second expression ratio. Further, the computer-readable medium includes one or more instructions for calculating the probability that the combination of the discretized first expression ratio for each of the plurality of experiments and the discretized second expression ratio for each of the plurality of experiments arises by random chance. In addition, the computer-readable medium includes one or more instructions for determining whether the first gene and the second gene are related. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a block diagram of various steps for predicting gene relationships based on gene expression data according to an embodiment of the present invention. 
         FIG. 2  illustrates a control cell sample and n experimental cell samples. 
         FIG. 3  illustrates extraction of mRNA and reverse transcription of cDNA from control cell sample and experimental cell samples. 
         FIG. 4  illustrates measurement of gene expression data with microarray technology. 
         FIG. 5  is a block diagram of various steps for processing gene expression ratio between experimental sample and control sample. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     According to the present invention, techniques for predicting gene relationships and processing gene expression data are provided. 
       FIG. 1  is a block diagram for predicting gene relationships according to an embodiment of the present invention. The method includes eight steps: experimental sample processing step  110 , sample preparation step  120 , gene expression measurement step  130 , expression ratio calculation step  140 , discretization of expression ratio step  150 , co-expression table formation step  160 , probability calculation step  170 , gene relation prediction step  180 , and gene correlation database step  190 . 
     At experimental sample processing step  110 , experimental samples are processed in various studies. In each study, a control sample and N experimental samples are provided as shown in FIG.  2 . The control sample preferably has a cell type similar to those of experimental samples. Subsequently, the experimental samples are processed under various experimental conditions. For instance, experimental sample  1  undergoes experiment  1 , experimental sample  2  undergoes experiment  2 , . . . , experimental sample m undergoes experiment m, . . . , and experimental sample N undergoes experiment N. Experiments  1  through N may be performed under the same or different experimental conditions. 
     Next, at sample preparation step  120 , the control sample and the N experimental samples are prepared for gene expression measurement. For example, as shown in  FIG. 3 , the messenger RNA (mRNA) is isolated from the control sample and the N experimental samples, and subsequently converted to the reverse complementary DNA (cDNA) through reverse transcription. To allow comparison, the cDNA for the control sample and the experimental samples are labeled either directly with different flurochromes or indirectly with different epitopes. Alternatively, other detection schemes for cDNA may be employed, such as radioisotopes, electronic detection, or electrochemical detection. In addition, mRNA may be measured without conversion to cDNA. For example, mRNA may be directly measured at sample preparation step  120 . 
     Subsequently, at the gene expression measurement step  130 , the cDNAs for the control sample and the experimental samples are applied to a microarray, which contains DNA probes for genes. For example, as shown in  FIG. 4 , there are DNA probes for genes A, B, and C on a microarray  410 . If a cDNA sample possesses a nucleotide sequence complementary to that of a DNA probe, hybridization would occur. Upon hybridization, the fluorescent intensities at each spot of the microarray may be measured with a scanning system. Consequently, for each gene and for each experiment, the expression ratio r may be calculated in various means. For example, 
             r   =       log   2     ⁡     (         Signal   experiment     -     Background   experiment           Signal   control     -     Background   control         )               (     Equation   ⁢           ⁢   1     )             
 
where Signal experiment  is the fluorescent intensity for a gene, corresponding an experimental sample, and Signal control  is the fluorescent intensity for the control sample corresponding to the same gene. Background experiment  and Background control  are fluorescent intensities resulting from variability in overall fluorescence intensity and pixel intensity, with respect to experimental sample and control sample respectively.
 
     Afterwards, at the discretization of expression ratio step  150 , the expression ratio r is discretized into three states: the gene up regulated, not changed, and down regulated during the corresponding experiment. As shown in Table 1, if r is larger than a noise threshold T H , the gene was up regulated and is assigned a value of “+.” If r is larger than or equal to T L  and is smaller than or equal to T H , the gene was not changed and is assigned a value of 0. If r is smaller than T L , the gene was down regulated and is assigned a value of “−.” The noise thresholds T H  and T L  are used to filter out measurement errors associated with the experimental sample processing step  110 , the sample preparation step  120 , the gene expression measurement step  130 , and the expression ratio calculation step  140 . Even if a gene in an experiment sample has not changed in the corresponding experiment, the expression ratio r may not equal zero because Signal experiment −Background experiment  may not equal Signal control −Background control . For example, the different signal strength for the experiment sample and the control sample for the unchanged gene may result from property variations of the experiment sample and the control sample even if the two samples are of the same cell type. In addition, the microarray measurements may introduce noises related to fluorescent excitation and emission. Moreover, measurement noises for up-regulation and for down-regulation may be different, so magnitudes of T H  and T L  may differ. Furthermore, up-regulation and down-regulation of genes may be determined by other means even without computing an expression ratio. 
     
       
         
               
             
               
               
               
               
             
           
               
                 TABLE 1 
               
             
             
               
                   
               
               
                 Discretization of Expression Ratio 
               
             
          
           
               
                   
                 Expression Ratio 
                 Discretized State 
                 Discretized Value 
               
               
                   
                   
               
               
                   
                 r &gt; T H   
                 up-regulated 
                 + 
               
               
                   
                 T L  ≦ r ≦ T H   
                 not changed 
                 0 
               
               
                   
                 r &lt; T L   
                 down-regulated 
                 − 
               
               
                   
                   
               
             
          
         
       
     
     Next, the discretized expression ratios for various experiments and for various genes are organized into a summary table. As exemplified in Table 2, discretized expression ratios for ten experimental samples and three genes A, B, and C are tabulated. “X” indicates that the expression ratio is not available for the corresponding experiment and the corresponding gene. Unavailability of expression ratio may result from reasons such as lack of gene measurement on microarray or poor quality of measurement data. For example, for experimental sample #1, both genes A and B are up-regulated, but the expression ratio for gene C is not available. 
     
       
         
               
             
               
               
             
               
               
               
               
               
               
               
               
               
               
               
             
           
               
                 TABLE 2 
               
             
             
               
                   
               
               
                 Gene Expression Ratio Summary 
               
             
          
           
               
                   
                 Experiment 
               
             
          
           
               
                 Gene 
                 #1 
                 #2 
                 #3 
                 #4 
                 #5 
                 #6 
                 #7 
                 #8 
                 #9 
                 #10 
               
               
                   
               
               
                 A 
                 + 
                 0 
                 0 
                 − 
                 + 
                 0 
                 + 
                 + 
                 X 
                 X 
               
               
                 B 
                 + 
                 + 
                 0 
                 − 
                 + 
                 − 
                 + 
                 0 
                 X 
                 X 
               
               
                 C 
                 X 
                 X 
                 − 
                 + 
                 0 
                 0 
                 0 
                 + 
                 − 
                 0 
               
               
                   
               
             
          
         
       
     
     Next, at the co-expression table formation step  160 , a co-expression table for a pair of genes is obtained based on an expression ratio summary table. Any experiment that does not provide a gene expression ratio for both genes is ignored. For example, table 3, as derived from table 2, shows numbers of occurrences that genes A and B have various combinations. In three experiments, both genes A and B were both up regulated, in one experiment, gene A was up regulated while gene B was not changed, and in no experiment, gene A was up regulated and gene B was down regulated. Therefore, the total number of experiments in which gene B was up regulated is 4. 
     
       
         
               
             
               
               
             
               
               
               
               
               
               
             
           
               
                 TABLE 3 
               
             
             
               
                   
               
               
                 Co-expression Table for Genes A and B 
               
             
          
           
               
                   
                 Gene B 
               
             
          
           
               
                 Gene A  
                   
                 + 
                 0 
                 − 
                 Total 
               
               
                   
               
               
                   
                 + 
                 (++) = 3 
                 (+0) = 1 
                 (+−) = 0 
                 (+total) = 4 
               
               
                   
                 0 
                 (0+) = 1 
                 (00) = 1 
                 (0−) = 1 
                 (0 total) = 3 
               
               
                   
                 − 
                 (−+) = 0 
                 (−0) = 0 
                 (−−) = 1 
                 (−total) = 1 
               
               
                   
                 Total 
                 (total +) = 4 
                 (total 0) = 2 
                 (total −) = 2 
                 (total 
               
               
                   
                   
                   
                   
                   
                 total) = 8 
               
               
                   
               
             
          
         
       
     
     Given a co-expression table for a pair genes, there is the probability that the two genes are unrelated in their respective up-regulation, down-regulation, and no change in the same experiments. At the probability calculation step  170 , such probability P can be estimated with the Fisher Exact test, as described in Mehta, C. R., Patel, N. R. Algorithm 643, FEXACT: A FOTRAN Subroutine for Fisher&#39;s Exact Test on Unordered rxc Contingency Tables. ACM Transactions on Mathematical Software: 12(2), at 154-161 (June 1986), and in Clarkson, D. B., Fan, Y-A., and Joe, H. A Remark on Algorithm 643: FEXACT: An Algorithm for Performing Fisher&#39;s Exact Test in rxc Contingency Tables. ACM Transactions on Mathematical Software: 19(4), at 484-488 (December 1993). Alternatively, probability P may be calculated with approximated Fisher Exact test as described below. 
               P     A   ,   B       =           (     +   total     )     !     ⁢       (     0   ⁢   total     )     !     ⁢       (     -   total     )     !     ⁢       (     total   +     )     !     ⁢       (     total   ⁢           ⁢   0     )     !     ⁢       (     total   -     )     !             (   totaltotal   )     !     ⁢       (   ++   )     !     ⁢       (     +   0     )     !     ⁢       (     +   -     )     !     ⁢       (     0   +     )     !     ⁢       (   00   )     !     ⁢       (     0   -     )     !     ⁢       (     -   +     )     !     ⁢       (     -   0     )     !     ⁢       (   --   )     !                 (     Equation   ⁢           ⁢   2     )             
 
     The approximated Fisher Exact test requires less computation than the Fisher Exact test and hence improves performance of data processing. The magnitude of the probability P can be effectively represented with its base 10 logarithm p, as shown in equation 3.
 
p A,B =log P A,B   (Equation 3)
 
     Next at the gene relation prediction step  180 , the relatedness of two genes are assessed. If p A,B  is smaller than a correlation factor, it is unlikely that the up-regulation, down-regulation, and no change of the two genes are unrelated; therefore the pair of genes are significantly related. If p A,B  is larger than or equal to the correlation factor, the pair of genes are not significantly related. The magnitude of correlation factor may depend on sample sizes, number of experiments, acceptable rates of false correlations, and other factors. For example, the correlation factor may be set larger than −8 but smaller than −2. When the correlation factor equals −8, few false correlations may be incorporated but some true correlations may be missed. In contrast, when the correlation factor equals −2, some false correlations may be incorporated but few true correlations may be missed. As exemplified in tables 2 and 3, the pair of genes may be genes A and B, genes A and C, and genes B and C. 
     Next, at the gene correlation database step  190 , the predicted gene relationships are summarized and form a gene correlation database, from which all genes significantly related to a given gene may be obtained. Such database may improve ease of access to the predicted gene relations, but is nonetheless optional to practice the present invention. 
     According to another embodiment of the present invention, gene relationship may be predicted based on gene expression data measured from experiments with different control samples. As shown in  FIGS. 2-4 , p A,B  may be calculated from gene expression data measured from experiments with the same control sample. But many experiments are performed with different control samples, so p A,B  values corresponding to different control samples may be combined to yield combined p A,B(combined)  to take full advantage of available experimental data. p A,B(combined)  may be calculated in various ways, such as Fisher combined p-value. Fisher combined p-value first calculates χ 2  and degree of freedom (“DOF”), as follows: 
               χ   2     =       -   2     ⁢       ∑     i   =   1     G     ⁢           ⁢     ln   ⁢           ⁢     p     A   ,   B   ,   i                     (     Equation   ⁢           ⁢   4     )             
 DOF=2G  (Equation 5)
 
Where G is number of experimental groups and p A,B,i  is p A,B  for experimental group i. Different experimental groups have different control samples, whereas different experiments within the same experimental group have the same control sample. Next, Fisher combined p-value uses χ 2  and DOF to calculate p A,B(combined) . Alternatively, other methods may be used to calculate p A,B(combined) . p A,B(combined)  may be determined by analyzing reliability of p A,B,i  for each experimental group. For example, the smallest p A,B,i  among all experimental groups may be adopted as p A,B(combined) . After p A,B(combined)  is obtained, the relatedness of two genes are assessed. As described above, if p A,B(combined)  is smaller than a correlation factor, it is unlikely that the up-regulation, down-regulation, and no change of the two genes are unrelated; therefore the pair of genes are significantly related. If p A,B(combined)  is larger than or equal to the correlation factor, the pair of genes are not significantly related.
 
       FIG. 5  illustrates a block diagram according to yet another embodiment of the present invention. The method includes six steps: data gathering step  540 , discretization of expression ratio step  550 , co-expression table formation step  560 , probability calculation step  570 , gene relation prediction step  580 , and gene relation database formation step  590 . 
     At the data gathering step  540 , a computer  530  gathers gene expression ratio data for one or several genes from a database  510  or other sources such as a data file or spreadsheet. The data may be transferred through the Internet  520  or by other means such as DVD, CD, floppy diskette, data printed on paper, or other media. The gene expression ratios may have been obtained from one or several experiments under one or several studies. Next the expression ratios are discretized in the step  550  and tabulated in the step  560 , substantially similar to the steps  150  and  160  as discussed above. Consequently, the probability for a pair of genes to be unrelated is calculated and the relationship between the genes are assessed in the steps  570  and  580 , substantially similar to the steps  170  and  180  as discussed above. Finally, at the gene correlation database step  590 , the predicted gene relationships are summarized and form a gene correlation database, from which all genes significantly related to a given gene may be obtained. The gene correlation database may be stored on the computer  530  or another computer connected to the computer  530  over a network such as the Internet. The gene correlation database would improve ease of access to predicted gene relationships, but the step  590  is nonetheless optional for practicing the present invention. 
     To prove the principle and operation of the invention, we performed prediction of gene relationships for various genes. Prediction of relationships includes data gathering, discretization of expression ratio, co-expression table formation, probability calculation, gene relation prediction, and gene relation database formation. Embodiments of the present invention include code that directs a processor to perform the inventive steps as discussed above. An exemplary computer code has been implemented using Pentium PC and C++ computer language. This computer code is not intended to limit the scope of the claims herein. One of ordinary skill in the art would recognize other variations, modifications, and alternatives. 
     Part of the computer code is reproduced below. The computer code has three major components: 
     Declarations, initializations, and parsing command line; 
     Establishing connection with database and importing discretized expression ratios; 
     Computing p-values for gene pairs and write output file. 
     Declarations, Initializations, and Parsing Command Line 
     In this section, the computer code directs the processor to perform declarations, initializations, and parsing command line, as shown below. 
                                                                                         // Copyright 2002       // Michael van Dam, Rebecca Shafee, Jim Brody, Steve Quake.       // Department of Applied Physics       // California Institute of Technology       //       // Heavily modified from an original version by jpbrody and quake.       //       // make_pvalues:       // Compute pvalues for all pairs of UniGene clusters from their       // ‘onebit’ (discretized) expression data.       // Output file format:       //       // Line 1: make_pvalues (program name)       // Line 2: 3.0 (version number)       // Line 3: Hs (organism)       // Line 4: 150 (unigene build)       // Line 5: 2 (expt group number)       // Line 6: 28063 (number of clusters with non-null ‘onebit’)       // Line 7: 43 (number of experiments in exptgroup)       // Line 8: −5 (threshold/cutoff p-value)       //       // Then the data: pvalue clust1 clust2 N dot dist (TAB-delimited)       //       // Data only stored if pvalue &lt; cutoff.       // Version info       #define MAJOR_VERSION 3       #define MINOR_VERSION 0       // Include files       #include “common.h”       #include “local.h”       #include “mysql.h”       #include &lt;stdio.h&gt;       #include &lt;stdlib.h&gt;       #include &lt;math.h&gt;       #include &lt;string&gt;       // Compute the Fisher Exact Probability for a contingency table.       // http://www.netlib.or/toms/643       // Algorithm 643, Collected Algorithms from ACM.       // Published in Transactions on Mathematical Software 19(4): 484-488, 1993.       //       // fexact (nrow, ncol, table, ldtabl, expect, percnt, emin, prt, pre)       //       // nrow, ncol : integer       // ldtable : integer (leading dimension of ‘table’)       // table : double precision (nrow x ncol matrix)       // expect, percnt, emin, prt, pre : double precision       //       // Note, in Fortran, all arguments are ‘passed-by-reference’.       // Also, arrays stored in order such that indices should be reversed; here       // doesn&#39;t matter       extern “C” void fexact_(int &amp;nrow, int &amp;ncol, double *table, int &amp;ldtabl,       double &amp;expect, double &amp;percnt, double &amp;emin, double &amp;prt, double &amp;pre,       double *fact);       // The Fisher exact probability for a contingency table is the summation       // of pvalues for all less-likely tables having the same margin totals.       //       // A less accurate, but still usable calculation is simply to take the       // p-value of the given table (no summation). For low p-values this tends       // to be quite accurate, but for higher p-values can be quite inaccurate.       // However, if we are only interested in the small fraction of pairs       // which have low-pvalues, then this provides a MUCH faster calculation.       float p_single(int N, int np, int nz, int nm, int mp, int mz, int mm,                int kpp, int kpz, int kmz,                int kzp, int kzm, int kpm,           int kmp, int kmm, int kzz);            int status_interval = 100000;       double* f = NULL; // will store table of ln(n!)       // Error codes (see program for description)       #define ERROR_USAGE   −1       #define ERROR_DB_CONNECT  −100       #define ERROR_FOPEN  −200       // Initialize mysql variables       MYSQL* dbc = NULL;       MYSQL_RES *result = NULL;  // pointer to a result       MYSQL_ROW row;       char query [1024];       int main (int argc, char **argv) {                // -------------------           // Parse command line           // -------------------           int nargs = 4+1;           if (argc &lt; nargs) {                cout &lt;&lt; endl;           cout &lt;&lt;“Usage: make_pvalues &lt;organism&gt;&lt;exptgroup&gt;&lt;cutoff&gt;&lt;outputfile&gt;” &lt;&lt; endl;           cout &lt;&lt; endl;           cout &lt;&lt; “ &lt;organism&gt; is the two-letter organism code (e.g. Hs)” &lt;&lt; endl;           cout &lt;&lt; “ &lt;exptgroup&gt; is the expt group id” &lt;&lt; endl;           cout &lt;&lt; “ &lt;cutoff&gt; is the log of highest p-value to store” &lt;&lt; endl;           cout &lt;&lt; “ &lt;outputfile&gt; is the name of the results file” &lt;&lt; endl;           exit (ERROR_USAGE);                }           char* organism = argv[1];           int exptgroup;           sscanf (argv[2], “%d”, &amp;exptgroup);           float cutoff;           sscanf (argv[3], “%f”, &amp;cutoff);           char* outputfile = argv[4];           // Summarize run parameters           cout &lt;&lt; “make_pvalues v” &lt;&lt; MAJOR_VERSION &lt;&lt; “.” &lt;&lt; MINOR_VERSION &lt;&lt; endl;           cout &lt;&lt; “Organism: ” &lt;&lt; organism &lt;&lt; endl;           cout &lt;&lt; “Exptgroup: ” &lt;&lt; exptgroup &lt;&lt; endl;           cout &lt;&lt; “Cutoff:  ” &lt;&lt; cutoff &lt;&lt; endl;           cout &lt;&lt; “Output file: ” &lt;&lt; outputfile &lt;&lt; endl;           int build = 150; // hardwired for now           // Figure out MySQL table names           char* table_onebit = new char [1024];           sprintf (table_onebit, “%s%d_onebit”, organism, build);           char* table_exptgroups = “exptgroups”;                        
Establishing Connection with Database and Importing Discretized Expression Ratios
 
     In this section, the computer code directs the processor to establish connection with database containing discretized expression ratios and import such ratios, as shown below. 
                                                                                                                                                                                 // ------------------------------------------       // Establish connection with database (MySQL)       // ------------------------------------------       cout &lt;&lt; “Connecting to database........”;       // Attempt connection       MYSQL *dbh = mysql_init (NULL);       dbc = mysql_real_connect (dbh, DB_HOST, DB_USER, DB_PASS, PVAL_DB,           DB_PORT, DB_SOCK, DB_FLAG);       if (mysql_errno(dbc)) {                cout &lt;&lt; endl;           cout &lt;&lt; “ERROR: couldn&#39;t connect to database server” &lt;&lt; endl;           cout &lt;&lt; “MYSQL_ERROR: ” &lt;&lt; mysql_error(dbc) &lt;&lt; endl;           exit (ERROR_DB_CONNECT);            }       cout &lt;&lt; “Successful” &lt;&lt; endl;       // --------------------------------------------       // Read in the onebit strings from the database       // --------------------------------------------       // first get the number of items so we know how big to make       // the arrays       ulong n_clusters = 0;       ulong n_expts = 0;       // SOMETIMES THE VALUE CAN BE NULL!!!!!       // COUNT counts only non-NULL records... we need to make sure       // to retrieve only non-NULL records later...       sprintf(query, “SELECT COUNT(group%d) FROM %s”, exptgroup, table_onebit);       mysql_query (dbc,query);       result = mysql_store_result (dbc);       row = mysql_fetch_row (result);       sscanf (row[0], “%lu”, &amp;n_clusters);       mysql_free_result(result);       cout &lt;&lt; “Found ” &lt;&lt; n_clusters &lt;&lt; “ clusters” &lt;&lt; endl;       sprintf(query, “SELECT nexpts FROM %s WHERE id=%d”, table_exptgroups,       exptgroup);       mysql_query (dbc,query);       result = mysql_store_result (dbc);       row = mysql_fetch_row (result);       sscanf (row[0], “%lu”, &amp;n_expts);       mysql_free_result(result);       cout &lt;&lt; “Found ” &lt;&lt; n_expts &lt;&lt; “ experiments” &lt;&lt; endl;       // Initialize arrays       cout &lt;&lt; “Allocating memory for arrays...”;       ulong* clusterid = new ulong [n_clusters];       char** onebitstring = new char* [n_clusters];       for (int i=0; i &lt; n_clusters; i++) {                onebitstring[i] = new char [n_expts];           if (onebitstring[i] == NULL) {                cout &lt;&lt; “ERROR: could not allocate for onebitstring[” &lt;&lt; i  &lt;&lt;“]” &lt;&lt; endl;                }            }       cout &lt;&lt; “Done” &lt;&lt; endl;       cout &lt;&lt; “Precomputing table of ln(n!)... ”;       f  = new double [n_expts+1];       f[0] = 0.0;       f[1] = 0.0;       f[2] = log(2.0);       for (int i=3; i &lt;= n_expts; i += 2) {                f[i] = f[i−1] + log(i);           int j = i+1;           if (j &lt;= n_expts) {                f[j] = f[i] + f[2] + f[j/2] − f[j/2−1];                }            }       cout &lt;&lt; “Done” &lt;&lt; endl;       #define MINP 1000       ulong* histo = new ulong [MINP];       int i;       for (i=0; i&lt;MINP; i++) histo[i] = 0;       // Read database into arrays       cout &lt;&lt; “Reading database into arrays... ”;       sprintf (query, “SELECT cluster,group%d FROM %s WHERE group%d IS NOT NULL”,       exptgroup, table_onebit, exptgroup);       mysql_query (dbc,query);       result = mysql_store_result (dbc);       ulong temp_cluster;       for (int i=0; i &lt; n_clusters; i++) {                row = mysql_fetch_row (result);           sscanf (row[0], “%lu”, &amp;temp_cluster);           clusterid[i] = temp_cluster;           if(row[1]) {                //sscanf (row[1], “%s”, temp_string);           strncpy (onebitstring[i], row[1], n_expts);                }           else {                onebitstring[i] = NULL;                }            }       mysql_free_result(result);       cout &lt;&lt; “Done” &lt;&lt; endl;                    
Computing p-Values for Gene Pairs and Writing Output File.
 
     In this section, the computer code directs the processor to compute p-values for gene pairs and write an output file. More specifically, the processor forms a co-expression table, calculate probabilities for gene correlations, and predicts gene relations as shown below. 
     
       
         
               
             
               
               
             
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
             
               
               
             
               
               
             
               
             
           
               
                   
               
             
             
               
                 // ----------------------------------------------------- 
               
               
                 // Compute p-values for gene pairs and write output file 
               
               
                 // ----------------------------------------------------- 
               
               
                 // --------------------------- 
               
               
                 // Initialization for FEXACT 
               
               
                 // 
               
               
                 // We must pass actual variables to the FORTRAN subroutine since 
               
               
                 // everything is passed by reference, never by value. 
               
               
                 // --------------------------- 
               
               
                 int nrow = 3; 
               
               
                 int ncol = 3; 
               
               
                 int ldtabl = 3; 
               
               
                 // Cochran conditions... if satisfied do chi-squared instead of 
               
               
                 // Fisher (faster) 
               
               
                 double expect = 5.0; 
               
               
                 double percnt = 80.0; 
               
               
                 double emin = 1.0; 
               
               
                 double prt; // output 
               
               
                 double pre; // output 
               
               
                 double* table = new double [9]; 
               
               
                 float p; 
               
               
                 int dotproduct; 
               
               
                 // Create output file 
               
               
                 cout &lt;&lt; “Creating output file...”; 
               
               
                 FILE* ofile = fopen (outputfile, “w”); 
               
               
                 if (ofile == NULL) { 
               
             
          
           
               
                   
                 cout &lt;&lt; endl; 
               
               
                   
                 cout &lt;&lt; “ERROR: could not create output file”; 
               
               
                   
                 exit (ERROR_FOPEN); 
               
             
          
           
               
                 } 
               
               
                 cout &lt;&lt; “Done” &lt;&lt; endl; 
               
               
                 // Write headers to output file 
               
               
                 fprintf (ofile, “make_pvalues\n”); 
               
               
                 fprintf (ofile, “%d.%d\n”, MAJOR_VERSION, MINOR_VERSION); 
               
               
                 fprintf (ofile, “%s\n”, organism); 
               
               
                 fprintf (ofile, “%d\n”, build); 
               
               
                 fprintf (ofile, “%d\n”, exptgroup); 
               
               
                 fprintf (ofile, “%d\n”, n_clusters); 
               
               
                 fprintf (ofile, “%d\n”, n_expts); 
               
               
                 fprintf (ofile, “%f\n”, cutoff); 
               
               
                 // Main program loop 
               
               
                 ulong pairs = 0; 
               
               
                 ulong pairs_nodata = 0; 
               
               
                 cout &lt;&lt; “Writing output file...” &lt;&lt; endl; 
               
               
                 for (int i=0; i &lt; n_clusters-1; i++) { 
               
               
                  for (int j=i+1; j &lt; n_clusters; j++) { 
               
             
          
           
               
                   
                 pairs++; 
               
               
                   
                 // These need to be reset every time 
               
               
                   
                 int N=0; 
               
               
                   
                 int np=0, nz=0, nm=0; 
               
               
                   
                 int mp=0, mz=0, mm=0; 
               
               
                   
                 int kpp=0, kpz=0, kpm=0, kzp=0, kzz=0, kzm=0, kmp=0, kmz=0, kmm=0; 
               
               
                   
                 int dotproduct = 0; 
               
               
                   
                 int distsquared = 0; 
               
               
                   
                 if (onebitstring[i] &amp;&amp; onebitstring[j]) { 
               
               
                   
                   for (int x=0; x &lt; n_expts; x++) { 
               
             
          
           
               
                   
                 //char* chars = substr(onebitstring[i],x,1). substr(onebitstring[j],x,1); 
               
               
                   
                 char char1 = onebitstring[i][x]; 
               
               
                   
                 char char2 = onebitstring[j][x]; 
               
               
                   
                 if (char1 == ‘.’ || char2 == ‘.’) { 
               
               
                   
                  // do nothing 
               
               
                   
                 } 
               
               
                   
                 else if (char1 == ‘+’ &amp;&amp; char2 == ‘+’) {kpp++; dotproduct ++; } 
               
               
                   
                 else if (char1 == ‘+’ &amp;&amp; char2 == ‘0’) {kpz++; distsquared += 1; } 
               
               
                   
                 else if (char1 == ‘+’ &amp;&amp; char2 == ‘−’) {kpm++; dotproduct −−; distsquared += 4; } 
               
               
                   
                 else if (char1 == ‘0’ &amp;&amp; char2 == ‘+’) {kzp++; distsquared += 1; } 
               
               
                   
                 else if (char1 == ‘0’ &amp;&amp; char2 == ‘0’) {kzz++; } 
               
               
                   
                 else if (char1 == ‘0’ &amp;&amp; char2 == ‘−’) {kzm++; distsquared += 1; } 
               
               
                   
                 else if (char1 == ‘−’ &amp;&amp; char2 == ‘+’) {kmp++; dotproduct −−; distsquared += 4; } 
               
               
                   
                 else if (char1 == ‘−’ &amp;&amp; char2 == ‘0’) {kmz++; distsquared += 1; } 
               
               
                   
                 else if (char1 == ‘−’ &amp;&amp; char2 ==‘−’) {kmm++; dotproduct ++; } 
               
               
                   
                 else { 
               
               
                   
                  // error 
               
               
                   
                 } 
               
               
                   
                 // Count statistics 
               
               
                   
                 // The coincidences are indicated by k[mzp][mzp] 
               
               
                   
                 // (m-minus,z-zero,p-plus) and the counts by m[mp]and n[mp]. 
               
               
                   
                 //     j 
               
               
                   
                 //      
               
               
                   
                 //    kpp kpz kpm np 
               
               
                   
                 //   i  kzp kzz kzm nz 
               
               
                   
                 //    kmp kmz kmm nm 
               
               
                   
                 // 
               
               
                   
                 //    mp mz mm N 
               
               
                   
                 } // endfor: x 
               
             
          
           
               
                   
                 } // endif: onebitstrings both defined 
               
               
                   
                 else { 
               
               
                   
                  pairs_nodata ++; 
               
               
                   
                  continue; // skip iteration if no data 
               
               
                   
                 } 
               
               
                   
                 np = kpp + kpz + kpm; 
               
               
                   
                 nz = kzp + kzz + kzm; 
               
               
                   
                 nm = kmp + kmz + kmm; 
               
               
                   
                 mp = kpp + kzp + kmp; 
               
               
                   
                 mz = kpz + kzz + kmz; 
               
               
                   
                 mm = kpm + kzm + kmm; 
               
               
                   
                 N = np + nz + nm; 
               
               
                   
                 if (N == 0) { 
               
               
                   
                  // No data  - jump to next iteration 
               
               
                   
                  // fexact terminates on this condition 
               
               
                   
                  pairs_nodata ++; 
               
               
                   
                  continue; 
               
               
                   
                 } 
               
               
                   
                 else if (N == 1) { 
               
               
                   
                  // fexact crashes... therefore set p manually 
               
               
                   
                  p = 0.0; 
               
               
                   
                 } 
               
               
                   
                 else { // N &gt; 1 
               
               
                   
                 // Load values into table for fexact... 
               
               
                   
                 table[0] = kpp; 
               
               
                   
                 table[1] = kzp; 
               
               
                   
                 table[2] = kmp; 
               
               
                   
                 table[3] = kpz; 
               
               
                   
                 table[4] = kzz; 
               
               
                   
                 table[5] = kmz; 
               
               
                   
                 table[6] = kpm; 
               
               
                   
                 table[7] = kzm; 
               
               
                   
                 table[8] = kmm; 
               
               
                   
                 fexact_(nrow, ncol, table, ldtabl, expect, percnt, emin, prt, pre, f); 
               
               
                   
                 p = log10(pre); 
               
               
                   
                 // float p_old = p_single (N,np,nz,nm,mp,mz,mm, 
               
               
                   
                 //   kpp,kpz,kmz,kzp,kzm,kpm,kmp,kmm,kzz)/2.30258; 
               
               
                   
                 // NOTE: ln(10)=2.30258 
               
               
                   
                 } // endif 
               
               
                   
                 // Output to file if meets cutoff 
               
               
                   
                 if (p &lt;= cutoff) { 
               
             
          
           
               
                   
                 fprintf (ofile, “%f\t%lu\t%lu\t%d\t%d\t%d\n”, p, clusterid[i], clusterid[j], N, dotproduct, 
               
             
          
           
               
                   
                 distsquared); 
               
               
                   
                 } 
               
               
                   
                 if (!(pairs % status_interval)) { 
               
             
          
           
               
                   
                 cout &lt;&lt; “Done ”&lt;&lt; pairs &lt;&lt;“ pairs (“ &lt;&lt; pairs_nodata &lt;&lt;” no data)” &lt;&lt;endl; 
               
             
          
           
               
                   
                 } 
               
               
                   
                 // Update histogram of p-values 
               
               
                   
                 if (p &gt; 0) { 
               
             
          
           
               
                   
                 if (p &gt; 1.0e−6) printf (“p %f&gt; 0, adding to 0:−1 bin\n”, p); 
               
               
                   
                 p = 0; 
               
             
          
           
               
                   
                 } 
               
               
                   
                 if (p &lt;= −MINP) { 
               
             
          
           
               
                   
                 printf (“p %f&lt; −MINP, adding to MINP−1:MINP bin\n”, p); 
               
               
                   
                 p = −MINP+1; 
               
             
          
           
               
                   
                 } 
               
               
                   
                 histo[int(-p)]++; 
               
               
                   
                 } // endfor: i 
               
             
          
           
               
                 } // endfor: j 
               
               
                 fclose(ofile); 
               
               
                 cout &lt;&lt;“end”&lt;&lt;endl; 
               
               
                 cout &lt;&lt; endl; 
               
               
                 cout &lt;&lt; “Final total pairs:” &lt;&lt; pairs &lt;&lt; endl; 
               
               
                 cout &lt;&lt; “Final total nodata:” &lt;&lt; pairs_nodata &lt;&lt; endl; 
               
               
                 cout &lt;&lt; endl; 
               
               
                 cout &lt;&lt; “Writing histogram...” &lt;&lt; endl; 
               
               
                 cout &lt;&lt; “Bin Count” &lt;&lt; endl; 
               
               
                 for (i=0; i&lt;MINP; i++) { 
               
               
                 cout &lt;&lt; -i &lt;&lt; “:” &lt;&lt; -(i+1) &lt;&lt; “ ” &lt;&lt; histo[i] &lt;&lt; endl; 
               
               
                 } 
               
               
                 } 
               
               
                 float p_single(int N, int np, int nz, int nm, int mp, int mz, int mm, int kpp, int kpz, int kmz, int 
               
               
                 kzp, int kzm, int kpm,\ 
               
               
                    int kmp, int kmm, int kzz) { 
               
               
                 // Can probably squeeze out a little more efficiency by passing 
               
               
                 // fewer parameters (some can be deduced from the others as only 
               
               
                 // 9 are independent) 
               
               
                 return f[np] + f[nz] + f[nm] + f[mp] + f[mz] + f[mm] − f[N] 
               
             
          
           
               
                   
                 − f[kpp] − f[kpz] − f[kpm] 
               
               
                   
                 − f[kzp] − f[kzz] − f[kzm] 
               
             
          
           
               
                   
                 − f[kmp] − f[kmz] − f[kmm]; 
               
             
          
           
               
                 } 
               
               
                   
               
             
          
         
       
     
     It is understood the examples and embodiments described herein are for illustrative purposes only and that various modifications or changes in light thereof will be suggested to persons skilled in the art and are to be included within the spirit and purview of this application and scope of the appended claims.