Phased Whole Genome Genetic Risk In A Family Quartet
In an embodiment of the present invention, three novel human reference genome sequences were developed based on the most common population-specific DNA sequence (“major allele”). Methods were developed for their integration into interpretation pipelines for highthroughput whole genome sequencing.
Latest The Board of Trustees of the Leland Stanford Junior University Patents:
This application claims priority to U.S. Provisional Application No. 61/474,749 filed Apr. 13, 2011, which is hereby incorporated by reference in its entirety for all purposes.
STATEMENT OF GOVERNMENT SPONSORED SUPPORTThis invention was made with Government support under contracts HL083914 and OD004613 awarded by the National Institutes of Health. The Government has certain rights in this invention.
FIELD OF THE INVENTIONThe present invention generally relates to the field of genomics. More particularly, the present invention relates to the field of whole genome sequencing and analysis.
BACKGROUND OF THE INVENTIONWhole genome sequencing of related individuals provides opportunities for investigation of human recombination and compound heterozygous loci contributing to Mendelian disease traits as well as error control. The recent publication of low-coverage sequencing data from large numbers of unrelated individuals offers a broad catalog of genetic variation in three major population groups that is complementary to deep sequencing of related individuals. (Durbin, R. M., et al. A map of human genome variation from population-scale sequencing. Nature 467, 1061-1073 (2010). This and all other cited references are incorporated herein by reference in their entirety.) Recently, investigators used a family-sequencing approach to fine map recombination sites, and combined broad population genetic variation data with phased family variant data to identify putative compound heterozygous loci associated with the autosomal recessive Miller syndrome. (Roach, J. C., et al. Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328, 636-639 (2010).)
One of the challenges to the interpretation of massively parallel whole genome sequence data is the assembly and variant calling of sequence reads against the human reference genome. Although de novo assembly of genome sequences from raw sequence reads represents an alternative approach, computational limitations and the large amount of mapping information encoded in relatively invariant genomic regions make this an unattractive option presently. The NCBI human reference genome in current use (Pruitt, K. D., Tatusova, T. & Maglott, D. R. NCBI reference sequences (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res 35, D61-65 (2007).) is derived from DNA samples from a small number of anonymous Caucasian donors and, therefore, represents a small sampling of the broad array of human genetic variation. Additionally, this reference sequence contains both common and rare disease risk variants (Chen, R. & Butte, A. J. The reference human genome demonstrates high risk of type 1 diabetes and other disorders. Pac Symp Biocomput, 231-242 (2011).) and may bias interpretation of genetic variants aligned and called against the sequence.
With the advent of the human genome project and the draft sequences of human genomic DNA came promises to revolutionize personalized health care by tailoring risk modification, medications, and health surveillance to patients' individual genetic backgrounds. High-throughput whole genome sequencing (“Next generation sequencing”) has revolutionized the study of genetic variation and facilitated a precipitous drop in the cost per quantum of genetic variation data generated. This trend has thus far outstripped Moore's law by a factor of two, opening the door for population-wide genome sequencing. As such, technologies for interpretation of the massive amounts of genetic data produced with each genome sequence must advance in step.
One of the upstream challenges to interpretation of this genetic data is the reliance on a human reference genome sequence to 1) identify the genomic location of the billions of short (30 to 500 base pair) sequence reads produced in a massively parallel fashion by high-throughput sequencing, and 2) identify variation in other individuals from this “normal” sequence.
Another barrier to the realization of the goal of genome interpretation pipelines is the difficulty in assigning “phase,” or parental origin of genetic variants, in sequencing studies, because genotypes, not haplotypes, are given at each position. The assessment of compound heterozygous and multigenic disease risk, integration of sex-specific risk inheritance, and integration of genetic background in areas of high recombination, such as the Human Leukocyte Antigen Loci, will be crucial to understanding genome wide risk in related individuals, because these analyses are wholly reliant on phased haplotype data.
Another barrier to the realization of the goal of genome interpretation pipelines is the difficulty in assigning risk estimates to the millions of genetic variants present in any individual. These variants harbor between 50 and 100 potentially damaging mutations in genes associated with Mendelian (“single gene”) disorders, the majority of which are of uncertain significance.
Therefore, there is a need for improved methods for improved methods for analyzing genome sequence data.
SUMMARY OF THE INVENTIONIn an embodiment of the present invention, three novel human reference genome sequences were developed based on the most common population-specific DNA sequence (“major allele”). Methods were developed for their integration into interpretation pipelines for highthroughput whole genome sequencing.
To address bias inherent to the mapping, assembly, and interpretation of genetic variation data from high-throughput sequencing, an embodiment of the present invention uses three “major allele” reference sequences that capture most of the population-specific genetic variation present in three major population groups. Using the current human genome reference sequence as a scaffold, three ethnicitysspecific human major allele references are developed for the major HapMap population groups (of Caucasian (CEU), East Asian (CHB/JPT), and African subpopulation (Yoruba, YRI) descent) using estimated allele frequency data at 7,917,426, 10,903,690, and 6,253,467 genomic positions cataloged in the 1000 genomes project. In an embodiment, though relatively insensitive for very rare genetic variation, the low coverage pilot sequencing data (which comprises the majority of population-specific variation data) has a sensitivity for an alternative allele of >99% at allele frequencies>10% and thus has high sensitivity for detecting the major allele. At every position at which the major allele differed from the reference sequence base, the major allele was substituted for the reference base. This resulted in reference base changes at 1,543,755, 1,658,360, and 1,676,213 positions in the CEU, YRI, and CHB/JPT populations, respectively. There were 796,548 positions common to all three HapMap population groups at which the major allele differed from the NCBI reference base. Variation from the NCBI reference genomes was relatively uniform across chromosomal locations with the exception of loci in and near the Human Leukocyte Antigen loci on chromosome 6.
Methods for integration of these novel human reference sequences into the pipeline for interpretation of high-throughput sequencing from individuals are as follows for certain embodiments of the present invention: 1) short sequence reads from massively parallel sequencing of individuals are aligned and mapped to the population specific reference genome which most closely recapitulates self-identified ethnicity, and 2) genetic variation in sequenced individuals is cataloged by comparison of mapped sequence data from individuals to the population specific major allele at that genomic position.
Another embodiment of the present invention is a method for resolving long-range haplotype phase based on family pedigree data, inheritance state determination, and population linkage disequilibrium data. A method according to an embodiment of the present invention provides for the evaluation of genome wide risk using phased haplotype data.
An embodiment of the present invention first applies the concept of inheritance state to the allele assortments resulting from single nucleotide variants. In any nuclear family, there are 2″ possible inheritance states, where n is the number of meiosis in the pedigree. An embodiment of the present invention uses a Hidden Markov Model (HMM) in which the hidden states correspond to the 2″ inheritance states as well as two fixed error states: 1) the compression/CNV state in which hemizygous structural variants in the study genomes or reference genome result in uniform heterozygosity across the quartet; and 2) the MIE-rich state which contains a high number of impossible allele assortments likely due to sequencing or assembly errors. Transition probabilities for each of the non-error and two error inheritance states are set according to the expected number of state transitions and the total number of allele assortments in the pedigree, with the remaining transition probability allocated to self-self transitions. In an embodiment, the Viterbi algorithm is used to find the most likely state path given the observed allele assortments, resulting in the assignment of an inheritance state to each allele assortment in which one or more family members have an allele differing from the reference sequence.
State transitions in this path correspond to recombination events, with recombination window resolution given by the distance between informative allele assortments. An embodiment of the present invention uses a combination of pedigree data and statistical phasing based on inheritance state and, for uniformly heterozygous positions, population linkage disequilibrium data to determine long-range haplotypes for families. In an embodiment, the inheritance state of the surrounding variants is used to phase the remaining heterozygous positions.
This information is informative for positions at which each of three individuals in a father-mother-child trio is heterozygous for a non-reference allele and the sibling is homozygous for the reference or non-reference allele. For uniformly heterozygous positions in which the family information is not informative, pair-wise pre-computed population linkage disequilibrium data is used to assign the minor allele to the paternal or maternal chromosome scaffolds according to maximization of aggregate r2. In an embodiment, phasing is performed for each adult in contigs according to passage of allele contigs to one, both, or neither of the children. Phased haplotype data are used to find compound heterozygous and multigenic disease risk alleles using prioritization schemes. Parent of origin is used in a prior expectation in combinatorial likelihood ratio estimation for common disease risk.
An embodiment of the present invention for phasing leverages pedigree and population linkage disequilibrium data and rapidly provides haplotype phased data. In contrast to existing phasing algorithms, the methods of the present invention for providing long-range phased data provide opportunities for investigation of multigenic risk for disease incorporating gene loci at widely disparate locations on the same or even different chromosomes. Current risk prediction pipelines do not incorporate parent-of-origin into pipelines for genome annotation, despite evidence for specific effects of risk alleles for common disease, including diabetes type II, depending on parent of origin. The methods of the present invention for phasing and incorporation into risk prediction pipelines for personal genomes address this limitation.
Another embodiment of the present invention is a methodology for prioritizing variants relevant to inherited Mendelian (“single gene”) disease syndromes according to disease phenotype, gene, and variant level information.
In an embodiment, non-synonymous variants (SNVs, also called nsSNPs) are annotated using a combination of prediction algorithms and manual curation. In order to prioritize the evaluation and delivery of information relating to the rare and novel variants, a rating schema is applied based on phenotype-level information about the variants in Mendelian-disease associated genes and predicted or experimentally derived variant pathogenicity. In an embodiment, this rating schema is applied to a data source built on the Online Mendelian Inheritance in Man (OMIM) database, curated to include phenotypic information about only Mendelian diseases, and disease specific and private databases.
These and other embodiments can be more fully appreciated upon an understanding of the detailed description of the invention as disclosed below in conjunction with the attached figures.
The following drawings will be used to more fully describe embodiments of the present invention.
Among other things, the present invention relates to methods, techniques, and algorithms that are intended to be implemented in a digital computer system 100 such as generally shown in
Computer system 100 may include at least one central processing unit 102 but may include many processors or processing cores. Computer system 100 may further include memory 104 in different forms such as RAM, ROM, hard disk, optical drives, and removable drives that may further include drive controllers and other hardware. Auxiliary storage 112 may also be include that can be similar to memory 104 but may be more remotely incorporated such as in a distributed computer system with distributed memory capabilities.
Computer system 100 may further include at least one output device 108 such as a display unit, video hardware, or other peripherals (e.g., printer). At least one input device 106 may also be included in computer system 100 that may include a pointing device (e.g., mouse), a text input device (e.g., keyboard), or touch screen.
Communications interfaces 114 also form an important aspect of computer system 100 especially where computer system 100 is deployed as a distributed computer system. Computer interfaces 114 may include LAN network adapters, WAN network adapters, wireless interfaces, Bluetooth interfaces, modems and other networking interfaces as currently available and as may be developed in the future.
Computer system 100 may further include other components 116 that may be generally available components as well as specially developed components for implementation of the present invention. Importantly, computer system 100 incorporates various data buses 116 that are intended to allow for communication of the various components of computer system 100. Data buses 116 include, for example, input/output buses and bus controllers.
Indeed, the present invention is not limited to computer system 100 as known at the time of the invention. Instead, the present invention is intended to be deployed in future computer systems with more advanced technology that can make use of all aspects of the present invention. It is expected that computer technology will continue to advance but one of ordinary skill in the art will be able to take the present disclosure and implement the described teachings on the more advanced computers or other digital devices such as mobile telephones or “smart” televisions as they become available. Moreover, the present invention may be implemented on one or more distributed computers. Still further, the present invention may be implemented in various types of software languages including C, C++, and others. Also, one of ordinary skill in the art is familiar with compiling software source code into executable software that may be stored in various forms and in various media (e.g., magnetic, optical, solid state, etc.). One of ordinary skill in the art is familiar with the use of computers and software languages and, with an understanding of the present disclosure, will be able to implement the present teachings for use on a wide variety of computers.
The present disclosure provides a detailed explanation of the present invention with detailed explanations that allow one of ordinary skill in the art to implement the present invention into a computerized method. Certain of these and other details are not included in the present disclosure so as not to detract from the teachings presented herein but it is understood that one of ordinary skill in the at would be familiar with such details.
Overview of Methods According to Embodiments of the Present InventionTo be described first is an overview of methods according to embodiments of the present invention. Afterwards, further details of methods according to embodiments of the present invention will be described.
Ethics and DNA Sequence Generation
Genetic counseling of test subjects was performed at the time of DNA isolation, and additional counseling was performed concurrent with interpretation of variants with disease risk and pharmacological response implications (see steps 160 and 162 of
Inheritance State Determination and Recombination Mapping
The inheritance state analysis algorithm developed by Roach, et al (Roach, J. C., et al. Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328, 636-639 (2010)), was foundational to resolve contiguous blocks of single nucleotide variants into one of four Mendelian inheritance states using a Hidden Markov Model: paternal identical, in which each child receives the same allele from the father but different alleles from the mother; maternal identical, in which each child receives the same allele from the mother but different alleles from the father; identical; and nonidentical. Two additional non-Mendelian inheritance states were modeled (compression and Mendelian inheritance error (MIE) rich, described further below (see also Roach, J. C., et al. Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328, 636-639 (2010).)). The modeling of two additional error states allowed for identification of error-prone regions that are difficult to sequence or properly map and genotype. After excluding error prone regions that are potential sources of spurious recombination site inferences, the variant allele assortments were re-analyzed using only four Mendelian inheritance states, identifying meiotic crossover windows as intervals between SNVs defining the end and start, respectively, of contiguous inheritance state blocks.
Phasing
A combination of per-trio pedigree information, inheritance state information, and population linkage disequilibrium data, described in full further below, were used to provide long-range phasing of each of the four family members (See
The inheritance state was determined for all allele assortments in the variant call sets as described in the methods herein. As shown in
Immunogenotyping
An iterative, leave-one-out heuristic search was used as described further below for the nearest tag haplotype for common HLA types (de Bakker, P. I., et al. A high-resolution HLA and SNP haplotype map for disease association studies in the extended human MHC. Nat Genet 38, 1166-1172 (2006)) using phased variant data, assigning an HLA type to each chromosome for each study subject based on this nearest tag haplotype.
Rare Variant Prioritization
A combination of prediction algorithms were used based on characteristics of amino acid change and predicted protein structural and functional changes (SIFT, Polyphen2) (Adzhubei, I. A., et al. A method and server for predicting damaging missense mutations. Nat Methods 7, 248-249 (2010); Ng, P. C. & Henikoff, S. SIFT: Predicting amino acid changes that affect protein function. Nucleic Acids Res 31, 3812-3814 (2003)), and a novel MSA of 46 mammalian species, in which the evolutionary rate and time span was computed at each genomic position according to the method of Fitch (Fitch, W. Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology. Systematic Zoology 20, 406-416 (1971)) to provide genetic risk predictions about non-synonymous coding variants in Mendelian-disease associated genes. These variants were further manually annotated according to disease phenotype features and variant pathogenicity as described with regard to Table 5 (
Common Variant Risk Prediction
A manually curated database was developed of greater than 4000 publications investigating associations between 35,997 SNPs and 1,194 diseases or traits. A combinatorial approach was applied for point estimation of likelihood ratios of disease-SNP association, and generated composite likelihood ratios for groups of SNPs and associated diseases as described previously (Ashley, E. A., et al. Clinical assessment incorporating a personal genome. Lancet 375, 1525-1535 (2010).) from phased genetic variant data (described further below). In this analysis, included were disease-SNP associations replicated in greater than 2 GWAS with a total sample size of greater than 2000 individuals and only SNPs genotyped in the HapMap CEU population to provide a population-risk framework for interpreting composite likelihood ratios.
Pharmacogenomics
432 clinical annotations between 298 SNPs and drugs (Pharmacogenomics Knowledge Base, PharmGKB) were compiled. For all family members in the quartet, associations were evaluated between 248 phased SNPs, including 147 heterozygous loci, and 141 drugs (Tables 9-11 (
To be describe below are certain details relating to methods according to embodiments of the present invention for generating and analyzing synthetic major allele reference genomes. Although certain details will be provided, one of skill in the art will understand that many variations are possible without deviating from the teachings of the present invention.
Data Sources
The 1000 genomes project has provided an extensive catalog of genetic variation in human populations, allowing for fine-scale mapping of variation in several different ethnic groups. (Durbin, R. M., et al. A map of human genome variation from population-scale sequencing. Nature 467, 1061-1073 (2010).) Data regarding allele frequencies in each of the three major Haplotype Map (HapMap) populations was obtained from genotypes in the low coverage whole genome sequencing pilot (pilot 1) of the 1000 genomes project (Dec. 13, 2010 release). This compendium catalogs variation at 7,917,426, 10,903,690, and 6,253,467 sites in the CEU, YRI, and CHB/JPT populations, respectively, with sensitivity for an alternative allele of >99% at allele frequencies>10%.
Synthetic Major Allele Reference Sequence Generation
Genomic positions were first identified corresponding to NCBI reference genome sequence 37.1 coordinates at which the major allele (defined as estimated allele frequency>50%) differed from the reference sequence base. The major allele was then substituted for each of the three HapMap populations for the reference base to create three synthetic ethnicity-specific major allele references, specific to each of the CEU, YRI, and CHB/JPT populations. This resulted in a reference base change at 1,543,755, 1,658,360, and 1,676,213 positions in the CEU, YRI, and CHB/JPT populations, respectively. As demonstrated in
Study Subjects
Clinical characteristics of study subjects are described in graphical form in the pedigree in
Clinical characteristics of study subjects are described in graphical form in the pedigree shown in
The study was approved by the Stanford University Institutional Review Board and all study subjects attended genetic counseling and provided informed written consent (or assent, in the case of the children). This consent process occurred at two points in time: before the sequencing was performed (overseen by Illumina, Inc., and conducted with a clinical geneticist) and before this clinical interpretation was performed (conducted with a genetic counselor and research assistant). Pedigree and genotyping results were discussed in a genetic counseling session in the context of information that may be obtained in a clinical interpretation of genome sequence data and the personal and family risks and benefits that may arise in obtaining this information. (Ormond, K. E., et al. Challenges in the clinical application of whole-genome sequencing. Lancet 375, 1749-1751 (2010).)
Genome Sequencing, Mapping, and Base Calling.Genomic DNA Isolation and Sequencing Library Preparation
Peripheral blood was obtained from study subjects and sample DNA was directly isolated according to standard protocols. Genomic DNA was fragmented by sonication and sequencing libraries were prepared by clonal amplification, end repair, and sequencing adapter ligation according to standard Illumina protocol by Illumina, Inc. (San Diego, Calif.). (Bentley, D. R., et al. Accurate whole human genome sequencing using reversible terminator chemistry. Nature 456, 53-59 (2008).)
Sequence Generation
Illumina Inc. performed all sequencing reactions on the GA II instrument using massively parallel reversible terminator (sequencing by synthesis) chemistry to generate 75 base pair paired-end reads. (Ashley, E. A., et al. Clinical assessment incorporating a personal genome. Lancet 375, 1525-1535 (2010).) See step 150 of
Sequence Alignment and Mapping
Paired-end short reads were aligned to NCBI reference genome build 37.1 (obtained from the UCSC Golden Path genome browser) using the Burrows-Wheeler Aligner (BWA) software version 0.5.8a as shown in step 152 of
Base Calling and Variant Calling
The samtools multi-sample pileup tool was used with default settings to compute likelihoods of observed base data at each covered position given each possible underlying genotype. BCFtools was used to apply prior probabilities of each genotype and perform base calling against both the NCBI reference genome version 37.1 and the CEU major allele reference genome sequence created as described above. Single nucleotide variants were identified at an average distance of 699 base pairs when compared with the NCBI reference and 809 base pairs when compared with the CEU major allele reference. Short indels were called concurrently with single nucleotide variants using the samtools multi-sample pileup tool. As described in
Shown in
Genetic Variant Quality Control
Both orthogonal genotyping technology and bioinformatics tools were used to perform variant quality control and exclude likely spurious genotype calls. For confirmatory testing of common variants and quality score calibration, genomic DNA isolated from saliva from all four study subjects was genotyped using a customized array built on the Illumina HumanHap 550K+Genotyping BeadChip (23 and ME, Inc., Mountain View, Calif.). This genotyping array contains probesets corresponding to approximately 578,000 single nucleotide variants, including ˜30,000 variants unique to the 23 and ME array implementation. The discordance rate between array-based genotyping and sequencing variants was calculated according to mapping quality, base depth, and genotype quality. As described above, the base quality score was filtered at the stage of initial BWA mapping. All discordant calls were excluded from risk interpretation of genetic variants, as neither technology (array based genotyping or whole genome sequencing) was considered the gold standard. This information was also subsequently used to choose quality score cutoffs as follows: variants were retained if the mean mapping quality was greater than 40, the average base coverage depth was greater than 10, and the average and minimum genotype qualities were greater than 45 and 15, respectively.
The information provided by family-based sequencing for quality control of variant calls was also leveraged. In a family quartet, of the 81 possible genotype combinations in a bialleleic state, 52 correspond to Mendelian inheritance errors (MIEs), or allele assortments that are impossible under Mendel's laws. A very small subset of these allele assortments will be due to germ-line or somatic de novo mutation events, gene conversions, or hemizygous structural variation. However, the overwhelming majority of these allele assortments result from sequencing errors. Therefore, the identification and sequestering of these variants can greatly reduce the genotyping error rate. Similarly, the identification of regions of the genome that are prone to sequencing errors or errors in mapping and consensus assembly due to discordant structural variation between the reference and sample sequence can greatly reduce the error rate. (Roach, J. C., et al. Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328, 636-639 (2010).) Lastly, the identification of allele assortments discordant with the neighboring inheritance state (state consistency errors, SCEs) allows for identification of sequencing errors. All MIEs and SCEs were excluded from subsequent annotation as well as all variants in error prone regions.
Error Rate Estimation
Genotyping error rate was estimated using three methods: 1) estimation of the MIE rate per base sequenced, 2) estimation of the MIE and SCE rate in the 24% of the genome in which the children were identical by descent, and 3) estimation of the discordant rate between the 23 and ME genotyping and the whole genome sequence call. All three methods yielded approximately the same error rate at each stage of variant quality control. As demonstrated in
Family Inheritance State Analysis
The concept of inheritance state developed by Roach, et al, was applied to the allele assortments resulting from single nucleotide variants in each of the four family members. (Kong, A., et al. Sequence variants in the RNF212 gene associate with genome-wide recombination rate. Science 319, 1398-1401 (2008).) This nuclear family of four has 4 possible inheritance states: maternal identical in which the children each inherit the same allele from the mother; paternal identical in which the children each inherit the same allele from the father; identical in which the children are identical by descent and inherit the same allele from both parents; and nonidentical in which the children inherit different alleles from both parents. Two algorithms were used to determine inheritance state for neighboring SNVs. The first heuristic algorithm binned allele assortments into 100 kb pair regions based on chromosomal position and assigned an inheritance state according to the total number of SNPs in the bin consistent with that inheritance state.
The second algorithm was based on a Hidden Markov Model (HMM) in which the hidden states correspond to the four inheritance states described above and two error states first described by Roach, et al. (Roach, J. C., et al. Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328, 636-639 (2010).) These two states were the compression/CNV state, in which hemizygous structural variants in the study genomes or reference genome result in uniform heterozygosity across the quartet, and the MIE-rich state, which contains a high number of impossible allele assortments likely due to sequencing or assembly errors. The emission probabilities for allele assortments consistent with each inheritance state were set equal to one another and the total probability of emitting an inconsistent allele assortment was set to 0.005. For the CNV/compression state, the emission probability for uniform heterozygosity was set to 0.66. The MIE-rich state was modeled to emit an MIE 33% of the time and a consistent allele assortment 67% of the time, with equal probability weight for each consistent allele assortment. Transition probabilities for each of the four non-error and two-error inheritance states were set according to the expected number of state transitions and the total number of allele assortments in the quartet with the remaining transition probability allocated to self-self transitions. Manipulation of these transition probabilities within four orders of magnitude did not qualitatively change the resulting inheritance state determination. The Viterbi algorithm was used to find the most likely state path given the observed allele assortments, resulting in the assignment of an inheritance state to each allele assortment in which one or more family members had an allele differing from the reference sequence. State transitions in this path correspond to recombination events, with recombination window resolution given by the distance between informative allele assortments.
As CNV/compression regions and MIE-rich regions are potential sources of spurious state switches and, therefore, incorrectly inferred recombination events. The identification of these regions can improve recombination inference accuracy. After excluding SNVs in these regions, a four state HMM was developed and the Viterbi path again found, resulting in an improved median recombination window resolution 963 base pairs. To identify enrichment for these recombination windows within known hotspots, a quantitative trait associated with PRDM9 allele status, calculated was the number of recombinations in which a maximum recombination rate of >10 cM/Mbp was observed. A Monte Carlo simulation was employed with 10,000 replicates to estimate the hotspot enrichment for 106 randomly placed recombination windows of width equal to that observed in the quartet, finding that 4.1% of these random windows were in hotspots. Given that 52 recombination windows in hotspots in the quartet were found, this corresponds to a p value for hotspot enrichment of 2.0×10−73.
Phasing
Haplotype phase is important to understanding genetic risk in patients with and without disease phenotypes. Resolution of genotype information provided by whole genome sequencing into phased haplotypes has long proved difficult. A variety of statistical phasing techniques exist for inferring haplotype phase in unrelated individuals that seek to minimize recombination events and/or maximize the likelihood of heterozygous positions in high linkage disequilibrium assorting together on contigs. (Williams, A. L., Housman, D. E., Rinard, M. C. & Gifford, D. K. Rapid haplotype inference for nuclear families. Genome Biol 11, R108 (2010); Marchini, J., Howie, B., Myers, S., McVean, G. & Donnelly, P. A new multipoint method for genome-wide association studies by imputation of genotypes. Nat Genet 39, 906-913 (2007); Kruglyak, L., Daly, M. J., Reeve-Daly, M. P. & Lander, E. S. Parametric and nonparametric linkage analysis: a unified multipoint approach. Am J Hum Genet 58, 1347-1363 (1996); Donnelly, K. P. The probability that related individuals share some section of genome identical by descent. Theor Popul Biol 23, 34-63 (1983); Abecasis, G. R., Cherry, S. S., Cookson, W. O. & Cardon, L. R. Merlin—rapid analysis of dense genetic maps using sparse gene flow trees. Nat Genet 30, 97-101 (2002).) These techniques do not provide long-range phasing; however, precluding assessment of multigenic contribution to disease phenotypes or assessment of parental contribution to risk profiles. Father-mother-child trio sequencing provides information on long-range haplotype phase but only provides definitive phasing information at ˜80% of heterozygous positions, as uniformly heterozygous positions are not informative. Sequencing an additional sibling allows for precise identification of meiotic crossovers, as well as a framework for understanding the inheritance state of contiguous polymorphic markers. A combination of pedigree data and statistical phasing based on inheritance state and, for uniformly heterozygous positions, population linkage disequilibrium data was used to determine long-range haplotypes for the family quartet (heuristic described in
where n is the number of heterozygous loci on h, and ri2 is the value for r2 between the minor allele at 1 and the minor allele at heterozygous loci i on h. The minor allele was then placed at 1 on haplotype scaffold hm or hp according to Lmax=max(Lm,Lp).
This methodology allowed for phasing 34% of the remaining uniformly heterozygous positions. Finally, phasing was performed for each adult in contigs according to passage of allele contigs to one, both, or neither of the children. Phasing of adult contigs across recombination sites was not attempted due to largely uninformative linkage disequilibrium structure across recombination sites, which fell into areas of high population-averaged recombination rates. This combination of pedigree and population-linkage disequilibrium-based phasing resulted in phase resolution of 97.9% of heterozygous positions; family information alone informed phasing of 96.8% of positions.
Immunogenotyping
Phased genetic variant information simplifies greatly the combinatorial problem of resolving clinically relevant haplotypes in genomic locations with high recombination rates, in which traditional population-based statistical methods for haplotype determination are most problematic. The Human Leukocyte Antigen Loci (HLA) are examples of such regions. A combination of phased haplotype information was used for each family member and known tag haplotypes (de Bakker, P. I., et al. A high-resolution HLA and SNP haplotype map for disease association studies in the extended human MHC. Nat Genet 38, 1166-1172 (2006)) to estimate HLA type based on genotype for all four family members. A leave-one-out iterative search was performed using the phased haplotype information from each family member for the nearest common tag haplotype for HLA type, assigning the HLA type for each chromosome (paternal or maternal origin for the children, and contigs transmitted to one, both, or neither of the children in adults). This resulted in HLA types for each haploid chromosome 6 as displayed in
Principle Components Analysis of Ancestry
To determine the ancestral origins for the quartet, principal components analysis (PCA) using the maternal and paternal genotypes and a subset of individuals of European ancestry was used from the Population Reference Sample (POPRES) data set. (Nelson, M. R., et al. The Population Reference Sample, POPRES: a resource for population, disease, and pharmacological genetics research. Am J Hum Genet 83, 347-358 (2008).) Briefly, maternal and paternal genotypes were combined with the genotypes from 10 individuals from Eastern/Southeastern Europe, 190 individuals from central Europe, 78 individuals from Northern/Northeastern Europe, 167 individuals from Northwestern Europe, 238 individuals from Southern Europe, 99 individuals from Southeastern Europe, 272 individuals from Southwestern Europe and 403 individuals from Western Europe (EuropeESE, EuropeC, EuropeNNE, EuropeNW, EuropeS, EuropeSE, EuropeSW, and EuropeW respectively in
Definitions and Heuristic
The operational definition of a rare variant was a variant with an allele frequency of less than 5%. A hierarchical search was used for allele frequencies that first considered ethnicity-specific allele frequencies and then population-wide allele frequencies as follows:
-
- 1. Ethnicity-specific allele frequencies from the dbSNP 132 database
- 2. Ethnicity specific allele frequencies from the 1000 genomes pilot 1 data
- 3. Allele frequencies from the dbSNP 132 full database
- 4. Allele frequencies from the 1000 genomes full project data Dec. 13, 2010 release
Alleles with an assigned rsid but no published allele frequency were considered rare but with no frequency data, and alleles not found in any published database were considered novel. This resulted in sets of 351,555 and 354,074 rare or novel variants from the HG19 and CEU reference variant call sets, respectively.
The overall heuristic for searching novel and rare variants for significant disease associations is presented in
Shown in
Nonsynonymous Coding Variants
Non-synonymous variants (nsSNPs) were annotated using a combination of prediction algorithms and manual curation. All rare and novel variants were annotated using the Sorting Intolerant from Tolerate (SIFT) algorithm (see step 1226 of
Because the prediction accuracies of SIFT and PolyPhen have been shown to depend heavily on the evolutionary anatomy of genomic position at which an nsSNP occurs (Kumar, S., et al. Positional conservation and amino acids shape the correct diagnosis and population frequencies of benign and damaging personal amino acid mutations. Genome Res 19, 1562-1569 (2009)), nsSNPs were further annotated using position-specific evolutionary features derived from a 46 species multiple sequence alignment (MSA) of vertebrate genomes obtained from the UCSC Genome Browser (http://genome.ucsc.edu/; Accessed Nov. 30, 2010). For each genomic position, the evolutionary rate (mutations per billion years) was computed (see step 1232 of
In order to prioritize the evaluation and delivery of information relating to the rare and novel variants, a rating schema was developed and applied that was based on phenotype-level information about the variants in Mendelian-disease associated genes and predicted or experimentally derived variant pathogenicity. This rating schema is summarized in Table 5 (
The children were compound heterozygous for variants in four disease-related genes: BLM (mother, son and daughter), associated with Bloom syndrome, MLH3 (son), associated with familial non-polyposis colon cancer, SLC4A5 (son), associated with proximal renal tubular acidosis, and COG7 (daughter), associated with type IIe congenital disorder of glycosylation. There were two instances of homozygosity for rare/novel variants in disease-related genes in the daughter (KRT8, associated with monilethrix, and ASAH1, associated with Farber lipogranulomatosis) and one such instance in the son (PLEC1, associated with epidermolysis bullosa simplex, Table 7 (
Synonymous Coding Variant Risk Prediction
Apart from amino-acid substitutions, there are several ways that synonymous single nucleotide polymorphisms (sSNPs) can affect a gene and its resulting protein products. Alteration of splice sites can modify how a gene is spliced and result in important changes in the resulting mRNAs; most of these alterations result in premature mRNA degradation. Creation of spurious splice sites may affect the resulting protein sequence. Other factors that affect protein production and structure include mRNA decay rates and mRNA structural motifs surrounding important regulatory sites (such as 5′ and 3′ UTRs). Finally, codon usage bias can have a direct effect on protein elongation and translational kinetics, a consequence of the correlation between codon usage frequency and tRNA availability. This leaves three main mechanisms that can be computationally explored to detect putative phenotypic changes provoked by sSNPs (See
Three models for the association between a synonymous SNPs gene function were developed.
Aberrant splicing is a phenomenon that has been linked to synonymous mutations in various studies. Creation and disruption of 5′ donor splice sites and exonic splice site enhancers through synonymous alterations have been reported to be part of the etiology of diseases such as type 1 neurofibromatosis, multiple sclerosis, and phenylketonuria. (Chamary, J. V., Parmley, J. L. & Hurst, L. D. Hearing silence: non-neutral evolution at synonymous sites in mammals. Nat Rev Genet 7, 98-108 (2006); Kimchi-Sarfaty, C., et al. A “silent” polymorphism in the MDR1 gene changes substrate specificity. Science 315, 525-528 (2007).) Splice site prediction algorithms exist and are primarily used for genome-wide gene detection of splice sites. However, they can also be used to detect putative disruption or creation of splicing sites in a simplistic fashion by comparing predictions when applying the algorithm to reference and the variant DNA sequences. Using this criteria, the maximum entropy splice site detection algorithm (Yeo, G. & Burge, C. B. Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals. J Comput Biol 11, 377-394 (2004)) was applied to the flanking sequence of a sSNP with and without the polymorphic substitution. Predictions resulting in a positive odds ratio for the reference sequence but in a negative odds ratio for the sequence with the polymorphism are flagged as putative splice site disruptions. Conversely, a combination of a negative prediction for the reference sequence and a positive score for the SNP-affected sequence is reported as putative creation of a splice site.
Several mRNA structural factors are associated with important effects on phenotype. Secondary structure may directly affect mRNA decay rates as well as confer protection from premature degradation. Furthermore, highly structured UTRs can prevent regulatory molecules, such as microRNAs, from performing proper regulatory functions. Thus, investigating the effects of sSNPs in mRNA structure becomes a great pivotal point to indirectly study putative changes in the resulting protein. A small set of articles have already laid ground on the case, by analyzing the influence of sSNPs in mRNA secondary structure and its effects on mRNA stability and decay.
RNA secondary structure prediction is a classical problem in computational biology and there are many methods that give reasonable estimates. Most of them report the resulting free energy, AG, of the predicted secondary structure, thereby giving a thermodynamic measure of structure. Algorithms for detecting non-coding RNAs use free energy along with other heuristics to detect putative biologically active transcripts. (Rivas, E. & Eddy, S. R. Secondary structure alone is generally not statistically significant for the detection of noncoding RNAs. Bioinformatics 16, 583-605 (2000).) In particular, these algorithms attempt to find a ‘structural signal’ in a certain window of nucleotides while scanning a genome. One approach to do this is by performing free energy calculations for randomized samples of the same size and monomeric or dimeric conformations than that of the current window. A Z-score is then given to the window, defined as:
where G(seq) is the free energy of the RNA sequence seq, Gμ (seq,S) is the average free energy of the sequences of the sample set S that have the same length and monomeric (or dimeric, if desired) conformation than seq, and Gσ (seq,S) is the standard deviation of the free energies. There has been evidence demonstrating that secondary structure by itself does not give a strong signal from random sequences with the same monomer or even dimer conformations. (Rivas, E. & Eddy, S. R. Secondary structure alone is generally not statistically significant for the detection of noncoding RNAs. Bioinformatics 16, 583-605 (2000).) This is reasonable since permutation of nucleotides is a far more benign alteration than deletion, insertion, or replacement. To express this in the Z-score, the definition of the sample set is modified to a set of random sequences of the same length of the window but not necessarily with the same n-meric conformation.
To apply the Z-score notion to probe if a change in secondary structure occurs with a SNP, the structural significance of the subsequence flanking the SNP was assessed. This was done by taking two windows: the flanking window, Wf, and the sampling window, Ws. The flanking window is the sequence that contains the SNP position in its midpoint. The sampling window is a subsequence of the flanking window and also contains the SNP position. The sampling is then performed from the set S(Wf,Ws) of sequences with length of the flanking window that vary only in the sampling window. Finally the Z-score, as defined previously in equation 3, is taken using this sample set:
This is done using the ViennaRNA folding package. (Hofacker, I. L. Vienna RNA secondary structure server. Nucleic Acids Res 31, 3429-3431 (2003).) The Z-score of the reference sequence is then compared with the Z-score of the sequence containing the SNP substitution and obtain a difference score. This score expresses the difference between structural importance of the sequence in the sampling window in the reference and SNP-containing sequence.
Two genes that code for the same protein using synonymous codons do not necessarily give the same result. This is mainly due to the fact that tRNA iso-acceptors do not have equal abundance in the cell. Even though this statement was confirmed in vitro several years ago, only recently has such a situation been observed as occurring in vivo. The demonstration that codon usage bias can alter translational kinetics opens an interesting new venue to search for relations between phenotype alterations and sSNPs. Codon usage bias analysis is not new, and has been fairly well studied since the beginning of the genomic era. (Eyre-Walker, A. C. An analysis of codon usage in mammals: selection or mutation bias? J Mol Evol 33, 442-449 (1991); Ikemura, T. Codon usage and tRNA content in unicellular and multicellular organisms. Mol Biol Evol 2, 13-34 (1985).) Several results confirm that, in some organisms, codon usage is also related with position, since it is not rare to see codons with similar relative frequency (relative frequency is the frequency of a codon occurring in a genome with respect to codons that code for the same amino-acid while absolute frequency is the frequency of codon occurrence with respect to the set of all codons) cluster together in particular sites. (Eyre-Walker, A. C. An analysis of codon usage in mammals: selection or mutation bias? J Mol Evol 33, 442-449 (1991); Zhang, G. & Ignatova, Z. Generic algorithm to predict the speed of translational elongation: implications for protein biogenesis. PLoS One 4, e5036 (2009).) This has led to the speculation that codon choice is directed by evolution, given that there could be selection constraints acting in some aspects of translational kinetics, such as protein elongation. Following this conceptualization, changes in codon bias via a clustering criterion are assessed.
Given an exon sequence seq, a set of pairs is first produced:
for all possible n in seq, where n is the nth codon in the sequence given the ith open reading frame, N is the total number of codons in the sequence, and reln is the relative frequency of the nth codon. The k-means clustering algorithm is then applied to Ci(seq) for each open reading frame (ORF) with a given k. This is performed with both the reference and SNP-modified sequence, SNPseq. Finally, for all ORFs, the resulting centroids are compared between both sequences and the sum of their distances is computed, taking the minimum of these values. In other words, the final codon usage score is:
CU=mini(Σdist(Ck,i(seq),Ck,i(SNPseq)) (4)
where Ck,i is the set of k centroids in the ith ORF.
Noncoding Variant Risk Prediction
A search was conduced for rare and novel variants in noncoding regions associated with introns, 3′ and 5′ UTRs, and miRNA target regions in 3′ UTRs of genes associated with Mendelian disorders as well as pre-miRNA and mature miRNA sequences targeting genes associated with Mendelian disorders. miRNA target regions and sequence coordinates were obtained from the miRbase database. As a supplement to the maximum entropy splice site disruption algorithm described above, a search was also conducted for known splice donor and acceptor sites for rare and novel variants. Rare and novel variants in pre-miRNA, mature miRNA, miRNA target regions, and splice sites were annotated as putative loss of function variants.
Structural Variant Risk Prediction
Indels were annotated based on their rarity (with allele frequencies derived from the 1000 genomes pilot 1 data), association with coding regions, whether they disrupted a splice site, and whether they were predicted to cause a frame-shift in an open reading frame. As allele frequencies for indels are less reliable than for single nucleotide variants, novel frameshift indels in coding regions or novel indels in splice sites of genes associated with OMIM-curated diseases were considered to be loss of function variants. There were 27 such variants when compared to the HG19 reference genome and 29 such variants when compared to the CEU major allele reference genome.
Common Genetic Variant Risk PredictionQuantitative Disease-SNP Association Database
As described previously, quantitative human disease-SNP associations were manually curated from the full text, figures, tables, and supplemental materials of 4,022 human genetics papers. (Ashley, E. A., et al. Clinical assessment incorporating a personal genome. Lancet 375, 1525-1535 (2010).) Also, more than 100 features were recorded from each paper, including the disease name (e.g. coronary artery disease), specific phenotype (e.g. acute coronary syndrome in coronary artery disease), study population (e.g. Finnish individuals), case and control population (e.g. 2,508 subjects with coronary artery disease proven by angiography), gender distribution, genotyping technology, major/minor risk alleles, odds ratio, 95% confidence interval of the odds ratio, published p-value, and genetic model. Studies on similar diseases were categorized and mapped to the Concept Unique Identifiers (CUI) in the Unified Medical Language System (UMLS). (Bodenreider, O. The Unified Medical Language System (UMLS): integrating biomedical terminology. Nucleic Acids Res 32, D267-270 (2004)) For each study, the frequency of each genotype and allele in the case and control populations was recorded. Strand ambiguities were resolved with an automatic strand detection algorithm described previously.
Calculation of Predicted Personal Genetic Risk for 28 Common Diseases
For each of 28 common diseases, all SNPs that had been significantly associated with the disease with a p-value of ≦10−6 in two or more Genome-Wide Association Studies (GWAS) with a total sample size of 2,000 or more subjects were identified. Genetic risk was estimated using a likelihood ratio for each SNP defined by the relative frequency of the individual's genotype in the diseased vs. healthy control populations (e.g., given an allele “A”, LR=Pr(A|diseased)/Pr(A|control)). The LR incorporates both the sensitivity and specificity of the test and provides a direct estimate of how much a test result will change the odds of having a disease. (Morgan, A. A., Chen, R. & Butte, A. J. Likelihood ratios for genome medicine. Genome Med 2, 30 (2010).) Studies with diseased subjects in the control group and studies on non-Caucasian populations were excluded. For each SNP, the LRs from multiple studies were averaged with a weight of the square root of the sample size to give higher confidence to studies with larger sample size. After removing SNPs in linkage disequilibrium (R2≧0.3 in the corresponding population group), each locus was assumed as an independent genetic test and LRs were multiplied to report the summarized score or predicted genetic risk. Pre- and post-test estimates of disease risk were calculated (see step 160 of
Calculation of Relative Population Based Disease Risk
To evaluate the relative population based disease risk, the personal genetic risk was calculated on 58 diseases for 174 CEU individuals from the HapMap project version II and III (see
Family Differential Risk and Parental Contribution to Common Disease Risk
The disease risk for each of the family members was calculated according to phased SNP genotypes as described above; this information is displayed in
Pharmacogenomics Knowledge Base
Genome-wide pharmacogenomic associations were annotated using the Pharmacogenomics Knowledge Base (www.pharmgkb.org), an online pharmacogenomics resource containing manually created annotations on a large collection of pharmacogenomics literature. Over 1400 drug-variant-phenotype relationships were curated from the literature and these relationships were used to create clinical annotations for 298 known variants (see step 162 of
All primary annotations referring to the same drug-variant-phenotype association were combined into one summary for each clinical annotation in the database. For example, the relationship between the TPMT*3B variant (rs1800460) and adverse reactions to purine analogs has been studied for years and published multiple times. These combined clinical annotations were then linked to the PMIDs of citations that were used to create the annotation, creating a written summary for each possible genotype. Using the same TPMT example, the AA genotype contains 2 copies of the *3B variant and is associated with significantly increased risk of side effects due to decreased enzyme levels; the AG genotype contains 1 copy of the *3B variant and is associated with slightly increased risk of side effects due to moderately decreased enzyme levels; and the GG genotype contains no copies of the *3B variant and is not associated with increased risk of side effects. The level of risk for any given genotype is written with respect to the risk for other genotypes. Thus, the AA genotype is associated with increased risk of side effects as compared to the AG and GG genotypes, but not necessarily at increased risk of side effects for patients on the drug in general. Risk relative to other genotypes is used because the incidence of efficacy or adverse events for any given drug has usually not been quantified. Also, the typical genotype for any given population can be difficult to determine. Although many groups use HapMap frequencies to calculate population major alleles or typical genotypes, the HapMap populations are small and ethnically very specific. Their frequencies do not necessarily represent larger population frequencies (nor were they meant for such purpose). Therefore, no attempt is made to state what the risk of any given drug reaction is compared to “normal,” since normal is not defined, but rather the risk is reported as compared to other possible genotypes.
A “Strength of evidence” score is then assigned to each clinical annotation and this score is used as a measure of confidence in the association. This strength of evidence score is based on several criteria, including independent replication of the association, population size of the studies and p-value of the association (after correction for multiple hypothesis testing, if applicable). A score of high confidence was assigned if the association has been replicated in populations of at least 1000 cases and 1000 controls of the same ethnicity with p-values>0.05 after correction. A score of medium confidence required p-values<0.05 after correction and at least one population studied of at least 100 people, with or without further replication. Finally, a low confidence score was assigned to all clinical annotations that did not meet the criteria for high or medium confidence scores, including those annotations where all evidence of the association is based on in vitro studies or pharmacokinetic data alone.
Variant Level Annotation
Drug-variant phenotypes can be very complicated and specific to certain patient populations. For simplicity and illustration purposes, clinical annotations were binned into the following groups in an embodiment of the present invention: Drug(s) More Likely to Work, Drug(s) Less Likely to Work, Drug(s) More Likely to Cause Side Effect, Drug(s) Less Likely to Cause Side Effect, Drug Dose(s) Easy to Predict, Drug Dose(s) Difficult to Predict, Drug Dose(s) Above Average, Drug Dose(s) Below Average, No Pharmacogenomic Action and/or Phenotype Unknown and/or Phenotype Not Applicable (example: ovarian cancer risk for males). An example annotation is given in Table 8 (
Development of Ethnicity Specific Major Allele References
In an embodiment of the present invention, three ethnicity-specific major allele references were developed for the CEU, YRI, and CHB/JPT HapMap population groups using estimated allele frequency data at 7,917,426, 10,903,690, and 6,253,467 positions cataloged in the 1000 genomes project. Though relatively insensitive for very rare genetic variation, the low coverage pilot sequencing data (which comprises the majority of population-specific variation data) has a sensitivity for an alternative allele of >99% at allele frequencies>10% and, thus, has high sensitivity for detecting the major allele. In an embodiment, substitution of the ethnicity-specific major allele for the reference base resulted in single base substitutions at 1,543,755, 1,658,360, and 1,676,213 positions in the CEU, YRI, and CHB/JPT populations, respectively. There were 796,548 positions common to all three HapMap population groups at which the major allele differed from the NCBI reference base (see
Study Subjects and Genome Sequence Generation
Clinical characteristics of study subjects are described in graphical form in the pedigree shown in
The Illumina sequencing platform was used to provide 39.3× average coverage of 92% of known chromosomal positions in all four family members (See
Inheritance State Analysis Identifies >90% of Sequencing Errors
Sequencing family quartets allows for precise identification of meiotic crossover sites from boundaries between inheritance states and superior error control. (Roach, J. C., et al. Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328, 636-639 (2010).) Contiguous blocks of single nucleotide variants were resolved into one of four Mendelian inheritance states or two error states. Using this methodology, 3.1% of variant positions were identified as associated with error prone regions (see
Using a combination of these methods and quality score calibration with orthogonal genotyping technology, 94% of genotyping errors (see
Consistent with PRDM9 allelic status, approximately half of all recombinations in each parent occurred in hotspots as shown in
RNF212 Haplotypes Predict Nearly Equal Male and Female Recombination Frequency
Boundaries between contiguous inheritance state blocks defined 55 maternal and 51 paternal recombination events across the quartet at a median resolution of 963 base pairs. A parallel heuristic analysis of recombination sites confirmed the observation of nearly equal paternal and maternal recombination frequency. Two SNPs (rs3796619 and rs1670533), separated by 17 kb in the gene RNF212, form a haplotype that is significantly associated with genome wide recombination rate in a sex-specific manner, such that the haplotype associated with the highest recombination rate in males is associated with a low female recombination rate. (Kong, A., et al. Sequence variants in the RNF212 gene associate with genome-wide recombination rate. Science 319, 1398-1401 (2008).) Fine scale recombination mapping (see step 158 of
PRDM9 Allelic Status Predicts High Hotspot Usage
Several loci in the vicinity of the PRDM9 gene are associated with recombination hotspot usage. (Baudat, F., et al. PRDM9 is a major determinant of meiotic recombination hotspots in humans and mice. Science 327, 836-840 (2010); 8. (Myers, S., et al. Drive against hotspot motifs in primates implicates the PRDM9 gene in meiotic recombination. Science 327, 876-879 (2010); Parvanov, E. D., Petkov, P. M. & Paigen, K. Prdm9 controls activation of mammalian recombination hotspots. Science 327, 835 (2010); Kong, A., et al. Fine-scale recombination rate differences between sexes, populations and individuals. Nature 467, 1099-1103 (2010).) The heritability for hotspot use explained by this locus is among the highest of all described quantitative trait loci and is determined by both single nucleotide substitutions in and near PRDM9 and the number of zinc finger α-helices in exon 12 of the gene. (Kong, A., et al. Fine-scale recombination rate differences between sexes, populations and individuals. Nature 467, 1099-1103 (2010).) The use of whole genome sequencing allows for fine mapping of these sites and investigation of the relationship between PRDM9 haplotypes and hotspot usage. Of the SNPs near PRDM9, rs2914276 is most significantly associated with hotspot usage heritability. In this quartet, both parents are homozygous for the rs2914276-A allele that is associated with high hotspot usage as well as the number of zinc finger repeats in PRDM9. (Kong, A., et al. Fine-scale recombination rate differences between sexes, populations and individuals. Nature 467, 1099-1103 (2010).) Zinc finger repeat number ranges between 12 and 16 in population studies; both parents carry 13 zinc finger repeats in the PRDM9. Consistent with previous population-wide observations for individuals with this combination of rs2914276 alleles and number of zinc finger repeats, it was found that 25 of 51 paternal recombination windows (49%) and 27 of 55 maternal recombination windows (49%) were in hotspots (defined by maximum recombination rate of >10 cM/Mbp), while only ˜4 (4.1%) would be expected by chance alone (p=2.0×10−73 for hotspot enrichment according to Monte Carlo permutation).
Prior Population Mutation Rate Estimates are Biased Upwards by the Reference Sequence
After excluding variants in sequencing-error prone regions, 4,302,405 positions were identified at which at least one family member differed from the NCBI reference sequence and 3,733,299 positions at which at least one family member differed from the CEU major allele reference sequence (See
After short read mapping and local realignment, variants were called against the NCBI reference genome 37.1 and the CEU major allele reference. Likely spurious variant calls were first filtered by mapping quality, read depth and genotyping quality. The inheritance state for all allele assortments was determined by HMM and error prone regions (compression regions and Mendelian inheritance error rich (MIE)-rich regions, which represent likely sequencing errors) were identified and excluded. 606,757 fewer variants were identified when compared the CEU major allele reference than the NCBI reference genome 37.1 (HG19 reference) (See
Phased Whole Genome Sequence Data Reveal HLA Types
Human leukocyte antigen (HLA) group typing from whole genome sequence data has previously been challenging due to the high recombination frequency on chromosome 6 that complicates phasing, as well as the extreme genetic diversity in the HLA loci. This approach that incorporates fine-mapping of recombination events with long-range phasing from polymorphic markers specific to the ethnic background of the sequenced family (e.g., variants called against the CEU major allele reference) and an iterative heuristic, described herein, for identifying the nearest HLA tag haplotype (de Bakker, P. I., et al. A high-resolution HLA and SNP haplotype map for disease association studies in the extended human MHC. Nat Genet 38, 1166-1172 (2006)) allowed for chromosome-specific HLA typing from whole genome sequence data as well as determination of the parental origin of HLA types (
Rare Variant Analysis Identifies Multi-Genic Risk for Familial Thrombophilia
It has been estimated from population sequencing data that apparently healthy individuals harbor between 50 and 100 putative loss of function variants in genes associated with Mendelian diseases and traits. (Durbin, R. M., et al. A map of human genome variation from population-scale sequencing. Nature 467, 1061-1073 (2010).) Many of these variants have little or no significance either because they result in subtle phenotypes or have no biological effect. Thus, prioritization of putative loss of function variants remains a significant challenge.
To further characterize the potential adverse effects of non-synonymous single nucleotide variants (nsSNVs), a multiple sequence alignment (MSA) of 46 mammalian genomes was implemented. For coding variants of unknown significance, the mammalian evolutionary rate is proportional to the fraction of selectively neutral alleles (Kimura, M. Evolutionary rate at the molecular level. Nature 217, 624-626 (1968)) and can therefore serve as a prior expectation in determining the likelihood that an observed nsSNV is deleterious. This MSA was used in combination with existing algorithms (Adzhubei, I. A., et al. A method and server for predicting damaging missense mutations. Nat Methods 7, 248-249 (2010); Ng, P. C. & Henikoff, S. SIFT: Predicting amino acid changes that affect protein function. Nucleic Acids Res 31, 3812-3814 (2003)) to provide genetic risk predictions for phased putative loss of function variants for the family quartet (see heuristic in
Of 354,074 rare or novel variants compared with the CEU major allele reference sequence, 200 non-synonymous variants were identified in coding regions, one nonsense variant, one SNV in the conserved 3′ splice acceptor dinucleotides, and 27 novel frame-shifting indels in genes associated with Mendelian diseases or traits. Consistent with prior observations and the extremely conserved regulatory role for miRNAs, no rare or novel SNVs were found in mature miRNA sequence regions or miRNA target regions in 3′ UTRs. There were four compound heterozygous variants in disease-related genes and three homozygous variants in disease-related genes (See Table 7 (
The father has experienced two venous thromboemboli, the second of which occurred on systemic anticoagulation. One variant in the gene F5 encoding the coagulation factor V confers activated protein C resistance and increased risk for thrombophilia. (Koster, T., et al. Venous thrombosis due to poor anticoagulant response to activated protein C: Leiden Thrombophilia Study. Lancet 342, 1503-1506 (1993); Ridker, P. M., et al. Mutation in the gene coding for coagulation factor V and the risk of myocardial infarction, stroke, and venous thrombosis in apparently healthy men. N Engl J Med 332, 912-917 (1995).) Another variant (the thermolabile C677T variant) in the gene MTHFR encoding methylenetetrahydrofolate reductase, predisposes heterozygous carriers to hyper-homocysteinemia and may have a synergistic effect on risk for recurrent venous thromboembolism (Ridker, P. M., et al. Interrelation of hyperhomocyst(e)inemia, factor V Leiden, and risk of future venous thromboembolism. Circulation 95, 1777-1782 (1997); Margaglione, M., et al. The methylenetetrahydrofolate reductase TT677 genotype is associated with venous thrombosis independently of the coexistence of the FV Leiden and the prothrombin A20210 mutation. Thromb Haemost 79, 907-911 (1998)), though this association has more recently been questioned. (Bezemer, I. D., Doggen, C. J., Vos, H. L. & Rosendaal, F. R. No association between the common MTHFR 677C->T polymorphism and venous thrombosis: results from the MEGA study. Arch Intern Med 167, 497-501 (2007).)
Follow-up serological analysis demonstrated the father's serum homocysteine concentration was 11.5 μmol/L. A homozygous variant in F2, a gene known to be associated with hereditary thrombophilia, was able to be excluded based on its high evolutionary rate in multiple sequence alignment (See Table 6 (
Synonymous SNP Prediction Identifies a Putative Loss of Function Variant in ATP6V0A4
Association between synonymous SNPs (sSNPs) and disease has recently been described. (Macaya, D., et al. A synonymous mutation in TCOF1 causes Treacher Collins syndrome due to mis-splicing of a constitutive exon. Am J Med Genet A 149A, 1624-1627 (2009).) sSNPs may affect gene product function in several ways, including codon usage bias, mRNA decay rates, and splice site creation and/or disruption (See
The results of a Z score method for predicting free energy change conferred by synonymous single nucleotide variants identifies loci in the coding region of ATP6V0A4 associated with a significant change in mRNA stability are shown in
Common Variant Risk Prediction Identifies Risk for Obesity and Psoriasis
Results from Genome Wide Association Studies (GWAS) provide a rich data source for assessment of common disease risk in individuals. To provide a population risk framework for genetic risk predictions for this family quartet, ancestral origins were first localized based on principal components analysis of common SNP data in each parent and the Population Reference Sample (POPRES) dataset (Nelson, M. R., et al. The Population Reference Sample, POPRES: a resource for population, disease, and pharmacological genetics research. Am J Hum Genet 83, 347-358 (2008)), demonstrating North/Northeastern and Western European ancestral origins for maternal and paternal lineages, respectively (see
Composite likelihood ratios were then calculated for 28 common diseases for 174 ethnically-concordant HapMap CEU individuals, and provided percentile scores for each study subject's composite LR for each disease studied (see
All four family members were at high risk for psoriasis, with the mother and daughter at highest risk (98th and 79th percentiles, respectively), consistent with clinical manifestations of psoriasis in both the father and daughter. It was also found that both parents were predisposed to obesity, while both children had low genetic risk for obesity. This approach was extended to incorporate haplotype phased SNP data in assessment of genetic risk, providing estimates of disease risk according to parental risk haplotypes (see
Pharmacogenomic Variant Annotation Predicts Empiric Warfarin Dosing
The Pharmacogenomics Knowledge Base, PharmGKB, was used to evaluate associations between phased variants and 141 drugs (See Tables 9-11 (
It was also found that the mother and daughter are both homozygous for an ultra-rapid metabolism allele at CYP2C19, which encodes a key metabolizer of the pro-drug clopidogrel, an antiplatelet agent used in the prevention and therapy of cardiovascular disease. Because the metabolic activity of CYP2C19 is directly correlated with the antiplatelet activity of clopidogrel, there is a higher bleeding risk associated with clopidogrel use in the mother and daughter. This finding has significant implications for the daughter who has multigenic risk for thrombophilia and may require anticoagulant therapy should she develop thrombosis; concomitant use of clopidogrel in this setting may contribute further to bleeding risk associated with systemic anticoagulation. Full details of pharmacogenetic variants in other key metabolic enzymes and associated pharmacokinetic and pharmacodynamic profiles are displayed in Table 4 (
Here, phased genetic risk assessment in a family quartet using whole genome sequencing and a novel ethnicity-specific major allele reference is discussed. At approximately 1.6 million genomic positions, it was found that the NCBI reference sequence differed from the major allele in each of the three HapMap populations. New estimates were provided for the population mutation rate which illustrate the upward bias in estimates based on comparison with the NCBI human reference genome sequence and show that use of the reference sequence for variant calling may result in spurious genotype calls at medically relevant SNP loci.
Family quartet sequencing allows for providing a fine map of meiotic crossover sites to sub-kb resolution for the first time, finding that meiotic crossovers happen with nearly equal frequency in male and female parents in this family quartet. This is in contrast to previous observations in a family quartet and to observations of greater recombination frequency in females than in males in mice. (Roach, J. C., et al. Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328, 636-639 (2010); Paigen, K., et al. The recombinational anatomy of a mouse chromosome. PLoS Genet 4, e1000119 (2008); Petkov, P. M., Broman, K. W., Szatkiewicz, J. P. & Paigen, K. Crossover interference underlies sex differences in recombination rates. Trends Genet 23, 539-542 (2007).) Estimates of human sex-specific recombination rates have demonstrated significant variation in recombination number, particularly in females. (Kong, A., et al. Fine-scale recombination rate differences between sexes, populations and individuals. Nature 467, 1099-1103 (2010); Broman, K. W., Murray, J. C., Sheffield, V. C., White, R. L. & Weber, J. L. Comprehensive human genetic maps: individual and sex-specific variation in recombination. Am J Hum Genet 63, 861-869 (1998); Kujovich, J. L. Factor V Leiden thrombophilia. Genet Med 13, 1-16 (2010).) Notably, the RNF212 haplotypes revealed by long range phasing in this quartet are associated with lower than average female and average male recombination rates (Kong, A., et al. Sequence variants in the RNF212 gene associate with genome-wide recombination rate. Science 319, 1398-1401 (2008)), potentially explaining the described observations.
In this quartet, the power of sequencing nuclear families was leveraged for identification of >90% of sequencing errors, providing unprecedented accuracy of sequence information used for genetic risk interpretation and identification of compound heterozygous and multigenic disease risk. This approach, first applied to a family quartet in which two family members had Miller syndrome and ciliary dyskinesia (Roach, J. C., et al. Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328, 636-639 (2010)) was extended into a tool for phasing genetic variants.
Algorithms for multiple sequence alignment for novel and rare nonsynonymous variant risk prediction, and functional prediction for the effects of synonymous SNPs were developed and applied. In doing so, multigenic risk for thrombophilia in the father and daughter was identified, consistent with a history of recurrent venous thromboembolism in the father despite systemic anticoagulation. This multigenic risk for thrombophilia is more consistent with the father's clinical history of recurrent thromboembolism on systemic anticoagulation than monogenic risk conferred by heterozygous factor V Leiden alone, which was previously identified as the putative etiological factor in clinical assessment. (Kujovich, J. L. Factor V Leiden thrombophilia. Genet Med 13, 1-16 (2010).) Furthermore, multigenic risk for thrombophilia identified in the daughter prior to first venous thrombosis has significant clinical implications in terms of risk mitigation.
Phased variant information allows for determination of parental contribution to common disease risk. Furthermore, several common SNPs have been associated with disease risk, most notably diabetes mellitus 2, in a sex specific manner, such that maternal origin confers direct association with risk of disease and paternal origin confers indirect association with disease risk. (Kong, A., et al. Parental origin of sequence variants associated with complex diseases. Nature 462, 868-874 (2009).) Additionally, phasing variant data demonstrated the origin of discordant parent-child risks for common disease. The presently described methods for determining parental origin of risk alleles allows for precise risk assessment utilizing this information as a prior, in contrast to existing un-phased genotype data.
The ethnicity-specific, family-based approaches to interpretation of genetic variation presented here represent the next generation of genetic risk assessment using whole genome sequencing.
The current iteration of the human population-specific reference sequences includes only genetic variation data at single nucleotide polymorphisms, or substitutions at individual base pairs. Structural variation (stretches of DNA repeats, deletions, inversions, and translocations) comprises a large degree of human genetic variation. In a further embodiment of the present invention the major allele reference genomes are revised to incorporate known common structural variants specific to the population groups described above. As the quantification of genetic variation becomes more precise, mapping and variant calling algorithms can be developed that consider the full distribution of DNA bases at variant positions in the human genome, as opposed to simply the major allele as would be known to those of ordinary skill in the art.
As genetic variation data is cataloged on other population groups (e.g., Latin American, South Asian, Pacific Islander), additional major allele reference sequences can be developed for these population groups. They can then be implemented as described above for ethnicity-aware genome sequence interpretation.
It should be appreciated by those skilled in the art that the specific embodiments disclosed above may be readily utilized as a basis for modifying or designing other image processing algorithms or systems. It should also be appreciated by those skilled in the art that such modifications do not depart from the scope of the invention as set forth in the appended claims.
Claims
1. A method for interpreting genetic data, comprising:
- receiving a reference genome sequence associated with a sub-population;
- receiving at least one short sequence read from at least one individual;
- receiving information that the at least one individual is a member of the sub-population;
- aligning the received at least one short sequence read to the human major allele reference; and
- determining at least one difference between the at least one short sequence read and the reference genome sequence associated with the sub-population.
2. The method of claim 1, wherein the sub-population is an ethnicity.
3. The method of claim 1, further comprising mapping the received short sequence reads to the human major allele reference.
4. The method of claim 1, wherein determining differences between the short sequence reads and the reference genome sequence associated with the sub-population is performed for predetermined genomic positions.
5. The method of claim 1, further comprising identifying at least one genetic variation responsive to the determined at least one difference.
6. The method of claim 1, further comprising determining at least one genetic risk for at least one disease responsive to the determined at least one difference.
7. The method of claim 1, further comprising determining the reference genome sequence associated with the sub-population using estimated allele frequency information.
8. The method of claim 1, further comprising determining a diagnosis for the individual responsive to the determined at least one difference.
9. The method of claim 1, further comprising determining a drug for the individual for treatment of effects from the determined at least one difference.
10. The method of claim 1, further comprising determining a prognosis for the individual from the determined at least one difference.
11. A computer-readable medium including instructions that, when executed by a processing unit, cause the processing unit to interpret genetic data, by performing the steps of:
- receiving a reference genome sequence associated with a sub-population;
- receiving at least one short sequence read from at least one individual;
- receiving information that the at least one individual is a member of the sub-population;
- aligning the received at least one short sequence read to the human major allele reference; and
- determining at least one difference between the at least one short sequence read and the reference genome sequence associated with the sub-population.
12. The computer-readable medium of claim 11, wherein the sub-population is an ethnicity.
13. The computer-readable medium of claim 11, further comprising mapping the received short sequence reads to the human major allele reference.
14. The computer-readable medium of claim 11, wherein determining differences between the short sequence reads and the reference genome sequence associated with the sub-population is performed for predetermined genomic positions.
15. The computer-readable medium of claim 11, further comprising identifying at least one genetic variation responsive to the determined at least one difference.
16. The computer-readable medium of claim 11, further comprising determining at least one genetic risk for at least one disease responsive to the determined at least one difference.
17. The computer-readable medium of claim 11, further comprising determining the reference genome sequence associated with the sub-population using estimated allele frequency information.
18. The computer-readable medium of claim 11, further comprising determining a diagnosis for the individual responsive to the determined at least one difference.
19. The computer-readable medium of claim 11, further comprising determining a drug for the individual for treatment of effects from the determined at least one difference.
20. The computer-readable medium of claim 11, further comprising determining a prognosis for the individual from the determined at least one difference.
21. A computing device comprising:
- a data bus;
- a memory unit coupled to the data bus;
- a processing unit coupled to the data bus and configured to receive a reference genome sequence associated with a sub-population; receive at least one short sequence read from at least one individual; receive information that the at least one individual is a member of the sub-population; align the received at least one short sequence read to the human major allele reference; and determine at least one difference between the at least one short sequence read and the reference genome sequence associated with the sub-population.
Type: Application
Filed: Apr 13, 2012
Publication Date: Dec 24, 2015
Patent Grant number: 10127346
Applicant: The Board of Trustees of the Leland Stanford Junior University (Palo Alto, CA)
Inventors: Frederick Dewey (Redwood City, CA), Euan Ashley (Menlo Park, CA), Matthew Wheeler (Sunnyvale, CA), Michael Snyder (Stanford, CA), Carlos Bustamante (Emerald Hills, CA)
Application Number: 13/445,923
