Structured statistical modelling of gene expression data

1. Windsor, October 2004 Structured statistical modelling of gene expression data Peter Green (Bristol) Sylvia Richardson, Alex Lewin, Anne-Mette Hein (Imperial) with Clare Marshall, Natalia Bochkina (Imperial) Graeme Ambler (Bristol) Tim Aitman and Helen Causton (Hammersmith) BGX

2. Statistical modelling and biology • Extracting the message from microarray data needs statistical as well as biological understanding • Statistical modelling – in contrast to data analysis – gives a framework for formally organising assumptions about signal and noise • Our models are structured, reflecting data generation process: ‘highly structured stochastic systems’

3. Background and 3 studies • Hierarchical modelling • A fully Bayesian gene expression index (BGX) • Differential expression and array effects • Two-way clustering

4. Part 1 • Hierarchical modelling • A fully Bayesian gene expression index (BGX) • Differential expression and array effects • Two-way clustering

5. Gene expression using Affymetrix chips * * * * * Zoom Image of Hybridised Array Hybridised Spot Single stranded, labeled RNA sample Oligonucleotide element 20µm Millions of copies of a specific oligonucleotide sequence element Expressed genes Approx. ½ million different complementary oligonucleotides Non-expressedgenes Slide courtesy of Affymetrix 1.28cm Image of Hybridised Array

6. condition/treatment biological array manufacture imaging technical within/between array variation gene-specific variability Variation and uncertainty Gene expression data (e.g. Affymetrix) is the result of multiple sources of variability Structured statistical modelling allows considering all uncertainty at once

7. Advantages of avoiding plug-in approach Uncertainties propagated throughout model Realistic estimates of variability Avoid bias The price you pay – computational costs Intricate implementation Longer run times (but far less than experimental protocol!) Costs and benefits of this approach

8. Part 2 • Hierarchical modelling • A fully Bayesian gene expression index (BGX) • Differential expression and array effects • Two-way clustering

9. A fully Bayesian Gene eXpression indexfor Affymetrix GeneChip arraysAnne-Mette HeinSylvia Richardson, Helen Causton, Graeme Ambler, Peter Green Gene specific variability (probe) PM MM PM MM PM MM PM MM BGX Gene index

10. Single array model: motivation Key observations: Conclusions: PMs and MMs both increase with spike-in concentration (MMs slower than PMs) MMs bind fraction of signal Multiplicative (and additive) error; transformation needed Spread of PMs increase with level Considerable variability in PM (and MM) response within a probe set Varying reliability in gene expression estimation for different genes Estimate gene expression measure from PMs and MMs on log scale Probe effects approximately additive on log-scale

11. Model assumptions and key biological parameters • The intensity for the PM measurement for probe (reporter) j and gene g is due to binding • of labelled fragments that perfectly match the oligos in the spot (the true signal Sgj) • of labelled fragments that do not perfectly match these oligos (the non-specific hybridisation Hgj) • The intensity of the corresponding MM measurement is caused • by a binding fraction Φ of the true signal Sgj • by non-specific hybridisation Hgj

12. PMgj N( Sgj + Hgj , τ2) MMgj N(Φ Sgj + Hgj,τ2) Background noise, additive fraction Non-specific hybridisation signal log(Sgj+1)  TN (μg ,ξg2) log(Hgj+1)  TN(λ, η2) j=1,…,J Array-wide distribution Gene specific error terms: exchangeable Gene expression index (BGX): qg=median(TN (μg ,ξg2)) “Pools” information over probes j=1,…,J “Empirical Bayes” log(ξg2)N(a, b2) Priors: “vague”t2 ~ G(10-3, 10-3) F~ B(1,1), mg ~ U(0,15) h2 ~ G(10-3, 10-3),l ~ N(0,103) BGX single array model g=1,…,G (thousands), j=1,…,J (11-20)

13. Markov chain Monte Carlo (MCMC) computation • Fitting of Bayesian models hugely facilitated by advent of these simulation methods • Produce a large sample of values of all unknowns,  from posterior given data • Easy to set up for hierarchical models • BUT can be slow to run (for many variables!) • and can fail to converge reliably

14. Sample in place of a distribution - 1D

15. Sample in place of a distribution - 2D

16. Single array model performance • Data set : varying concentrations (geneLogic): • 14 samples of cRNA from acute myeloid leukemia (AML) tumor cell line • In sample k: each of 11 genes spiked in at concentration ck: • sample k: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 • conc. (pM): 0.0 0.5 0.75 1.0 1.5 2.0 3.0 5.0 12.5 25 50 75 100 150 • Each sample hybridised to an array • Consider subset consisting of 500 normal genes + 11 spike-ins

17. Signal & expression indices 10 arrays: gene 1 spiked-in at increasing concentrations Lines: 95% credibility intervals for log(Sgj+1) Curves: posterior for signal `true signal`/ expression index BGX increases with concentration

18. Non-specific hybridisation 10 arrays: gene 1 spiked-in at increasing concentrations Lines: 95% credibility intervals for log(Hgj+1) Curves: posterior for signal Non-specific hybridisation does not increase with concentration

19. Comparison with other expression measures 11 genes spiked in at 13 (increasing) concentrations BGX index qgincreases with concentration ….. … except for gene 7 (incorrectly spiked-in??) Indication of smooth & sustained increase over a wider range of concentrations

