Bayesian hierarchical modelling of genomic data Sylvia Richardson

Bayesian hierarchical modelling of genomic data Sylvia Richardson Centre for Biostatistics Imperial College, London In collaboration with Natalia Bochkina, Anne Mette Hein, Alex Lewin (St Mary’s) Helen Causton and Tim Aitman (Hammersmith) Peter Green (Bristol) Philippe Broët (INSERM, Paris) BBSRC Exploiting Genomics grant 1

Outline • Introduction • A fully Bayesian gene expression index (BGX) • Differential expression and array effects • Mixture models • Discussion 2

Part 1 Introduction • Recent developments in genomics have led to techniques – Capable of interrogating the genome at different levels – Aiming to capture one or several stages of the biological process DNA m. RNA protein phenotype 3

Fundamental process Protein-encoding genes are transcribed into m. RNA (messenger), and the m. RNA is translated to make proteins DNA -> m. RNA -> protein 4 Pictures from http: //www. emc. maricopa. edu/faculty/farabee/BIOBK/Bio. Book. TOC. html

What are gene expression data ? DNA Microarrays are used to measure the relative abundance of m. RNA, providing information on gene expression in a particular cell, under particular conditions The fundamental principle used to measure the expression is that of hybridisation between a sample and probes: – Known sequences of single-stranded DNA representing genes are immobilised on microarray –Tissue sample (with unknown concentration of RNA) fluorescently labelled – Sample hybridised to array – Array scanned to measure amount of RNA present for each sequence gene expression measure The expression level of ten of thousands of probes are measured on a single microarray ! gene expression profile 5

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

Example of within vs between strains gene variability • 7 cross-bred strains of mice that differ only by a small portion of chromosome 1 • Strains have different phenotypes related to immunological disorders • For each line, 9 animals used to obtain 3 pooled RNA extracts from spleen 7 x 3 samples Excellent experimental design to minimise “biological variability between replicate animals” Aim: to tease out differences between expression profiles of the 7 lines of mice and relate these to locations on chromosome 1 7

Biological variability is large ! Ratio within/total Total variance calculated over the 21 samples Average (over the 7 groups) of within strain variance 8 calculated from the 3 pooled samples

1000 genes most variable between strains: hierarchical clustering recovers the cross-bred lines structure Random set of 1000 genes 9

Common characteristics of genomics data sets • High dimensional data (ten of thousands of genes) and few samples • Many sources of variability (low signal/noise ratio) Common issues • • Pre-processing and data reduction Multiple testing Need to borrow information Importance to include prior biological knowledge 10

Part 2 • Introduction • A fully Bayesian gene expression index (BGX) – Single array model – Multiple array model • Differential expression and array effects • Mixture models • Discussion 11

A fully Bayesian Gene e. Xpression index for Affymetrix Gene. Chip arrays Anne Mette Hein SR, Helen Causton, Graeme Ambler, Peter Green PM MM Background correction Gene specific variability (probe) PM MM Raw intensities Gene index BGX 12

Affymetrix Gene. Chips: Zoom Image of Hybridised Array * * * Hybridised Spot Each gene g represented by probe set: (J: 11 -20) Expressed PM Perfect match: PMg 1, …, PMg. J Mis-match: Non-expressed PM Image of Hybridised Array Slide courtesy of Affymetrix MMg 1, …, MMg. J expression measure for gene g 13

Commonly used methods for estimation expression levels from Gene. Chips MAS 5: • uses PM and MMs. Imputes IMs from MMs to obtain all PMMMs positive • gene expression measure : estimate obtained by applying Tukey Biweight to the set of log(PM-MM) values in the probe set Characteristics: positive, robust, noisy at low levels RMA: • uses PMs only. • Fits an model with additive gene and probe effects to logscale background corrected PMs using median polish 14 Characteristics: positive, robust, attenuated signal detection

Mean (left) and Empirical standard deviation (right) over 7 conditions (arrays) for 45000 genes estimated by 2 different methods for quantifying gene expression mean Variability across conditions is conditioned by the 15 choice of summary measure ! Beware of filtering

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 16

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

