Computational and statistical approaches for genotype imputation, haplotype reconstruction and analysis of genome variation

PhD Thesis Computational and statistical approaches for genotype imputation, haplotype reconstruction and analysis of genome variation PhD Thesis Computational and statistical approaches for genotype imputation, haplotype reconstruction and analysis of genome variation Nathan Medina Rodríguez Las Palmas de Gran Canaria November, 2015 ULPGC Anexo I D. PEDRO PÉREZ CARBALLO SECRETARIO DEL INSTITUTO UNIVERSITARIO DE MICROELECTRÓNICA APLICADA DE LA UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA, CERTIFICA, Que el Consejo de Doctores del Instituto en su sesión de fecha 23 de Noviembre de 2015 tomó el acuerdo de dar el consentimiento para su tramitación, a la tesis doctoral titulada “Computational and statistical approaches for genotype imputation, haplotype reconstruction and analysis of genome variation” presentada por el doctorando D. Nathan Medina Rodríguez y dirigida por los doctores Ángelo Santana del Pino, Ana María Wägner Fahlin y José María Quinteiro González. Y para que así conste, y a efectos de lo previsto en el Artº 6 del Reglamento para la elaboración, defensa, tribunal y evaluación de tesis doctorales de la Universidad de Las Palmas de Gran Canaria, firmo la presente en Las Palmas de Gran Canaria, a 23 de Noviembre de 2015. Anexo II Departamento/Instituto/Facultad: Instituto Universitario de Microelectrónica Aplicada Programa de doctorado: Tecnologías de Telecomunicación Título de la Tesis: Computational and statistical approaches for genotype imputation, haplotype reconstruction and analysis of genome variation Tesis Doctoral presentada por: Nathan Medina Rodríguez Dirigida por el Dr.: Ángelo Santana del Pino Dirigida por la Dra.: Ana María Wägner Fahlin Dirigida por el Dr.: José María Quinteiro González Director Directora Director Doctorando Las Palmas de Gran Canaria, a 19 de Noviembre de 2015 EN CONFORMIDAD CON LOS REQUERIMIENTOS SOLICITADOS PARA LA OBTENCIÓN DEL GRADO DE DOCTOR POR LA UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA IUMA – TECNOLOGÍAS DE LA INFORMACIÓN DEPT. DE MATEMÁTICAS – GRUPO DE ESTADÍSTICA Computational and statistical approaches for genotype imputation, haplotype reconstruction and analysis of genome variation Autor: Nathan Medina Rodríguez Director: Dr. Ángelo Santana Directora: Dra. Ana Mª Wägner Director: Dr. José Mª Quinteiro Las Palmas de Gran Canaria, A 19 de Noviembre de 2015 A THESIS IN CONFORMITY WITH THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA IUMA – INFORMATION AND COMMUNICATION SYSTEMS DEPARTMENT OF MATHEMATICS – GROUP OF STATISTICS Computational and statistical approaches for genotype imputation, haplotype reconstruction and analysis of genome variation Autor: Nathan Medina-Rodriguez Supervisor: Angelo Santana, PhD Supervisor: Ana M. Wägner, MD, PhD Supervisor: Jose M. Quinteiro, PhD Las Palmas de Gran Canaria, November 19th, 2015 i To my wife and family. ii Acknowledgments I would like to thank my PhD supervisors, who with their support, trust, help, and patience have facilitated the realization of this work. First of all, I want to thank my mentor Angelo Santana, from whom I learned so much. He addressed my thesis tirelessly and unbreakable, guiding me in each of my steps. Thanks for your time, for your support, for your sincere friendship, for your knowledge; in short, thank you for everything. To Ana M Wagner, for her help and for always being present when I needed her, advising me and facilitating me all necessary to conduct my research tasks. To Jose M Quinteiro, for allowing me to join his group as a doctoral student and for his timely comments on the design of our research. In addition, I would also like to thank the Department of Mathematics, especially to the Group of Statistics at the Universidad de Las Palmas, and to the Institute of Applied Microelectronics for facilitating me the development of my thesis and for helping with all administrative tasks. I also extend special thanks to Dr. Torben Hansen of The Novo Nordisk Foundation Center for Basic Metabolic Research at the University of Copenhagen, for making me feel at home. Thanks to my family for supporting me in every one of my journeys. Thanks to Elizabeth for being my most faithful companion. Finally, my sincere gratitude to all the studies’ participants, since without their contribution and willingness our research would not have been possible. Part of this work has been supported by the European Foundation for the Study of Diabetes with an Albert Renold Travel Fellowship grant. Nathan Medina-Rodriguez Las Palmas de Gran Canaria November 2015 Contents i Introduction 17 1 Introduction 19 1.1. Background....................... 19 1.2. Motivation and Original Contributions . . . . . . . . 20 1.3. Document Structure . . . . . . . . . . . . . . . . . . 21 ii Basic Concepts and State of the Art 23 2 Concepts of Human Genetics 25 2.1. Molecular Genetics Terminology . . . . . . . . . . . . 25 2.2. Transmission of Genetic Information . . . . . . . . . . 31 2.3. Human Genetic Diversity . . . . . . . . . . . . . . . . 33 3 Concepts of Statistical Genetics 39 3.1. Statistical Genetics . . . . . . . . . . . . . . . . . . . 39 3.2. Population Studies . . . . . . . . . . . . . . . . . . . 40 3.3. Family-based Studies . . . . . . . . . . . . . . . . . . 50 4 Concepts of Computational Genomics 55 4.1. Computational Genomics . . . . . . . . . . . . . . . . 55 4.2. Computational Biology vs. Bioinformatics . . . . . . . 55 4.3. Machine Learning in Biology . . . . . . . . . . . . . . 56 4.4. Applied Computational Techniques . . . . . . . . . . 58 4.5. International Genetic Databases . . . . . . . . . . . . 63 5 State of the Art 69 5.1. Introduction....................... 69 5.2. Genotype Imputation . . . . . . . . . . . . . . . . . . 70 1 2CONTENTS 5.3. Haplotype Reconstruction . . . . . . . . . . . . . . . 76 5.4. Analysis of Genome Variation . . . . . . . . . . . . . 83 iii Approaches for Population Data 87 6 Quality Control 89 6.1. Introduction....................... 89 6.2. Marker Quality Measures . . . . . . . . . . . . . . . . 92 6.3. Sample Quality Measures . . . . . . . . . . . . . . . . 95 7 Alignment and Phasing 103 7.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . 103 7.2. Alignment........................107 7.3. Phasing .........................109 8 Imputation 113 8.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . 113 8.2. Imputation .......................115 9 GWA Testing of Imputed Data 121 9.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . 121 9.2. Association Method Description . . . . . . . . . . . . 123 9.3. Association testing of imputed data . . . . . . . . . . . 126 9.4. Results using HapMap as Reference Panel . . . . . . . 129 9.5. Results using 1000 Genomes as Reference Panel . . . . 130 9.6. Selection of Significant Results . . . . . . . . . . . . . 133 iv Approaches for Family-based Data 141 10 alleHap Package: Description 143 10.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . 143 10.2. Theoretical Description . . . . . . . . . . . . . . . . . 143 10.3. Practical Description . . . . . . . . . . . . . . . . . . 147 11 alleHap Package: Performance 171 11.1. Computing Times . . . . . . . . . . . . . . . . . . . . 171 11.2. Genotype Imputation Rates . . . . . . . . . . . . . . 175 11.3. Reconstructed Haplotypes . . . . . . . . . . . . . . . 182 CONTENTS 3 12 alleHap Package: Applications 191 12.1. alleHap into T1DGC database . . . . . . . . . . . . . 191 12.2. Comparison of the distribution of risk haplotypes between the Canary Islands and the rest of Spain . . . . . 201 12.3. Comparison of the distribution of risk haplotypes between Spain and the rest of Europe . . . . . . . . . . . 202 v Conclusion 205 13 Main conclusions 207 vi Appendix 211 A GWA Tutorial 213 A.1. Quality Control . . . . . . . . . . . . . . . . . . . . . 214 A.2. Pre-processing . . . . . . . . . . . . . . . . . . . . . . 218 A.3.Imputation .......................220 A.4. GWA Analysis . . . . . . . . . . . . . . . . . . . . . 222 A.5. Data Representation . . . . . . . . . . . . . . . . . . 224 B alleHap Manual 227 B.1. Input Format . . . . . . . . . . . . . . . . . . . . . . 227 B.2. Data Simulation . . . . . . . . . . . . . . . . . . . . . 228 B.3.Workflow ........................231 vii Resumen en Español 243 14 Introducción 245 15 Objetivos 249 16 Planteamiento y Metodología 251 16.1. Planteamiento . . . . . . . . . . . . . . . . . . . . . . 251 16.2. Metodología . . . . . . . . . . . . . . . . . . . . . . . 251 17 Resultados 275 17.1. Resultados del análisis de bases de datos poblacionales . 275 17.2. Resultados del análisis de bases de datos familiares . . . 287 4CONTENTS 18 Conclusiones 301 Bibliografía 305 List of Figures 2.1. From cell to gene . . . . . . . . . . . . . . . . . . . . . . 25 2.2. Karyogram of human chromosomes . . . . . . . . . . . . 27 2.3. DNA structural diagram . . . . . . . . . . . . . . . . . . 28 2.4. Alleles in chromosomes . . . . . . . . . . . . . . . . . . . 30 2.5. Meiosis two stage scheme . . . . . . . . . . . . . . . . . . 32 2.6. Meiotic recombination . . . . . . . . . . . . . . . . . . . 33 2.7. DNA Structural Variation . . . . . . . . . . . . . . . . . 34 2.8. SNPdiagram ........................ 35 2.9. MHC-HLA Complex . . . . . . . . . . . . . . . . . . . 36 2.10. HLA Alleles Nomenclature . . . . . . . . . . . . . . . . . 38 3.1. Hardy-Weinberg proportions for two alleles . . . . . . . . 42 3.2. Linkage and Linkage Disequilibrium . . . . . . . . . . . . 47 3.3. Indirect Association . . . . . . . . . . . . . . . . . . . . . 50 4.1. Machine Learning Topics . . . . . . . . . . . . . . . . . . 57 4.2. HMM Architecture . . . . . . . . . . . . . . . . . . . . . 59 4.3. HMM for Haplotypic Data . . . . . . . . . . . . . . . . . 60 4.4. HMMexample....................... 61 4.5. Overview of the EM algorithm . . . . . . . . . . . . . . . 62 4.6. SNPs, haplotypes and tag SNPs. . . . . . . . . . . . . . . 64 4.7. Pedigree structures into the T1DGC . . . . . . . . . . . . 67 5.1. Genotype imputation overview . . . . . . . . . . . . . . . 70 5.2. Association example using imputed data . . . . . . . . . . 72 5.3. IMPUTE2 standard imputation scenario . . . . . . . . . . 74 5.4. Imputation scenarios . . . . . . . . . . . . . . . . . . . . 78 5.5. Phasing Methods Comparison . . . . . . . . . . . . . . . 79 5.6. Haplotype correction example using DuoHMM . . . . . . 80 6 List of Figures 7 5.7. GWA published reports . . . . . . . . . . . . . . . . . . . 83 5.8. Published Genome-Wide Associations . . . . . . . . . . . 84 6.1. Flowchart of the Quality Control process . . . . . . . . . 90 6.2. Genotyping Efficiency . . . . . . . . . . . . . . . . . . . 93 6.3. QQ plot of HW control p-values . . . . . . . . . . . . . . 95 6.4. Sample missingness . . . . . . . . . . . . . . . . . . . . . 97 6.5. Example of relatedness networks . . . . . . . . . . . . . . 99 6.6. Relatedness among samples . . . . . . . . . . . . . . . . . 99 6.7. IBD Histograms . . . . . . . . . . . . . . . . . . . . . . 99 6.8. Heterozygosity Histograms before SQC . . . . . . . . . . 100 6.9. Heterozygosity Histograms after SQC . . . . . . . . . . . 101 7.1. SHAPEIT method example . . . . . . . . . . . . . . . . 111 8.1. Standard Imputation Scenario . . . . . . . . . . . . . . . 119 9.1. Manhattan plot using HapMap . . . . . . . . . . . . . . . 129 9.2. QQ plot using HapMap . . . . . . . . . . . . . . . . . . 129 9.3. Manhattan plot using 1000G . . . . . . . . . . . . . . . . 131 9.4. QQ plot using 1000G . . . . . . . . . . . . . . . . . . . . 131 9.5. Manhattan plot using 1000G . . . . . . . . . . . . . . . . 132 9.6. QQ plot using 1000G . . . . . . . . . . . . . . . . . . . . 132 9.7. First region of Typed SNPs . . . . . . . . . . . . . . . . . 136 9.8. First region of Typed/Imputed SNPs . . . . . . . . . . . . 136 9.9. Second region of Typed SNPs . . . . . . . . . . . . . . . 137 9.10. Second region of Typed/Imputed SNPs . . . . . . . . . . 137 9.11. Third region of Typed SNPs . . . . . . . . . . . . . . . . 138 9.12. Third region of Typed/Imputed SNPs . . . . . . . . . . . 138 10.1. Alleles, haplotypes and pedigree scheme . . . . . . . . . . 144 10.2. Package Description . . . . . . . . . . . . . . . . . . . . . 148 11.1. Computing Times per Families . . . . . . . . . . . . . . . 172 11.2. Computing Times per Markers . . . . . . . . . . . . . . . 173 11.3. Computing Times per Alleles . . . . . . . . . . . . . . . . 174 11.4. Initial vs. Final Imputation rates . . . . . . . . . . . . . . 176 11.5. Imputation rates vs. missing genotypes -1/4- . . . . . . . 176 11.6. Imputation rates vs. missing genotypes -2/4- . . . . . . . 177 11.7. Imputation rates vs. missing genotypes -3/4- . . . . . . . 178 11.8. Imputation rates vs. missing genotypes -4/4- . . . . . . . 179 8List of Figures 11.9. Imputation rates vs. number of alleles per marker . . . . . 181 11.10.Imputation rates vs. number of markers . . . . . . . . . . 182 11.11.Reconstructed haplotypes vs. missing genotypes -2/4- . . . 183 11.12.Reconstructed haplotypes vs. missing genotypes -3/4- . . . 184 11.13.Reconstructed haplotypes vs. missing genotypes -4/4- . . . 185 11.14.Reconstructed haplotypes vs. number of alleles per marker 186 11.15.Reconstructed haplotypes vs. number of markers . . . . . 187 11.16.Reconstructed vs. number of markers . . . . . . . . . . . . 188 17.1. Eficiencia de Genotipado . . . . . . . . . . . . . . . . . . 277 17.2. Gráfico cuantil-cuantil de los p-valores de sujetos controles 278 17.3. Tasa de pérdidas por individuo . . . . . . . . . . . . . . . 279 17.4. Ejemplo de redes de parentesco . . . . . . . . . . . . . . . 281 17.5. Parentesco entre los sujetos del estudio . . . . . . . . . . . 281 17.6. Histogramas de heterocigosidad antes del control de calidaddemuestras.......................281 17.7. Histogramas de heterocigosidad después del control de calidad de muestras . . . . . . . . . . . . . . . . . . . . . . 282 17.8. Gráfico Manhattan usando 1000G y MINIMAC3 . . . . 283 17.9. Gráfico Cuantil-Cuantil usando 1000G y MINIMAC3 . . 284 17.10.Gráfico Manhattan usando 1000G y IMPUTE2 . . . . . . 284 17.11.Gráfico Cuantil-Cuantil usando 1000G e IMPUTE2 . . . 285 17.12.Primera región con SNPs significativos . . . . . . . . . . . 286 17.13.Segunda región con SNPs significativos . . . . . . . . . . 286 17.14.Tercera región con SNPs significativos . . . . . . . . . . . 287 i PART Introduction 17 Chapter 1 Introduction 1.1. Background Genotype imputation and haplotype reconstruction have achieved an important role in Genome-Wide Association Studies (GWAS) during recent years. Estimation methods are frequently used to infer either missing genotypes as well as haplotypes from databases containing related or unrelated subjects. The majority of these analyzes have been developed using several statistical methods [1] which can impute genotypes as well as perform haplotype phasing (also known as haplotype estimation) of the corresponding genomic regions. Currently, algorithms do not carry out genotype imputation or haplotype reconstruction using deterministic techniques on pedigree databases, despite thefact thatcomputational inferencebyprobabilisticmodels may cause some incorrect results. These methods are usually focused on population data. In the case of pedigree data, families typically are comprised by duos (parent-child) or trios (parents-child) [2], whereas those studies focused on more than two offspring (for each line of descent) are uncommon. On the other hand, certain genomic regions are very stable against recombination but at the same time, they may be highly polymorphic. For this reason, in some well-studied regions, such as HLA loci [3] in the extended human Major Histocompatibility Complex (MHC) [4], an alphanumeric nomenclature is needed to facilitate later analysis. At this juncture, the available typing techniques usually are not able to determine the allele phase and, therefore, the constitution of the appropri19 . Introduction ate haplotypes is not possible. Although some computational methods have been evaluated for the reconstruction of haplotypes [5], none of them is capable to perform haplotype phasing or genotype imputation of missing data without using reference panels. Finally, although there is a growing number of bioinformatic/biostatistical tools for processing genomic databases, the documentation related to each of them is often somewhat confusing. Therefore, a need exists for clarification and simplification in the documentation relating to processes that include quality control, imputation of missing values and statistical analysis of association in genetic/genomic databases. 1.2. Motivation and Original Contributions The motivation for the realization of this work came from the necessity of organizing and applying different biostatistical and bioinformatics methods to solve several problems posed by diverse research groups in the field of endocrinology at the Complejo Hospitalario Universitario Insular Materno Infantil of Las Palmas de Gran Canaria. These problems were related to diabetes genetics and used different kind of data. On one hand, family-type genetic data (genetic information of parents and children in a number of families) were selected and, in other, population data came from a case-control study. The data analysis process required knowledge not only of statistical and computing methods but also of the basics of human genetics and the recent developments in methodologies for treating genetic data. For that reason, this document includes an informative part that intends to summarize these concepts and ideas, previous to the development of our original contributions, which can be synthesized as: 1) Identification of genetic variants associated with advanced diabetic nephropathy in a Type 2 Diabetes (T2D) population from the Gran Canaria Island. 2) GWAS Tutorial:Quality Control, Imputation, Analysis of population data. 3) Identification of haplotype associations in the international Type 1 Diabetes Genetics Consortium (T1DGC) pedigree database. 4) Development of the R package alleHap. 20 1.3. Document Structure 5) alleHap Manual:Allele Imputation and Haplotype Reconstruction from Pedigree Databases. Together with previous milestones, the author of this PhD thesis was also co-author of several publications in peer-reviewed journals as well as international proceedings/conferences. 1.3. Document Structure This document has been structured in 6 parts, each one containing, at least, three chapters. A brief description of each one is listed as follows: Part i will explain the motivation for this PhD dissertation, the original contributions of its author and the document structure. Since this thesis covers diverse research areas, Part ii intends to establish the state of the art, as well as to specify and clarify some basic concepts that may be useful for those who are not familiarized with genetics, biostatistics, and/or bioinformatics. Part iii will cover all necessary topics to develop the identification of genetic variants associated with a phenotype (disease) in a population. To achieve this purpose, quality control, alignment, haplotype phasing, genotype imputation and association analysis of genomic data were implemented. Part iv will comprise a description of the alleHap package, an analysis of its performance and the study of its application in the T1DGC database. Part v will present the main conclusions of previous chapters. Part vi will consist of two tutorials/manuals, one for the management of Genome-Wide Association (GWA) data and other for proper utilization of the alleHap package. As completion of this PhD dissertation, Part vii will include a summary of this thesis in Spanish, containing proposed goals, methodology, original contributions and final conclusions. 21 ii PART Basic Concepts and State of the Art 23 Chapter 2 Concepts of Human Genetics 2.1. Molecular Genetics Terminology Some basic terms used in human genetics are important to define before going further in this dissertation. From cell to gene, essential concepts of molecular genetics are explained and clarified in this chapter. Some of these concepts are represented in Figure 2.1. Figure 2.1: Cell, Chromosome, DNA and Gene representation, adapted from [6]. 25 . Concepts of Human Genetics Figure 2.5: Meiosis includes two nuclear divisions. The four daughter cells resulting from meiosis are haploid and genetically distinct. The daughter cells resulting from mitosis are diploid and identical to the parent cell. Meiosis two stage scheme, adapted from [32]. are separated, and the cell division results in two haploid cells. If this were the full description of meiosis, each of the 22 autosomes in a gamete would be an exact copy of one of the two parental homologous chromatids. In fact, a process called meiotic recombination mixes the genetic material of the homologous chromatids during meiosis, so that each chromosome present in the gamete has contributions from both parents [32]. 2.2.1. Meiotic recombination In meiosis I, homologous chromatids pair up and form physical connections called chiasmata (singular, chiasma). Chiasmata are essential for correct chromosome alignment and segregation and thus are thought to perform a role in meiosis I. Each chromosome arm normally forms at least one chiasma [9]. Meioticrecombinationorcrossing-overoccurs at the chiasmata. Specific enzymes break the DNA strands and repair the break in a way that swaps material from one chromosome with material from another (see Figure 2.6). The most important consequence of meiotic recombination is that gametes receive contributions from both homologs of a chromosome pair (thus from both grandparents) [9]. Two alleles that are linked (on the same chromosome) in the parent may or may not be linked in the offspring. As a result, a single person could theoretically produce an almost infinite number of genetically different gametes. 32 2.3. Human Genetic Diversity Figure 2.6: Effects of crossing over: the blue chromosome came from the individual’s father and the red chromosome came from the individual’s mother. Meiotic recombination, adapted from [32]. 2.3. Human Genetic Diversity Human genetic diversity is considerably lower than other species, including our nearest evolutionary relative, the chimpanzee. Genetic diversity is a function of a population’s ”age”, i.e. the amount of time during which mutations accumulate to generate diversity and its size. Our genetic homogeneity implies that anatomically modern humans arose relatively recently (200000 years ago) and that our population size was quite small at one time (10000 breeding individuals) [33]. According to Baker [34], the human genetic diversity can be estimated as 0.1-0.5%. Taking into account that humans have approximately 3 billion (3 ×109) base pairs in a haploid cell, it can be said that any pair of humans differs by approximately 3 to 15 million base pairs. These differences contain much useful information about the evolutionary history of our species [33]. Genetic variations mainly include mutations and polymorphisms, described in subsection 2.3.1. These DNA variations can be single base pair changes, deletions, insertions, inversions, translocations, changes in the number of copies of a given DNA sequence, or even duplications of whole chromosome sections, see Figure 2.7. 33 . Concepts of Human Genetics Figure 2.7: Genome structural variation encompasses polymorphic rearrangements 50 base pairs to hundreds of kilobases in size and affects about 0.5 percent of the genome of a given individual. Genetic Variations, adapted from [34]. 2.3.1. Polymorphism vs. Mutation DNAsequencevariations aresometimesdescribedasmutationsand sometimes as polymorphisms. It is important to clarify what is the difference between these terms and how they are applied to the human genome. The term polymorphism (a term that comes from the Greek words poly: ”many” and morphe: ”form”) is generally restricted to those variations that are relatively common (present in more than 1 percent of individuals) and that usually do not have highly deleterious consequences. Highly deleteriousrarevariantsareoftenreferredtoasmutations(present in less than 1 percent of the population) [35]. Therefore, to be classified as a polymorphism the least common allele must have a frequency of 1% or more in the population. If the frequency is lower than this, the allele is regarded as a mutation. The above definitions cannot be applied rigorously. Thus, within a population, a variant can be characterized by its lower allelic frequency (minor allele frequency) which it is simply the lowest frequency corresponding to the two alleles of the at a given locus. Given the variations among human populations, the least frequent allele of a certain locus within a population may be the most prevalent in another, namely a mutation in one population can become a polymorphism in another if it confers an advantage and increases in frequency. A good example is the allele of sickle-cell disease. In Caucasian populations, this is a rare 34 2.3. Human Genetic Diversity sequence variant of the beta-globin gene that causes a severely debilitating blood disorder. In certain parts of Africa, however, the same allele is polymorphic because it confers resistance to the blood-borne parasite that causes malaria [36]. 2.3.2. SNP: Single Nucleotide Polymorphism ThemostcommontypeofvariationinthehumangenomeistheSimple Nucleotide Polymorphism (SNP), which is a change in one base pair at a particular location of the genome. SNPs can be divided into two types, depending on the base substitution [37]: Transitions, the substitution of one purine for another (A ↔G) or one pyrimidine for another (C ↔T), are the most common type of SNP. Transversions, in which a purine is replaced by a pyrimidine, or vice versa, are less common. Figure 2.8: DNA molecule 1 differs from DNA molecule 2 at the same location. SNP representation, adapted from [38]. 35 . Concepts of Human Genetics SNPs are found more frequently in those DNA regions containing genes and may serve as biomarkers to identify what genes are associated with specific diseases. When SNPs occur within the gene, or in a regulatory region near a gene, they can significantly affect gene function [39]. The frequency of a particular SNP tends to remain stable in the population. Unlike the other, rarer kinds of variations, many SNPs occur in genes and in the surrounding regions of the genome that control their expression. The effect of a single polymorphism in a gene may not be large –perhaps influencing the activity of the encoded protein in a subtle way– but even subtle effects can influence susceptibility to common diseases [36]. 2.3.3. Highly Polymorphic Regions: MHC - HLA Complex Certain DNA regions of the human genome have a high polymorphism rate. From all these stands the MHC that contains the most diverse genes known in vertebrates. These highly polymorphic genes encode cell surface receptors that play a central role in distinguishing self/own from foreign proteins. The polymorphisms of MHC genes has been maintained by natural selection over long periods of evolutionary time [40]. The MHC is also known as HLA in humans [41]. This region encompasses 7.6×106bases on chromosome 6p21 and is the most gene dense region within the human genome encoding 252 loci [42] including several key immune response genes [43]. The region can be subdivided into Class I, Class II and Class III regions (see Figure 2.9). Figure 2.9: MHC complex (HLA region) in Human Chromosome 6. MHC-HLA Complex [44]. 36 2.3. Human Genetic Diversity The enormous polymorphism of MHC alleles is difficult to explain because natural selection should eliminate all but the most diseaseprotective allele. For example, a particular human MHC allele confers resistance to malaria in Africa. This process of directional selection should lead to all of the susceptible alleles becoming extinct, leaving only the most resistant allele in the population (fixation). Therefore, there must be some other evolutionary pressure that maintains MHC diversity [40]. Specific HLA alleles are associated with susceptibility and resistance to autoimmune and infectious diseases. The disparity between donor and recipient HLA-A, B, C, DR, DQA, DQB and DPA and DPB alleles impacts the outcome of both bone marrow and solid organ transplantation [44]. 2.3.3.1. HLA Allele nomenclature Early in their study, it was recognized that the genes encoding the HLA molecules were highly polymorphic and that there was a need for a systematic nomenclature. The HLA complex contains more than 220 genes of diverse function. Many of the genes encode proteins of the immune system [45]. The naming of new HLA genes and allele sequences and their quality control is the responsibility of theWHO Nomenclature Committeefor Factors of the HLA System. In 2010, a new HLA nomenclature system was adopted. The main drive for the change was that the old system could no longer accommodate the increasing number of HLA alleles that were being described, due to the fact that HLA complex is the most polymorphic region of the entire human genome with close to 9000 different HLA alleles characterized thus far [44]. The list continues to expand rapidly as increasing numbers of new alleles continue to be identified. The list containing the most updated alleles is in IMGT/HLA database. As is depicted in Figure 2.10, have been adopted colons ‘:’ as separators between pairs of digits. HLA-A*02010102L therefore became HLA-A*02:01:01:02L.Thepairsofdigitsseparatedbycolonsareknown as Fields. The first and second digits of the old nomenclature form the 1st Field of the new nomenclature. The third and fourth digits of the old nomenclature form the 2nd Field of the new nomenclature. To help reduce confusion in adopting the new nomenclature, the leading ‘0’ in alleles 1-9 of each allele group was kept [46]. 37 . Concepts of Human Genetics Figure 2.10: Each HLA allele name has a unique number corresponding to up to four sets of digits separated by colons. The length of the allele designation is dependent on the sequence of the allele and that of its nearest relative. HLA Alleles Nomenclature [46]. 38 Chapter 3 Concepts of Statistical Genetics 3.1. Statistical Genetics Statistical genetics has experienced exponential growth during last years, with an increasing number of people involved in this area. It can be affirmed that principles of statistical genetics synthesize previous statistical models and methods applied to genome study. Statistical methodologies behind certain computational software packages may help researchers to choose the optimal tool to analyze their data and also interpret their studies better [47]. Primarily, contemporary statistical genetics focuses on the development and implementation of data analysis methodologies that may facilitate the identification of genes and genetic variations that influence phenotypic expression and disease susceptibility [48]. As it is possible that a single SNP or a set of a few SNPs contribute significantly to a disease, the fact that there are millions of SNPs makes it particularly difficult to detect their real causal association with a specific disease. Furthermore, most conditions for current public health research, are complex and multifactorial, since they usually are composed of many genes and environmental factors, as well as their corresponding expression or interactions [49]. There are two general designs in genetic association studies:familybased designs that use pedigrees and population-based studies that employ unrelated individuals. The recruitment of unrelated individuals is easier than the recruitment of families, but they are subject to bias in the presence of population stratification (individuals with a significantly 39 . Concepts of Statistical Genetics different genetic ancestry and phenotype from the rest of the study). As a compromise between linkage studies (which usually involve large families where the disease affects individuals in several generations) and population-based association studies, family-based association designs can have similar power as population-based designs and are more robust in the presence of population stratification [50]. 3.2. Population Studies Throughout the history of population genetics, statistical models haveplayeda significantrolein explaining theeffectsof geneticdiversity of organisms. Such models today are crucial in the development of statistical tools for analyzing molecular biology data. Regardless of assumptions about the genetic model of a trait, or technologies used to assess genetic variation, no genetic study will have meaningful results without a thoughtful approach to characterize the phenotype of interest. When embarking on a genetic study, the initial focus should be on identifying precisely what genetic variation influences [51]. The following subsections will provide a brief introduction to some concepts of population genetics. 3.2.1. Case-Control and Quantitative Designs There are two primary classes of phenotypes: categorical (often binary case/control) or quantitative. From the statistical point of view, quantitative traits are preferred because they improve power to detect a genetic effect, and often have a more interpretable outcome [51]. Anexampleofquantitativetrait ischolesterol levels, which are strong predictors of heart disease. The analysis of such levels is very useful for clinical practice since they are precise and ubiquitous measurements that are easy to obtain [51]. Genetic variants that influence these levels have a clear interpretation (e.g. a unitary level change per allele), having, therefore, an easily measurable effect on the quantitative trait. Other disease traits do not have well-established quantitative measures. In these circumstances, individuals are usually classified as either affected or unaffected (case or control), i.e. a binary categorical variable [51]. Although quantitative outcomes are preferred, they are not always required for a successful study. 40 3.2. Population Studies 3.2.2. Population Association Analysis In recent years, GWAS have become the most used tool for the identification of loci associated with complex traits. Through this method, association between a trait of interest and genetic polymorphisms is studied using individuals’ samples typed for millions of SNPs [52], allowing for the discovering of numerous statistical associations between genomic variants and quantitative traits or complex diseases. When genotypes are collected, and a well-defined phenotype has been selected for a population study, the statistical analysis of genetic datacan begin[51]. Suchanalysis can be performed eitherwith seriesof single-locus tests (examining each SNP independently for association with phenotype) or by multi-loci tests (considering interactions among different genetic variants throughout the genome). 3.2.2.1. Single-locus Analysis For both quantitative and dichotomous trait analysis, there are numerous ways to handle genotype data for association tests. The choice of the model can have implications for the statistical power of the test, as the degrees of freedom may change depending on genotype classes (i.e.: homozygous/heterozygous, dominant/recessive, etc.). Allelic association tests, therefore, examine the association between one allele of the SNP and a single or multiple phenotypes. Statistical analysis of possible genotype independence at a single locus is simple because the absolute frequencies observed in both of them will be presented in a double entry table. In such a table, independence will be studied by an independent χ2test (see next paragraph) using Yates correction for continuity (or Yates’ χ2test) [53] if the expected cell frequencies are less than 5. 3.2.2.1.1. Hardy-Weinberg Equilibrium Hardy-WeinbergEquilibrium(HWE) (also knownasHardy-Weinberg principle) states that alleleandgenotype frequenciesin a populationwill remain constant from generation to generation in the absence of other evolutionary influences [54]. These influences include mate choice, mutation, selection, genetic drift, gene flow and meiotic drive. HWE denotes independence of alleles at a single site between two homologous chromosomes [55]. 41 . Concepts of Statistical Genetics 3.2.2.2.1.1. Measures of Linkage Disequilibrium To explain the principal measuresof LD, firstly we have to consider a distribution of alleles for nindividuals across two loci. Then, assuming that the two loci are independent of each other, i.e. they are in linkage equilibrium, the presence of an allele at one locus should not influence the particular allele observed at the second locus. Supposing Aand aas possible alleles at Locus 1, and Band bas possible alleles at Locus2, the marginal probabilities of the alleles A, a,B,bwill be pA,pa,pB, and pb, respectively. Since each individual carries two homologous chromosomes, there will be a total of N= 2n homologs across the nsubjects in a population. Locus 2 Locus 1 A a AN(pApB+D)N(pApb−D) aN(papB−D)N(papb+D) Table 3.2: Observed allele distributions under LD, adapted from [55]. On the other hand, if the two loci are associated, the expected values will have be deviated a quantity. Such deviation is commonly symbolized by the scalar Dand its amount is represented in Table 3.2. We can express Dregarding the joint probability of Aand B, and the product of the individual allele probabilities as follows: D=pAB−pApB(3.13) the value that estimates the Disequilibrium being the scalar D′: D′=|D| Dmax (3.14) where Dmax represents the upper bound on Dand is given by: Dmax ={min(pApb, papB)ifD > 0 min(pApB, papb)ifD < 0(3.15) Note that D′will be a value with range: 0≤D′≤1, so that D′ values close to 0 will suggest linkage equilibrium (almost no association), 48 3.2. Population Studies while values close to 1 will indicate high levels of LD and depending on the study, a possible genetic association with a disease. Another measure that is also a function of the scalar Dis the quantity r2. This measure is based on Pearson’s χ2test of no association between the rows and columns of table 3.2. Specifically, r2can be defined as: r2=χ2 1 N(3.16) Where, if r2is written in terms of the scalar D, results: r2=D2 pApBpapb (3.17) The difference between Dand r2rests in the type of adjustment madetothescalarD. In bothcases, this adjustmentinvolves themarginal allele frequencies since the value of Dwill depend on these. Investigators commonly use r2[55], due to its simpler relationship to the usual Pearson’s χ2test. 3.2.2.2.2. Indirect Association The presence of LD creates two possible positive outcomes from a genetic association study. In the first one, the SNP influencing a biological system that ultimately leads to the phenotype is directly genotyped in the study and found to be statistically associated with the trait. The above is referred to as a direct association, and the genotyped SNP is sometimes referred to as the functional SNP. The second possibility is that the influential SNP is not directly typed, but instead, a tag SNP in high LD with the influential SNP is typed and statistically associated with the phenotype (see Figure 3.3). This is referred to as an indirect association [60]. Because of these two possibilities, a significant SNP association from a GWAS should not be assumed as the causal variant and may require additional studies to map the precise location of the influential SNP [51]. 49 . Concepts of Statistical Genetics Figure 3.3: Genotyped SNPs often lie in a region of high linkage disequilibrium with an influential allele. The genotyped SNP will be statistically associated with the disease as a surrogate for the disease SNP through an indirect association. Indirect Association [51]. 3.3. Family-based Studies As was mentioned above, Family-based association studies have several advantages in comparison to population-based association studies. Firstly, the need for match cases to controls in population studies may lead to selection bias and confounding effects if gene frequencies differ between case and control populations. By choosing controls from the same families as the cases, the confounding effects will be substantially reduced [61]. Furthermore, as a family-based association can only be detected when the linkage is present, the identification of such association confirms, despite some difficulties in interpreting results, that truly associated markers are physically close to the causal genetic variant, supporting the fact that phenotypes have been inherited [61]. On the other hand, family-based studies have some disadvantages that mainly arise from practical matters of recruitment and cost [62]. This is the main reason why population controls are commonly used for large-scale genome association studies. 3.3.1. Family-based Association Tests Family-Based Association Test (FBAT) include several ways for detecting associations between specific markers and quantitative phenotypes or diseases. In these studies, nuclear families consisting at least of two parents and a number of full siblings are widely used, but extended pedigrees may also be used for testing association [61], often improv50 3.3. Family-based Studies ing expected results. Such association methods can be broadly classified into two groups: nonparametric methods (based on allele counting) and parametric methods (based on the likelihood function). For the simplest family-based association design (two parents and one affected offspring) previous methods result in similar test statistics, i.e. their extensions on more complex situations vary considerably [50]. Typicalstatisticaltestsusegeneralpedigreesandmaydetect genotypeenvironment interactions. The subsections that follow briefly describe some of the most used ones. 3.3.1.1. Transmission/Disequilibrium Test The simplest family-based design for association studies is the case– parent design, in which an affected offspring and both parents are genotyped at bi-allelic markers. A method called TDT for evaluating this was proposed by Spielman and Ewens [63]. In this, the alleles transmitted from parents to the affected offspring and the alleles of not transmitted can be determined based on the observed genotype data [50]. Thus, a two by two transmission/nontransmission table for a bi-allelic marker with alleles Aand Bfrom tcase–parent trios can be constructed: Non-transmitted Transmitted A B Total AnAA nAB NA· BnBA nBB NB· Total N·AN·B4t Table 3.3: Summary of the transmission/nontransmission allele counting for a biallelic marker. TDT allele counting, adapted from [50]. In Table 3.3, nrepresents the number of parents who have genotype A,B and transmit allele Ato the affected offspring. Another approach for the above data would be: TDT =(nAB −nBA )2 nAB +nBA (3.18) whichcomparesthenumberof Aalleles transmittedto the offspring from theirs parents and the number of Aalleles not transmitted. 51 . Concepts of Statistical Genetics Also, it is important to note that the TDT has two key benefits: 1) It only assumes Mendel’s first law of inheritance. The specification of the disease model and the distribution of the disease in the general population are not required and will not affect its validity. Thus, the TDT is not only robust to population stratification but also robust to any misspecification of the disease model and the distribution of the disease. 2) The TDT test statistic has an asymptotic chi-square distribution with one degree of freedom if either θ= 12 or δ= 0, where θand δare the recombination fraction and the linkage disequilibrium, respectively [50]. 3.3.1.2. Sib Transmission/Disequilibrium Test The TDT requires marker genotypes for affected individuals and their parents. However, for some diseases data from parents may be difficult or impossible to obtain. A method called sib TDT or Sib Transmission/Disequilibrium Test (S-TDT) for describing this was implemented by Spielman and Ewens [64]. This overcomes the previous problem by use of marker data from unaffected sibs instead of from parents, thus allowing application of the TDT to sibships without parental data. Namely, the S-TDT method uses sibships consisting of at least one affected and one unaffected sib. In the S-TDT method the number of variant alleles in affected sibs is counted, calculating its mean and variance for each family under the assumption that the proportion of variant alleles is the same in affected sibs as it is in unaffected sibs. These counts are summed over the set of families to form a z score. Several other sibling-based methods have been suggested [65, 66, 67] though the S-TDT remains most closely related to the more recent approaches. Some families will be suitable to be analysed only by the TDT, and others could be described by the S-TDT. The work by Spielman et al. [64] also explains how all the data may be used jointly in one overall TDT-type procedure that tests for linkage in the presence of association. These extensions of the TDT could be interesting for the study of diseases associated with aging [64]. 3.3.1.3. Pedigree Disequilibrium Test The Pedigree Disequilibrium Test (PDT) combines the principles of the TDT and S-TDT into a test for general pedigrees [68]. It splits a 52 3.3. Family-based Studies pedigree into a list of all case-parent trios and discordant sib pairs with genotype data. For a trio j,XTjis defined as the number of transmissions of the variant allele minus the number of its non-transmissions. For a sib pair j,XSjis defined as the number of copies of the variant allele in the affected sib minus the number in the unaffected sib [61]. The measure of association Dfor the pedigree would be: D=1 NT+NS[NT ∑ j=1 XTj+ NS ∑ j=1 XSj](3.19) where NTand NSare the total number of trios and discordant sib pairs, respectively. D has expectation 0 in any pedigree. After computing this measure for each i= 1, ..., N pedigrees, the PDT statistic Tis then: T= N ∑ i=1 Di N ∑ i=1 D2 i (3.20) This gives a valid test of linkage or association in any pedigree structure, although some pedigrees are uninformative, notably the affected sib pair. The PDT has been adapted to test quantitative traits [69], haplotypes [70] and genotypes [71]. A very flexible approach for constructing unbiased tests is implemented in the software FBAT [72]. 3.3.2. Pedigree Structures with Missing Data Two general approaches have emerged to deal with the problems of missing family members: 1) Fitting a statistical model to the missing data and conducting an analysis that takes all of the possible completions into account. 2) Developing test statistics that are unbiased under the null hypothesis while using only the available data. This approach retains complete robustness to population stratification and can be readily applied to arbitrary pedigree structures [61]. For instance, a case of particular interest is the late-onset disease in missing parents [61], where unaffected siblings may be used as controls, as long as their relationship to the affected ones has been considered. 53 Chapter 4 Concepts of Computational Genomics 4.1. Computational Genomics Computational genomics refers to the use of computational and statistical analysis to decipher biology from genome sequences and related data [73]. Computational genomics focuses on understanding the human genome, and more generally the principles of how DNA controls the biology of any species at the molecular level. With the current abundance of massive biological datasets, computational studies have become one of the most important means to biological discovery [74]. During thepast few years, therehave beenenormous advances ingenomics and molecular biology, which carry the promise of understanding the functioning of whole genomes in a systematic manner [75]. The challenge of interpreting the vast amounts of biological data has led to the development of new tools in the fields of computational biology and bioinformatics, and opened new connections to areas such as chemometrics, exploratory data analysis, statistics, machine learning, and graph theory. 4.2. Computational Biology vs. Bioinformatics Computational biology can be defined as the study of biology using computational techniques. Its primary goal is to learn new biology, i.e. knowledge about living systems from a more scientific perspective. 55 . Concepts of Computational Genomics Bioinformatics is more focused on the creation of useful tools, algorithms, and software to solve problems relating to biological data. It could also be said that bioinformatics usually manages biological problems and data from an engineering point of view. Mathematical (statistical) and computational techniques are continuously being developed, analyzed, and improved to offer greater utility in biological arenas. Computational biology and bioinformatics employ powerful tools from computer science, applied mathematics, and statistics to solve those biological problems. Majoravenuesofresearchincludegenomeassembly–piecing together a large collection of short DNA sequences–, gene finding –locating patches of DNA that have biological function–, genomic sequence alignment – ordering of multiple sequences to elucidate their parallel structure–, protein structure prediction –determining the 3D structure of proteins from their chemical makeup–, and phylogenetic analysis –studying and modeling evolutionary relationships between species– [12]. 4.3. Machine Learning in Biology The term Machine Learning [76] (also known as Data Mining) is a research area of computer science which refers to a set of topics dealing with the creation and evaluation of algorithms that facilitate pattern recognition, classification, and prediction, based on models derived from existing data. Its development in recent years is due to the advances in data analysis research, growth in the database industry and the resulting market needs for methods that are capable of extracting valuable knowledge from large databases [77]. In biology-related research areas, machine learning appears as one of the main drivers of progress, where most of the targets of interest deal with complex structured objects: sequences, 2D, and 3D structures or interaction networks. At the same time, bioinformatics and systems biology have already induced new significant developments of general interest in machine learning, for example in the context of learning with structured data, graph inference, semi-supervised learning, system identification, and novel combinations of optimization and learning algorithms [78]. Molecular biology in particular and more generally all biomedical sciences are undergoing a genuine revolution as a result of the emergence and growing impact of a series of new disciplines sharing the 56 4.3. Machine Learning in Biology -omics suffix in their name. These include in particular genomics, transcriptomics, proteomics, and metabolomics, devoted respectively to the examination of the entire systems of genes, transcripts, proteins and metabolites present in a given cell or tissue type [79]. Figure 4.1 shows a scheme of some of these biological domains where computational methods are applied for knowledge extraction from biological data. Figure 4.1: Classification of some topics where machine learning methods are applied. Machine learning topics in Biology, adapted from [80]. 4.3.1. Machine Learning Approaches and Paradigms As machine learning is mainly concerned with the discovery of models, patterns and other regularities in data, its approaches can be firstly categorized as follows: Symbolic approaches, including inductive learning of symbolic descriptions, such as rule learning, decision trees or logical representations. Statistical approaches, including Statisticalorpattern-recognition methods as k-Nearest Neighbors (k-NN) or instance-based learning such as Bayesian classifiers, neural network learning, and Support Vector Machine (SVM). 57 . Concepts of Computational Genomics Figure 4.6: a) SNPs. A short stretch of DNA from four versions of the same chromosome region in different people. b) Haplotypes. Observed genotypes for 20 SNPs that extend across 6000 bases of DNA. c) Tag SNPs. Genotyping just the three tag SNPs out of the 20 SNPs is sufficient to identify these four haplotypes uniquely. SNPs, haplotypes, and tag SNPs [95]. The DNA samples of HapMap Project in Phase I,II and III came from a total of 1184 individuals from the following 11 populations: ASW: African ancestry in Southwest (USA), CEU: Utah residents with Northern and Western European ancestry (USA), CHB: Han Chinese in Beijing (China), CHD: Chinese in Metropolitan Denver (USA), GIH: Gujarati Indians in Houston (USA), JPT: Japanese in Tokyo (Japan), LWK: Luhya in Webuye (Kenya), MXL: Mexican ancestry in Los Angeles (USA), MKK: Maasai in Kinyawa (Kenya), TSI: Toscani in Italy and YRI: Yoruba in Ibadan (Nigeria) [96]. The development of the HapMap has enabled geneticists and other specialists to take the advantage of how SNPs and other genetic variants are organized on the same chromosome [95]. The HapMap project has also helped researchers to find functional regions that influence human health outcomes as well as responses to therapeutic drugs and environmental factors. The project itself has not identified such regions directly. Instead, HapMap has provided a tool that can be used in both population-based and family-based disease association studies [97]. 64 4.5. International Genetic Databases 4.5.2. 1000 Genomes Project The 1000 Genomes Project [98] has sequenced the genomes of more than 1000 people, to provide a comprehensive resource on human genetic variation. Its primary goal was to find most genetic variants that have frequencies of at least 1 percent in the populations studied [98]. The samples used in this project were mostly anonymous and had no associated medical or phenotypic data. Despite this fact, the genetic variation data produced by the 1000G Project has been used by researchers to study many diseases, in sets of case and control samples thatwerecarefullyphenotyped [98]. Extended informationabout using the project data is available in 1000 Genomes Project website. This resource, which captured up to 98 percent of accessible SNPs at a frequency of 1 percent in related populations, has enabled numerous analysis of common and low-frequency variants in subjects from diverse, including admixed, populations. The individual samples collected were grouped in five geographic areas according to the regions of the ancestries [99]: I. EastAsianAncestry(EAS): Chinese Dai in Xishuangbanna (CDX), Han Chinese in Bejing (CHB), Vietnamese Kinh in Ho Chi Minh City (KHV), Southern Han Chinese (CHS). II. South Asian Ancestry (SAS): Bengali in Bangladesh (BEB), Gujarati Indian in Houston (GIH), Indian Telugu in the UK (ITU), Pakistani Punjabi in Lahore (PJL), Sri Lankan Tamil in the UK (STU). III. AfricanAncestry(AFR): AfricanAncestryinSouthwestUS(ASW), African Caribbean in Barbados (ACB), Esan in Nigeria (ESN), Gambian in Western Division (GWD), Kenyan Luhya in Webuye (LWK), Mende in Sierra Leone (MSL), Nigerian Yoruba in Ibadan (YRI). IV. EuropeanAncestry(EUR): BritishinEnglandandScotland(GBR), Finnish in Finland (FIN), Iberian populations in Spain (IBS), Toscani in Italia (TSI), Utah residents with Northern and Western European ancestry (CEU). V. AmericasAncestry(AMR): ColombianinMedellin(CLM),Mexican Ancestry in Los Angeles (MXL), Peruvian in Lima (PEL), Puerto Rican in Puerto Rico (PUR). 65 . Concepts of Computational Genomics In this project 2500 samples at 4X coverage have been sequenced. The first set of samples for sequencing included 1167 samples and were collected from 13 populations during 2010 and early 2011. The second set included 633 samples that was collected from 7 populations in early 2011. The third set, consisting of 700 samples, was collected for sequencing in late 2011 [100]. Full details of the samples are shown in Table 4.1. Population Group Pilot Samples Set  Samples Set  Samples Set  Samples Total EAS 185 286 515 504 523 SAS 0 0 494 489 494 AFR 208 246 669 661 691 EUR 160 379 505 503 514 AMR 0 181 352 347 355 Total 553 1092 2535 2504 2577 Table 4.1: Summary of the sequenced samples the 1000 Genomes project. 1000 Genomes Samples, adapted from [99]. Finally, the international 1000 Genomes Project has provided a validated haplotype map of more than 84 million SNPs, 3.1 million short insertions and deletions, 42.279 biallelic deletions, 6.025 biallelic duplications, 2.929 mCNVs (multiallelic copy-number variants), 786 inversions, 168 Nuclear Mitochondrial Insertion (NUMT)s, and 16.631 mobile element insertions of 2.504 individuals from 26 populations. Despite the great genetic diversity, most variants (86% of 84.7 million) are restricted to a single continental group, particularly among sub-Saharan populations in Africa [100]. 4.5.3. T1DGC Project The T1DGC [101] is an international, multicenter program organized to identify genes and their alleles that determine an individual’s risk for Type 1 Diabetes (T1D). The program had two primary goals: I. Identifying genomic regions and candidate genes whose variants modify an individual’s risk of T1D and help explain the clustering of the disease in families. 66 4.5. International Genetic Databases II. Making research data available to and establish resources that can be used by the research community. TheT1DGC assembledaresource ofaffected sib-pairfamilies, parentchild trios, and case-control collections with banks of DNA, serum, plasma, and cell lines [101]. In addition to T1DGC-recruited Affected Sib-Pair (ASP) families, the T1DGC recruited trio families from ethnic groups with a lower prevalence of T1D. The T1DGC also welcomed the inclusion of earlier ascertained case-control collections [102]. Research with T1DGC data has included genome-wide linkage scans, evaluation of the human MHC, examination of published candidate genes for T1D, and analysis of autoimmune disease genes and those affecting β-cell function in type 2 diabetes [101]. All the information on T1DGC can be accessed at the T1DGC website. Figure 4.7 shows the criteria of pedigree structure for inclusion of families. The maximal included pedigree structure includes five affected and two unaffected siblings in affected sibling pair families; no additional siblings were collected in trio families. All recruited family members were typed for Class I and II HLA loci [101]. Figure 4.7: Dark fill represents a family member with type 1 diabetes; no fill, unaffected; and crosshatch may be either. The dotted line indicates the minimum inclusion criteria for family recruitment into the T1DGC collection. Pedigree structures into the T1DGC, adapted from [103]. Eligibility criteria [102] included the following: 1. Siblings with a diagnosis of T1D. 2. Diagnosis before 35 years of age. 3. Use of insulin within 6 months of diagnosis. 4. Continuous use of insulin (without stopping for 6 months or more). 5. Informed consentforbloodcollection, geneticanalysis, andexam. 67 . Concepts of Computational Genomics In addition, trio families were collected in selected populations with a low prevalence of the disease throughout the Asia-Pacific, European and North American Networks. The required family structure was an affected child and both biological parents. Eligibility criteria are the same as listed above. Cases and controls were collected throughout the Asia-Pacific, European and North American Networks in selected lowprevalence populations. Outcome measures included the establishment of resources for research into the genetic origins of T1D and identification of genomic regions and genes whose variants contribute to an individual’s risk of T1D. Phenotype and genotype data from study participantshas beenwidely used in researchstudies concerningthe genetic origins of T1D risk in families and the general population [102]. 68 Chapter 5 State of the Art 5.1. Introduction The knowledge about human genetic variation has been growing exponentially over the last decade. Collaborative efforts such as the International HapMap [97] and 1000 Genomes [98] projects have contributed to increase the rate of discovery about human genetic diversity. Genome population-based studies usually employ DNA microarrays to produce genotype information for a set of individuals. If those studies where collected with different array types, some markers may not be assayed at the same genomic positions. However, recent computationaladvancesenableresearcherstousealgorithms tofill inor impute genotypes (see Figure 5.1) at the markers that are not common among the genotyping arrays [104]. Given the genotypes of a sample of individuals from a population, using haplotype pre-phasing to infer firstly the haplotypes of the sample (using haplotype sharing information within the sample) allows to build a phased reference panel which is used to estimate later missing markers of the sample. Therefore, genotype imputation (together with haplotype pre-phasing) techniques usually increases the sample size at each marker, boosting thus the study’s power to detect fine-map associations and facilitating the combination of results across different studies using meta-analysis [85]. This chapter comprises the state of the art of several different statistical methods for genotype imputation,haplotype reconstruction, and analysis of genome variation. 69 . State of the Art Figure 5.1: SNPs 1–9 form three blocks of high LD, indicated by the red diamonds between the SNPs. Imputation overview, adapted from [104]. 5.2. Genotype Imputation Genotype imputation can be fundamental in the analysis of GWA. The most common approach works by finding haplotype segments that are shared between study individuals, which are typically genotyped on a commercial array with 300000–2500000 SNPs [105]. This process can be carried out across the whole genome as part of a GWAS or in a more focused region as part of a fine-mapping study. The goal is to predict the genotypes at the SNPs that are not directly genotyped in the study sample [85]. These in-silico genotypes can then be used to boost the number of SNPs which can be tested for association, thus increasing the power of the study, improving fine-mapping of causal variants and facilitating meta-analysis. 70 5.2. Genotype Imputation On the other hand, many existing genotype imputation methods require substantial computing power to run using large reference datasets. This problem may be aggravated by the fact that reference panels are regularly improved and expanded, in order to investigators can re-impute their samples multiple times over the course of a study [105]. 5.2.1. Uses of Genotype Imputation The uses of genotype imputation are several and can range from the imputation of untyped variation to meta-analysis. Next subsections enumerate and describe them briefly. 5.2.1.1. Imputation of Untyped Variation Imputation of SNPs that have not been typed in either the haplotype reference panel or the study sample is also possible. Some methods do this via inference of the genealogy between study sample haplotypes [106, 107] while others aim to identify haplotype effects more directly [108]. These methods can lead to a boost in power, especially when the causal variant is rare, or where there is local heterogeneity in the signal of association [109]. 5.2.1.2. Imputation of Non-SNP Variation The general idea of imputation is readily extended to other types of genetic variation such as Copy Number Variant (CNV)s and classical HLA alleles [110]. The imputation of large numbers of small insertions and deletions (indels) which will be discovered from recent international sequencing based projects (as the 1000 Genomes project [100]) are likely to be widely adopted in GWAS studies [116]. 5.2.1.3. Boosting Power Imputation can lead to a notable boost in the power of a GWAS. Spencer et al. [111] illustrate how imputation could produce a power improvement of 5 to 10% if the density of the chips is close to that of a hypothetical complete chip consisting of all SNPs. Other simulations have shown that the biggest benefit occurs for rare SNPs that are harder to tag [112]. 71 . State of the Art 5.2.1.4. Fine-mapping Imputation provides a much higher resolution view of an associated region than would be seen by just considering genotyped SNPs (see Figure 5.2) and increases the chance that a causal SNP can be directly identified. When imputedSNPsproducelargersignalsthananyofthegenotyped SNPs they can become better candidates for replication in new samples. Imputation methods can also help elucidate when multiple variants or allelic heterogeneity occurs in a region of interest [113, 109]. Figure 5.2: Association of genetic variants near LDLR with LDL-cholesterol levels using imputed data. An example of genetic variants association using imputed data [114]. 5.2.1.5. Meta-analysis Imputationhasbeenwidelyusedtofacilitatemeta-analysisofGWAS from different cohorts that may have been genotyped using different genotyping chips, i.e., different sets of SNPs. Imputation effectively enlarges and equates the set of SNPs in each study for which genotypes are available for testing. A useful practical guide is provided by [115]. Results from cohort-specific GWAS are combined using fixed effects models rather than combining the raw data from all studies and then carrying out one association test. 72 5.2. Genotype Imputation 5.2.2. Genotype Imputation Methods Severalmethodshavebeenproposedforgenotypeimputation. These methods can provide a boost in imputation accuracy, mostly at rarer SNPs and those SNPs that are not well tagged by a small number of flanking SNPs. Most imputation methods are based on HMM, a very useful class of statistical model readily applicable in genetics. HMM model is used for relating an observed process across the genome to an underlying, unobserved process of interest [116]. Thus, the most commonly used programs for genotype imputation are: 5.2.2.1. IMPUTE v1 IMPUTE1 [112] is based on an extension of the HMM models originally developed as part of importance sampling schemes for simulating coalescent trees [117, 118] and for modeling linkage disequilibrium and estimating recombination rates [119]. The method is based on the HMM of each individual’s vector of genotypes conditional upon a reference set of haplotypes and an estimate of the fine-scale recombination map across the region. Exact marginal probability distributions for the missing genotypes are obtained using the forward-backward algorithm for HMMs [120]. 5.2.2.2. IMPUTE v2 IMPUTE2 [121] takes a different and more flexible approach. SNPs are first divided into two sets: a T set that is typed in both study sample and reference panel, and a U set that is untyped in the study sample but typed in the reference panel. As is depicted in Figure 5.3 the algorithm estimates haplotypes at SNPs in T (blue) using IMPUTE1 and then imputing alleles at SNPs in U (green) conditional upon the current estimated haplotypes. Since imputation performance is driven by accurate matching of haplotypes, the method focuses on accurate haplotype estimation at the SNPs in T using as many individuals as possible [85]. 5.2.2.3. MACH MACH [122] uses an HMM model very similar to that used by HOT-SPOTTER[119] andIMPUTE. Usingthis method, genotypephasing can be performed and, therefore, used for imputation. The method works by successively updating the phase of each individual’s genotype 73 . State of the Art 5.3.3. Haplotype Reconstruction within Pedigrees Existingapproachesfor haplotype reconstructioncan be categorised according to the type of cohort each method is designed to phase, and the level of relatedness between the individuals in such cohort [142]. Much of the recent literature is devoted to phasing nominally unrelated (or distantly related) individuals. Currently, the most accurate methods use HMMs to model local haplotype sharing between individuals [143, 140] and take advantage of LD. Some of these methods can also handle mother-father-child trios and parent-child duos [144, 2]. For more complex pedigrees there are several comprehensive pedigree analysis software packages [145, 146]. As a general strategy for phasing cohorts with any level of implicit or explicit relatedness between individuals, a novel HMM method has been recently proposed by O’Connell et al. [142]. To utilize this program (called DuoHMM), firstly SHAPEIT2 has to be run ignoring all explicit family information, to combine later those haplotypes with any family information to infer the inheritance pattern. This method has the advantage that allows the detection of recombination events, genotyping errors and the correction of switch errors (see Figure 5.6). Figure 5.6: Paths for a three father-child duos from a nuclear family on chromosome 10. Two possible IBD states are shown using the colors light and dark red. Haplotype correction example using the DuoHMM method [142]. In Figure 5.6, the left panel shows the path prior to any corrections, and the right panel after a minimum recombinant correction is applied. The second and third sibling initially had a transition at around 25 mb, which is more likely a recombination event in the first child hence the parental haplotypes are switched at this point. The panel on the right 80 5.3. Haplotype Reconstruction has the corrected haplotypes, the number of recombination events required to explain the observed data has been reduced [142]. 5.3.4. Haplotype Reconstruction programs Numerous software packages are tackling the haplotype reconstruction problem. The following tables group them in three categories: most used, R packages and other related programs. Package Name Last Ver. Main aims Samples Type Ref. BEAGLE 4.1 Imputation, phasing and analysis. Unrelated and nuclear families. [108] fastPHASE 1.2 Imputation and phasing. Unrelated [126] IMPUTE 2.3.2 Imputation and phasing. Unrelated [121] MACH 1.0.18 Imputation and phasing. Unrelated [122] Minimac3 1.0.13 Imputation and phasing. Unrelated [124] SHAPEIT 2.20 Alignment and phasing. Unrelated and nuclear families. [141] Table 5.2: Most used programs for genotype imputation and haplotype phasing. There are almost no R packages related to haplotype reconstruction, inference, phasing or assembly. The ones listed in the CRAN repository are shown in the following table: Package Name Last Ver. Main aims Samples Type Ref. haplo.ccs 1.3.1 Relative risk estimation Unrelated case-control [147] haplo.stats 1.7.1 Inference and analysis Unrelated [148] hsphase 2.0.1 Imputation and phasing Half-sib families [149] Table 5.3: R packages related to haplotype reconstruction. 81 . State of the Art Other related but less-used programs encompass from genotype calling, haplotype assembly to LD mapping or haplotype association analysis. Program Name Last Ver. Main aims Samples Type Arlequin 3.5.2.2 Phasing and analysis Unrelated DuoHMM 0.1.7 Phasing Complex pedigrees HapCUT 0.6 Assembly Unrelated HAPI-UR 1.01 Phasing. Unrelated HaploBlock 1.2 Phasing and LD mapping Unrelated HaploRec 2.3 Reconstrucion Unrelated HapFerret - Phasing Unrelated HapSeq 2 Calling and phasing Unrelated HARSH 0.21 Phasing. Unrelated MERLIN 1.1.2 Phasing, LD mapping and analysis General pedigrees PLINK 1.07 Imputation, phasing and analysis Unrelated and pedigrees PedPhase 3.0 Phasing Pedigrees S-MIG++ 1.0.0 Phasing and LD mapping Unrelated SpeedHap - Phasing Unrelated ReFHap 1.0.0 Phasing. Unrelated WinHAP 2.0 Phasing Unrelated Table 5.4: Other programs related to haplotype reconstruction. 82 5.4. Analysis of Genome Variation 5.4. Analysis of Genome Variation Over the last decade, there has been an increasing number of highprofile GWAS for a variety of different human diseases (see Figure 5.7). These studies have revealed hundreds of disease-associated loci and have provided insights into the study of complex traits. All these kind of studies produced a test statistic and a p-value indicating how significant the statistical association between a given SNP and a phenotype is, i.e. how likely a specific phenotype (disease) may have occurred by chance. Figure 5.7: Published GWA reports between the beginning of 2005 and end of 2013. GWA published reports, adapted from [150]. Usually, GWAS require stringent significance levels (p≤5×10−8) to overcome the multiple testing problem incurred when testing SNPs throughout the genome. However, in studies with a short-medium sample sizes (less than 1000 subjects), a more lax threshold has to be utilized in order to detect associations (in spite of losing statistical power). In these cases, GWAS may be insufficient to explain the functional or causal variants since the identified associations discovered commonly have small ORs (< 1.5), which suggests that these effects are not such significant. Recently, the National Human Genome Research Institute (NHGRI) in collaboration with European Bioinformatics Institute (EBI) have published a catalogue of GWAS that provides a publicly available manually curated collection of published GWAS assaying at least 100000 SNPs 83 . State of the Art and all SNP-trait associations with p≤1×10−5. Up to December 2013, this catalogueincluded almost2000curated publicationsformore than 12000 SNPs [150]. A screen-shot of this interactive catalogue grouping 17 GWA trait categories is shown in Figure 5.8. Figure 5.8: Published GWAs at p≤5×10−8for 17 trait categories. Published Genome-Wide Associations, adapted from [151]. 5.4.1. GWA Testing using Imputed Data The probabilistic nature of imputed SNPs means that association analysis for these SNPs requires some care, regardless of what software or reference sets are used to generate the imputed data [85, 152]. While genotype platforms usually produce exact genotype calls (i.e. each individual is assigned genotype AA, AB, or BB –coded as 0, 1 or 2–), imputation programs generate probabilities for each of the three possible genotypes [152]. Using just those imputed genotypes with posterior probability above some threshold (or using the best guess genotype) is a reasonable method of comparing the accuracy across methods, but it is not recommended when carrying out association tests at imputed SNPs. Removing genotypes in this way can lead to both false positives and loss 84 5.4. Analysis of Genome Variation of power [85]. Luckily, many GWA packages can analyze the genotype dosages (the expected number of copies of a specified allele, from 0 to 2) that are produced by imputation programs [152]. 5.4.1.1. GWA Testing programs Several GWA testing tools can be applied to conduct analysis of imputed SNPs using the corresponding posterior probabilities. These can provide additional insights beyond what is provided by testing on typed tagging SNPs only. For this reason, numerous stand-alone packages as BEAGLE [125],BIMBAM [113],MACH2dat [122],ProbABEL [153], PLINK [130],SNPMStat [131] and SNPTEST [154] have been proposed. Despite the amount of GWAS performed using imputed data, only a few reviews [85, 155] carry out a comprehensive comparison of specific GWA testing programs for imputed data. Since the performance of these programs are affected by a variety of genetic factors, Pei et al. [155] investigated through a comprehensive comparison of these methods the effects, for example, of LD, Minor Allele Frequency (MAF) of untyped causal SNPs, and imputation accuracy rate. Hence, as part of these results, Table 5.5 and Table 5.6 were generated. Table 5.5 shows Type-I error rates of various imputation-based association methods for the causal SNP at the significant level of 5%. Quantitative Trait Qualitative Trait Low LD Mid LD High LD Low LD Mid LD High LD SNPTEST 5.0 5.0 4.8 5.0 5.0 4.8 SNPTEST-BG 5.0 5.0 5.1 5.0 5.1 5.1 MACH2qtl/dat 5.0 5.0 5.0 4.9 4.9 5.0 BIMBAM 5.0 4.8 5.0 5.0 4.8 5.0 BEAGLE - - - 4.9 5.1 4.9 PLINK - - - 4.4 4.3 4.4 ProbABEL 5.1 5.1 5.0 5.1 5.2 5.1 SNPMStat - - - 7.0 6.0 4.7 Table 5.5: GWA imputation-based tools comparison: type-I error rates [155]. From Table 5.5, it can be seen that, when testing association at the imputed potential causal SNP, all tools had type-I error rates close 85 . State of the Art to 5%. When testing for an entire genomic region, all programs but SNPMStat continue to have reasonable error rates, whereas SNPMStat had an inflated type-I error rate under low LD level. However, when testing for the whole region under high LD level, all methods were conservative. On the other hand, Table 5.6 depicts the accuracy of the GWA imputation-based methods for testing a whole genomic region under the significant level of 5%. Quantitative Trait Qualitative Trait Low LD Mid LD High LD Low LD Mid LD High LD SNPTEST-BG 49.6 52.8 65.5 50.1 52.3 62.2 SNPTEST 49.9 54.2 67.2 50.1 53.2 62.4 MACH2qtl/2dat 50.4 54.0 66.3 50.0 53.7 62.7 BEAGLE - - - 50.3 51.3 61.9 PLINK - - - 50.3 51.2 56.9 ProbABEL 50.7 53.9 66.6 50.1 52.7 62.9 SNPMStat - - - 50.5 51.2 59.0 Table 5.6: GWA imputation-based tools comparison: accuracy [155]. For both quantitative and qualitative traits, MACH2qtl/dat,ProbABEL and SNPTEST had the best performance under most situations, followed by SNPTEST-BG.BEAGLE program had similar performance to SNPTEST-BG under high LD level, but was inferior under medium LD level. SNPMStat and PLINK had the lowest power. As BIMBAM estimated p-value through permutation with 1000 replicates, its output had a resolution 1×10−3, which did not reach the significant level (2×10−4) with Bonferroni correction. For this reason, BIMBAM was not included in the analysis. Note that for all developed analysis, SNPTEST-BG utilized the best-guess genotype method, while SNPTEST considered the uncertainty and took the posterior probability into analyses [155]. 86 iii PART Approaches for Population Data This part of the thesis is motivated by the problem of identification of genetic variants associated with advanced diabetic nephropathy in a T2D population from the Gran Canaria island. For this purpose, we analysed a dataset from a case-control study where more than 2.5 millionof SNPs weregenotyped. Theobjective was to identify which SNPs could be associated with the time since diagnosis of diabetic nephropathy. Both cases and controls were individuals with several years of T2D progress. Cases were selected after nephropathy development and controls were individuals of similar characteristics (age, sex, history of diabetes, etc.) but not affected by nephropathy. Once these SNPs were identified, a new sample of subjects will be genotyped in those regions where such SNPs were located to confirm the association and identify relatedgenestothedisease. Thispart will describediverseprocessesthat have to be carried out to this aim: quality control, alignment and phasing, imputation, GWAS testing, as well as it will be explained obtained results and selection of candidate SNPs. These tasks have required an extensive use of complex bioinformatics resources. In Appendix A, we have included a tutorial that synthesizes the different steps of the analysis which we expect that can be useful to others facing the same problems. 87 Chapter 6 Quality Control 6.1. Introduction Quality control (QC) is a critical stage previous to applying association tests to SNP data. Therefore, the need for performing QC lies on several reasons. Firstly, since hundreds of thousands or million of genotypes may be properly generated, sometimes genotyping errors appear (in a small proportion), which if unidentified, may lead to spurious GWAS results [156] or in false positive association errors. Secondly, the quality of genotype calling usually is lower in practice than in the manufacturer’s panels. This is due to real samples normally are not such carefully prepared as the ones used to benchmark the panels. Finally, some developments such as custom SNPs panels, separately-genotyped reference control cohorts, and combined analyses of separate studies, may all increase the need for a QC phase [157]. In summary, QC can be differentiated by two aspects: issues related togenotypecalls (fromgenotypingchips) and downstreamissues [156]. This part of the dissertation will exclusively focus the last (i.e. those procedures to be applied once genotype calling is already performed), covering both subject-based and variant-based quality measures. These measures will include QC on individuals (samples) and variants (markers). A flowchart overview of the entire QC process is depicted in Figure 6.1. Each QC aspect is detailed along the corresponding section of this chapter. 89 . Quality Control To deepen in the sample missingness analysis, a threshold has to be specified through a balance between the maximization of genotyping efficiency and the minimization of the number of samples to remove. Thus, a sensible limit has to be placed where there is a qualitative change in data loss (as is depicted in Figure 6.4). The *.imiss file, generated by the –missing option of PLINK [130] (as mentioned in section 6.2.1), shows the missingness rate for each individual sample as follows: FID IID MISS_PHENO N_MISS N_GENO F_MISS 4 4 N 64370 2389506 0.02694 12 12 N 25846 2389506 0.01082 13 13 N 23048 2391739 0.009637 14 14 N 33146 2391739 0.01386 17 17 N 116259 2391739 0.04861 22 22 N 29320 2391739 0.01226 24 24 N 97750 2389506 0.04091 . . Where the column names are: FID (family identifier), IID (individual identifier), MISS_PHENO (missing phenotype?), N_MISS (number of missing SNPs), N_GENO (number of non-obligatory missing genotypes) and F_MISS (proportion of missing SNPs). Figure 6.4 shows the cumulative non-missingness distribution, in order to investigate the proportion of failed SNPs per sample, according to the F_MISS values of the *.imiss file. Based on Figure 6.4, from a call rate >91%, all individuals should have good genotyping quality. In the case of being more stringent (regarding call rate per sample), we could choose a threshold >97% using the PLINK [130] option –mind and a value of 0.03. In that case, we should eliminate 11 of 110 samples, i.e. the 10% of the total number of individuals in our study, which we thought highly not inadvisable. 6.3.2. Gender Mismatches This QC step is performed to check that the gender of individuals (samples) matches with the corresponding number of X chromosomes. Those subjects where the X-chromosome data disagrees with the reported gender are defined as problematic subjects. So, a PROBLEM arises if the two sexes do not match, or if the SNP data or pedigree data are ambiguous with regard to sex. A male call is made if F is more than 0.8 and a female call is made if F is less than 0.2 [130]. 96 6.3. Sample Quality Measures Figure 6.4: Individual Call Rate cumulative distribution. The–check-sex optionof PLINK [130]will generatea *.sexcheck file containing the gender mismatches as follows (except the Explanation column): FID IID PEDSEX SNPSEX STATUS F Explanation 35 35 1 1 OK 0.99 Male 41 41 1 1 OK 0.993 Male 45 45 2 2 OK 0.03151 Female 53 53 1 1 OK 0.9905 Male 55 55 1 1 OK 0.9938 Male 57 57 1 1 OK 0.9892 Male 60 60 2 0 PROBLEM 0.4347 Likely a female 61 61 2 2 OK 0.004779 Female . . Where field names are: FID (family identifier), IID (individual identifier), PEDSEX (sex as determined in pedigree file: 1=male, 2=female and 0=unknown), SNPSEX (sex as predicted based on genetic data -stored in X chromosome-), STATUS (displays ”PROBLEM” or ”OK” for each individual) and F(the actual X chromosome inbreeding -homozygosityestimate). 97 . Quality Control If there are loads of mismatches, it can be assumed that all Sample Identifiers could have become scrambled in some way, but in our case, there is only one mismatch, so we have assumed that the majority of sample labels were allocated correctly between our clinical and genetic data. 6.3.3. Population Stratification Population outliers (also known as ethnic outliers) occur when the study samples comprise multiple groups of individuals who differ systematically in both genetic ancestry and phenotype. Spurious apparent associations in admixed populations may be due to differences in ancestry rather than a true association of alleles to disease, leading to both false positives or false negatives [62]. So, although this is an important step of QC for population data, in our study it was not necessary to take this into account since the subjects were collected from the same population. 6.3.4. Individual Relatedness Individual relatedness occurs when pairs or groups of subjects are more closely related to each other than the population average, thus indicating they are close family members [157]. Those individuals induce a correlation structure which may cause mistaken associations, i.e. they can introduce false positive or false negative results. Furthermore, this is particularly problematic if the number or degree of cryptic relatedness differs between cases and controls. The –genome command of PLINK [130] generates a *.genome file which can be used to estimate genetic relatedness for all pairs of samples in the dataset. The PI_HAT column of *.genome file stores the Identityby-Descent (IBD) estimates of the individuals. A PI_HAT value close to 1 indicates a sample duplicate or monozygotic twins, a value close to 0.5 indicates 1st degree relatives(full sibs, parent-offspring), a valueclose to 0.25 indicates 2nd degree relatives (half-sibs, uncle/aunt-nephew/niece, grandparent-grandchild), and a value close to 0.125 indicates 3rd degree relatives (cousins, etc.) [161]. The applied criterion in our analysis to decide which sample would be dropped was to remove the one with the greater proportion of missing SNP data (extracted from the *.imiss file, described in Section 6.3.1). So, according to Figure 6.6, two individual samples had to be removed. 98 6.3. Sample Quality Measures Figure 6.5: Example of more complex relatedness networks [161]. Figure 6.6: Relatedness networks in our study, where two cryptic relatednesses emerge. Next figure comprises the IBD estimates histograms before and after the Variant Quality Control (VQC) and before Sample Quality Control (SQC). Figure 6.7: Comparison between IBD estimates histograms. 6.3.5. Heterozygosity Rate The heterozygosity rate (H) is the proportion of heterozygous genotypes for a given individual [158]. This proportion is predictable from Hardy-Weinberg expectations and the MAF at each SNP according to the Wright’s Inbreeding Coefficient (F), which becomes one minus the ob99 . Quality Control served number of heterozygotes in a population divided by its expected number of heterozygotes at HWE, i.e.: F= 1 −O.HET E.HET (6.2) Positive F indicates an excess of homozygotes (low heterozygosity), negative F indicates an excess of heterozygotes (high heterozygosity) [157]. High heterozygosity can also indicate sample contamination (i.e., a mixture of two or more DNAs, leading to more apparent heterozygotes. Low heterozygosity can indicate membership in a different population (the Wahlund effect [162]) or indeed could indicate inbreeding. Thereby, deviation from expected heterozygosity may indicate either low genotyping quality (efficiency) or relatedness of individuals [157], which can justify their later removal. TheFparametercan be estimated withthe–het optionofPLINK [130] (using the standard MLE for one locus, and then using a method-ofmoments). This method is an unbiased procedure for combining information across loci that involve separate summations of the number of observed and expected homozygous genotypes at each locus [161]. Thereupon, if we represent the H and F values in the same figure we can compare if there are two individuals with unusually high F (indicating either a genotyping problem or that they come from a different population), and thus should be dropped. Figure 6.8: Histogram of Heterozygosity H, and an inversely related value F (before Sample-QC). 100 6.3. Sample Quality Measures The heterozygosity histograms for before and after the Sample-QC in Figure 6.8 and Figure 6.9, respectively. In both figures some outliers can be noted, exceeding in three times the Standard Deviation (SD), which is the usual limit used to analyse heterozygosity [158]. Although such H,F values were out of the limits, we decided not remove their corresponding subjects, since this would unbalance the number of casecontrol samples. Figure 6.9: Histogram of Heterozygosity H, and an inversely related value F (after Sample-QC). 101 Chapter 7 Alignment and Phasing 7.1. Introduction This chapter will detail the alignment and phasing of data once it has been passed the QC stage (previously explained in Chapter 6). The main purpose of aligning and phasing the study genotypes is to get later a faster imputation from a large reference panel of haplotypes such as HapMap [95] or 1000 Genomes [98] projects. The procedure to implement at this stage is composed of two steps: 1) alignment of study samples with the reference panel and 2) estimation of haplotypes from genotype data (also known as phasing). 7.1.1. Data Format Resulting data of QC stage have to be divided as many parts as chromosomes, to process later each part in a manageable file by the program SHAPEIT2 [141] (which will be used to perform both alignment and phasing). Data formats used by SHAPEIT2 [141] are quite different, and they may be classified into three types: Input Data file (described in subsection 6.1.1 of Chapter 6), Genetic Map file (genetic map per chromosome) and Reference Panel files (haplotype text file format used by IMPUTE2 [121]). 103 . Alignment and Phasing 7.1.1.1. Genetic Map File Description Each genetic map file (per chromosome) should have the following structure: Position Combined_rate Genetic_Map 10906723 0.6231151531 0.025416867094949 10906915 0.4976528486 0.0255124164418802 10906989 0.4965798052 0.025549163347465 10907208 0.4953800875 0.0256576515866275 10913973 0.442217683 0.0286492542121225 10916916 0.4231906694 0.0298947043521667 . . Where the column names are: Position (physical position in base pairs ’bp’), Combined_rate (recombination rate in centiMorgans per Megabase ’cM/Mb’), and Genetic_Map (genetic position in centiMorgans ’cM’). 7.1.1.2. Reference Panel Files Description The reference panel files used by SHAPEIT2 [141] have the default text file format that can handle IMPUTE2 [121], namely the *.sample, *.legend and *.hap formats. To describe them shortly, we can consider four unrelated individuals from EUR and AMR groups, specifically from three different populations, for which two haplotypes per individual are available at three markers (SNP1,SNP2, and SNP3): Haplotype 1 Haplotype 2 Sample Group Popul SNP_1 SNP_2 SNP_3 SNP_1 SNP_2 SNP_3 Indiv1 EUR IBS A T A A C T Indiv2 EUR IBS G C T A T A Indiv3 EUR TSI A T T A C T Indiv4 AMR MXN G T T G C T . . From the previous description, we can construct a very large file per chromosome, but this would be inefficient and computationally costly. For this reason, the adaptation into SAMPLE/LEGEND/HAP files is necessary, since SHAPEIT2 [141] program need to reduce the dimensions of each file to decrease the computing time. So, a brief explanation of each format is the following: 104 7.1. Introduction ▶The SAMPLE file describes the non-genetic information per individual sample: Sample Population Group Sex Indiv1 IBS EUR male Indiv2 IBS EUR female Indiv3 TSI EUR female Indiv4 MXN AMR female . . Where the column names represent: sample (individual identifier), Population (population identifier), Group (group inside the population), Sex (1=male, 2=female, 0=unknown). ▶The LEGEND file describes the genetic information per single marker: id position a0 a1 SNP_1 10906723 A G SNP_2 10906915 T C SNP_3 10906989 A T . . Where the column names are: id (marker identifier), position (marker position), a0 (main allele) and a1 (alternate allele). ▶The HAP file contains the haplotypes of the reference panel in binary format, where 0 stands for the main allele and 1 represents the alternate allele: 00100011 01100101 01101111 . . Where each line corresponds to a SNP, storing the allele pairs for all individuals of the dataset. For example, the allele pair ”0_1” means that the first haplotype carries the main allele while the second carries the alternate allele. Haplotypes are given in the same order than in the *.sample file. 105 Chapter 8 Imputation 8.1. Introduction The main purpose of the imputation process is estimating unobserved genotypes in study samples by the extrapolation of the allelic correlations from a reference panel. This task must necessarily include a dataset composed of the study samples (genotyped at a subset of SNPs) and a reference panel (genotyped at a denser set of SNPs). It can be carried out by different programs, such as MaCH [122],BIMBAM [113], BEAGLE [108],MINIMAC3 [124] or IMPUTE2 [121]. After analyzing the state of the art (see Chapter 5), we have decided to use the last two, because of their good relationship performance/precision. In both cases, a list of the imputed alleles (with the corresponding probability distribution on them (in *.gen and *.vcf formats) can be obtained as output. In this chapter, data format and methods implemented by the corresponding imputation programs will be described. 8.1.1. Data format Dataformatsused throughout the imputation stagemay be classified into three types: Input Data file (described in subsection 6.1.1 of Chapter 6), Genetic Map file (genetic map per chromosome) and Reference Panel files (haplotype text file format used by IMPUTE2 [121]). 113 . Imputation 8.1.1.1. Imputed Files Description The default imputed files of IMPUTE2 [121] and MINIMAC3 [124] programs are *.gen format and the *.dose.vcf format, respectively. The latter is a modified version of the *.vcf format previously described in Section 7.1.1.3. A brief description of both is the following: ▶The GEN file contains the SNPs information together with the haplotypes of the phased dataset. In this format, the SNPs are clearly identifiable from the study samples (respect to those from the referencepanel). So, being thestudy SNPslabelledas SNP_01, SNP_02, SNP_03, etc. and the panel SNPs as rs0031, rs0032, rs0033, etc., we would get the following genotypes for 2 individuals: Panel SNPs Study SNPs rs0031: AA AA SNP_01: AA AA rs0032: GG GC SNP_02: GG GG rs0033: CC CC SNP_03: GT GC rs0034: CC TT SNP_04: CC GG . . .. . . The corresponding *.gen file would be: 21 rs0031 10979170 A G 1 0 0 1 0 0 21 rs0032 10979323 G C 1 0 0 0 1 0 21 SNP_01 10979896 T C 1 0 0 1 0 0 21 rs0033 10979913 C T 0.960 0.040 0 0.986 0.014 0 21 rs0034 10980199 C T 1 0 0 0 0 1 . . 21 SNP_02 14622336 G A 1 0 0 1 0 0 21 rs0191 14622684 C T 0 0.984 0.016 0 0.986 0.014 21 rs0192 14622862 G A 1 0 0 0 1 0 . . Where each line corresponds to a single SNP with the following columns: chromosome number,marker identifier,marker position, main allele,alternate allele and the rest of columns store the three probabilities for each individual corresponding to the genotypes AA, AB and BB respectively. ▶The default format of MINIMAC3 [124] is the VCF file although it can output files in both *.vcf and *.dose formats, being the last one the usual MINIMAC [105] output format. In order to facilitate the posterior statistical analysis, the probabilities of the im114 8.2. Imputation puted alleles by setting the –format command with the GT,GP labels, where each one means: ·GT: Estimated most likely genotype. ·GP: Estimated posterior genotype probabilities. Then, thisformat isthesamethattheonedescribed inSection 7.1.1.3, but including posterior genotype probabilities (GP) in the individual columns, namely those columns would look like as follows: . . FORMAT Ind1 Ind2 GT:GP 0|0:1.000 ,0.000 ,0.000 1|0:1.000 ,0.000 ,0.000 GT:GP 0|1:0.002 ,0.000 ,0.998 1|0:0.992 ,0.000 ,0.008 GT:GP 0|1:0.000 ,0.000 ,1.000 1|0:1.000 ,0.000 ,0.000 GT:GP 0|0:1.000 ,0.000 ,0.000 0|0:1.000 ,0.000 ,0.000 GT:GP 2|2:0.410 ,0.550 ,0.040 2|1:0.882 ,0.110 ,0.008 . . Where Ind1, Ind2 columns now include the three probabilities for each individual (one per genotype ’AA, AB or BB’). 8.2. Imputation Once the SNP alleles of the study samples have been phased (previously described in Chapter 7), the imputation stage starts. If a haplotype composed by kmarkers is identified in a study sample, normally it will be possible to find several haplotypes of length K > k in the reference panel that match with the study haplotype in the kconsidered makers. Genotypes present in the remaining K−kmarkers will be imputed to the study haplotype. Obviously, there is no single allocation, but it must be conducted in a probabilistic way, according to the observed frequencies in the reference panel in K−kimputed markers. Most HMM-based imputation methods infer the missing genotypes for each sample independently, depending on the reference panel. In this way, the reference data is used to integrate the phase uncertainty in each individual’s genotype. All imputation methods normally run Markov Chain algorithms, such asMarkovChain MonteCarlo (MCMC) or HMM, which are essentially based on this two steps: 1. Imputation of any sporadic missing genotype at the study SNPs, once all the observed genotypes have been phased. 2. Use of the reference panel to impute missing alleles at untyped SNPs. 115 . Imputation Imputation programs repeat the previous steps over multiple iterations, averaging the imputation probabilities inferred for each missing genotype and also integrating the uncertainty phase. 8.2.1. Imputation Method Description The imputation method proposed by Marchini et al. [112] employs a HMM network to compare the set of genotypes for each study sample to the reference panel in order to solve the untyped SNPs. This HMM method contains hidden states as haplotype backgrounds or SNPs loci (pointsintime). ThemovementbetweenthesestatesisaMarkovChain, being the probability of belonging to state j(where j={1..k}at time idepends on time i−1). The kstates are the kpossible haplotype backgrounds of the reference panel, being the probability of moving between states based on recombination rates. In a state jat point i, the emission probabilities of the HMM will determine what an outcome could be (an allele, in this case). A detailed explanation of how such imputation is performed starts with the description of the following variables: ◦H={H1, . . . , HN}is the reference haplotypes set, which is formed by Ndistinct haplotypes. ◦Lrepresent the number of loci placed in the SNPs of interest. ◦Hi= (Hi1, . . . , HiL)is ith haplotype, being Hij ∈ {0,1}. ◦Kare the number of individuals of the study. ◦Giis the genotype of the ith subject, where Gi= (Gi1, Gi2, . . . , GiL)with Gij ∈ {0,1,2, missing}. ◦G={G1, . . . , GK}represents the set of genotypes of Ksubjects in the study. For each subject, the corresponding genotype Giis distributed as {GiO, GiM }being GiO the observed genotyped and GiM the missing genotype (not observed in the study samples, but available in reference haplotypes). We can then consider GO={G1O, G2O, . . . , GKO}as the observed genotype set for Ksubjects and GM={G1M, G2M, . . . , GKM }as the missing genotype set. To carry out the GMimputation, we must first calculate the probability distribution of GM. Obviously, this will depend on the reference haplotypes and the observed genotypes per individual. Therefore, we have: 116 8.2. Imputation P(GM=gM|GO=gO, H ) = P(G1M=g1M, . . . , GKM =gKM |G1O=g1O,...GKO =gKO, H ) = K ∏ i=1 P(GiM =giM |GiO =giO, H )(8.1) Where we have supposed that the study subjects are independent each other. Thus, from the basic rules of conditional probability, we know that the joint probability is the product of the probabilities: P(A|B∩C) = P(A∩B∩C) P(B∩C)= =P(A∩B∩C) P(C)·P(C) P(B∩C)=P(A∩B|C) P(B|C)(8.2) Being P(A|B∩C)proportional to P(A∩B|C), which can be written as P(A|B∩C)∝P(A∩B|C). According to our case, this rule results: P(GM=gM|GO=gO, H )∝ K ∏ i=1 P(GiM =giM , GiO =giO |H) = K ∏ i=1 P(Gi=gi|H)(8.3) Now, for calculating the probabilities P(Gi=gi|H), a HMM model is used, where the hidden states are (Z(1) i, Z(2) i)haplotype pairs formed from the Hset, namely Z(j) i=(Z(j) i1, Z(j) i2, . . . , Z(j) iL )and Z(j) il ∈ {1, ..., N}. Furthermore, the Z(j) il value indicates from which haplotype comes the j(where j={1,2}) of llocus. Notethat previous definitionmeansthatZ(j) imaynotexactlymatch with an Hhaplotype, but consist of linked fragments from such haplotypes (namely, Z(j) iwould become a new haplotype that is obtained by recombining those found in H). The HMM model, therefore, implies that: 117 . Imputation P(Gi|H) = ∑ Z(1) i,Z(2) i P(GiZ(1) i, Z(2) i, H )·P(Z(1) i, Z(2) i|H)(8.4) For each locus, the hidden states can then be considered as the haplotype pairs of Hset which are copied at that locus to form later the Gi vector of genotypes. The P(Z(1) i, Z(2) i|H)term defines the prior probability of changing Z(1) iand Z(2) ialong the ith subject. Thus, starting from the first locus, the (Z(1) i1, Z(2) i1)initial state would have the probability: P(Z(1) i1, Z(2) i1|H)=1 N2(8.5) Where the transition probabilities would be: P({Z(1) il , Z(2) il }→{Z(1) i(l+1), Z(2) i(l+1)}|H)=                (e− ρl N+1−e−ρl N N)2 No recombination (e− ρl N+1−e−ρl N N)(1−e−ρl N N)Recombination in a haplotype (1−e−ρl N N)Recombination in both haplotypes (8.6) The ρterm in the previous expression depends on a estimation of the recombination rate throughout the genome. Such rates have been collected from HapMap project data in order to build a recombination map of the human genome [112]. Based on the initial state probability and such transition probabilities, the P(Z(1) i, Z(2) i|H)probability could calculated as: P(Z(1) i, Z(2) i|H)= P(Z(1) i1, Z(2) i1|H)L−1 ∏ i=1 P({Z(1) il , Z(2) il }→{Z(1) i(l+1), Z(2) i(l+1)}|H) (8.7) The P(GiZ(1) i, Z(2) i, H )term models how the observed genotypes are close, but not exactly the same as those that would result from 118 8.2. Imputation Hhaplotypes that are being copied. Somehow this expression mimics the effects of mutation. Assuming that mutations occur independently along loci, we have: P(GiZ(1) i, Z(2) i, H )= L ∏ l=1 P(Gil Z(1) il , Z(2) il , H )(8.8) The P(Gil Z(1) il , Z(2) il , H )probabilities depend on the mutation rate are further described in [112]. A graphical description of the typical set up for an imputation study isdepicted in Figure 8.1. The bluedata arestudy genotypes, and the aim is to estimate genotypes up to the density of the reference haplotypes (in black). In some studies, however, the aim may be to impute only up to the density of some reference genotypes (green data). It is possible to have only haplotype reference data, genotype reference data or both. Sporadic missing data (codedas’?’, can also be imputed both in reference data and inference data genotypes) [167]. Figure8.1: Schematicdrawingof a standard imputationscenario. Example of a standard imputation, adapted from [167]. 119 . Imputation 8.2.2. IMPUTE and MACH/MINIMAC programs The program IMPUTE [112] was created to determine the probability distribution of missing genotypes conditional upon a set of known haplotypes and an estimated fine-scale recombination map (calculating the historical mutation through a theoretical result from population genetics theory). Later, a new version called IMPUTE2 [121] was released. This handles larger population genetic datasets and may use additional reference panels (or study samples) typed on different chips. On the other hand, we have MACH/MINIMAC programs. The first of these set of programs is MACH [122], which imputes unobserved genotypes for case-control studies, employing a reference dataset and a HMM approach. The second set of algorithms, called MINIMAC, are more computationally efficient implementations of the MACH program, which were initially developed by [105], and later by [123] (MINIMAC2) and [124] (MINIMAC3). Both IMPUTE and MACH/MINIMAC programs have ”best guess” option for imputed genotypes (i.e. those with probability of 0.50 or greater), which is the simplest way to carry out the imputation. IMPUTE2 [121] performs a single imputation step rather than running MCMC in this situation, and the main change on the command line is to replace the –g file (unphased genotypes) with a –known_haps_g file. The program imputes the untyped alleles in each phased haplotype, to later convert the corresponding allelic probabilities into diploid genotype probabilities [168]. Otherwise, MACH/MINIMAC programs output a file (*.dose file) that provides an estimated dosage of each imputed genotype along with the dosage of observed genotypes. The dosage is a range between 0 and 2, where 0 represents no copies of the SNP reference allele, and 2 represents two copies of the reference allele. Therefore, while observed genotypes have discrete values included in [0 1 2], imputed genotypes are an estimate of the number of copies of the reference allele, represented by a decimal falling anywhere between 0 and 2. MACH [122] also outputs a file of posterior probabilities and a file with quality scoring for each SNP. The ”best-guess” genotypes can be analyzed just as observed genotypes, with any standard statistical genetics software package, caution should be used in interpretation of results due to the uncertainty of the imputed data. The dosage file is particularly well-suited for analysis methods that take this into account, such as logistic and linear regression [167]. 120 Chapter 9 GWA Testing of Imputed Data 9.1. Introduction This chapter will describe the final stage of Part iii of the dissertation: thestatisticalanalysis (i.e. associationtests) ofpreviously imputed data. During lastyears, imputed SNPshave beenfrequently usedin GWAS analyzes, although, for a properly implementation of association tests, the uncertainty of such imputed data has to be taken into account. We have chosen SNPTEST2 [154] as the program for analyzing single marker associations in GWAS, since it has implemented different methods for dealing with imputed SNPs, including the analysis of binary (case-control) phenotypes, single, and multiple quantitative phenotypes and Bayesian or Frequentist tests. 9.1.1. Data format Input data formats used by the SNPTEST2 are the same as described in Section 8.1.1.1 (namely SAMPLE,GEN and VCF files). The-summary_stats_only optionofSNPTEST2 calculatesthedata summaries associated with each SNP (i.e. genotype counts, allele frequencies, missing dataproportionsandodds ratios), leadingto an output file as follows: id rsid chr pos allele \_A allele \_B index average \_maximum\ _posterior \ _call info . . 21:10979170 rs0031 21 10979170 A G 1 1 1 121 . GWA Testing of Imputed Data This expression shows that when there is uncertainty in some SNP, the likelihood function depends on P(Gij =k|H, G, θ) terms form that correspond to the probability distribution of the unobserved genotypes (conditioned by the reference haplotypes and by the model parameter values). As in the case of known genotypes, there are several ways to test the association between SNP and disease, and in this context we can consider: a) Using the likelihood ratio test. In this case, it is necessary to find the maximum likelihood estimator of θfor which the Newton-Raphson method can be used identically as the above, although now considering the new U∗(θ)score and I∗(θ)information functions (either use an EM algorithm). More details are shown in the supplementary information of [85]. b) Using a Score test. As we have seen for the estimate with known genotypes, it is not necessary to estimate β1, since in this test it is assumed that the null hypothesis (β1= 0) is true. The parameter theta is thus reduced to (β0,0). As with known genotypes β0=log pi 1−pi=log N1/N 1−N1/N= log N N−N1=log N1 N0, so β0is a fixed value and therefore P(Gij =k|H, G, θ) = P(Gij =k|H, G) = pijk. The supplementary material of Marchini and Bryan [85] the Score test is explicitly described in this case. c) Using a Wald test. The statistic is SW ALD =θ2 I∗(ˆ θ), which is distributed according to χ2 1when H0is true. Further explanations of statistical theory for missing data problems can be found in [112] and [85]. 128 9.4. Results using HapMap as Reference Panel 9.4. Results using HapMap as Reference Panel This section will depict the resulting log-p-values of additive test using HapMap as the reference panel. We have chosen Manhattan and QQ plots in order to perform the most representative illustrations of data. Figure 9.1: Manhattan plot of resulting log-p-values of additive test using HapMap as reference panel and IMPUTE2 as imputation software. Figure 9.2: QQ plot of the resulting p-values of the additive test using HapMap as reference panel and IMPUTE2 as imputation software. 129 . GWA Testing of Imputed Data TheManhattanplotofFigure9.1 representsalllog-p-valuesthroughout the Y-axis, being possible a clearly distinction among the corresponding values for each chromosome. It can be noted numerous ”gaps” per chromosome, which are due to the low density of HapMap panel. On the other hand, the QQ plot of Figure 9.2 shows the log-pvalues, but differentiating between observed and expected, although with no distinction among chromosomes. It also can be noted that the observed log-p-values are quite similar to the expected ones. Therefore, we can assert that the imputation of study data with the HapMap reference panel do not lead to significant results. 9.5. Results using 1000 Genomes as Reference Panel This section will illustrate the resulting log-p-values of additive tests using 1000 Genomes as the reference panel. We have chosen Manhattan and QQ plots in order to perform the most representative illustrations of data. We have analysed the imputation results of both MINIMAC3 [124] and IMPUTE2 [121] programs in order to detect if imputation method could have influence the results. 9.5.1. Manhattan and QQ plots from MINIMAC3 imputation Inthissubsection, onceperformed theimputationwithMINIMAC3, the resulting log-p-values of additive test (using 1000 Genomes as reference panel) are depicted in Figures 9.3. Such figure shows more density than results from HapMap, leading to an improvement of the significant log-p-values in all tests contrasted. It can be noted that in certain genome regions these values surpass a threshold of 10−5. As well, it has been illustrated the observed and expected log-pvalues of the additive test (using 1000 Genomes as reference panel) in Figure 9.4. Such QQ plot illustrates that a greater quantity of expected logp-values differ from the observed ones (in comparison with depicted values in QQ plot of Section 9.4). From a threshold of 10−5, it can be noted how numerous observed log-p-values walk away from the expected ones. 130 9.5. Results using 1000 Genomes as Reference Panel Figure 9.3: Manhattan plot of resulting log-p-values of the additive test using 1000 Genomes as reference panel and MINIMAC3 as imputation software. Figure 9.4: QQ plot of the resulting p-values of the additive test using 1000G as reference panel and MINIMAC3 as imputation software. 131 . GWA Testing of Imputed Data 9.5.2. Manhattan and QQ plots from IMPUTE2 imputation In this subsection, once performed the imputation with IMPUTE2, the resulting log-p-values of the additive test (using 1000 Genomes as reference panel) is depicted in Figure 9.5. Figure 9.5: Manhattan plot of resulting log-p-values of the additive test using 1000G as reference panel and IMPUTE2 as imputation software. Figure 9.6: QQ plot of the resulting p-values of the additive test using 1000G as reference panel and IMPUTE2 as imputation software. 132 9.6. Selection of Significant Results Figure 9.5 shows more density of values than results using HapMap as reference panel (see Section 9.4), leading to even more significant logp-values in comparison with MINIMAC3 results. The observed and expected log-p-values of the additive test (using 1000 Genomes as reference panel) has also been represented in Figure 9.6. In such QQ plot depicts how observed log-p-values move away of expected values from lower threshold than in Section 9.4. So, the total amount of significant values (comparing them from a given threshold) is greater using IMPUTE2 as imputation software than with MINIMAC3 program. 9.6. Selection of Significant Results The corresponding association analysis was performed using all frequentist tests (i.e. additive, dominant, general, recessive and heterozygous). We ultimately chose the additive results with and without adjusting by the ”retinopathy” co-variable (i.e., if subjects were affected by the retinopathy disease or not). To filter and select the most significant results, we have decided to group them in regions containing, at least, three significant SNPs (with p-values smaller than 10−5). Furthermore, we have selected those INFO values greater than 0.99. In Tables 9.5 and 9.6 above results’ selection are shown. In both tables, each line stores the information of a genomic region whereChr column isthe chromosome, Regis the selectedregion, rsID is the original SNP name, rsIDnew is the updated SNP name (by the Build 142 of the SNP locations for Homo sapiens), Position is the genomic position where each SNP is located (in base-pairs), alA is the main allele, alB is the alternate allele, Pvalue is the resulting p-value, testINFO is the imputation quality, Beta is the coefficient β, SE is the standard error, GENEsymbol is the symbol of the gene to which the SNPs belong. From Table 9.5 it can be seen that the three lowest p-values are in regions 8.1, 10.2 and 11.1, as well in Table 9.6 but with the difference SNP that, in such table the p-values adjusted by the retinopathy co-variable are better. From each region we selected one SNP, being the chosen: rs4841106, rs2358658, and rs35649357. Figures 9.7, 9.9 and 9.11 show the typed SNPs. Representations of the same SNPs but with imputed information are arranged in Figures 9.8, 9.10 and 9.12. 133 . GWA Testing of Imputed Data Chr Reg Typed rsID rsIDnew Position alA alB Pvalue testINFO Beta SE GENEsymbol 1 1.2 No rs3850619 rs3850619 199109114 G C 5.4161×10−60.994034 -1.37199 0.336038 1 1.2 Yes rs10800611 rs10800611 199115592 T C 5.24046×10−61 -1.35607 0.331267 3 3.1 No rs73060069 rs73060069 28948053 G A 5.29236×10−60.998066 -2.05171 0.53348 3 3.1 Yes rs11710772 rs11710772 28960287 C T 5.29709×10−61 -2.0502 0.53266 4 4.2 No rs1490674 rs1490674 155330282 T C 4.32929×10−60.99295 1.48872 0.361206 DCHS2 4 4.2 No rs990185 rs990185 155330882 G C 4.19076×10−60.994467 1.4906 0.361129 DCHS2 4 4.2 No rs6818438 rs6818438 155332511 G A 3.93608×10−60.997257 1.49442 0.361039 DCHS2 6 6.1 No rs78261651 rs78261651 85799282 A C 3.54827×10−60.999626 -8.4999 16.1244 6 6.1 Yes kgp4641172 rs78335147 85802481 C A 3.55061×10−61 -8.50257 16.1523 6 6.1 Yes kgp4434658 rs55897238 85804399 C T 3.55061×10−61 -8.50257 16.1523 8 8.1 Yes rs4841106 rs4841106 9044856 A G 1.49343×10−61 -1.53146 0.360835 LOC101929128 8 8.1 No rs4841107 rs4841107 9044897 C G 1.31739×10−60.997214 -1.54226 0.361905 LOC101929128 10 10.2 No rs2884507 rs2884507 20305489 A T 1.25442×10−60.997256 1.54494 0.361962 PLXDC2 10 10.2 Yes kgp517295 rs2358658 20307083 G A 4.57521×10−61 1.46649 0.358001 PLXDC2 11 11.1 Yes kgp8554447 rs35649357 49552399 T G 1.44891×10−61 -8.28856 14.3315 18 18.1 No rs1784763 rs1784763 9067924 G A 5.46945×10−60.999983 8.33617 15.4622 19 19.1 No rs4254447 rs4254447 24477702 T C 4.17053×10−61 1.89291 0.473502 19 19.1 No rs10414165 rs10414165 24478348 T C 4.17053×10−61 1.89291 0.473502 Table 9.5: Most significant SNPs from additive test (without adjusting by any co-variable). 134 9.6. Selection of Significant Results Chr Reg Typed rsID rsIDnew Position alA alB Pvalue testINFO Beta SE GENEsymbol 3 3.1 No rs73060069 rs73060069 28948053 G A 3.67826×10−60.997356 -2.12155 0.545348 3 3.1 Yes rs11710772 rs11710772 28960287 C T 3.68872×10−61 -2.11937 0.544212 8 8.1 Yes rs4841106 rs4841106 9044856 A G 1.36741×10−61 -1.54889 0.364094 LOC101929128 8 8.1 No rs4841107 rs4841107 9044897 C G 1.42376×10−60.996237 -1.54557 0.363528 LOC101929128 8 8.3 Yes rs17289507 rs17289507 131493326 C T 3.93318×10−61 -8.65919 16.7397 9 9.1 No rs9696078 rs9696078 92503499 G A 5.99158×10−60.99762 1.3629 0.335209 10 10.1 No rs2884507 rs2884507 20305489 A T 1.00068×10−60.998724 1.56999 0.36486 PLXDC2 10 10.1 Yes kgp517295 rs2358658 20307083 G A 3.2376×10−61 1.50004 0.360917 PLXDC2 11 11.1 Yes kgp8554447 rs35649357 49552399 T G 1.30259×10−61 -8.33769 14.2969 18 18.1 No rs1784763 rs1784763 9067924 G A 2.74509×10−60.999994 8.49673 15.3799 19 19.1 No rs62116441 rs62116441 24189406 C T 4.1583×10−60.993668 1.92528 0.470914 19 19.2 No rs4254447 rs4254447 24477702 T C 4.96705×10−61 1.87717 0.473332 19 19.2 No rs10414165 rs10414165 24478348 T C 4.96705×10−61 1.87717 0.473332 19 19.3 Yes rs7259375 rs7259375 43885085 T C 4.78323×10−61 1.61444 0.392688 19 19.3 No rs8109116 rs8109116 43901799 G C 2.18859×10−61 -1.68328 0.397624 TEX101 19 19.3 Yes rs8109124 rs8109124 43901831 G A 2.18859×10−61 -1.68328 0.397624 TEX101 Table 9.6: Most significant SNPs from additive test (adjusted by the retinopathy co-variable). 135 . GWA Testing of Imputed Data 0 2 4 6 8 −log10(p−value) 0 20 40 60 80 100 Recombination rate (cM/Mb) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ●● ● ● ● ●● ●● ● ● ● ● ● ● ●●● ●● ● ● ● ● ● ●●● ● ●●●● ● ● ● ● ● ● ● ●● ● ● rs4841106 0.2 0.4 0.6 0.8 r2 PPP1R3B 9 9.02 9.04 9.06 9.08 Position on chr8 (Mb) Typed Figure9.7: First regionof Typed SNPs, wherethers4841106 SNP is located at the center of the image. 0 2 4 6 8 −log10(p−value) 0 20 40 60 80 100 Recombination rate (cM/Mb) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ●● ● ● ●● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● rs4841105 0.2 0.4 0.6 0.8 r2 PPP1R3B 9 9.02 9.04 9.06 9.08 Position on chr8 (Mb) P1000G Typed Imputed Figure 9.8: First region of typed and imputed SNPs, where the rs4841106 SNP is located at the center of the image. 136 9.6. Selection of Significant Results 0 2 4 6 8 −log10(p−value) 0 20 40 60 80 100 Recombination rate (cM/Mb) ● ● ● ● ● ● ● ● ●● ● ●● ● ● ●●● ● ● ●● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ●● ● ● ●● ● ● ● ● ●●● ●● ● ●●● ● ●● rs2358658 0.2 0.4 0.6 0.8 r2 PLXDC2 20.26 20.28 20.3 20.32 20.34 Position on chr10 (Mb) Typed Figure 9.9: Second region of Typed SNPs, where the rs2358658 SNP is located at the center of the image. 0 2 4 6 8 −log10(p−value) 0 20 40 60 80 100 Recombination rate (cM/Mb) ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●● ● ● ● ●● ● ● ●●●●● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●●●● ● ● ● ● ● ● ● ● ● rs2358657 0.2 0.4 0.6 0.8 r2 PLXDC2 20.26 20.28 20.3 20.32 20.34 Position on chr10 (Mb) P1000G Typed Imputed Figure 9.10: Second region of typed and imputed SNPs, where the rs2358658 SNP is located at the center of the image. 137 B. alleHap Manual Mk1_1 Mk1_2 Mk2_1 Mk2_2 Mk3_1 Mk3_2 Mk4_1 Mk4_2 Mk5_1 Mk5_2 Mk6_1 Mk6_2 CCAGTTCTTTGT CCAGTTCTTTGT C C <NA> <NA> T T C C T T <NA> <NA> CCAGTTCTTTGT TCGACTCCCCTT TCGGCTTCTCTG TCGGCTCTCTT<NA> <NA> <NA> G G <NA> <NA> C C C C T <NA> > # Reconstructed haplotypes > fams2List[’haplotypes’] famID indID patID matID sex phen hap1 hap2 FAM01 1 0 0 1 2 C?TCT? C?TTT? FAM01 2 0 0 2 1 C?TCT? C?TTT? FAM01 3 1 2 2 1 C?TCT? C?TCT? FAM01 4 1 2 1 1 C?T?T? C?T?T? FAM02 1 0 0 1 1 ?G?CCT ?A?CCT FAM02 2 0 0 2 1 ?G?TT? ?G?CC? FAM02 3 1 2 2 1 ?G?CCT ?G?TT? FAM02 4 1 2 2 2 ?G?CCT ?G?CC? Example B.3.7 Haplotype reconstruction of a family containing parental and offspring missing data from a PED file. > ## PED file path > family3path <- file.path(find.package(”alleHap”),”examples”,”example3 .ped”) > ## Loading of the ped file placed in previous path > family3Alls <- alleLoader(family3path,dataSummary=FALSE) > ## Haplotype reconstruction of previous loaded data > family3List <- alleHaplotyper(family3Alls,dataSummary=FALSE) > # Original data > family3Alls famID indID patID matID sex phen Mk1_1 Mk1_2 Mk2_1 Mk2_2 Mk3_1 Mk3_2 1 1 0 0 2 1 C T <NA> <NA> <NA> <NA> 1 2 0 0 2 1 <NA> <NA> <NA> <NA> <NA> <NA> 131212CCAGAT 1 4 1 2 1 2 C T A C <NA> <NA> 151212CTAGCT 161212CTAGCT 1 7 1 2 2 1 <NA> <NA> <NA> <NA> C G 240 B.3. Workflow Mk4_1 Mk4_2 Mk5_1 Mk5_2 Mk6_1 Mk6_2 Mk7_1 Mk7_2 Mk8_1 Mk8_2 <NA> <NA> <NA> <NA> A C <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> AAGTACACAG ATCGCCCTCG <NA> <NA> G T A C A C A A AAGTACACAA A T C G <NA> <NA> C T <NA> <NA> > # Re-imputed alleles > family3List[’reImputedAlls’] famID indID patID matID sex phen Mk1_1 Mk1_2 Mk2_1 Mk2_2 Mk3_1 Mk3_2 110021CTGCTG 120021CTAAAC 131212CCGATA 141212TCCAGA 151212CTGATC 161212CTGATC 171221TTCAGC Mk4_1 Mk4_2 Mk5_1 Mk5_2 Mk6_1 Mk6_2 Mk7_1 Mk7_2 Mk8_1 Mk8_2 ATTCACATAC AAGGCCCCGA AATGACACAG TACGCCTCCG AATGACACAA AATGACACAA TACGCCTCCA > # Reconstructed haplotypes > family3List[’haplotypes’] famID indID patID matID sex phen hap1 hap2 1 1 0 0 2 1 CGTATAAA TCGTCCTC 1 2 0 0 2 1 CAAAGCCG TACAGCCA 1 3 1 2 1 2 CGTATAAA CAAAGCCG 1 4 1 2 1 2 TCGTCCTC CAAAGCCG 1 5 1 2 1 2 CGTATAAA TACAGCCA 1 6 1 2 1 2 CGTATAAA TACAGCCA 1 7 1 2 2 1 TCGTCCTC TACAGCCA 241 vii PART Resumen en Español 243 Capítulo 14 Introducción La motivación inicial para la realización de esta tesis doctoral surgió de la necesidad de organizar y aplicar diferentes métodos bioestadísticos y/o bioinformáticos para la resolución de varios problemas planteados por el grupo de investigación en diabetes y endocrinología del Complejo Hospitalario Universitario Insular Materno Infantil de Las Palmas de Gran Canaria. Los problemas en cuestión estaban relacionados con la genética de la diabetes, para lo que se debían manejar datos genéticos de diferente naturaleza. Para algunas cuestiones debíamos utilizar datos de un estudio de múltiples familias (con datos genéticos de padres e hijos para un varios miles de familias). Otros de los problemas a estudiar requerían utilizar datos genéticos de un estudio de casos y controles. El proceso de análisis de estos conjuntos de datos requirió no sólo del conocimiento de métodos estadísticos e informáticos, sino también de fundamentos de la genética humana y de metodologías desarrolladas recientemente para el tratamiento avanzado de bases de datos genéticas. La realización de estudios de asociación estadística entre datos genéticos/genómicos y la presencia de enfermedades (o fenotipos) en los individuos de una población plantea retos importantes. Por una parte los derivados el enorme número de datos que se deben manejar (fácilmente varios millones de variables, correspondientes a distintos marcadores observados a lo largo del genoma); por otro lado, los derivados de las propias técnicas de obtención de la información genética, que obligan a realizar controles de calidad sobre los datos, y a tener en cuenta características subyacentes a la composición de las poblaciones y su acervo genético. Además en muchos casos la detección de una asociación en245 . Introducción tre un marcador y una enfermedad (con características hereditarias) no implica necesariamente que dicho marcador esté involucrado en el proceso causal de la enfermedad; puede ser simplemente que ese marcador se encuentre en cierta posición del genoma ligado, por proximidad y por la fortaleza de los enlaces químicos, al verdadero gen causal. Otro problema frecuente es la pérdida de información en el proceso de genotipado. Este es un proceso químico que no está exento de fallos, por lo que es habitual que la base de datos que reúne la información genética de los individuos presente muchos valores perdidos. Los valores perdidos pueden, en muchos casos, determinarse mediante el adecuado uso de técnicas de imputación, que acuden a paneles de referencia públicamente disponibles –en la práctica auténticos mapas del contenido de nuestro genoma–, para determinar cuáles son los alelos faltantes en los marcadores con valores perdidos. Asimismo, aún cuando se hayan observado millones de marcadores, es posible que la componente genética que se asocia a la enfermedad esté en alguna sección no directamente observada del genoma. Los métodos de imputación pueden utilizarse también para rellenar (imputar) las secciones del genoma del individuo que no han sido directamente observadas. Estos valores imputados pueden utilizarse también para evaluar posibles asociaciones con la enfermedad. Citemos, por último, que en muchas ocasiones la relación con la enfermedad tiene que ver no con marcadores más o menos dispersos a lo largo del genoma, sino con los haplotipos de los que dichos marcadores forman parte. Es por ello que resulta de gran interés práctico disponer de herramientas que permitan determinar los haplotipos que dan lugar a los genotipos observados. Losmétodos(algoritmos) que se utilizanhabitualmenteparalaimputación de genotipos o la reconstrucción de haplotipos usan técnicas probabilistas para inferir el contenido de los marcadores perdidos en las bases de datos genéticas. Cuando los individuos del estudio son independientes entre sí (no comparten herencia genética más allá de la que corresponde a pertenecer a la misma especie) la única fuente de información para la imputación son los paneles de referencia, y cualquier imputación que se haga será probabilística por naturaleza: en dicho panel puede haber varias secciones distintas compatibles con la sección concreta del genoma que se quiere imputar, por lo que el valor imputado se elige como el más probable entre los posibles. Sin embargo, en el caso particular de las bases de datos familiares, cuando varios miem246 bros de la misma familia están disponibles para un estudio 1es posible combinar la información disponible en padres, madres e hijos para determinar a partir de unos el contenido faltante en otros, teniendo en cuenta que la constitución genética de los hijos ha sido necesariamente heredada de los padres; de esta forma en estos casos la imputación de muchos valores perdidos (incluso todos), así como la identificación de los haplotipos compartidos, puede realizarse de manera determinista, disminuyendo las falsas asociaciones que pueden derivarse de los métodos de imputación probabilista. Ciertas regiones genómicas son muy estables frente a la recombinación, al mismo tiempo que pueden llegar a ser altamente polimórficas, presentando una gran variabilidad. Una de esas regiones, que ha sido bastante bien estudiada, es la denominada Complejo Mayor de Histocompatibilidad (MHC). En humanos, los genes MHC conforman los denominados antígenos leucocitarios humanos (sistema HLA)[3]. Específicamente para esta región, debido a su alta tasa polimórfica, ha surgido incluso la necesidad de la creación una nomenclatura alfanumérica para facilitar los posteriores análisis. Muchos de los métodos actuales por lo general no son capaces de identificar correctamente los haplotipos en esta región, dado que la alta variabilidad hace que los métodos probabilistas que han de utilizarse cuando no se dispone de datos familiares tengan una alta tasa de error [5]. Cuando hay datos familiares disponibles la identificación de haplotipos (y la imputación de valores perdidos) es más precisa, si bien son necesarios procedimientos eficientes para ello. Por último, cabría destacar también que aunque existe un número elevadoycrecientedeherramientasbioestadísticas/bioinformáticas para el procesamiento y análisis de bases de datos genómicas, la documentación relacionada con cada una de ellas es a menudo un tanto confusa y se encuentra dispersa. Particularmente para aquellos usuarios que quieren utilizar dichas herramientas por primera vez, suele resultar tediosa y complicada la puesta a puesto de los algoritmos necesarios para analizar los datos. Por lo tanto, existe una necesidad de clarificación y simplificación de la documentación relativa a los procesos que incluyen el control de calidad, la imputación de valores perdidos y/o el análisis estadístico de asociación en las bases de datos genómicas. 1normalmente duos formados por el padre o la madre y un hijo, tríos, con los dos padre y un hijo, o familias con dos hijos. Los estudios con más de dos hijos por familia son poco comunes 247 Capítulo 15 Objetivos Los principales objetivos de esta tesis han sido: 1) Estudio exhaustivo del estado actual de los métodos bioestadísticos/bioinformáticos que se aplican en el manejo de bases de datos genéticas/genómicas. 2) Aprendizajedelastécnicas y metodologíasnecesarias paraelanálisis de datos genómicos, especificamente para las tareas de control de calidad, imputación, haplotipado, y análisis de asociación con datos poblacionales. 3) Identificación de nuevas variantes genéticas asociadas a la nefropatía diabética avanzada en una población de pacientes con diabetes tipo 2 de la isla de Gran Canaria. 4) Desarrollo de un paquete informático en el lenguaje R para la imputación de genotipos perdidos y/o reconstrucción de haplotipos en bases de datos familiares. 5) Identificación deposibles asociacioneshaplotípicas para datos familiares, concretamente para los sujetos de la base de datos internacional T1DGC. 249 . Planteamiento y Metodología Pasa los datos de cada familia secuencialmente a la función famImputer función,queseencargaderealizarlaimputación marcador a marcador llamando repetidamente a la función mkrImpputer, Devuelve un conjunto de datos con el mismo formato y dimensiones que el conjunto de datos de entrada, con los valores imputados en aquellos alelos donde ha sido posible llevar a cabo la imputación. Paso 3: Opcionalmente, se mostrará un breve resumen del proceso de imputación: número de alelos imputados, incidencias detectadas (número de marcadores cancelados debido a los problemas detectados en el proceso de control de calidad), tasa de imputación (cociente de los alelos imputados sobre el número inicial de alelos perdidos) desaparecidos originalmente) y tiempo de proceso. Paso 4: Si los datos de la familia han sido leídos desde un archivo externo, se crea un nuevo archivo con el mismo nombre, pero con extensión imputed.ped, que contiene el conjunto de datos devuelto por famsImputer en formato PED. Paso 5: El conjunto de datos resultante de la aplicación de famsImputer se devuelve como un dataframe de R. 16.2.3. alleHaplotyper Esta función tiene como objetivo identificar los haplotipos que han dado lugar a los genotipos observados en un conjunto de marcadores, en un número arbitrario de familias, cuando no hay recombinación. Consideraremos también –como de hecho sucede en la práctica– que inicialmente cuando hay datos faltantes en un marcador, faltan siempre los dos alelos (esto es, no puede haber un marcador en un sujeto con un único alelo perdido; o están perdidos los dos alelos o no está perdido ninguno). Ahora bien, si quien tiene los alelos perdidos es un hijo y uno de sus progenitores es homocigoto, por ejemplo GG, entonces una G habrá sido imputada en el hijo por alleImputer, con lo que el hijo presentará el genotipo G-NA. Lo mismo ocurresi el marcador completamente perdido ocurre en el padre y hay un hijo homocigoto en ese marcador. Cuando se observan Kmarcadores se utilizará la siguiente notación para describir los alelos en el i−ésimo sujeto de la familia (i= 1 es el padre, i= 2 la madre, i > 2son los descendientes): 256 16.2. Metodología Ai=[A11iA12i. . . A1Ki A21iA22i. . . A2Ki ](16.1) Cadacolumnakdeestamatrizrepresentaun marcador, siendo (A1ki, A2ki)el par de alelos identificados en dicho marcador. Cada alelo (o ambos) puede estar perdido, en cuyo caso se denotaría como NA. Asociada a esta matrix se define otra matrix de identificadores de herencia (IDSi) del siguiente modo: IDSi=[IDS11iIDS12i. . . IDS1Ki IDS21iIDS22i. . . IDS2Ki ](16.2) donde para h= 1,2, se podría afirmar que: IDShki =     0si el alelo Ahki no pertence al haplotipo h, o está perdido 1si el alelo Ahki pertenece al haplotipo h De esta forma, si todos los términos de la matriz IDSison 0, la fase (haplotipo al que pertenece) de cada alelo es desconocida. A su vez, cuando todos los términos son iguales a 1, los alelos están alineados en el haplotipo al que pertenecen y las filas de la matriz Aipueden leerse directamente como los haplotipos del i−ésimo miembro de la familia. Cuando se leen los datos genotípicos de una familia, inicialmente las matrices IDSitienen todos sus valores idénticamente iguales a 0 para todos los miembros de la familia ya que se desconoce la fase de los genotipos. El objetivo de la función alleHaplotyper es ordenar los alelos Ahki en cada marcador de cada individuo, de tal manera que las matrices IDSicontengan tantos valores iguales a 1 como sea posible. Cuando la fila hde IDSiestá completamente (parcialmente) rellenada con unos, la correspondiente fila hde la matriz de alelos Aiestá completamente (parcialmente) ordenada, con todos sus alelos en fase. Para lograr este objetivo, el algoritmo comienza considerando sólo los hijos, tratando de ordenar los alelos en cada marcador de tal manera que el alelo en la primera fila de la matriz Aisea el heredado del padre, y el alelo de la segunda fila sea el heredado de la madre. De esta forma, si todos los marcadores pudieran ser ordenados de esta manera, la primera fila de la matriz Aisería el haplotipo heredado del padre y la segunda, el 257 . Planteamiento y Metodología heredado de la madre. Una vez que estos haplotipos han sido identificados en los hijos, se pueden identificar fácilmente en los padres. Lo que complica esta idea y hace difícil su aplicación directa es el hecho de que en algunos casos ambos padres y un hijo comparten el mismo genotipo (digamos GT para los tres sujetos), y por lo tanto no es posible conocer qué alelo se ha heredado de qué progenitor. También puede haber alelos perdidos en los padres o hijos, lo que impide determinar la procedencia de los alelos en algunos marcadores. En particular, si ambos padres tienen todos los alelos perdidos en un marcador, es imposible determinar la procedencia de los alelos de ese marcador en los hijos, al menos si hay menos de tres hijos en la familia. Como veremos en la sección 16.2.3.1 cuando la familia tiene tres o más hijos, si no hay alelos perdidos en al menos tres hijos, es posible identificar loa haplotipos incluso cuando los progenitores está completamente perdidos. Además, en la sección 16.2.3.1 se muestra que en el caso particular de tener sólo dos hijos, si los alelos parentales están disponibles en algunos marcadores y perdidos completamente en otros marcadores, es posible bajo ciertas condiciones identificar los haplotipos que contienen los marcadores perdidos en los padres. En las siguientes secciones consideramos cuatro escenarios para el procedimiento de identificación de haplotipos. En el primero no hay marcadores completamente perdidos en los padres, mientras que los hijos pueden no tener alelos perdidos en ningún marcador, o tener alelos perdidos completa o parcialmente en algunos marcadores; en el segundo escenario consideraremos el caso de que todos los marcadores en los padres están completamente perdidos y hay por lo menos tres hijos en la familia sin alelos perdidos; el tercer escenario es una mezcla de los dos anteriores: algunos marcadores tienen padres con alelos totalmente perdidos y por lo menos tres hijos completos (hijos sin alelos perdidos); algunos marcadores tienen alelos perdidos en los padres, y en algunos hijos; y algunos marcadores no tienen alelos perdidos en padres ni en hijos. Por último, en el cuarto escenario mostramos las condiciones en que es posible identificar parcialmente los haplotipos a partir de solo dos hijos cuando los padres tienen algunos marcadores completamente perdidos. 16.2.3.1. alleHaplotyper: Escenarios Escenario 1: No hay marcadores con alelos completamenteperdidos en los padres 258 16.2. Metodología El algoritmo para identificar el haplotipo al que pertenece cada alelo en este escenario es como sigue. Com los datos genotípicos de una familia: 1. Hacer i= 3,k= 1 (recuérdese que los miembros de la familia se indexan de forma que i= 1 el el padre, i= 2 la madre and i= 3,4, . . . los hijos). 2. Dado el marcador kel el i-ésimo miembro de la familia (i≥3, por lo que solo se consideran los hijos), comprobar si es posible determinar inequívocamente para cada marcador kqué alelo ha sido heredado del padre y qué alelo de la madre. Esto puede hacerse directamente en los siguientes casos: a) Si el niño es homocigoto en ese marcador; los dos alelos son iguales, por lo que es trivial asignar una copia a cada padre. b) Si el niño tiene al menos un alelo que está presente sólo en uno de los padres (y como máximo sólo hay un alelo perdido en uno de los progenitores); ese alelo se asigna a ese progenitor y el otro alelo al otro progenitor. c) Si el niño tiene un alelo perdido, el otro alelo ha debido ser imputado desde un progenitor homocigoto, por lo que la procedencia de ese alelo corresponde a dicho progenitor. d) Si un padre es homocigoto, el alelo se ha transmitido necesariamente a todos sus hijos, por lo que este alelo en cada hijo es asignado a ese progenitor, incluso si el otro progenitor tiene alelos perdidos. Si los dos alelos del hijo están presentes en ambos padres (por ejemplo, el niño tiene alelos TG y ambos padres tienen también TG), no es posible determinar la procedencia de los alelos. 3. Si la procedenciadelosalelossehadeterminadodemanera inequívoca, colocar el alelo heredado del padre en la primera fila de la matriz Aiy el alelo heredado de la madre en la segunda fila. Hacer IDS1ki =IDS2ki = 1 4. Repetir los pasos 2 y 3 para todos los marcadores en todos los hijos de la familia. 259 . Planteamiento y Metodología 5. Calcular las sumas por filas de las matrices IDS para cada hijo. Sea c1el hijo con el mayor valor de dicha suma en la primer fila de su matriz IDSi, y sea c2el hijo con el mayor valor en la suma de la segunda fila. Asimismo, sea m11, m12, . . . , m1l1el conjunto de marcadores con IDS=1 en el hijo c1, y m21, m22, . . . , m2l2el conjunto de marcadores con IDS=1 en el hijo c2.Estos marcadores tienen ya sus alelos debidamente colocados en sus haplotipos correspondientes. 6. En el padre, ordenar los alelos en los marcadores m11, m12, . . . , m1l1de tal forma que la primera fila de la matriz A1(matriz de alelos del padre) sea igual a la primera fila de la matriz Ac1(matriz de alelos del hijo c1) en las columnas m11, m12, . . . , m1l1. En la primera fila de IDS1hacerIDS1k1= 1 parak∈ {m11, m12, . . . , m1l1}. En la segunda fila de IDS1,para k∈ {m11, m12, . . . , m1l1}hacer IDS2k1= 1 si el alelo A2k1no está perdido e IDS2k1= 0 si el alelo A2k1está perdido. 7. En la madre, ordenar los alelos de los marcadores m21, m22, . . . , m2l2de tal forma que la primera fila de la matriz A2(matriz de alelos de la madre) sea igual a la segunda fila de la matriz Ac2 (matriz de alelos del hijo c2) en las columnas m21, m22, . . . , m2l2 . En la primera fila de IDS2hacer IDS1k2= 1 para k∈ {m21, m22, . . . , m2l2}. En la segunda fila de IDS2,para k∈ {m21, m22, . . . , m2l2}hacer IDS2k2= 1 si el alelo A2k2no está perdido y IDS2k2= 0 si el alelo A2k2está perdido. 8. Si todos los valores en todas las matrices IDSison iguales a 1, PARAR. Todos los genotipos están identificados. En otro caso proceder a actualizar de forma iterativa las matrices AieIDSi siguiendo el procedimiento descrito a continuación hasta que no hay ningún cambio en estas matrices entre dos iteraciones sucesivas. El objetivo de este procedimiento es esencialmente la localización de aquellos alelos que ya están correctamente colocados en los hijos, pero no en los padres y viceversa, y trasladar dicha información de unos a otros: a) Crearlamatriz idHaps,condoscolumnas y tantasfilascomo hijos en la familia (n) 260 16.2. Metodología      h11 h12 h21 h22 . . .. . . hn1hn2      (16.3) donde, siendo j= 1 el padre y j= 2 la madre: hij =                    0si el haplotipo jdel hijo ies el primer haplotipo en el progenitor j 1si el haplotipo jdel hijo ies el segundo haplotipo en el progenitor j 2si la procedencia del haplotipo jen el hijo i no puede ser determinada Esta matriz se crea inicialmente comparando los términos en las matrices de alelos Aicon IDS=1 entre padres e hijos. En el paso inicial es esperable que no todos los haplotipos puedan ser unívocamente determinados, por lo que esta matriz podría contener algunos ceros. b) Comenzando en el hijo 1 y procediendo hasta el hijo n, para cada hijo i: 1) Para la fila j(j= 1,2) en la matriz IDSicalcular el vector de diferencias Dj= (Dj1, Dj2, . . . , DjK )entre tal fila y la fila hij en la matriz IDS del j−esimo progenitor (IDSj): Unvalor de −1en la posiciónkde Dindica que para el marcador jel correspondiente haplotipo del progenitor tiene sus alelos correctamente identificados, pero no así el hijo; Un valor de 1indica que es el hijo el que tiene los alelos correctamente identificados (colocados en sus haplotipos correspondientes), pero no así el padre. Un valor de cero indica que el alelo está correctamente identificado en padre e hijo, o que no está identificado en ninguno. 261 . Planteamiento y Metodología 2) Para aquellos marcadores donde Djk = 0: Si Djk =−1: actualizar el marcador ken Aiand IDSia partir de la matriz de alelos Ajdel padre . Esta actualización cosiste en comprobar si el alelo en la posición jk de la matriz Ai(hijo) coincide con el alelo en la posición kde la fila hij en Aj (progenitor). Seaacel alelos en el hijo y apel alelo en el progenitor. Sea también bcel otro alelo en ese marcador en el hijo. Si Djk = 1 seguir el mismo procedimiento que para Djk =−1pero intercambiando los papeles de progenitor e hijo. 3) Si, como resultado del paso anterior, las matrices IDS de los padres cambian, y hay algún cero en la matriz idHaps, entonces esta matriz debe revisarse para comprobar si los haplotipos no identificados aún en hijos pueden ser ahora emparejados con alguno de los haplotipos de los padres, actualizando idHaps en consecuencia . Un marcador se considera no informativo si es heterocigoto con los mismos dos alelos en todos los miembros de la familia (ya que en tales condiciones es imposible determinar qué alelo viene del padrey qué alelo de la madre). Como resultado del algoritmo anterior, todos los marcadores informativos habrán identificado cuál de sus alelos procede del padre y cuál de la madre. Nótese también qu een el paso 2 de este último procedimiento, algunos alelos no imputados previamente por alleImputer pueden quedar imputados en padres y/o hijos al miosmo tiempo que se realiza el proceso de identificación de haplotipos. s. Escenario 2: Alelos completamente perdidos en los progenitores En este escenario surgen nuevas dificultades: los algoritmos anteriores se basan en realizar constantes comparaciones entre los alelos presentes en los padres y en sus hijos para lograr el objetivo final de identificar los haplotipos; cuando los padres tienen todos sus alelos perdidos tales comparaciones sencillamente no pueden realizarse. Ahora bien, si hay al menos tres descendientes genotípicamente diferentes entre sí, entonces 262 16.2. Metodología sí que es posible determinar los haplotipos en los hijos, aunque no será posible saber cual es el haplotipo que proviene del padre y cual proviene de la madre. Para entender cómo funciona el procedimiento, hay que tener en cuenta que las combinaciones de haplotipos en los padres para las cuales puede haber al menos tres posibles hijos genotípicamente diferentes en la familia son las que se muestran en la tabla 16.1. Supondremos inicialmente que no hay alelos perdidos en los hijos. Caso I Caso II Caso III Progenitor 1 A/B A/B A/B Progenitor 2 A/B A/C C/D A|A A|A A|C Posible A|B A|B A|D Descendencia B|B A|C B|C B|C B|D Tabla 16.1: Configuraciones haplotípicas cuando son posibles al menos tres hijos genotípicamente diferentes. Enlatabla16.1puedeobservarsefácilmentequeparacualquiercombinación de tres hijos en cualquiera de los casos anteriores, siempre hay al menos un haplotipo común compartido entre dos de ellos. La idea para la identificación de los haplotipos sin conocer los alelos de los padres comienza con la identificación del haplotipo compartido entre dos hijos. Una vez que se identifica este haplotipo, automáticamente quedan también identificados los correspondientes haplotipos complementarias. Si los dos hijos seleccionados son heterocigotos, habríamos identificado el haplotipo común y los dos complementarios, en total tres haplotipos. Si uno de los dos hijos fuese homocigoto, habíamos terminado con encontrando dos haplotipos diferentes (el haplotipo repetido que posee el hijo homocigoto y el complementario en el otro hijo). Volviendo nuevamente la tabla 16.1, observamos también que cualquiera que sea el grupo de dos hijos con el que comience el algoritmo, en el tercer hijo debe estar presente al menos uno de los haplotipos identificados en aquellos dos. De esta manera, el algoritmo para encontrar los haplotipos a partir de tres hijos sin alelos perdidos, y sin conocer a los padres, procede como sigue: 263 . Planteamiento y Metodología 1. Dados tres hijos genotípicamente distintos, 1, 2 y 3, para cada marcador k= 1,2, . . . , K encontrar los alelos comunes entre los hijos 1 y 2, 1 y 3 y 2 y 3. En aquellos en que existan tales alelos comunes existan, construir los correspondientes haplotipos , así como sus haplotipos complementarios. Denotemos por H(m) ij un conjunto de haplotipos encontrado de esta manera a partir de los hijos iyj. Nótese que, en muchos casos, dependiendo del número de alelos en cada marcador, puede haber más de un conjunto de dichos haplotipos derivados de los genotipos de los hijos iy j. También algunas veces, cuando la pareja de hijos seleccionada tenga los haplotipos A/CyB/D, no habrá ningún haplotipo común a ambos hijos. En cualquier caso, si todos los H(m) ij están vacíos es que ha habido un error de genotipado o una recombinación. En ambos casos se genera una incidencia (se informa de la situación al usuario creando una entrada en una lista de incidencias) y el algoritmo se detiene. 2. En otro caso, para cada H(m) ij =Ø determinar si al menos uno de los haplotipos del conjunto H(m) ij está también presente en el tercer hijo: Si no se cumple esta condición, se crea una incidencia y el algoritmo se detiene. Si hay más de un conjunto de haplotipos que cumplen esta condición (esto puede ocurrir dependiendo de la configuración alélica, que podría ser compatible con varias estructuras de haplotipos distintas) entonces no hay una solución única y algoritmo se detiene sin haber identificado los haplotipos. Si sólo hay un conjunto de haplotipos Hij que cumpla esta condición, ir al paso 3. 3. Si hay más de 3 hijos en la familia, identificar el conjunto de haplotipos que han sido encontrado en el paso 2 con un juego de tres haplotipos compatible con alguno de los casos que aparecen en la tabla 16.1 y determinar los haplotipos en el progenitor 1 y en el progenitor 2 (no será posible saber quien es el padre y quien la madre). Emparejar estos haplotipos parentales con el resto de los hijos, determinando qué haplotipos particulares están presentes en cada hijo. En caso de que algún hijo tenga alelos perdidos, imputarlos cuando sea posible. 264 16.2. Metodología 4. PARAR y devolver, para cada hijo, la pareja de haplotipos que se ha identificado. Esto significa devolver las matrices Aide cada hijo, con los alelos de cada haplotipo en una fila, y poner unos en aquellas posiciones de la matriz IDSipara aquellos alelos que han sido ubicados en su haplotipo correspondiente. Para que este algoritmo pueda funcionar es necesario que por lo menos tres hijos tienen todos sus alelos completos, sin valores perdidos. El resto de los hijos (en el caso de que la familia tenga más de tres hijos) pueden tener alelos perdidos, que podrían resultar imputados en el paso 3 del algoritmo en caso de que sus alelos no perdidos fuesen compatibles con una única posible pareja de haplotipos parentales; si los genotipos de los hijos con alelos perdidos encajan con varias posibles parejas distintas de haplotipos de los padres, entonces no es posible determinar qué haplotipos concretos son los que portan dichos hijos. Escenario 3: Combinación de los dos escenarios anteriores El proceso de construcción de haplotipos en este escenario es obviamente más complicado que en los otros dos casos, y el algoritmo para la identificación de haplotipos es una mezcal de los dos anteriores. Llamaremos algoritmo 1 al algoritmo descrito en el primer escenario, y algoritmo 2 al segundo. 1. Si no hay ningún marcador con alelos perdidos en ambos progenitores, aplicar el algoritmo 1. 2. Contar el número de unos en las matrices IDSipara todos los miembros de la familia. Sea IDSNr ese número. 3. Si hay marcadores con alelos perdidos en ambos padres, localizar combinaciones de marcadores que tengan al menos tres hijos completamente genotipados sin valores perdidos, y que incluyan al menos un marcador con padres completamente perdidos. Sea S={S1, S2, . . . , Sr}el conjunto de tales combinaciones. Los Sitales que Si⊆∪j=iSjson eliminados de S. 4. Si S=Ø PARAR. Si no, ordenar los conjuntos Sien orden decreciente de su tamaño. Para i= 1 hasta r: 265 . Planteamiento y Metodología Esta tabla puede implementarse fácilmente en forma de un algoritimo mediante una cadena de condiciones Si-Entonces. 16.2.3.2. alleHaplotyper: Implementación El núcleo de la función alleHaplotyper es la función famHaplotyper. Esta función: 1. Recibe como datos de entrada la matriz de datos imputados devueltos por alleImputer para una familia. 2. Aplica los algoritmos descritos en el escenario 3 anterior (téngase en cuenta que este algoritmo se adapta también a los escenarios 1 y 2) o en el escenario 4, de acuerdo a la disponibilidad de hijos y de información genotípica. 3. Devuelve: a) Una matriz igual a la matriz de entrada, pero con los nuevos alelos imputados b) Unamatriz conlasmismasdimensionesquelaanteriorllena de ceros y unos. El valor cero indica un alelo que no está en fase y el 1 que sí lo está. c) Una matriz con dos columnas que se corresponden con los haplotipos encontrados en cada miembro de la familia. Como función auxiliar, alleHaplotyper incluye también la función famsHaplotyper, encargada de aplicar famHaplotyper secuencialmente a todas las familias en el dataframe. De este modo, el funcionamiento de alleHaplotyper puede reducirse al siguiente algoritmo: 1. Llamar a la función alleImputer para leer los datos familiares e imputar marcador a marcador todos aquellos alelos que sea posible. 2. Llamar a la función famsImputer. Esta función: Identifica todas las familias en el dataframe. Pasa secuencialmente los datos de cada familia a la función famHaplotyper, que lleva a cabo la identificación de haplotipos aplicando los algoritmos descritos en la sección anterior. 272 16.2. Metodología Devuelve una lista que contiene el conjunto de datos original, los datos genotípicos imputados por alleImputer, los imputados por alleHaplotyper, la matriz IDS con ceros y unos, y los haplotipos (completos o parciales) hallados en todos los sujetos. 3. Opcionalmente, muestra un breve resumen del proceso de identificación de haplotipos, que contiene la tasa de imputación final conseguida tras el haplotipado, la proporción de alelos en fase, las proporciones de haplotipos completos, parciales y perdidos, y el tiempo empleado en todo el proceso . El paquete alleHap , tal como se encuentra en CRAN, dispone de una vignette que explica su funcionamiento y contiene numerosos ejemplos que permiten ver parte de la casuística a la que nos hemos enfrentado en su desarrollo y que se ha contemplado en los algoritmos anteriores. 273 Capítulo 17 Resultados 17.1. Resultados del análisis de bases de datos poblacionales Con respecto al análisis de datos poblacionales, el trabajo que abordamos consistión en identificar las variantes genéticas con mayor significación asociadas a la nefropatía diabética avanzada en la diabetes tipo 2 (T2D) en la población de la isla de Gran Canaria. Para ello dispusimos de una base de datos en la que se disponía del genotipo de algo más de 4 millones de marcadores en 110 sujetos –todos afectados con diabetes tipo 2– clasificados en casos o controles (55 sujetos en cada grupo) según que tuviesen una nefropatía avanzada o no. El objetivo final era determinar si hay condiciones genéticas diferentes entre ambos grupos que pudieran estar relacionadas con el avance de la nefropatía. Esta es la fase preliminar de un estudio más amplio; los marcadores que puedan detectarse en esta criba deberán ser posteriormente observados en otra muestra independiente que permita confirmar -o nola validez de la asociación detectada. El análisis de estos datos requirió la realización de varias fases: control de calidad por muestras y marcadores (eliminando los marcadores y sujetos que no superan los requisitos de calidad), imputación de valores perdidos y marcadores adyacentes a los observados, y ejecución del análisis estadístico de asociación con la base de datos resultante. Debemos señalar que el bajo tamaño de nuestra muestra dificulta la detección de posibles asociaciones, por lo que en la fase de control de calidad hemos relajado alguna de las condiciones sobre la selección de individuos para 275 . Resultados no reducir aún más el tamaño muestral. Describimos a continuación brevemente las distintas fases de este proceso. 17.1.1. Resultados del control de calidad Los resultados del control de calidad obtenidos incluyen medidas de calidad por muestras (individuos) y por variantes (marcadores). 17.1.1.1. Medidas para el control de calidad de marcadores Las medidas para el control de calidad de marcadores que hemos considerado son las siguientes: eficiencia de genotipado (proporción de genotipos perdidos por marcador), MAF (frecuencia del los alelos alternativos) and HWE (frecuencias alélicas (genotípicas) que permanecen constantes en una población de una generación a la siguiente). 17.1.1.1.1. Eficiencia de genotipado La figura 17.1 se ha generado para investigar y/o analizar la eficiencia de genotipado (también denominada como SNP coverage). Para la evaluación de la eficiencia o tasas de genotipado, en nuestro caso, y de acuerdo con la figura 17.1, hemos utilizado 0.1 como límite máximo aceptable de pérdidas. Con este umbral, han sido conservados solo aquellos SNPs con menos de un 10% de tasa de pérdidas (o más de un 90% de eficiencia de genotipado). Finalmente,despuésdelanálisisdelosgenotiposperdidos,seobtuvo que la tasa/eficiencia de genotipado para todos los individuos fue 0.97. 17.1.1.1.2. Frecuencia de alelos alternativos Para el evaluar las frecuencia de los alelos alternativos (MAF) hemos elegido un umbral dependiendo del número de de sujetos (n), si bien no la relación habitual MAF = 10/n, puesto el número de sujetos disponibles en el estudio era bajo, sino que hemos seleccionado el umbral resultante de la siguiente expresión: MAF =1 n×2=1 110 ×2= 0,0045 (17.1) 276 17.1. Resultados del análisis de bases de datos poblacionales Figura 17.1: Distribución acumulativa de SNPs Con este umbral de frecuencias alélicas alternativas, todos las variantes alélicas homocigotas (SNPs homozigotos) han sido excluidas de nuestro conjunto de datos, ya que durante análisis posteriores (específicamente en los test de asociación) estos valores no aportan información por ser iguales en casos y cotroles. Como consecuencia de este del control de calidad, se eliminaron un total de 484556 SNPs. 17.1.1.1.3. Equilibrio de Hardy-Weinberg El principio de Hardy-Weinberg establece que la variación genética en una población se mantendrá constante de una generación a la siguiente, en ausencia de factores perturbadores. Ello se traduce en que la frecuencia relativa con que se observan los genotipos debe ser igual al producto de las frecuencias relativas de los alelos que los forman. Cuando un marcador no se encuentra en equilibrio de Hardy-Weinberg puede ser confundido con un marcador asociado a la enfermedad, por lo que los marcadores fuera de esta condición deben eliminarse del estudio. Para evaluar ladesviaciónde los sujetos controles(de nuestroestudio caso-control), hemos generado un gráfico Cuantil-Cuantil (QQ) para 277 . Resultados poder apreciar si existe desviación en cada SNP. Figura 17.2: Gráfico cuantil-cuantil de sujetos controles de aquellos p-valores en Equilibrio de Hardy Weinberg. Del gráfico anterior se puede apreciar cómo existe un SNP con un p-valor desviado considerablemente del equilibrio de Hardy Weinberg. Dicho SNP (denominado rs114833138) está localizado en la posición 186204334 of the chromosome 4. Esta desviación tan grande podría indicar un error de genotipado. 17.1.1.2. Medidas para el control de calidad de muestras Las medidas para el control de calidad de muestras (sujetos) que hemos considerado son: tasa de pérdidas (proporción de genotipos perdidos por sujeto), discordancias de género (comprobación de la concordancia entre el género que se puede extraer del análisis de los genotipos y el de la identificaión del individuo), estratificación de la población (individuos con un origen genético significativamente diferente del resto de la muestra de estudio), tasa de heterocigosidad (proporción de genotipos heterocigóticos para un individuo dado) y parentesco entre individuos (comprueba si los sujetos son familiares cercanos). 278 17.1. Resultados del análisis de bases de datos poblacionales 17.1.1.2.1. Tasa de pérdidas Para el estudio de la tasa de pérdidas por individuo, establecimos un umbral en donde se detetaba un cambio cualitativo en la tasa de pérdida de datos. Figura 17.3: Distribución acumulativa del Call Rate (1 - tasa de pérdidas) por individuo. La figura 17.3 muestra la distribución acumulativa del Call Rate (1 - tasa de pérdida) de los genotipos de los individuos, con el fin de investigar la proporción de SNPs genotipados perdidos. Basándose en la figura 17.3, a partir de un call rate >91%, todas los individups tienen una buena calidad de genotipado. Siendo más rigurosos (en términos de call rate por muestra), y eligiendo un umbral >97%, tendríamos que eliminar 11 de 110 muestras de nuestro estudio, es decir, el 10% del número total de individuos genotipados, con lo que decidimos no usar dicho umbral tan estricto. 17.1.1.2.2. Discordancias de género Este paso del control de calidad se realiza para comprobar que el sexo declarado de los individuos coincide con el determinado por su número de cromosomas X. 279 . Resultados Si existe un número elevado de discordancias de género, se puede suponer que todos los identificadores de ejemplo podrían haberse mezclado de alguna manera. En nuestro caso, sólo había una discordancia, por lo que asumimos que la mayoría de los identificadores de los sujetos fueron asignados correctamente entre datos clínicos y los genéticos. 17.1.1.2.3. Estratificación de la población Para el estudio de la estratificación de la población (o detección de valores atípicos étnicos) cabe destacar que, cuando las muestras de estudio comprenden múltiples grupos de individuos que difieren sistemáticamente en tanto ascendencia genética como en fenotipos, suelen aparecer asociaciones espúrias entre las poblaciones mezcladas. Estas pueden deberse a diferencias en la ascendencia y no a una verdadera asociación genética con la enfermedad, lo que conllevaría tanto a falsos positivos como a falsos negativos [62]. Por tanto, aunque este es un paso importante en el control de calidad de datos poblacionales en muestras muy grandes donde cabe esperar mezclas étnicas, en nuestro caso decidimos que no era necesario tenerlo en cuenta, puesto que todos los sujetos genotipados pertenecían a la misma población. 17.1.1.2.4. Parentesco entre individuos El parentesco entre individuos ocurre cuando parejas o grupos de sujetos están más estrechamente relacionados entre sí que la media de la población, lo que indica que son familiares cercanos [157]. Tales individuos suelen presentar correlaciones que pueden provocar asociaciones erróneas (falsos positivos o falsos negativos). Así, de acuerdo a la figura 17.5, se puede apreciar que claramente existen dos parejas de sujetos que son familiares entre sí, con lo que uno de cada pareja tuvo que eliminarse. El criterio aplicado para decidir qué muestra eliminar, fue seleccionar aquel con la mayor proporción de datos (SNPs) perdidos. 17.1.1.2.5. Tasa de heterocigosidad Para analizar la tasa de heterozigosidad (H) hemos representado en la misma figura sus correspondientes valores con los de coeficiente de consanguinidad de Wright (F), donde una Fpositiva indicaría un exceso de 280 17.1. Resultados del análisis de bases de datos poblacionales Figura 17.4: Ejemplo de una red de parentesco más compleja (aparecen varias relaciones de parentesco) [161]. Figura 17.5: Red de parentesco entre sujetos de nuestro estudio (aparecen dos relaciones de parentesco). homocigotos (baja heterocigosidad), y una Fnegativa indicaría un exceso de heterocigotos (alta heterocigosidad) [157]. La presencia de una Finusualmente alta en un individuo podría indicar que ha habido una problema de genotipado o que la muestra provenía de una población diferente, y por lo tanto debería ser eliminada. Mediante la representación de la figura 17.6 se puede identificar la existencia de individuos con una inusual tasa de heterocigosidad. Figura 17.6: Histogramas de heterocigosidad H, y valores inversamente proporcionales F(antes del control de calidad de muestras). 281 . Resultados 4 años, del 3,1% para los niños de 5 a 9 años de edad, y un 2,4% para los de 10 a 14 años [170]. Se espera que de 2005 a 2020 se duplque el número de niños que debutan con menos de 5 años y que el número de los que debutan antes de los 15 años aumente en un 70% [171]. Hay grandes diferencias entre la diabetes T1D con inicio en la infancia y la T1D en la que el debut se produce cuando el individuo es adulto. El inicio de T1D en la infancia se asocia con cetosis/cetoacidosis más frecuentes, con una severa descompensación metabólica, con una mala función de las células beta residuales, fuerte autoinmunidad humoral contra las células de los islotes y la insulina, una mayor frecuencia de infecciones, una duración más corta de los síntomas y más independencia de los mecanismos de activación estacional, todo ello en mayor medida que en aquellos sujetos cuyo debut en T1D se produce a la edad adulta, lo que apunta a una forma más agresiva de la diabetes, [172],[173], [189],[174]. De acuerdo con algunos estudios [175], la mayor frecuencia de la enfermedad en hombres es otra característica aún no explicada de la diabetes tipo 1 en adultos los jóvenes. La región HLA es el determinante genético más importante de la susceptibilidad a la diabetes tipo 1 [176]. De acuerdo con Gillespie et al [177], en el Reino Unido más de 90% de los niños con diabetes tipo 1 portan haplotipos HLA de clase II: DRB1 * 03-DQB1 * 02: 01 (DR3DQ2) y / o DRB1 * 04-DQB1 * 03: 02 (DR4-DQ8), y el diplotipo de mayor riesgo, DR3-DQ2 / DR4-DQ8, está presente en el 50% de los casos de diabetes de inicio muy precoz. Los autoanticuerpos asociados a la diabetes pueden ser utilizados como marcadores de T1D para sujetos jóvenes con mayor susceptibilidad genética a la enfermedad, y se ha encontrado asociación entre una edad precoz de aparición de T1D con la presencia de ciertos haplotipos HLA de alto riesgo, que se encuentran con mayor frecuencia en niños T1D diagnosticados antes de los 5 años de edad que en los diagnosticados cuando son mayores [178, 177]. Los estudios en parejas de gemelos sugieren que gran parte de la variabilidad de la edad de inicio T1D está determinada genéticamente [179]. La edad de inicio puede ser considerada como un indicador de la susceptibilidad genética, estando un inicio más temprano de la enfermedad relacionado con una componente genética más fuerte y por lo tantoconunmayorriesgopara losfamiliaresenprimergrado[177, 181]. Ahora bien, el incremento que se registra en la incidencia de T1D en la población joven no puede ser exclusivamente debido a cambios en el acervo genético de la población,sinoquemás bien sugiereunainfluencia 288 17.2. Resultados del análisis de bases de datos familiares temprana del medio ambiente, por ejemplo modificaciones epigenéticas queseproducenyadurante elperiodoperinatal. Variosestudioshan analizado la influencia de factores maternales en el riesgo de diabetes y han mostrado la existencia de una asociación entre la edad de la madre en el parto y un mayor riesgo que T1D se inicie en la infancia; asimismo, el orden de nacimiento ha mostrado también cierta asociación con una disminución significativa en el riesgo de la enfermedad (Sumnik et al. 2004, Cardwell et al. 2005, Haynes et al. 2007) (Bingley et al. 2000). El debut en la infancia (pero no en la edad adulta) del padre parece asociarse con la edad de debut de los hijos, mientras que solo las madres con una edad de debut anterior a los 10 años parece afectar a la edad de debut de sus hijos, pero no de sus hijas [183, 184]. Otros estudios indican que el hecho de que la madre debute con T1D antes o durante el embarazo no afecta al riesgo de diabetes del hijo de manera distinta a madres con edad de debut adulta [184]. Cabe destacar también que para el estudio de la genética y la patogénesis de la diabetes tipo 1, se ha desarrollado un importatne un esfuerzo internacional denominado T1DGC. Este proyecto ha sido constituido con miles de familias afectadas por la T1D, incluidas de todas partes del mundo. La colección de datos que se ha reunido representa un recurso extraordinario, no sólo por los datos genéticos, sino también por la información clínica asociada. La base de datos (de enero 2009) contiene información de 14494 sujetos en 3275 familias, de las cuales 2849 contienen al menos dos hermanos afectados con T1D y 426 contienen sólo un hijo afectado. Hay un total de 6271 hijos afectados y 1673 no afectados en esta base de datos. Entre los progenitores, 194 padres y 85 madres se están afectados por T1D. La edad de inicio está disponible para todos los hijos afectados y para algunos de los padres afectados (concretamente 130 padres y 67 madres). Los sujetos en la base de datos T1DGC han sido reclutados en cuatro regiones: Asia-Pacífico (561 familias, 2289 sujetos), Europa (con exclusión del Reino Unido, 1221 familias, 5502 sujetos), América del Norte (1330 familias, 5967 sujetos) y Reino Unido (163 familias, 736 sujetos) La base de datos contiene información de los alelos de varios marcadores en el complejo mayor de histocompatibilidad humano HLA, en particular HLA-A, HLA-B, HLA-CW, HLA-DPA, HLA-DPB, HLA-DQA, HLA-DQB, HLA-DRB, así como CTLA4 y el gen de la insulina-HPH SNPs. Los alelos en estos marcadores están completos para 12370 sujetos en la base de datos: 2215 padres y 2651 madres, 289 . Resultados así como 6005 hijos afectados y 1499 hermanos no afectados (los 2124 sujetos restantes tienen los genotipos completamente perdidos, normalmente padres o hermanos que no aportaron muestras para el genotipado, aunque se reunió toda o parte de su información clínica). Hemos utilizado nuestro paquete alleHap para identificar los haplotipos HLA de este conjunto de datos. Los haplotipos que comprenden los marcadores DRB-DQA-DQB se sabe que están relacionados con el riesgo de T1D. Entre los 12.370 individuos con genotipo completo en estos marcadores, fue posible obtener los haplotipos también completos en 11.095 (89,7%). 493 (4%) fueron parcialmente haplotipados y en 782 (6,3%) no pudo hallarse ninguno de sus haplotipos de forma unívoca. Entre los 2124 sujetos con genotipos completamente perdidos, 849 (40%) pudieron haplotiparse completamente, 53 (2.5%) fueron parcialmente haplotipados y en 1222 (57.5%) no pudo obtenerse ningún haplotipo. El objetivo de nuestro estudio fue estudiar los factores maternos asociados con el inicio precoz e infantil de la T1D, que pudieran utilizarse como predictores de esta forma de la enfermedad en la descendencia utilizando para ello el conjunto de datos disponible en T1DGC el 1 de octubre 2009. 17.2.1.2. Distribución del número de haplotipos de Alto Riesgo en la base de datos T1DGC La tabla 17.1 muestra la distribución de frecuencias del número de haplotipos de alto riesgo DR3-DQ2 y DR4-DQ8, dependiendo de si los sujetos tienen la enfermedad, o no. Como puede verse, globalmente el 60% de los sujetos afectados es protador de dos haplotipos de riesgo frente a sólo el 35,3% en los no afectados. En el caso del Reino Unido la frecuencia de portadores de dos haplotipos de riesgo entre los afectados de T1D duplica a la frecuencia observada entre los no afectados. En todo caso, debe tenerse en cuenta, a la hora de interpretar estos datos, que la base de datos T1DGC no constituye una muestra aleatoria de las poblaciones estudiadas, sino una meustra compuesta por familias que tienen al menos dos hijos T1D. Por tanto los sujetos no afectados son padres y hermanos de sujetos afectados, con los que comparten su genética, por lo que es esperable una alta frecuencia de estos haplotipos incluso entre las personas no afectadas. 290 17.2. Resultados del análisis de bases de datos familiares Número de Haplotipos de Riesgo T1D 0 1 2 Datos globales: No 1544 (19.7%) 3519 (45%) 2761 (35.3%) Sí 620 (9.5%) 1997 (30.5%) 3933 (60%) Asia-Pacífico: No 311 (23.2%) 517 (38.6%) 511 (38.2%) Sí 115 (12.2%) 281 (29.8%) 546 (58%) Europa: No 544 (18.8%) 1291 (44.7%) 1051 (36.4%) Sí 233 (9%) 794 (30.8%) 1551 (60.2%) Norte América: No 654 (20.2%) 1511 (46.7%) 1068 (33%) Sí 254 (9.5%) 841 (31.6%) 1569 (58.9%) Reino Unido: No 35 (9.6%) 200 (54.6%) 131 (35.8%) Sí 18 (4.9%) 81 (22.1%) 267 (73%) Tabla 17.1: Distribución del número haplotipos de riesgo DR3-DQ2 y DR4-DQ8 para sujetos afecto y no afectos en la base de datos T1DGC. Se muestran datos regionales y globales. Cuando se considera la edad de inicio de los sujetos, la tabla 17.2 muestra la distribución de frecuencias del número de haplotipos de altoriesgo DRB-DQA-DQB (en particular DR3-DQ2 y DQ8 DQ4) en sujetos de la base de datos T1DGC, a nivel mundial y por regiones. Puede observarse que más del 93% de los sujetos con un inicio de la diabetes tipo 1 antes de la edad de 5 años tienen al menos un haplotipo de riesgo, y más de 61 % tienen dos haplotipos de riesgo. Por lo tanto, el número de haplotipos de alto riesgo se relaciona no sólo con la presencia de T1D, sino también con una edad más temprana de inicio de la enfermedad. 17.2.1.3. Factores maternos asociados con la aparición temprana y en la infancia de T1D Para evaluar si hay algún tipo de efecto materno en la edad de aparición de T1D más allá de lo que puede explicarse por la presencia de los haplotipos de riesgo consideramos las siguientes variables: 291 . Resultados Número de Haplotipos de Riesgo Edad de debut 0 1 2 Datos globales: [0,5) 102 (7.1%) 451 (31.4%) 885 (61.5%) [5,10) 169 (9.3%) 585 (32.2%) 1065 (58.5%) [10,15) 180 (10.6%) 552 (32.6%) 963 (56.8%) 15 o más 164 (10.8 %) 402 (26.5%) 950 (62.7%) Sin T1D 1544 (19.7%) 3519 (45%) 2761 (35.3%) Asia-Pacifico: [0,5) 16 (7.8%) 77 (37.4%) 113 (54.9%) [5,10) 29 (11.2%) 80 (31 %) 149 (57.8%) [10,15) 29 (11.4%) 76 (29.9%) 149 (58.7%) 15 o más 41 (19.4 %) 48 (22.7%) 122 (57.8%) Sin T1D 311 (23.2%) 517 (38.6%) 511 (38.2 %) Europa: [0,5) 33 (7.4%) 146 (32.7%) 267 (59.9 %) [5,10) 58 (8.8%) 229 (34.8%) 371 (56.4 %) [10,15) 63 (9.9%) 214 (33.8%) 357 (56.3 %) 15 o más 77 (9.5%) 201 (24.9%) 529 (65.6%) Sin T1D 544 (18.8%) 1291 (44.7%) 1051 (36.4%) Norte América: [0,5) 49 (7.2%) 207 (30.4%) 424 (62.4 %) [5,10) 76 (9.7%) 253 (32.3%) 455 (58%) [10,15) 82 (11.7%) 236 (33.7%) 382 (54.6%) 15 o más 44 (9.4%) 142 (30.3%) 282 (60.3%) Sin T1D 654 (20.2%) 1511 (46.7%) 1068 (33%) Reino Unido: [0,5) 4 (3.8%) 21 (19.8%) 81 (76.4%) [5,10) 6 (5%) 23 (19.3%) 90 (75.6%) [10,15) 6 (5.6%) 26 (24.3%) 75 (70.1%) 15 o más 2 (6.7%) 11 (36.7%) 17 (56.7%) Sin T1D 35 (9.6%) 200 (54.6 %) 131 (35.8%) Tabla 17.2: Distribución del número haplotipos de riesgo DR3-DQ2 y DR4-DQ8 dependiendo de la edad de debut en la base de datos. Se muestran datos regionales y globales. 292 17.2. Resultados del análisis de bases de datos familiares Estado de la Madre (tiene T1D/ no tiene T1D). La edad de la madre en el momento del nacimiento del hijo. Para madres afectadas con T1D: • La edad de inicio de la madre. • Si el inicio dela diabetes materna seproduceantes o después del nacimiento del niño afectado. • Número de años desde el diagnóstico de la T1D materna hasta el momento del nacimiento del niño afectado (si procede). Con el fin de evitar efectos indeseados del tamaño de la familia (es decir, el sesgo en favor de factorespresentes en las familias más grandes), se incluyeron en el análisis sólo los dos primeros hermanos afectados en cada familia (2849 familias). El género del sujeto, la positividad de anticuerpos, el número de enfermedades autoinmunes asociadas, AAID, el número de haplotipos HLA de riesgo, y los genotipos INS y CTLA4 fueron incluidos en el modelo como variables independientes y se analizaron en todas las familias. 17.2.1.4. Factores maternos, considerando todas las madres de la muestra Considerando el modelo lineal: Y=β0+ p ∑ i=1 βiXi+ε dónde: La variable dependiente Yes la edad de inicio de la T1D (en los dos primeros hijos afectados de cada familia). Las variables independientes Xison: •BirthAgeMother: la edad de la madre en el momento del nacimiento del niño. •NRiskHaps: número de haplotipos de riesgo DR3-DQ2 y DR4-DQ8 del individuo. •R_gad65 yr_ia2: TAG y IA2 positividad de anticuerpos. •Género: Masculino o Femenino 293 . Resultados •Ins_hph1: genotipo del gen de la insulina HPH1. El genotipo de referencia es AA y el modelo analiza los efectos de TA y TT en comparación con AA. •CTLA4: genotipo CTLA4. El genotipo de referencia es AA, con respecto al cual se analizan los efectos de AG y GG. •AIDn: número de enfermedades autoinmunes. •T1DM: Variable indicadora de si la madre tiene T1D. •T1DF: Variable indicadora de si el padre tiene T1D. La estimación del modelo se muestra en la tabla 17.3. Se puede apreciar que las variables ins_hph1, CTLA4 yAIDn no son significativas. El reajuste de este modelo sin estas variables produce el resultado se muestra en la tabla 17.4. La diferencia entre ambos modelos (diferencia en suma residual de cuadrados) no es significativa (p =0.2304). Estimate Std. Error t value Pr(>|t|) (Intercept) 24.0789 0.7601 31.68 0.0000 birthAgeMother -0.2265 0.0195 -11.60 0.0000 nRiskHaps -0.9039 0.1491 -6.06 0.0000 r_gad65 -3.3821 0.1992 -16.98 0.0000 r_ia2 -0.6928 0.1991 -3.48 0.0005 gender.female -0.7993 0.1981 -4.03 0.0001 ins_hph1.TA 0.0840 0.2332 0.36 0.7186 ins_hph1.TT 0.9910 0.5824 1.70 0.0889 ctla4.AG -0.3306 0.2194 -1.51 0.1319 ctla4.GG -0.1993 0.2828 -0.70 0.4811 AIDn 0.3251 0.2557 1.27 0.2038 T1DM.Yes -2.1153 0.6293 -3.36 0.0008 T1DF.Yes -0.7408 0.3918 -1.89 0.0587 Tabla 17.3: Estimación del modelo lineal para la edad de debut del hijo/hija. Como podemos ver, después de ajustar por el número de haplotipos de riesgo, la positividad de anticuerpos GAD y IA2 y el género del sujeto, el efecto de variables maternas (edad de la madre al nacer el hijo y presencia de diabetes tipo 1 en la madre) es aún apreciable. Incluso es perceptible un ligero efecto de la presencia de T1D en el padre. De hecho, la edad de inicio se reduce, como se esperaba, con el aumento en el 294 17.2. Resultados del análisis de bases de datos familiares Estimate Std. Error t value Pr(>|t|) (Intercept) 23.9721 0.7379 32.49 0.0000 birthAgeMother -0.2272 0.0195 -11.66 0.0000 nRiskHaps -0.9025 0.1488 -6.06 0.0000 r_gad65 -3.3876 0.1991 -17.02 0.0000 r_ia2 -0.6813 0.1987 -3.43 0.0006 gender.female -0.7632 0.1961 -3.89 0.0001 T1DM.Yes -2.1128 0.6292 -3.36 0.0008 T1DF.Yes -0.7393 0.3917 -1.89 0.0592 Tabla 17.4: Estimación del modelo lineal para la edad de debut del hijo/hija excluyendo variables predictivas no significativas. número de haplotipos de alto riesgo DR3-DQ2 y DR4-DQ8; menores edades de inicio también se asocian a la positividad GAD e IA2; y las niñas tienden a debutar antes que los niños. Teniendo en cuenta estas variables, la mayor edad de la madre en el parto se asocia con un inicio más temprano de la enfermedad en el hijo. Cuando la madre (y tal vez el padre) tienen diabetes tipo 1, también se detecta una tendencia a una reducción en la edad de debut en el hijo, lo que significa que probablemente hay otros factores genéticos implicados en la edad de inicio. Estimate Std. Error t value Pr(>|t|) (Intercept) 24.0912 0.7290 33.05 0.0000 birthAgeMother -0.2262 0.0192 -11.77 0.0000 nRiskHaps -1.0420 0.1523 -6.84 0.0000 r_gad65 -3.3601 0.1968 -17.07 0.0000 r_ia2 -0.7213 0.1968 -3.67 0.0003 gender.female -0.7982 0.1935 -4.12 0.0000 T1DM.Yes -2.0788 0.6267 -3.32 0.0009 T1DF.Yes -0.9063 0.3854 -2.35 0.0187 HLA.ACwB.A1-B8 0.6052 0.2392 2.53 0.0114 HLA.ACwB.A24-B39 -3.4473 0.9700 -3.55 0.0004 Tabla 17.5: Estimación del modelo lineal para la edad de debut del hijo/hija incluyendo haplotipos HLA A-Cw-B. Nuestro paquete permite la exploración del efecto de otros haplotipos posibles en la edad del sujeto de inicio. Por ejemplo, al considerar haplotipos HLA clase I en los loci A-CW-B, algunos estudios indi295 . Resultados can [185] que el haplotipo A1-B8 (HLA-A*0101-Cw*0701-B*0801) puede asociarse con DR3-DQ2 y modificar el riesgo asociado a este haplotipo. Hemos utilizado alleHap para explorar los haplotipos en esta región y hemos encontrado que A1-B8 es un haplotipo relativamente frecuente, presente en 1104 sujetos en la base de datos. AsimismotambiénhemosobservadoqueelhaplotipoA24-B39(HLAA*2402-CW*0702-B*3906) parece estar asociado a una menor edad de inicio. Cuando estos haplotipos se incluyen en el modelo anterior se obtiene la estimación que se muestra en la tabla 17.5, en la que se aprecia un efecto significativo de ambos haplotipos. 17.2.1.5. Factores maternos considerando sólo las madres con T1D Cuando sólo se consideran las madres con diabetes tipo 1, pueden introducirse en el modelo los efectos de la edad de inicio de T1D de la madre, o el tiempo transcurrido desde el diagnóstico de la madre hasta el momento del parto (años de evolución de la enfermedad). Como el número de madres con diabetes tipo 1 en la base de datos es bajo (n = 67) no se puede esperar gran resolución por parte del modelo. De hecho, si tenemos en cuenta la mismas variables que antes, llegamos a los resultados de la tabla 17.6, donde la única variable significativa resulta ser el número de haplotipos de riesgo. Estimate Std. Error t value Pr(>|t|) (Intercept) 18.8559 4.2940 4.39 0.0000 birthAgeMother -0.2264 0.1254 -1.81 0.0741 nRiskHaps -3.8879 0.8426 -4.61 0.0000 r_gad65 0.9741 1.2350 0.79 0.4321 r_ia2 -0.3782 1.1688 -0.32 0.7469 sexfemale 0.1404 1.1648 0.12 0.9043 inshph1TA 2.0843 1.3626 1.53 0.1293 inshph1TT 1.2318 2.9629 0.42 0.6785 ctla4AG -1.0317 1.3931 -0.74 0.4607 ctla4GG -0.0495 1.7396 -0.03 0.9774 numEnfAuto 0.8142 1.6676 0.49 0.6265 T1DFYes -1.9170 1.8754 -1.02 0.3092 Tabla 17.6: Estimación del modelo lineal para la edad de debut del hijo/hija usando sólo datos de familias en la que la madre tiene T1D. 296 17.2. Resultados del análisis de bases de datos familiares Si reajustamos el modelo dejando sólo el número de haplotipos de riesgo y la edad de la madre en el parto (tabla 17.7) vemos que el efecto de esta variable sigue siendo significativo. Estimate Std. Error t value Pr(>|t|) (Intercept) 19.7262 3.2209 6.12 0.0000 nRiskHaps -3.8667 0.7033 -5.50 0.0000 birthAgeMother -0.2342 0.1103 -2.12 0.0358 Tabla 17.7: Estimación del modelo lineal para la edad de debut del hijo/hija usando sólo datos familiares procedentes de familias en las que la madre tenda T1D, considerando sólo el número de haplotipos de riesgo en el sujeto y la edad de debut de la madre en la infancia. Introduciendo ahora las variables OnsetM (que especifica la edad de inicio de T1D en la madre) y motherEvolTime (tiempo desde el diagnóstico de la diabetes tipo 1 en la madre hasta el nacimiento del niño) se obtienen los resultados mostrados en la tabla 17.8. Teniendo en cuenta el significado de estas variable, cabe esperar fuerte colinealidad entre ellas (al fin y al cabo, la edad de la madre en el parto es igual a su edad de debut más el número de años de evolución hasta el parto), lo que produce que no se detecte significación en ninguna. Estimate Std. Error t value Pr(>|t|) (Intercept) 14.7636 3.8348 3.85 0.0002 nRiskHaps -3.7556 0.7067 -5.31 0.0000 birthAgeMother -0.3558 0.2180 -1.63 0.1054 onsetM 0.3191 0.2319 1.38 0.1715 motherEvolTime 0.2285 0.2580 0.89 0.3776 Tabla 17.8: Estimación del modelo lineal para la edad de debut del hijo/hija usando sólo datos de familias en las madre tenga T1D, considerando el número de haplotipos de riesgo en el sujetos, la edad de debut de la madre en la infancia, la edad de debut, y el tiempo de evolución de T1D (in years). Debido a esta colinealidad, dado que la edad de la madre en el parto es una variable que ha resultado significativa en la muestra global, tendría sentido incluir en el modelo sólo la edad de debut de la madre o sólo el número de años de evolución con T1D. Tras probar ambos modelos, 297 . Conclusiones b) La edad de inicio de los hijos es menor en promedio para las madres afectadas por T1D. c) Para las madres afectadas por T1D, se detecta cierta asociación negativa entre el tiempo de evolución de la enfermedad y la edad de inicio de la enfermedad en los hijos diabéticos (aunque hay que tener cautela con este resultado debido a los posibles factores de confusión). X. Al comparar la frecuencia de los haplotipos de riesgo en la muestra de las Islas Canarias frente a muestras de la España peninsular, no se han detectado diferencias significativas. Lo mismo ocurre al comparar la muestra española con el resto de la muestra europea. En cualquier caso, los resultados no son concluyentes debido al bajo número de familias canarias en la muestra. 304 Bibliografía [1] Sharon R Browning and Brian L Browning. Haplotype phasing: existing methods and new developments. Nature Reviews Genetics, 12(10):703–714, 2011. [2] Brian L Browning and Sharon R Browning. A unified approach to genotype imputation and haplotype-phase inference for large data sets of trios and unrelated individuals. The American Journal of Human Genetics, 84(2):210–223, 2009. [3] S. J. Mack et al. Common and well-documented hla alleles: 2012 update to the cwd catalogue. Tissue Antigens, 81(4):194–203, 2013. [4] Paul IW de Bakker et al. A high-resolution hla and snp haplotype map for disease association studies in the extended human mhc. Nature genetics, 38(10):1166– 1172, 2006. [5] EC Castelli, CT Mendes-Junior, LC Veiga-Castelli, NF Pereira, ML Petzl-Erler, and EA Donadi. Evaluation of computational methods for the reconstruction of hla haplotypes. Tissue Antigens, 76(6):459–466, 2010. [6] Virtual Medical Centre. Genetic DNA. © Virtual Medical Centre, 2010. [7] Harvey F Lodish, Arnold Berk, S Lawrence Zipursky, Paul Matsudaira, David Baltimore, James Darnell, et al. Molecular cell biology, volume 4. WH Freeman New York, 4 edition, 2007. [8] Neil A. Campbell, Brad Williamson, and Robin J. Heyden. Biology: Exploring Life. Pearson Prentice Hall, Boston, Massachusetts, 2006. [9] Duncan C. Thomas. Statistical methods in genetic epidemiology. Oxford University Press, 2004. [10] RERF. Japan-US Research Foundation. Characteristics of chromosome groups: Karyotyping. © Radiation Effects Research Foundation, 2007. [11] Departament of Biology. Chromatid Definition. © Biology Online, 2008. [12] PBworks Online Team Collaboration. Online Computational Biology Textbook. © Radiation Effects Research Foundation, 2007. 305 Bibliografía [13] Scitable. DNA Is a Structure That Encodes Biological Information. © Nature Education, 2009. [14] Wikibooks Principles of Biochemistry. Nucleic acid I: DNA and its nucleotides. Wikibooks, 2011. [15] Genetics Home Reference. Base pair. Lister Hill National Center for Biomedical CommunicationsU.S. National Library of Medicine, 2007. [16] Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, and Peter Walters. Molecular Biology of the Cell. New York and London: Garland Science, fourth edition, 2002. [17] Departament of Biology. DNA Structure. © Penn State University, 2004. [18] International Human Genome Sequencing Consortium et al. Finishing the euchromatic sequence of the human genome. Nature, 431(7011):931–945, 2004. [19] Howard Gest. Evolution of knowledge encapsulated in scientific definitions. Perspectives in biology and medicine, 44(4):556–564, 2001. [20] Indira Rajagopal. Genome Organization. Oregon State University, 2009. [21] Suzanne Clancy and William Brown. Translation: DNA to mRNA to protein, volume 1. Nature Education, 2008. [22] Robert C Elston, Jaya M Satagopan, and Shuying Sun. Genetic terminology. In Statistical Human Genetics, pages 1–9. Springer, 2012. [23] Anne Cronin and Mary Beth Mandich. Human development and performance throughout the lifespan. Cengage Learning, 2015. [24] Bruce Alberts, Dennis Bray, Karen Hopkin, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, and Peter Walter. Essential cell biology. Garland Science, 2013. [25] Talking Glossary of Genetic Terms. Locus. National Human Genome Research Institute, 2009. [26] N Malats and F Calafell. Basic glossary on genetic epidemiology. Journal of epidemiology and community health, 57(7):480–482, 2003. [27] E.P. Muljadi. Human genetics. © Google Books, 2012. [28] Daniel L. Hartl. Essential genetics: A genomics perspective. Jones & Bartlett Publishers, 5 edition, 2014. [29] Wikipedia. Allele - Dominant and recessive alleles. © Wikipedia, 2013. [30] GCSE Bitesize. Recessive and dominant alleles. © British Broadcasting Corporation - BBC, 2007. [31] Scitable. Haplotype / Haplotypes. © Nature Education, 2009. 306 Bibliografía [32] Jeffrey Mahr. Operating a Sex Machine - Meiosis. OpenStax CNX, © Rice University, 2015. [33] Lynn B. Jorde. Genetic Variation and Human Evolution. The American Society of Human Genetics, Inc., 2003. [34] Monya Baker. Structural variation: the genome’s hidden architecture. Nature methods, 9(2):133–137, 2012. [35] RGH Cotton and CR Scriver. Proof of “disease causing” mutation. Human mutation, 12(1):1–3, 1998. [36] Richard Twyman. A variable genome. Wellcome Trust, 2003. [37] Jocelyn E Krebs, Benjamin Lewin, Elliott S Goldstein, and Stephen T Kilpatrick. Lewin’s essential genes. Jones & Bartlett Publishers, 2013. [38] Sirius Genomics. What is a Single Nucleotide Polymorphism? Sirius Genomics Inc., 2013. [39] Genetics Home Reference. What are single nucleotide polymorphisms (SNPs)? Lister Hill National Center for Biomedical CommunicationsU.S. National Library of Medicine, 2007. [40] Dustin J. Penn. Major Histocompatibility Complex (MHC). Macmillan Publishers Ltd., 2002. [41] SCL Gough and MJ Simmonds. The hla region and autoimmune disease: associations and mechanisms of action. Current genomics, 8(7):453, 2007. [42] Roger Horton, Laurens Wilming, Vikki Rand, Ruth C Lovering, Elspeth A Bruford, Varsha K Khodiyar, et al. Gene map of the extended human mhc. Nature Reviews Genetics, 5(12):889–899, 2004. [43] AJ Mungall, SA Palmer, SK Sims, CA Edwards, JL Ashurst, L Wilming, MC Jones, R Horton, SE Hunt, CE Scott, et al. The dna sequence and analysis of human chromosome 6. Nature, 425(6960):805–811, 2003. [44] HLA Complex. HLA Complex. Scisco Genetics Inc., 2013. [45] WHO Committee. Nomenclature for Factors of the HLA System. © Anthony Nolan Research Institute, 2010. [46] BhadranBose,David WJohnson,and Scott B Campbell. TransplantationAntigens and Histocompatibility Matching. INTECH Open Access Publisher, 2013. [47] Shizhong Xu. Principles of statistical genomics. Springer, 2013. [48] Nicholas J Schork, Tiffany A Greenwood, and David L Braff. Statistical genetics concepts and approaches in schizophrenia and related neuropsychiatric research. Schizophrenia bulletin, 33(1):95–104, 2007. 307 Bibliografía [49] W Maxwell Cowan, Kathy L Kopnisky, and Steven E Hyman. The human genome project and its impact on psychiatry. Annual Review of Neuroscience, 25(1): 1–50, 2002. [50] Shili Lin and Hongyu Zhao. Handbook on Analyzing Human Genetic Data. Springer, 2010. [51] William S Bush and Jason H Moore. Chapter 11: Genome-wide association studies. PLoS Comput Biol, 8(12):e1002822, 2012. [52] Lucia A Hindorff, Praveen Sethupathy, Heather A Junkins, Erin M Ramos, Jayashri P Mehta, Francis S Collins, and Teri A Manolio. Potential etiologic and functional implications of genome-wide association loci for human diseases and traits. Proceedings of the National Academy of Sciences, 106(23):9362–9367, 2009. [53] Frank Yates. Contingency tables involving small numbers and the χ2 test. Supplement to the Journal of the Royal Statistical Society, pages 217–235, 1934. [54] Curt Stern. The hardy-weinberg law. Science, 97(2510):137–138, 1943. [55] Andrea S. Foulkes. Applied statistical genetics with R: for population-based association studies. Springer Science & Business Media, 2009. [56] Stephen Turner, Loren L Armstrong, Yuki Bradford, Christopher S Carlson, Dana C Crawford, Andrew T Crenshaw, Mariza Andrade, Kimberly F Doheny, Jonathan L Haines, Geoffrey Hayes, et al. Quality control procedures for genomewide association studies. Current protocols in human genetics, pages 1–19, 2011. [57] C Andrews. The hardy-weinberg principle. Nature Education Knowledge, 1(8):65, 2010. [58] Anuj Gupta. Classification of complex uci datasets using machine learning and evolutionary algorithms. IJSTR, 4(5):85–94, 2015. [59] David A Freedman. Statistical models: theory and practice. cambridge university press, 2009. [60] Joel N Hirschhorn and Mark J Daly. Genome-wide association studies for common diseases and complex traits. Nature Reviews Genetics, 6(2):95–108, 2005. [61] David J Balding, Martin Bishop, and Chris Cannings. Handbook of statistical genetics, volume 1. John Wiley & Sons, 2008. [62] Lon R Cardon and Lyle J Palmer. Population stratification and spurious allelic association. The Lancet, 361(9357):598–604, 2003. [63] Richard S Spielman, Ralph E McGinnis, and Warren J Ewens. Transmission test for linkage disequilibrium: the insulin gene region and insulin-dependent diabetes mellitus (iddm). American journal of human genetics, 52(3):506, 1993. [64] Richard S Spielman and Warren J Ewens. A sibship test for linkage in the presence of association: the sib transmission/disequilibrium test. The American Journal of Human Genetics, 62(2):450–458, 1998. 308 Bibliografía [65] D Curtis. Use of siblings as controls in case-control association studies. Annals of Human genetics, 61(4):319–333, 1997. [66] Michael Boehnke and Carl D Langefeld. Genetic association mapping based on discordant sib pairs: the discordant-alleles test. The American Journal of Human Genetics, 62(4):950–961, 1998. [67] Steve Horvath and Nan M Laird. A discordant-sibship test for disequilibrium and linkage: no need for parental data. The American Journal of Human Genetics, 63(6):1886–1897, 1998. [68] Eden R Martin, Stephanie A Monks, Liling L Warren, and Norman L Kaplan. A test for linkage and association in general pedigrees: the pedigree disequilibrium test. The American Journal of Human Genetics, 67(1):146–154, 2000. [69] SA Monks and NL Kaplan. Removing the sampling restrictions from familybased tests of association for a quantitative-trait locus. The American Journal of Human Genetics, 66(2):576–592, 2000. [70] Frank Dudbridge. Pedigree disequilibrium tests for multilocus haplotypes. Genetic epidemiology, 25(2):115–121, 2003. [71] ER Martin, MP Bass, JR Gilbert, MA Pericak-Vance, and ER Hauser. Genotypebased association test for general pedigrees: The genotype-pdt. Genetic epidemiology, 25(3):203–213, 2003. [72] Stephen L Lake, Deborah Blacker, and Nan M Laird. Family-based tests of association in the presence of linkage. The American Journal of Human Genetics, 67(6): 1515–1525, 2000. [73] Eugene V Koonin. Computational genomics. Current Biology, 11(5):R155–R158, 2001. [74] BioEECS. Computational Genomics and Proteomics. Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, 2012. [75] Robert Gentleman, Vincent Carey, Wolfgang Huber, Rafael Irizarry, and Sandrine Dudoit. Bioinformatics and computational biology solutions using R and Bioconductor. Springer Science & Business Media, 2006. [76] Adi L Tarca, Vincent J Carey, Xue-wen Chen, Roberto Romero, and Sorin Drăghici. Machine learning and its applications to biology. PLoS Comput Biol, 3(6): e116, 06 2007. [77] Johannes Fürnkranz, Dragan Gamberger, et al. Foundations of rule learning. Springer Science & Business Media, 2012. [78] MLSB14. Workshop on Machine Learning for Systems Biology. European Conference on Computational Biology, 2014. [79] Sašo Džeroski, Simon Rogers, and Guido Sanguinetti. Machine learning in systems biology. In Proceedings of The Fourth International Workshop, 2010. 309 Bibliografía [80] Pedro Larrañaga, Borja Calvo, Roberto Santana, Concha Bielza, Josu Galdiano, Iñaki Inza, José A Lozano, Rubén Armañanzas, Guzmán Santafé, Aritz Pérez, et al. Machine learning in bioinformatics. Briefings in bioinformatics, 7(1):86–112, 2006. [81] Jason Weston, Christina Leslie, Eugene Ie, Dengyong Zhou, Andre Elisseeff, and William Stafford Noble. Semi-supervised protein classification using cluster kernels. Bioinformatics, 21(15):3241–3247, 2005. [82] Pierre Baldi, Yves Chauvin, Tim Hunkapiller, and Marcella A McClure. Hidden markov models of biological primary sequence information. Proceedings of the National Academy of Sciences, 91(3):1059–1063, 1994. [83] Lawrence Rabiner. First Hand: The Hidden Markov Model. IEEE Global History Network, 2012. [84] Michael Nothnagel. Genotype Imputation. University of Kiel, 2010. [85] Jonathan Marchini and Bryan Howie. Genotype imputation for genome-wide association studies. Nature Reviews Genetics, 11(7):499–511, 2010. [86] Ion Mandoiu and Alexander Zelikovsky. Bioinformatics algorithms: techniques and applications, volume 3. John Wiley & Sons, 2008. [87] Jeff A Bilmes et al. A gentle tutorial of the em algorithm and its application to parameter estimation for gaussian mixture and hidden markov models. International Computer Science Institute, 4(510):126, 1998. [88] Juan Manuel Górriz, Elmar W Lang, and Javier Ramírez. Recent advances in biomedical signal processing. Bentham Science Publishers, 2011. [89] Longbing Cao, Yong Feng, and Jiang Zhong. Advanced Data Mining and Applications: 6th International Conference, ADMA 2010, Chongqing, China, November 19-21, 2010, Proceedings, volume 6440. Springer, 2010. [90] Todd K Moon. The expectation-maximization algorithm. Signal processing magazine, IEEE, 13(6):47–60, 1996. [91] Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the royal statistical society. Series B (methodological), pages 1–38, 1977. [92] Radford M Neal and Geoffrey E Hinton. A view of the em algorithm that justifies incremental, sparse, and other variants. In Learning in graphical models, pages 355– 368. Springer, 1998. [93] Laurent Excoffier and Montgomery Slatkin. Maximum-likelihood estimation of molecular haplotype frequencies in a diploid population. Molecular biology and evolution, 12(5):921–927, 1995. [94] Montgomery Slatkin and Laurent Excoffier. Testing for linkage disequilibrium in genotypic data using the expectation-maximization algorithm. Heredity, 76: 377–383, 1996. 310 Bibliografía [95] Richard A Gibbs, John W Belmont, Paul Hardenbol, Thomas D Willis, Fuli Yu, Huanming Yang, Lan-Yang Ch’ang, Wei Huang, Bin Liu, Yan Shen, et al. The international hapmap project. Nature, 426(6968):789–796, 2003. [96] Gudmundur A Thorisson, Albert V Smith, Lalitha Krishnan, and Lincoln D Stein. The international hapmap project web site. Genome research, 15(11):1592– 1593, 2005. [97] International HapMap Consortium et al. A haplotype map of the human genome. Nature, 437(7063):1299–1320, 2005. [98] 1000 Genomes Project Consortium et al. An integrated map of genetic variation from 1,092 human genomes. Nature, 491(7422):56–65, 2012. [99] 1000 Genomes Project Consortium. 1000 Genomes - A Deep Catalog of Human Genetic Variation. © 1000 Genomes, 2012. [100] Peter H Sudmant, Tobias Rausch, Eugene J Gardner, Robert E Handsaker, Alexej Abyzov, John Huddleston, Yan Zhang, Kai Ye, Goo Jun, Markus Hsi-Yang Fritz, et al. An integrated map of structural variation in 2,504 human genomes. Nature, 526(7571):75–81, 2015. [101] Stephen S Rich, Patrick Concannon, Henry Erlich, Cecile Julier, Grant Morahan, Jorn Nerup, Flemming Pociot, and John A Todd. The type 1 diabetes genetics consortium. Annals of the New York Academy of Sciences, 1079(1):1–8, 2006. [102] NIDDK Central Repository. Type 1 Diabetes Genetics Consortium. © The National Institute of Diabetes and Digestive and Kidney Diseases, 2010. [103] Josyf C Mychaleckyj, Janelle A Noble, Priscilla V Moonsamy, Joyce A Carlson, Michael D Varney, Jeff Post, Wolfgang Helmberg, June J Pierce, Persia Bonella, Anna Lisa Fear, et al. Hla genotyping in the international type 1 diabetes genetics consortium. Clinical Trials, 7(1 suppl):S75–S87, 2010. [104] Technical Note: DNA Analysis. Imputation Estimates Genotypes at Un-Genotyped Loci. © Illumina, Inc., 2013. [105] Bryan Howie,ChristianFuchsberger,MatthewStephens, JonathanMarchini, and Gonçalo R Abecasis. Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. Nature genetics, 44(8):955–959, 2012. [106] Sebastian Zöllner and Jonathan K Pritchard. Coalescent-based association mapping and fine mapping of complex trait loci. Genetics, 169(2):1071–1092, 2005. [107] Mark J Minichiello and Richard Durbin. Mapping trait loci by use of inferred ancestral recombination graphs. The American Journal of Human Genetics, 79(5): 910–922, 2006. [108] Brian L Browning and Sharon R Browning. Efficient multilocus association testing for whole genome association studies using localized haplotype clustering. Genetic epidemiology, 31(5):365–375, 2007. 311 Bibliografía [109] Zhan Su, Niall Cardin, Wellcome Trust Case Control Consortium, Peter Donnelly, Jonathan Marchini, et al. A bayesian method for detecting and characterizing allelic heterogeneity and boosting signals in genome-wide association studies. Statistical Science, pages 430–450, 2009. [110] Stephen Leslie, Peter Donnelly, and Gil McVean. A statistical method for predicting classical hla alleles from snp data. The American Journal of Human Genetics, 82(1):48–56, 2008. [111] Chris C Spencer, Zhan Su, Peter Donnelly, and Jonathan Marchini. Designing genome-wide association studies: sample size, power, imputation, and the choice of genotyping chip. PLoS Genet, 5(5):e1000477, 2009. [112] Jonathan Marchini, Bryan Howie, Simon Myers, Gil McVean, and Peter Donnelly. A new multipoint method for genome-wide association studies by imputation of genotypes. Nature genetics, 39(7):906–913, 2007. [113] Bertrand Servin and Matthew Stephens. Imputation-based analysis of association studies: candidate regions and quantitative traits. PLoS Genetics, 3(7):e114, 2007. [114] Yun Li, Cristen Willer, Serena Sanna, and Gonçalo Abecasis. Genotype imputation. Annual review of genomics and human genetics, 10:387, 2009. [115] Paul IW de Bakker, Manuel AR Ferreira, Xiaoming Jia, Benjamin M Neale, Soumya Raychaudhuri, and Benjamin F Voight. Practical aspects of imputationdriven meta-analysis of genome-wide association studies. Human molecular genetics, 17(R2):R122–R128, 2008. [116] Eleftheria Zeggini and Andrew Morris. Analysis of complex disease association studies: a practical guide. Academic Press, 2010. [117] Matthew Stephens and Peter Donnelly. Inference in molecular population genetics. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(4): 605–635, 2000. [118] Paul Fearnhead and Peter Donnelly. Estimating recombination rates from population genetic data. Genetics, 159(3):1299–1318, 2001. [119] Na Li and Matthew Stephens. Modeling linkage disequilibrium and identifying recombination hotspots using single-nucleotide polymorphism data. Genetics, 165 (4):2213–2233, 2003. [120] Lawrence R Rabiner. A tutorial on hidden markov models and selected applications in speech recognition. Proceedings of the IEEE, 77(2):257–286, 1989. [121] Bryan N Howie, Peter Donnelly, and Jonathan Marchini. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet, 5(6):e1000529, 2009. [122] Yun Li, Cristen J Willer, Jun Ding, Paul Scheet, and Gonçalo R Abecasis. Mach: using sequence and genotype data to estimate haplotypes and unobserved genotypes. Genetic epidemiology, 34(8):816–834, 2010. 312 Bibliografía [123] Christian Fuchsberger, Gonçalo R Abecasis, and David A Hinds. minimac2: faster genotype imputation. Bioinformatics, 31(5):782–784, 2015. [124] Sayantan Das. Minimac3. Center for Statistical Genetics, University of Michigan, 2015. [125] Sharon R Browning. Multilocus association mapping using variable-length markov chains. The American Journal of Human Genetics, 78(6):903–913, 2006. [126] Paul Scheet and Matthew Stephens. A fast and flexible statistical model for largescale population genotype data: applications to inferring missing genotypes and haplotypic phase. The American Journal of Human Genetics, 78(4):629–644, 2006. [127] Yongtao Guan and Matthew Stephens. Practical issues in imputation-based association mapping. PLoS Genet, 4(12):e1000279, 2008. [128] Gillian CL Johnson, Laura Esposito, Bryan J Barratt, Annabel N Smith, Joanne Heward, Gianfranco Di Genova, Hironori Ueda, Heather J Cordell, Iain A Eaves, Frank Dudbridge, et al. Haplotype tagging for the identification of common disease genes. Nature genetics, 29(2):233–237, 2001. [129] David M Evans, Lon R Cardon, and Andrew P Morris. Genotype prediction using a dense map of snps. Genetic epidemiology, 27(4):375–384, 2004. [130] Shaun Purcell, Benjamin Neale, Kathe Todd-Brown, Lori Thomas, Manuel AR Ferreira, David Bender, Julian Maller, Pamela Sklar, Paul IW De Bakker, Mark J Daly, et al. Plink: a tool set for whole-genome association and population-based linkage analyses. The American Journal of Human Genetics, 81(3):559–575, 2007. [131] DY Lin, Y Hu, and BE Huang. Simple and efficient analysis of disease association with missing genotype data. The American Journal of Human Genetics, 82(2):444– 452, 2008. [132] Frank Dudbridge. Likelihood-based association analysis for nuclear families and unrelated subjects with missing genotype data. Human heredity, 66(2):87–98, 2008. [133] Dan L Nicolae. Testing untyped alleles (tuna)—applications to genome-wide association studies. Genetic epidemiology, 30(8):718–727, 2006. [134] Autumn Laughbaum. Comparing BEAGLE, IMPUTE2, and Minimac Imputation Methods for Accuracy, Computation Time, and Memory Usage. Golden Helix, Inc., 2013. [135] Romeo Rizzi, Vineet Bafna, Sorin Istrail, and Giuseppe Lancia. Practical algorithms and fixed-parameter tractability for the single individual snp haplotyping problem. In Algorithms in Bioinformatics, pages 29–43. Springer, 2002. [136] Russell Schwartz et al. Theory and algorithms for the haplotype assembly problem. Communications in Information & Systems, 10(1):23–38, 2010. 313