Bayesian Joint Disease-Marker-Expression Analysis: Applied To Chronic Fatigue Syndrome

1. Bayesian Joint Disease-Marker-Expression Analysis: Applied To Chronic Fatigue Syndrome Madhuchhanda Bhattacharjee Department of Mathematics and Statistics University of Lancaster UK Jointly with Mikko J. Sillanpää Department of Mathematics and Statistics University of Helsinki Finland

2. Bayesian Joint Disease-Marker-Expression Analysis: Applied To Chronic Fatigue Syndrome Madhu Bhattacharjee

3. Phenotype-Marker-Expression Association • CFS-data makes it possible to study also joint phenotype-marker-expression association. • It also allows for screening context-specific marker and expression-effects with respect to considered subset of clinical variables. • Genetic association primarily focuses on identifying phenotype-marker association. • The considerable size of the gene-expression information available for this study necessitates subset selection. • Our approach of modeling to find either association with the genotype or analyzing the expression data differs with from the commonly known methods.

4. Phenotype data: Empiric variable • The variable “Empiric” was selected as a comprehensive summary of the disease phenotype. • Based on Reeves et al. (2005), variable empiric spans the space very similar to the first few principal components extracted from the original clinical variables. • Therefore this phenotype can be seen as a linear combination of clinical variables.

5. Phenotype data: Cluster variable • The other possible alternate was the “Cluster” variable which when compared shows a clear relationship with the Empiric-variable.

6. Phenotype data: Blood samples • Information on several aspect related to blood samples were also available. • However blood-sample information is rather too general without specific knowledge about how to use them for this particular disorder. • Unfortunately our lack of knowledge about CFS prevented us from using clinical information on blood data.

7. Phenotype data: Illness class • We selected 23 out of 84 illness class related variables based on the data availability (number of missing values) or for the ease of interpretation • The some of the variables that appeared to be "continuous" could have been very informative, unfortunately they had to be omitted due to lack of information on their actual scale or their distributional behavior. • Incidentally all of the 23 selected variables were of binary type, although that is not essential for most of the analyses presented here.

8. Marker data: Gene regions • We selected 9 of the 10 candidate genes (location information of MAOB was missing in the original data files provided). • The location information for the remaining nine gene-regions were also not clear.

9. Marker data: SNPs • In total the data includes 39 SNPs on these nine gene-regions.. • Unfortunately we were unable to locate location information on the SNPs also. • Our experience of working with this kind of data indicates that availability of accurate location information for both gene-regions and the SNPs within those could potentially increase inferential powers. • Distance/location information typically facilitates identifying and accounting for possible dependence in behavior of two closely placed “markers”.

10. Expression data: Arrays • Of the 177 arrays five were excluded due to non-availability of clinical data on these subjects. • The remaining 172 arrays included 8 replicate arrays on 8 subjects. • The 164 arrays pertaining to 164 individuals were used for further analysis after carrying out quality check of the data contained. • However it is possible to use unbalanced data in all our models.

11. Expression data: Missing entries • The intensity cut-off was set at 100 and all values below the cut-off were treated as missing. • Of the 20160 spots with data on more than 20 individuals were missing for any of the phenotype group (namely, NF, ISF, CFS and others) were eliminated. • This resulted in the selection of 9953 spots. • The data was also checked for positional information, intensity quality and annotation information.

12. Expression data: Checking • Quality of few of the arrays was found quite doubtful, however for present analysis they have not been excluded. • Position check of the spots confirmed predominantly more values missing from the top four blocks. • Also few images clearly showed poor quality around the top edges.

13. Expression data: Annotation • Spot annotations were obtained from the array manufacturer’s site. • Incidentally the ordering of annotations in the file thus obtained were slightly different compared to the order of reporting the intensities in the expression data files. • Annotations being crucial for correct interpretation of results of analysis, this mismatch unfortunately consumed some effort at the initial phase of data-exploration. • Using the manufacturer’s description of spots further annotations were obtained from public databases. • The location information in particular was not satisfactory.

14. Model Implementation • We implemented the model and performed parameter estimation using WinBUGS (Gilks et al. 1994, Spiegelhalter et al. 1999).