Inference: • Implemented in Win. Bugs and C • allows: - Joint estimation of parameters in full Bayesian framework • obtain: - posterior distributions of parameters (and functions of these) in model: mean 2. 5 -97. 5% credibility interval For each gene g: Log-scale true signals: Log-scale non-spec. hybr: BGX: gene expr: log(Sgj+1): j=1, …, J log(Hgj+1): j=1, …, J q g: 1 2 3 4 1. 75 2 2. 25 NB! A distribution 18

Computational issues • We found mixing slow for gene specific parameters (μg , ξg 2) and large autocorrelation • For low signal (bottom 25%) more variability of Sgj and Hgj , and less separation So less information on (μg , ξg 2) and longer runs are needed • For the full hierarchical model, the convergence of the hyperparameters for the distribution of ξg 2 was problematic • We studied sensitivity to a range of plausible values for those and implemented an “empirical Bayes” version of the model which was reproducible with sensible run length 19

Posterior mean of qg using a run of 30 000 versus those obtained from runs of 5 000, 10 000 and 20 000 sweeps Reproducibility is obtained with short runs for large expression values Longer runs are necessary for low expression values 20

Single array model performance: Data set : varying concentrations (gene. Logic): • 14 samples of c. RNA from acute myeloid leukemia (AML) tumor cell line • In sample k: conc. ck(p. M): each of 11 genes spiked in at concentration ck: 1 2 0. 0 0. 5 3 4 5 6 7 8 9 10 11 12 13 14 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 21

Single array model: examples of posterior distributions of BGX expression indices Each curve (truncated normal with median param. ) represents a gene Examples with data: o: log(PMgj-MMgj) j=1, …, Jg (at 0 if not defined) Mean +- 1 SD 22

Single array model performance: Comparison with other expression measures 11 genes spiked in at 13 (increasing) concentrations BGX index qg increases with concentration …. . … except for gene 7 (spiked-in? ? ) Indication of smooth & sustained increase over a wider range of concentrations 23

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

Integrated modelling of Affymetrix data Condition 2 Condition 1 PM MM PM MM Gene specific variability (probe) Gene index BGX Hierarchical model of replicate (biological) variability and array effect Distribution of expression index for gene g , condition 2 Distribution of differential expression parameter 25

BGX Multiple array model: conditions: c=1, …, C, replicates: r = 1, …, Rc PMgjcr N( Sgjcr+ Hgjcr , τcr 2) MMgjcr N(ΦSgjcr+ Hgjcr , τcr 2) log(Sgjcr+1) TN (μgc , ξ gc 2) Gene and condition specific BGX Background noise, additive Array specific log(Hgjcr+1) TN(λcr, ηcr 2) Array-specific distribution of non -specific hybridisation qgc=median(TN(μgc, ξ gc 2)) “Pools” information over replicate probe sets j = 1, …J, r = 1, …, Rc 26

Subset of Affy. U 133 A spike-in data set (Affy. Comp) Consider: • Six arrays, 1154 genes (every 20 th and 42 spike-ins) • Same c. RNA hybridised to all arrays EXCEPT for spike-ins: `1` `2` `3` … `12` `13` Spike-in genes: `14` 1 -3 4 -6 7 -9 … 34 -36 37 -39 40 -42 0. 0 0. 25 0. 50 … 128 256 512 Condition 2 (array 4 -6): 0. 25 0. 50 1. 00 … 256 512 0. 00 2 - Spike-in conc (p. M): Condition 1 (array 1 -3): Fold change: - 2 2 … 2 27

BGX: measure of uncertainty provided Posterior mean +- 1 SD credibility intervals diffg=bgxg, 1 - bgxg, 2 } Spike in 1113 -1154 above the blue line 28 Blue stars show RMA measure

Part 3 • Introduction • A fully Bayesian gene expression index (BGX) • Differential expression and array effects – Non linear array effects – Model checking • Mixture models • Discussion 29

Differential expression and array effects Alex Lewin SR, Natalia Bochkina, Tim Aitman 30