20. 95% credibility intervals for Bayesian gene expression index 11 spike-in genes at 13 different concentrations (data set A) Each colour corresponds to a different spike-in gene Gene 7 : broken red line Note how the variability is substantially larger for low expression level

21. Part 3 • Hierarchical modelling • A fully Bayesian gene expression index (BGX) • Differential expression and array effects • Two-way clustering

22. Bayesian modelling of differential gene expression, adjusting for array effectsAlex LewinSylvia Richardson, Natalia Bochkina,Clare Marshall, Anne Glazier, Tim Aitman • The spontaneously hypertensive rat (SHR): A model of human insulin resistance syndromes. • Deficiency in gene Cd36 found to be associated with insulin resistance in SHR • Following this, several animal models were developed where other relevant genes are knocked out comparison between knocked out and wildtype (normal) mice or rats. See poster!

23. Data set & biological question Microarray Data Data set A (MAS 5) ( 12000 genes on each array) 3 SHR compared with 3 transgenic rats Data set B (RMA) ( 22700 genes on each array) 8 wildtype (normal) mice compared with 8 knocked out mice Biological Question Find genes which are expressed differently in wildtype and knockout / transgenic mice

24. Condition 1 (3 replicates) Condition 2 (3 replicates) Exploratory analysis showing array effect

25. Differential expression model The quantity of interest is the difference between conditions for each gene: dg , g = 1, …,N Joint model for the 2 conditions : yg1r = g - ½ dg + 1r(g)+ g1r , r = 1, … R1 yg2r = g + ½ dg + 2r(g) + g2r , r = 1, … R2 where ygcr is log gene expression for gene g, condition c, replicate r g is overall gene effect cr() is array effect - a smooth function of  gcr is normally distributed error, with gene- and condition- specific variance

26. Differential expression model Joint modelling of array effects and differential expression: • Performs normalisation simultaneously with estimation • Gives fewer false positives Can work with any desired composite criterion for identifying ‘interesting’ genes, e.g. fold change and overall expression level

27. Data set A 3 wildtype mice compared to 3 knockout mice (U74A chip) MAS5 Criterion: Gene is of interest if |log fold change| > log(2) and log (overall expression) > 4 Plot of log fold change versus overall expression level Genes with pg,X > 0.5 (green) # 280 pg,X > 0.8 (red) # 46 The majority of the genes have very small pg,X : 90% of genes have pg,X < 0.2 pg,X = 0.49 Genes with low overall expression have a greater range of fold change than those with higher expression

28. Part 4 • Hierarchical modelling • A fully Bayesian gene expression index (BGX) • Differential expression and array effects • Two-way (gene by sample) clustering

29. Hierarchical clustering of samples The gene expression profiles cluster according to tissue of origin of the samples A subset of 1161 gene expression profiles, obtained in 60 different samples Red: more mRNA Green : less mRNA in the sample compared to a reference Ross et al, Nature Genetics, 2000

30. Non-model-based clustering • Many clustering algorithms have been developed and used for exploratory purposes • They rely on a measure of ‘distance’ (dissimilarity) between gene or sample profiles, e.g. Euclidean • Hierarchical clustering proceeds in an agglomerative manner: single profiles are joined to form groups using the distance metric, recursively • Good visual tool, but many arbitrary choices care in interpretation!

31. Model-based clustering • Build the cluster structure into the model, rather than estimating gene effects (say) first, and post-processing to seek clusters • Bayesian setting allows use of real prior information where it is exists (biological understanding of pathways, etc, previous experiments, …)

32. Additive ‘ANOVA’ models for (log-) gene expression The simplest model:gene + sample g=gene s=sample/condition Under standard conditions, the(least-squares) estimates of gene effectsare Themodelgenerates themethod, and in this case performs a simple form ofnormalisation

33. ... bring in mixture modelling … (single sample first!) g=gene Tg= unknown cluster to which gene g belongs This is a mixture model

34. … finally allow clusters to overlap – ‘Plaid’ model h denotes a ‘cluster’, ‘block’ or ‘layer’ – pathway? gh= 0 or 1 and sh= 0 or 1

35. ‘Plaid’ model samples genes

36. An early experiment : artificial raw data Artificial data from a very special case of the Plaid model: single sample s True H=3, b(h)=2.2, 3.4 and 4.7, N(0,2); 500 genes, some in each of 23=8 configurations ofgh  8 overlapping normal clusters

37. trueHwas 3 trueb(h)were 2.2, 3.4, 4.7

38. Human fibroblast data – Lemon et al (2002) • 18 samples split into 3 categories: serum starved, serum stimulated and a 50:50 mix of starved/stimulated. • We used the natural logarithm of Lemon et al.’s calculated LWF values as our measure of expression and subtracted gene and sample mean levels. • We then selected the 100 most variable genes across all 18 samples and used this 18×100 array as the input to our analysis.

39. Bayesian clustering • Hierarchical model allows us to learn about all unknowns simultaneously • In particular, this includes complete 2-way classification, gene by sample, with numerical uncertainties • We then construct visualisations of interesting aspects (marginal distributions) of this posterior

40. Bayesian clustering: samples

41. Bayesian clustering: genes

42. More details, papers and code • www.stats.bris.ac.uk/BGX/ • www.bgx.org.uk