15. Statistical Models and Estimation Disease-Marker Association mapping • Association mapping: variable selection in the model is based on indicator variables controlling inclusion / exclusion of the genetic effects from the model. • Due to lack of accurate location information the prior for indicator variables at SNP/marker level were assumed to be shared by all SNPs in a gene-region • Indicators for all nine gene-regions were modeled with independent Bernoulli variables with a user-specified shrinkage parameter (S). • The parameter S can be interpreted as the prior probability of selecting a candidate variable (that is, the corresponding indicator is one) in the model.

16. Statistical Models and Estimation Disease-Marker Association mapping Empiric variable Overall mean + + Indicator × Coefficients Error ˜ Strata level Gene-region level Strata × SNP × Allelic level Gene-region level shrinkage parameter SNP level variation parameter Stratifier: 23 Clinical variables Ref: Sillanpää & Bhattacharjee, Genetics, 2005 Sillanpää & Bhattacharjee, 2006, submitted

17. Statistical Models and Estimation Disease-Marker Association mapping Empiric variable Overall mean + + Indicator × Coefficients Error ˜ Strata level Gene-region level Strata × SNP × Allelic level SNP level Markov dependence model using distance information Gene-region level shrinkage parameter Strata × SNP level variation parameter Different levels of shrinkage Genotype instead of Allelic form Other choices of levels • No stratification • Reverse the role of Empiric and Clinical variables

18. Statistical Models and Estimation Disease-Marker Association mapping Onset-specific effects of SNPs in disease-marker association analysis. Sex-specific effects of SNPs in disease-marker association analysis.

19. Statistical Models and Estimation Expression data analysis • Normalization and expression analysis are done in one integrated model and not in two separate steps • Normalization was carried out on all 164 arrays • The normalization was done using the block-level-piecewise-linear-regression normalization method of Bhattacharjee et al. (2002, 2004). • Empiric variables was used as co-variate information • Differential expression was searched for based on expression variation between the two main groups of individuals of interest viz. NF and CFS groups.

20. Statistical Models and Estimation Expression data analysis Microarray data analysis is typically performed in several distinct steps. • Normalisation • Classification/Clustering/gene identification • Biological interpretation For each of these steps there are very many models available and also softwares to implement them But, most of the existing techniques ignore the effects of model choice in any of the steps on the following steps, of the numerous steps that are involved in a complete microarray data analysis

21. Normalising microarray data: An example Plots of raw data and piecewise linear regression transformation Parameter estimates under piecewise linear regression normalisation of data

22. Statistical Models and Estimation Expression data analysis • By carrying out a joint normalization and expression analysis we selected 21 genes for further analysis. • The selection was based on two aspects: • the similarity of their genomic positions (screened from the databases at the band-level as explained above in the annotations) to the nine candidate genes and • the expression difference they showed with respect to the disease phenotype (between CFS and NF groups). • Few of the expression values were marked as missing because of the poor quality of the particular expressions.

23. Statistical Models and Estimation Expression data analysis Some of the selected genes

24. Statistical Models and Estimation Handling of missing values • In the association analyses models, we used missing data model 2 of Sillanpää and Bhattacharjee (2005) to handle missing values in the genotype data. • In case there were values missing in the stratifying variables the augmentation was carried out using posterior frequency distribution resulting from Uniform-Bernoulli prior assumption on the respective distribution. • For expression analysis the missing values are augmented through the integrated model for normalization and differential analysis. • In the joint analysis of marker-expression data to handle missing observations in gene expressions we assumed normal prior with pre-specified mean and variance.

25. Statistical Models and Estimation Disease-Expression Association • Based on gene-expression data analysis 21 genes/spots were selected. • These were used as candidates to the disease-expression association model. • Each expression had own regression coefficient (with pre-specified prior variance). • Each expression had own indicator variable which controls inclusion / exclusion of the particular expression from the model. • The analyses were performed with shrinkage (S=1/10).

26. Statistical Models and Estimation Disease-Expression Association Empiric variable Overall mean + + Indicator × Coefficients × Expression Error ˜ No stratification using clinical variables Gene level shrinkage parameter Gene level Gene level Genes showing high association with “Empiric” variable when analyzed using expression data