Data Set and Biological question Biological Question Understand the mechanisms of insulin resistance Using animal models where key genes are knockout and comparison made between gene expression of wildtype (normal) and knockout mice Data set A (MAS 5) ( 12000 genes on each array) 3 wildtype mice compared with 3 mice with Cd 36 knocked out Data set B (RMA) ( 22700 genes on each array) 8 wildtype mice compared with 8 knocked out mice 31

Condition 2 Condition 1 Start with given point estimates of expression Hierarchical model of replicate Variability and array effect Posterior distribution Differential expression parameter (flat prior) Mixture modelling for classification 32

Exploratory analysis of array effect Mouse data set A Condition 1 (3 replicates) Needs ‘normalisation’ Spline curves shown Condition 2 (3 replicates) 33

Model for Differential Expression • Expression-level-dependent normalisation • Only few replicates per gene, so share information between genes to estimate variability of gene expression between the replicates • To select interesting genes: – Use posterior distribution of quantities of interest, function of, ranks …. – Use mixture prior on the differential expression parameter 34

Bayesian hierarchical model for differential expression Data: ygsr = log gene expression for gene g, replicate r g = gene effect δg = differential effect for gene g between 2 conditions r(g)s = array effect (expression-level dependent) gs 2 = gene variance • 1 st level yg 1 r N( g – ½ δg + r(g)1 , g 12), yg 2 r N( g + ½ δg + r(g)2 , g 22), Σr r(g)s = 0, r(g)s = function of g , parameters {a} and {b} • 2 nd level Priors for g , δg, coefficients {a} and {b} gs 2 lognormal (μs, τs) 35

Details of array effects (Normalization) • Piecewise polynomial with unknown break points: r(g)s = quadratic in g for ars(k-1) ≤ g ≤ ars(k) with coeff (brsk(1), brsk(2) ), k =1, … # breakpoints – Locations of break points not fixed – Must do sensitivity checks on # break points • Joint estimation of array effects and differential expression: In comparison to 2 step method – More accurate estimates of array effects – Lower percentage of false positive (simulation study) 36

Mouse Data set A For this data set, cubic fits well 3 replicate arrays (wildtype mouse data) Model: posterior means E( r(g)s | data) v. E( g | data) Data: ygsr - E( g | data) 37

Bayesian Model Checking • Check assumptions on gene variances, e. g. exchangeable variances, what distribution ? • Predict sample variance Sg 2 new (a chosen checking function) from the model specification (not using the data for this) • Compare predicted Sg 2 new with observed Sg 2 obs ‘Bayesian p-value’: Prob( Sg 2 new > Sg 2 obs ) • Distribution of p-values approx Uniform if model is ‘true’ (Marshall and Spiegelhalter, 2003) 38 • Easily implemented in MCMC algorithm

Data set A 39

Possible Statistics for Differential Expression δg ≈ log fold change δg* = δg / (σ2 g 1 / R 1 + σ2 g 2 / R 2 )½ (standardised difference) • We obtain the posterior distribution of all {δg} and/or { δ g* } • Can compute directly posterior probability of genes satisfying criterion X of interest: pg, X = Prob( g of “interest” | Criterion X, data) • Can compute the distributions of ranks 40

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

Experiment: 8 wildtype mice compared to 8 knockout mice RMA Criterion X: Gene is of interest if |log fold change| > log (1. 5) Plot of log fold change versus overall expression level Genes with pg, X > 0. 5 (green) # 292 pg, X > 0. 8 (red) # 139 The majority of the genes have very small pg, X : 97% of genes have pg, X < 0. 2 42

Posterior probabilities and log fold change Data set A : 3 replicates MAS 5 Data set B : 8 replicates RMA 43

Credibility intervals for ranks Data set B 100 genes with lowest rank (most under/ over expressed) Low rank, high uncertainty Low rank, low uncertainty 44

Using the posterior distribution of δg* (standardised difference) • Compute Probability ( | δg* | > 2 | data) Bayesian analogue of a t test ! • Order genes • Select genes such that Probability ( | δg* | > 2 | data) > cut-off ( in blue) By comparison, additional genes selected by a standard T test with p value < 5% are in red) 45

Part 4 • • Introduction A fully Bayesian gene expression index Differential expression and array effects Mixture models – Classification for differential expression – Bayesian estimate of False Discovery Rates – CGH arrays: models including information on clones spatial location on chromosome • Discussion 46

