STATISTICAL TOOLS FOR SYNTHESIZING LISTS OF DIFFERENTIALLY EXPRESSED

STATISTICAL TOOLS FOR SYNTHESIZING LISTS OF DIFFERENTIALLY EXPRESSED FEATURES IN MICROARRAY EXPERIMENTS Marta Blangiardo and Sylvia Richardson 1 1 Centre for Biostatistics, Imperial College, St Mary’s Campus, Norfolk Place London W 2 1 PG, UK. m. blangiardo@imperial. ac. uk SCOPE OF THE WORK Consider two different but related experiments, how to assess whethere are more differentially expressed genes in common than expected by chance? R(q*)=1. 0 RANKED LISTS q*=0. 05 Suppose we have two experiments, each reporting a measure (e. g. p-value, …) of differential expression on a probability scale: O 11(q*)=18 Small p value: MOST differentially expressed O 1+(q) Large p value: NOT differentially expressed Experiment A Experiment B p. A 1 p. A 2 … p. An p. B 1 p. B 2 … p. Bn All the computations have been performed in R and are available on BGX website (www. bgx. org. uk) O+1(q) SIMULATION We rank the genes according to the probability measures. For each cut off q we obtain a 2 X 2 table: We use three batches of simulations differing by level of association between experiments and percentage of DE genes. For every batch we simulate two lists of 2000 p-values (Allison et al. 2002) averaging the results over 100 simulations. Conditional Model Permutation Test Exp B Exp A DE DE O 11(q) DE O+1(q)-O 11(q) DE O 1+(q)-O 11(q) n-O 1+(q)- O+1(q)+O 11(q) O 1+(q) n- O+1(q) n O+1(q) The number of genes in common by chance is The number of genes observed in common is O 11(q) RATIO T(q*) r= 0 , 1. 1 DE = 10% q* Joint Model Bayesian Analysis O 11(q*) O 1+(q*) O+1(q*) MC p-value R(q*) 95% CI q* O 11(q*) O 1+(q*) O+1(q*) 0. 040 10 115 120 0. 550 1. 0 0. 050 [0. 4 -1. 5] 18 125 130 r= 0. 25 , 5. 7 0. 010 DE = 10% r= 0. 25 , 3. 0 0. 019 DE = 20% r= 0. 25 , 2. 5 0. 023 DE = 30% 6 49 50 0. 060 5. 0 0. 020 [2. 2 -10. 6] 8 59 59 11 82 82 0. 030 2. 9 0. 026 [1. 5 -4. 9] 17 105 106 21 125 126 0. 002 2. 4 0. 030 [1. 4 -3. 6] 28 148 150 NO association is declared when the two lists are not associated (MC p-value not significant, CI include 1) We propose to calculate the maximum of the observed to expected ratio: When there is a TRUE association: Conditional Model Joint Model • The ratio T(q*) decrease • q*, O 1+(q), O+1(q), and O 11(q) increase It is the maximal deviation from the underneath independence model. Increasing % of DE genes • MC p-value is more significant 1 • By using the maximum ratio, multiple testing issues for different list sizes are avoided • Returns a single list of O 11(q) genes for further biological investigation • The ratio R(q*) decreases • q*, O 1+(q), O+1(q), and O 11(q) increase • CI 95 are narrower R(q*) is always smaller than T(q*) and its q* is slightly bigger as it accounts for the additional variability APPLICATION: analysis of deleterious effect of mechanical ventilation on lung gene expression We re-analyse the experiment presented in Ma et al, 2005, investigating the deleterious effect of mechanical ventilation on lung gene expression through a model of mechanical ventilation-induced lung injury (VILI) on rodents (mice and rats). §We analyse separately the two dataset using Cyber-T (Baldi and Long, 2001) 0 0. 05 0. 1 § We use RESOURCERER to reconstruct the list of orthologs for the two species §We apply the methodology described to the lists of 2969 p-values (ortholog genes) PERMUTATION TEST Given a threshold q and fixed margins • 97 genes found in common between mice and rats But the distribution of T(q*) is not easily obtained since the tables are nested in each other. We take advantage of the empirical distribution for T(q*) obtained via permutations. • 15 genes in common with the original analysis (which highlighted 48 genes) We perform a Monte Carlo test of T under the null hypothesis of independence between the two experiments using permutations. This returns a Monte Carlo p-value. • Two enriched pathways with our methodology: Not associated Associated P-value 0. 8 P-value <0. 001 1) MAPK signalling activity. 6 out of the significant orthologs are involved in this KEGG pathway (Fgfr 1, Gadd 45 a, Hspa 8, Hspa 1 a, Il 1 b, Il 1 r 2) while only 4 were highlighted in the original one. T(q*)=1. 44 R(q*) = 1. 43 • The uncertainty of the margins is not taken into account q*=0. 01 q* = • The size of the list of genes in common can be vary small (typically when the total number of DE genes is small) and this can cause an instability in the estimate of T(q) MC p-value <0. 001 CI 95 = [1. 13 -1. 75] LIMITATIONS OF THE TEST O 11(q*) = 97 • We propose a Bayesian model treating also the margins as random variables 0. 01 2) Cytokine-Cytokine receptor interaction. 5 out of the significant orthologs are involved in this KEGG pathway (IL 6, Il 1 b, Il 1 r 2, CCL 2, Kit) while only 4 were highlighted in the original one. O 1+(q*) = 393 BAYESIAN MODEL Starting with the 2 x 2 table we specify a multinomial distribution of dimension 3 for the vector of joint frequencies: O+1(q*) = 886 DISCUSSION • This is a simple procedure to evaluate if two (or more) experiments are associated • The permutation test gives a first look under the model where the marginal frequencies are fixed and the vector of parameters q is modelled as non informative Dirichlet: q ~ Di(0. 05, 0. 05) The derived quantity of interest is the ratio of the probability that a gene is in common to the probability that a gene is in common by chance: • The Bayesian model permits to enlarge the scenario introducing variability on all the components • It is very flexible and adaptable for comparisons of several experiments at different levels (gene level, biological processes level) and for different problems (e. g. comparison between species , comparison between platforms ) REFERENCES Since the model is conjugated the posterior distribution for q is Dirichlet DECISION RULE We can obtain a sample from the posterior distribution of the derived quantities R(q) and calculate the credibility interval (CI) at 95% for each threshold q. We define q* as the value of the argument for which the median of R(q) attains its maximum value, only for the subset of credibility intervals which do not include 1: Then R(q*) is the ratio associated to q*. Allison et al. (2002), “A mixture model approach for the analysis of microarray gene expression data”, Computational Statistics And Data Analysis, 39, 1 -20. Baldi and Long, (2001) “A Bayesian framework for the analysis of microarray expression data: regularized t-test and statistical inferences of gene changes”, Bioinformatics, 17, 509 -519. Ma et al. , (2005) “Bioinformatics identification of novel early stress response genes in rodent models of lung injury”, Am J Physiol Lung Cell Mol Physiol 289(3), 468 -477. ACKNOWLEDGEMENTS We would like to thank Natalia Bochkina, Alex Lewin and Anne-Mette Hein for helpful discussions. This work has been supported by a Wellcome Trust Functional Genomics Development Initiative (FGDI) thematic award ``Biological Atlas of Insulin Resistance (BAIR)", PC 2910_DHCT