27. Statistical Models and Estimation Joint Disease-Marker-Expression Association • The 21 expressions were selected as explained before in expression analysis. • These 21 expressions were taken to the multilocus association model together with 39 SNPs to explain the disease phenotype. • Each expression and each gene region had own indicator variable which controls inclusion / exclusion of the particular expression or gene from the model. • Each expression had own regression coefficient (as above) and each SNP had two allelic effect coefficients (with gene-region specific common variance). • The analysis was performed with shrinkage (with own shrinkage parameters S(1)=1/9 for gene region and S(2)=1/10 for expressions).

28. Statistical Models and Estimation Joint Disease-Marker-Expression Association • Surprisingly SNPs continued to show no effects in this analysis also. • However gene-expressions showed association signals. • The associated genes were the same found above and two additional genes were also noticeable. • We also noted that the positions of several of the genes whose expressions showed some associations were located partly at the same gene regions than the markers showing some signals in stratified-disease-marker association analysis

29. Statistical Models and Estimation Joint Disease-Marker-Expression Association Genes showing high association with “Empiric” variable when analyzed using expression data in stratified analysis and gene-region close-by

30. Statistical Models and Estimation Joint Disease-Marker-Expression Association With stratification • This encouraged us to carryout an extension of the previous analysis using the clinical variables as stratifying factors. • We performed 23 different stratified analyses using the clinical variables. • Although implementation was simultaneous for all 23. • We allowed two overall mean parameters in the model to the both levels of the factor. • Moreover, we allowed factor-specific effects for both, SNPs (four coefficients with common variance) and gene-expressions (two coefficients independently), in the model. • These analyses were performed with shrinkage S(1) and S(2) as before.

31. Statistical Models and Estimation Joint Disease-Marker-Expression Association With stratification Indicator × Coefficients × Expression Empiric variable Overall mean Indicator × Coefficients + + Error + ˜ Strata level Gene-region level Strata × SNP × Allelic level Strata × Gene level Gene level Gene-region level shrinkage parameter SNP level variation parameter Gene level shrinkage parameter Stratification by 23 Clinical variables

32. Statistical Models and Estimation Joint Disease-Marker-Expression Association With stratification Indicator × Coefficients × Expression Empiric variable Overall mean Indicator × Coefficients + + Error + ˜ Strata level Gene-region level Strata × SNP × Allelic level Strata × Gene level Gene level SNP level Markov dependence model using distance information Gene-region level shrinkage parameter SNP level variation parameter Gene level shrinkage parameter Genotype instead of Allelic form Different levels of shrinkage Other choices of levels

33. Statistical Models and Estimation Joint Disease-Marker-Expression Association With stratification • Also here, SNPs did not show any clear association signals. • Only expressions were found to show some signs of the disease-association. • Here also the important expressions were partly overlapping in same genomic regions with the markers found in earlier stratified-disease-marker association analyses for that corresponding stratification analysis.

34. Summary of Models Clinical variables Empiric variables Other Clinical variables Marker data Expression data

35. Summary of Models Clinical variables Empiric variables Other Clinical variables Marker data Expression data

36. Summary of Models Clinical variables Empiric variables Other Clinical variables Marker data Expression data

37. Summary of Models Clinical variables Empiric variables Other Clinical variables Marker data Expression data

38. Summary of Models Clinical variables Empiric variables Other Clinical variables Marker data Expression data

39. Summary of Models Clinical variables Empiric variables Other Clinical variables Marker data Expression data

40. Discussion • To better understand why the markers did not show any effects in the joint disease-marker-expression association analysis, an additional genetical-genomics/e-QTL analysis was carried out. • This treats selected expressions as phenotypes in the phenotype-marker association • The model was similar to that in disease-marker association above and was done separately for all selected 21 expression-phenotypes. • However, we did not find any association signals (that would have indicated presence of in-cis effects) between 39 markers and 21 expressions. • It is still unclear if the modified model having SNP-specific indicators would have led to stronger conclusions.

41. Discussion • There would have been still room to do more extensive functional genomic analysis for the expression data. • Also, we have approached this problem without prior knowlegde of disease etiology. • We would like to emphasize that with input from experts of this particular disorder, the proposed models and methods of analyses could be easily modified / extended to reflect better knowledge and elicit newer dimensions of disease etiology.