Mixture and Bayesian estimation of false discovery rates Natalia Bochkina, Philippe Broët Alex Lewin, SR 47

Multiple Testing Problem • Gene lists can be built by computing separately a criteria for each gene and ranking • Thousands of genes are considered simultaneously • How to assess the performance of such lists ? Statistical Challenge Select interesting genes without including too many false positives in a gene list A gene is a false positive if it is included in the list when it is truly unmodified under the experimental set up Want an evaluation of the expected false discovery rate (FDR) 48

Bayesian Estimate of FDR • Step 1: Choose a gene specific parameter (e. g. δg ) or a gene statistic • Step 2: Model its prior (resp marginal) distribution using a mixture model -- with one component to model the unaffected genes (null hypothesis) e. g. point mass at 0 for δg -- other components to model (flexibly) the alternative • Step 3: Calculate the posterior probability for any gene of belonging to the unmodified component : pg 0 | data • Step 4: Evaluate FDR (and FNR) for any list assuming that all the gene classification are independent (Broët et al 2004) : Bayes FDR (list) | data = 1/card(list) Σg list p 49 g 0

Mixture framework for differential expression • To obtain a gene list, a commonly used method (cf Lönnstedt & Speed 2002, Newton 2003, Smyth 2003, …) is to define a mixture prior for δg : • H 0 δg = 0 point mass at 0 with probability p 0 δg ~ flexible 2 -sided distribution to model • H 1 pattern of differential expression Classify each gene following its posterior probabilities of not being in the null: 1 - pg 0 Use Bayes rule or fix the FDR to get a cutoff 50

Mixture prior for differential expression • In full Bayesian framework, introduce latent allocation variable zg to help computations • Joint estimation of all the mixture parameters (including p 0) avoids plugging-in of values (e. g. p 0) that are influential on the classification • Sensitivity to prior settings of the alternative distribution • Performance has been tested on simulated data sets Poster by Natalia Bochkina 51

Performance of the mixture prior yg 1 r = g - ½ δg + g 1 r , r = 1, … R 1 yg 2 r = g + ½ δg + g 2 r , r = 1, … R 2 (For simplification, we assume that the data has been pre normalised) Var( gsr ) = σ2 gs ~ IG(as, bs) δg ~ p 0δ 0 + p 1 G (1. 5, 1) + p 2 G (1. 5, 2) H 0 H 1 Dirichlet distribution for (p 0, p 1, p 2) Exponential hyper prior for 1 and 2 52

Estimation • Estimation of all parameters combines information from biological replicates and between condition contrasts • s 2 gs = 1/Rs Σr (ygsr - ygs. )2 , s = 1, 2 Within condition biological variability • 1/Rs Σr ygsr = ygs. , Average expression over replicates • ½(yg 1. + yg 2. ) Average expression over conditions • ½(yg 1. - yg 2. ) Between conditions contrast 53

DAG for the mixture model zg p a 1, b 1 1 , 2 δg g g = 1: G ½(yg 1. - yg 2. ) ½(yg 1. + yg 2. ) 2 g 1 s 2 g 1 2 g 2 s 2 g 2 a 2, b 2 54

Simulated data Plot of the true differences ygr ~ N(δg , σ2 g) (8 replicates) σ2 gs ~ IG(1. 5, 0. 05) δg ~ (-1)Bern(0. 5) G(2, 2), g=1: 200 δg = 0, g=201: 1000 Choice of simulation parameters inspired by estimates found in analyses of biological data sets 55

Bayes rule FDR (black) FNR (blue) as a function of 1 - pg 0 Observed and estimated FDR/FNR correspond well Important feature Post Prob (g H 1) = 1 - pg 0 56

Comparison of mixture classification and posterior probabilities for δg* (standardised differences) Post Prob (g H 1) 10 = 6% False positive In red, 200 genes with δg ≠ 0 31 = 4% False negative Probability ( | δg* | > 2 | data) 57

Wrongly classified by mixture: truly dif. expressed, truly not dif. expressed Classification errors are on the borderline: Confusion between size of fold change and biological variability 58

Another simulation 2628 data points Many points added on borderline: classification errors in red Can we improve estimation of within condition biological variability ? 59

DAG for the mixture model zg p a 1, b 1 1 , 2 δg g ½(yg 1. - yg 2. ) ½(yg 1. + yg 2. ) 2 g 1 2 g 2 a 2, b 2 g = 1: G s 2 g 1 s 2 g 2 The variance estimates are influenced by the mixture parameters Use only partial information from the replicates to estimate 2 gs and feed forward in the mixture ? 60

Mixture, full vs partial Classification altered for 57 points: In 46 data points with improved classification when ‘feed back from mixture is cut’ In 11 data points with changed but new incorrect classification 61 Work in progress

Mixture models in CGH arrays experiments • Philippe Broët, SR • Curie Institute oncology department CGH = Competitive Genomic Hybridization between fluorescein- labelled normal and pathologic samples to an array containing clones designed to cover certain areas of the genome 62

Aim: study genomic alterations Loss Tumor supressor gene Gain Oncogene In oncology, where carcinogenesis is associated with complex chromosomic alterations, CGH array can be used for detailed analysis of genomic changes in copy number (gains or loss of genetic information) in the tumor sample. Amplification of an oncogene or deletion of a tumor suppressor gene are considered as important mechanisms for tumorigenesis 63

Specificity of CGH array experiment A priori biological knowledge from conventional CGH : • Limited number of states for a sequence : - presence, - deletion, - gain(s) corresponding to different intensity ratios on the slide Mixture model to capture the underlying discrete states • Clones located contiguously on chromosomes are likely to carry alterations of the same type Use clone spatial location in the allocation model • Some CGH custom array experiments target restricted areas of the genome Large proportion of genomic alterations are expected 64

3 component mixture model with spatial allocation ygr N(θg , g 2) , normal versus tumoral change, clone g replicate measure r θg wg 0 N(μ 0 , 02) + wg 1 N(μ 1 , 12) + wg 2 N(μ 2 , 22) presence deletion gain μ 0 : known central estimate obtained from reference clones Introduce centred spatial autoregressive Markov random fields, {ug 0}, {ug 1}, {ug 2} with nearest neighbours along the chromosomes x x x Spatial neighbours of g g -1 g g+1 Define mixture proportions to depend on the chromosomic location via a logistic model: wgk = exp(ugk) / Σm exp(ugm) favours allocation of nearby clones to same component Work in progress 65

Curie Institute CGH platform Focus on Investigating deletion areas on chromosome 1 (tumour suppressor locus) Data on 190 clones μ 0 Presence ? Ref value μ 0 = - 0. 11 Deletion ? 66

Mixture model posterior probability p of clone being deleted Classification with cut-off at p ≥ 0. 8 Short arm 67

Summary Bayesian gene expression measure (BGX) Good range of resolution, provides credibility intervals Differential Expression-level-dependent normalisation Borrow information across genes for variance estimation Gene lists based on posterior probabilities or mixture classification False Discovery Rate Mixture gives good estimate of FDR and classifies Flexibility to incorporate a priori biological features, e. g. dependence on chromosomic location Future work Mixture prior on BGX index, with uncertainty propagated to mixture parameters, comparison of marginal and prior mixture approaches, clustering of profiles for more general experimental set-ups 68

Papers and technical reports: Hein AM. , Richardson S. , Causton H. , Ambler G. and Green P. (2004) BGX: a fully Bayesian gene expression index for Affymetrix Gene. Chip data (to appear in Biostatistics) Lewin A. , Richardson S. , Marshall C. , Glazier A. and Aitman T. (2003) Bayesian Modelling of Differential Gene Expression (under revision for Biometrics) Broët P. , Lewin A. , Richardson S. , Dalmasso C. and Magdelenat H. (2004) A mixture model-based strategy for selecting sets of genes in multiclass response microarray experiments. Bioinformatics 22, 2562 -2571. Broët, P. , Richardson, S. and Radvanyi, F. (2002) Bayesian Hierarchical Model for Identifying Changes in Gene Expression from Microarray Experiments. Journal of Computational Biology 9, 671 -683. Available at http : //www. bgx. org. uk/ Thanks 69