Systems and methods for exploiting relatedness in genomic data analysis
Patent Information
- Application Number
- JP2023076708
- Authority / Receiving Office
- JP · JP
- Patent Type
- Patents
- Current Assignee / Owner
- Priority Date
- 2017-09-07
- Filing Date
- 2023-05-08
- Publication Date
- 2025-01-06
- Estimated Expiration
- 2038-09-07
AI Technical Summary
【0129】 以上、本発明によれば、ゲノムデータ解析に近縁性を活用するシステム及び方法が提供できた。
Smart Images

Figure 0007608513000030 
Figure 0007608513000031 
Figure 0007608513000032
Abstract
Description
[Technical field]
[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims the benefit of U.S. Provisional Patent Application No. 62 / 555,597, filed September 7, 2017, the contents of which are incorporated herein by reference in their entirety. Additionally, a co-pending application entitled "System and Method for Predicting Relatedness in a Human Population," filed September 7, 2018, is also incorporated by reference in its entirety.
[0002] Field The present disclosure relates generally to methods and systems for the analysis of genomic data and to using relatedness in large population cohorts to link rare genetic variations to disease and disease susceptibility. More specifically, the present disclosure relates to systems and methods for establishing homoeologousity and phasing genetic variants as compound heterozygous mutations or de novo mutations. [Background technology]
[0003] Human pathologies are caused and influenced not only by environmental factors but also by genetic factors, therefore, understanding the genetic variation in human populations is important for understanding the etiology and progression of human diseases and for identifying novel drug targets for treating these diseases.
[0004] Genetic studies of health populations are particularly useful in this regard due to the availability of extensive health control data, which simplifies the study of how genetic variants contribute to human pathologies. In the past, such studies have typically been based on whole-genome genetic linkage analysis to map disease loci, which, once identified, could then be analyzed in more detail at the molecular level.
[0005] Over the last few years, the widespread availability of high-throughput DNA sequencing technologies has enabled the parallel sequencing of hundreds of thousands of human genomes. In theory, these data represent a powerful source of information that can be used to decipher the genetic basis of human disease. However, these ever-growing datasets require sustained innovation in bioinformatics tools and analysis pipelines to continue to efficiently process such very large datasets. Furthermore, the utility of relatedness and family structure in such large datasets, and the extent to which this can be exploited for variant identification and characterization, has not been fully realized or utilized. Summary of the Invention [Problem to be solved by the invention]
[0006] There is a need for improved bioinformatics tools for the analysis of large-scale genomic data, and this disclosure addresses this need. [Means for solving the problem]
[0007] In one aspect, the disclosure provides a method of phasing genetic variants in a population by leveraging population relatedness, comprising removing low quality sequence variants from a dataset of nucleic acid sequence samples obtained from a plurality of human subjects; establishing an ancestral superclass assignment for each of one or more of the samples; removing the low quality samples from the dataset; generating a first homoeologousness estimate for the subject within the ancestral superclass; generating a second homoeologousness estimate for the subject independent of the ancestral superclass of the subject; and, based on one or more of the second homoeologousness estimates, The method includes clustering subjects into a first first-degree family network, generating a third homoeologous estimate for the subjects within the first first-degree family network, integrating the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate, constructing a second first-degree family network for the subjects based on the combined homoeologous estimate, and phasing a variant as being a compound heterozygous mutation (CHM) or not a CHM according to the combined homoeologous estimate and the second first-degree family network, or identifying a variant as a de novo mutation (DNM) according to the combined homoeologous estimate and the second first-degree family network.
[0008] In some exemplary embodiments, integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to augment the first homoeologous estimate.
[0009] In some exemplary embodiments, phasing variants as compound heterozygous mutations (CHM) includes (1) phasing variants according to population allele frequency; (2) removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and phasing variants as compound heterozygous mutations (CHM). (3) removing single nucleotide polymorphisms (SNPs) with a QD of about 2 or less, or a read depth (DP) of less than about 5, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, and removing insertion or deletion (INDELs) with a QD of about 2 or less, or a DP of less than about 5, or an AB of about 10% or less, or a combination thereof; (4) selecting the remaining variants as potential compound heterozygous mutations (pCHMs), where one or more pairs of variants exist within the same sample and within the same gene; and (5) phasing the pCHMs as cis pCHMs or trans pCHMs, and then classifying the pCHMs phased as trans pCHMs as CHMs.
[0010] In some exemplary embodiments, phasing variants as compound heterozygous mutations includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of less than about 10, or an AB of about 20% or less, or a combination thereof.
[0011] In some exemplary embodiments, the method further comprises: (1) assigning scores to the CHMs according to functional effect priority; and (2) selecting the CHM with the highest functional effect priority score per gene per sample, such that when a human has multiple CHMs within the same gene, the CHM most likely to result in protein function inhibition is identified.
[0012] In some exemplary embodiments, phasing a variant as a de novo mutation includes (1) identifying the variant in the second first-degree family network and samples of the trio; (2) assigning genotype likelihood scores to the variant in parent samples and corresponding child samples in the trio, calculating the probability that the variant is a de novo mutation, and identifying the variant as a putative de novo mutation if the calculated probability is statistically significant; (3) identifying the variant in a child sample of the trio and identifying the variant as a putative de novo mutation if the variant is not present in any parent sample of the trio; (4) identifying a putative de novo mutation with a genotype quality (GQ) annotation of less than about 35 in the child sample, or a putative de novo mutation with an alternative allele count (AC) of 10 or more across samples from multiple human subjects, or a putative de novo mutation with a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or a putative de novo mutation with an allele balance (AB) of more than about 2% in either parent sample. (5) filtering the identified putative de novo mutations by removing putative de novo mutations that have novo mutations, or putative de novo mutations with an allele balance (AB) of less than about 15% in the child sample, or putative de novo mutations with an AB of more than about 90% in the child sample, or putative de novo mutations with alternative allele homozygosity in either parent sample, or a combination thereof; and (6) combining the filtered identified putative de novo mutations, thereby forming a putative de novo mutation dataset.
[0013] In some exemplary embodiments, the method further comprises classifying a putative de novo mutation in the putative de novo mutation dataset as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is greater than or equal to about 0.15 in the child samples and less than or equal to about 0.02 in each parent sample, and if the mapping quality of the putative de novo mutation is not less than about 40, and if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if the MAC of the putative de novo mutation is less than about 20 across all samples, and if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs.
[0014] In some exemplary embodiments, the method further comprises classifying the medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality annotation of about 90 or greater in the parent samples, and if the medium-confidence de novo mutation has a read depth of about 10 or greater in each parent sample, and if the medium-confidence de novo mutation has an alternative read depth of about 7 or greater in the child sample, and if the medium-confidence de novo mutation has a QD of greater than about 3 for a SNP, and if the medium-confidence de novo mutation has a QD of greater than about 5 for an INDEL.
[0015] In one aspect, the disclosure provides a method of identifying compound heterozygous mutations (CHM) in a population, comprising: identifying variants in DNA sequence samples from a plurality of human subjects; establishing an ancestral superclass assignment for the subjects based on the identified variants; generating a first homozygosity estimate for the subjects within the ancestral superclass; generating a second homozygosity estimate for the subjects independent of the ancestral superclass for the subjects; clustering the subjects into a first first-degree family network based on one or more of the second homozygosity estimates; generating a third homozygosity estimate for the subjects within the first first-degree family network; and integrating the first homoeologous estimate with the third homoeologous estimate to obtain a combined homoeologous estimate; constructing a second first-degree family network based on the combined homoeologous estimate; phasing variants in the sample according to population allele frequencies; classifying the phased variants as potential CHMs based on the presence of two or more variants in the same subject and gene; and phasing the potential CHMs as being in cis or in trans with another variant in the same subject and gene, and then classifying the potential CHMs phased as being in trans as CHMs.
[0016] In some exemplary embodiments, the method further comprises filtering the identified variants before an ancestral superclass assignment for the subject is established. In some exemplary embodiments, the method further includes filtering the identified variants before the first homoeologous estimate and the second homoeologous estimate for the subject are generated.
[0017] In some exemplary embodiments, filtering variants comprises removing variants with an alternative allele frequency of greater than about 10% across samples from multiple human subjects, or a p-value of greater than about 10. -6These include removing variants that are above 100% confidence interval (CI) or that violate Hardy-Weinberg equilibrium (HWE), or that are missing calls in more than about 5% of samples from multiple human subjects, or a combination of these.
[0018] In some exemplary embodiments, the method includes removing low quality samples after the identified variants have been filtered. In some exemplary embodiments, low quality samples are those with a D-stat greater than 0.12 or less than 75% read coverage at 20x, or both.
[0019] In some exemplary embodiments, integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to augment the first homoeologous estimate.
[0020] In some exemplary embodiments, homoeologous estimation involves the calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. In some exemplary embodiments, the method further comprises filtering the variants after they have been phased according to population allele frequencies.
[0021] In some exemplary embodiments, filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of about 5 or less, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of about 5 or less, or an AB of about 10% or less, or a combination thereof.
[0022] In some exemplary embodiments, phasing variants according to population allele frequency comprises dividing a DNA sequence sample of a human subject into genomic segments of approximately equal size and with substantial segment overlap and breakpoints in intergenic regions.
[0023] In some exemplary embodiments, the potential CHMs are phased based on trio data, or parent-offspring data, or full sibling data, or distant relative data, or a combination thereof, or are phased based on minor allele counts (MACs), or are phased based on population allele frequencies, or a combination thereof.
[0024] In some exemplary embodiments, the method further comprises scoring the CHMs according to a functional effect priority and selecting the CHMs with the highest functional effect priority score per gene per sample, thereby obtaining a set of medically important variants.
[0025] In some exemplary embodiments, the DNA sequence sample comprises an exome sequence. In some exemplary embodiments, the plurality of human subjects comprises more than 10K subjects.
[0026] In some exemplary embodiments, the D-stat of a low quality sample is determined by comparing the distribution of the actual allele balance of the sample to the distribution of the expected allele balance using the Kolmogorov-Smirnov (KS) test.
[0027] In some exemplary embodiments, filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of less than about 10, or an AB of about 20% or less, or a combination thereof.
[0028] In another aspect, the disclosure provides a non-transitory computer-implemented method for identifying compound heterozygous mutations (CHMs) in a population. In general, the non-transitory computer-implemented method includes using a data processor of a computing device to identify variants in DNA sequence samples from a plurality of human subjects, using the data processor to establish an ancestral superclass assignment for the subjects based on the identified variants, using the data processor to generate a first homozygosity estimate for the subjects within the ancestral superclass, using the data processor to generate a second homozygosity estimate for the subjects independent of the ancestral superclass for the subjects, using the data processor to cluster the subjects into a first first-degree family network based on one or more of the second homozygosity estimates, and using the data processor to generate a third homozygosity estimate for the subjects within the first first-degree family network. using a data processor to phase variants in the sample according to population allele frequencies; using a data processor to classify the phased variants as potential CHMs based on the presence of two or more variants in the same subject and gene; and using a data processor to phase the potential CHMs as being in cis or in trans with another variant in the same subject and gene and then classify the potential CHMs phased as being in trans as CHMs.
[0029] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to filter the identified variants before an ancestral superclass assignment for the subject is established.
[0030] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to filter the identified variants before a second homoeologous estimate for the subject is generated.
[0031] In some exemplary embodiments, filtering variants comprises removing variants with an alternative allele frequency of greater than about 10% across samples from multiple human subjects, or a p-value of greater than about 10. -6 These include removing variants that are above 100% confidence interval (CI) or that violate Hardy-Weinberg equilibrium (HWE), or that are missing calls in more than about 5% of samples from multiple human subjects, or a combination of these.
[0032] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to remove low quality samples after the identified variants have been filtered.
[0033] In some exemplary embodiments, low quality samples are those with a D-stat greater than 0.12 or less than 75% read coverage at 20x, or both.
[0034] In some exemplary embodiments, integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to augment the first homoeologous estimate.
[0035] In some exemplary embodiments, homoeologous estimation involves the calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. In some exemplary embodiments, the non-transitive computer-implemented method further comprises using a data processor to filter the variants after they have been phased according to population allele frequencies.
[0036] In some exemplary embodiments, filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of about 5 or less, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of about 5 or less, or an AB of about 10% or less, or a combination thereof.
[0037] In some exemplary embodiments, phasing variants according to population allele frequency comprises dividing a DNA sequence sample of a human subject into genomic segments of approximately equal size and with substantial segment overlap and breakpoints in intergenic regions.
[0038] In some exemplary embodiments, the potential CHMs are phased based on trio data, or parent-offspring data, or full sibling data, or distant relative data, or a combination thereof, or are phased based on minor allele counts (MACs), or are phased based on population allele frequencies, or a combination thereof.
[0039] In some exemplary embodiments, the non-transitive computer-implemented method further comprises using a data processor to score the CHMs according to a functional effect priority and to select the CHMs having the highest functional effect priority score per gene per sample, thereby obtaining a set of medically important variants.
[0040] In some exemplary embodiments, the DNA sequence sample comprises an exome sequence. In some exemplary embodiments, the plurality of human subjects comprises more than 10K subjects.
[0041] In some exemplary embodiments, the D-stat of a low quality sample is determined by comparing the distribution of the actual allele balance of the sample to the distribution of the expected allele balance using a KS test.
[0042] In some exemplary embodiments, filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of less than about 10, or an AB of about 20% or less, or a combination thereof.
[0043] In another aspect, the disclosure provides a system for performing the method and non-transitory computer-implemented method, the system generally includes a data processor, a memory coupled to the data processor, and a program stored in the memory, the program identifying variants in DNA sequence samples from a plurality of human subjects, establishing an ancestral superclass designation for the subjects based on the identified variants, generating a first homoeologousness estimate for the subjects within the ancestral superclass, generating a second homoeologousness estimate for the subjects independent of the ancestral superclass of the subjects, clustering the subjects into a first first-degree family network based on one or more of the second homoeologousness estimates, and clustering the subjects into a third homoeologousness estimate for the subjects within the first first-degree family network. The method includes instructions for generating an estimate, integrating the first homogologous estimate and the third homogologous estimate to obtain a combined homogologous estimate, constructing a second first-degree family network based on the combined homogologous estimate, phasing variants in the samples according to population allele frequencies, classifying the phased variants as potential CHMs based on the presence of two or more variants in the same subject and gene, and phasing the potential CHMs as being in cis or in trans with another variant in the same subject and gene, and then classifying the potential CHMs phased as being in trans as CHMs.
[0044] In some exemplary embodiments, the program includes instructions for filtering the identified variants before an ancestral superclass assignment for the subject is established. In some exemplary embodiments, the program includes instructions for filtering the identified variants before the first homoeologous estimate and the second homoeologous estimate for the subject are generated.
[0045] In some exemplary embodiments, filtering variants comprises removing variants with an alternative allele frequency of greater than about 10% across samples from multiple human subjects, or a p-value of greater than about 10. -6These include removing variants that are above 100% confidence interval (CI) or that violate Hardy-Weinberg equilibrium (HWE), or that are missing calls in more than about 5% of samples from multiple human subjects, or a combination of these.
[0046] In some exemplary embodiments, the program includes instructions for removing low quality samples after the identified variants have been filtered. In some exemplary embodiments, low quality samples are those with a D-stat greater than 0.12 or less than 75% read coverage at 20x, or both.
[0047] In some exemplary embodiments, integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to augment the first homoeologous estimate.
[0048] In some exemplary embodiments, homoeologous estimation involves the calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. In some exemplary embodiments, the program includes instructions for filtering the variants after they have been phased according to population allele frequencies.
[0049] In some exemplary embodiments, filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of about 5 or less, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of about 5 or less, or an AB of about 10% or less, or a combination thereof.
[0050] In some exemplary embodiments, phasing variants according to population allele frequency comprises dividing a DNA sequence sample of a human subject into genomic segments of approximately equal size and with substantial segment overlap and breakpoints in intergenic regions.
[0051] In some exemplary embodiments, the potential CHMs are phased based on trio data, or parent-offspring data, or full sibling data, or distant relative data, or a combination thereof, or are phased based on minor allele counts (MACs), or are phased based on population allele frequencies, or a combination thereof.
[0052] In some exemplary embodiments, the program includes instructions for scoring the CHMs according to functional effect priority and selecting the CHMs with the highest functional effect priority score per gene per sample, thereby obtaining a set of medically important variants.
[0053] In some exemplary embodiments, the DNA sequence sample comprises an exome sequence. In some exemplary embodiments, the plurality of human subjects comprises more than 10K subjects.
[0054] In some exemplary embodiments, the D-stat of a low quality sample is determined by comparing the distribution of the actual allele balance of the sample to the distribution of the expected allele balance using a KS test.
[0055] In some exemplary embodiments, filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of less than about 10, or an AB of about 20% or less, or a combination thereof.
[0056] In another aspect, the disclosure provides a method for identifying de novo mutations (DNMs) in a population. In general, the method includes identifying variants in DNA sequence samples from a plurality of human subjects, establishing an ancestral superclass designation for the subjects based on the identified variants, generating a first homoeologous estimate for the subjects within the ancestral superclass, generating a second homoeologous estimate for the subjects independent of the ancestral superclass for the subjects, clustering the subjects into a first first degree family network based on one or more of the second homoeologous estimates, generating a third homoeologous estimate for the subjects within the first first degree family network, merging the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate, constructing nuclear families based on the combined homoeologous estimate, identifying variants in the nuclear families, assigning genotype likelihood scores to variants in samples from parents and children of the constructed nuclear family trios to calculate the probability that the variant is a de novo mutation, and independently, simply identifying variants in the child sample that are not present in either parent sample of the trio and calculating the probability that the variant is a de novo mutation, and then combining both probabilities, thereby obtaining an estimated homoeologous estimate. This involves generating a dataset of novo mutations.
[0057] In some exemplary embodiments, the method further comprises filtering the identified variants before an ancestral superclass assignment for the subject is established. In some exemplary embodiments, the method further comprises filtering the identified variants before the second homoeologous estimate for the subject is generated.
[0058] In some exemplary embodiments, filtering variants comprises removing variants with an alternative allele frequency of greater than about 10% across samples from multiple human subjects, or a p-value of greater than about 10. -6 These include removing variants that are above 100% confidence interval (CI) or that violate Hardy-Weinberg equilibrium (HWE), or that are missing calls in more than about 5% of samples from multiple human subjects, or a combination of these.
[0059] In some exemplary embodiments, the method includes removing low quality samples after the identified variants have been filtered. In some exemplary embodiments, low quality samples are those with a D-stat greater than 0.12 or less than 75% read coverage at 20x, or both.
[0060] In some exemplary embodiments, integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to augment the first homoeologous estimate.
[0061] In some exemplary embodiments, homoeologous estimation involves the calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. In some exemplary embodiments, the genotype likelihood scores are based on DNA sequence samples from multiple human subjects in multiple nuclear families.
[0062] In some exemplary embodiments, the method further comprises filtering the variants after the probability of the variant being a de novo mutation is calculated based on the genotype likelihood score.
[0063] In some exemplary embodiments, the method further comprises filtering the variants after the probability of the variant being a de novo mutation is calculated based on simply identifying variants in the child sample that are not present in either parent sample.
[0064] In some exemplary embodiments, filtering the variants comprises removing variants with a genotype quality (GQ) annotation of less than about 35 in the child sample, or variants with an alternative allele count (AC) of 10 or more in the sample, or variants with a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or variants with an allele balance (AB) of more than about 2% in either parent sample, or variants with an allele balance (AB) of less than about 15% in the child sample, or variants with an AB of more than about 90% in the child sample, or variants with alternative allele homozygosity in either parent sample, or combinations thereof.
[0065] In some exemplary embodiments, the method further comprises annotating the variants with quality control metrics. In some exemplary embodiments, the method further includes filtering variants based on the sample BAM file data after putative de novo mutations are identified based simply on identifying variants in the child sample that are not present in either parent sample.
[0066] In some exemplary embodiments, the method further comprises classifying the putative de novo mutation as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is greater than or equal to about 0.15 in the child sample.
[0067] In some exemplary embodiments, the method further comprises classifying the putative de novo mutation as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is less than or equal to about 0.02 in each parent sample.
[0068] In some exemplary embodiments, the method further comprises classifying the putative de novo mutation as a medium confidence de novo mutation if the mapping quality of the putative de novo mutation is not less than about 40.
[0069] In some exemplary embodiments, the method further comprises classifying the putative de novo mutation as a medium confidence de novo mutation if the quality by depth (QD) value of the putative de novo mutation is not less than about 2.
[0070] In some exemplary embodiments, the method further comprises classifying the putative de novo mutation as a medium confidence de novo mutation if the MAC of the putative de novo mutation is less than about 20 across all samples.
[0071] In some exemplary embodiments, the method further comprises classifying the putative de novo mutation as a medium confidence de novo mutation if the soft clipped reads of the putative de novo mutation are about 3 or less at the variant site in a carrier of the putative de novo mutation.
[0072] In some exemplary embodiments, the method further includes classifying the putative de novo mutation as a medium confidence de novo mutation if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs.
[0073] In some exemplary embodiments, the method classifies the putative de novo mutation as a medium-confidence de novo mutation if the allele balance (AB) of the putative de novo mutation is greater than or equal to about 0.15 in the child sample and less than or equal to about 0.02 in each parent sample, and if the mapping quality (MQ) of the putative de novo mutation is not less than about 40, and if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if the minor allele count (MAC) of the putative de novo mutation is less than about 20 across samples, and if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs. It further includes classifying the mutation as a novo mutation.
[0074] In some exemplary embodiments, the method further comprises classifying the medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality (GQ) annotation of about 90 or more in the parent samples, and if the medium-confidence de novo mutation has a read depth (DP) of about 10 or more in each parent sample, and if the medium-confidence de novo mutation has an alternative DP of about 7 or more in the child sample, and if the medium-confidence de novo mutation has a SNP with a QD of greater than about 3, and if the medium-confidence de novo mutation has an INDEL with a QD of greater than about 5.
[0075] In some exemplary embodiments, the DNA sequence sample comprises an exome sequence. In some exemplary embodiments, the plurality of human subjects comprises more than 10K subjects.
[0076] In some exemplary embodiments, the D-stat of a low quality sample is determined by comparing the distribution of the actual allele balance of the sample to the distribution of the expected allele balance using a KS test.
[0077] In another aspect, the disclosure provides a non-transitive computer-implemented method for identifying de novo mutations (DNMs) in a population. In general, the non-transitory computer-implemented method includes using a data processor of a computing device to identify variants in DNA sequence samples from a plurality of human subjects; using the data processor to establish an ancestral superclass designation of the subjects based on the identified variants; using the data processor to generate a first homoeologous estimate for the subjects within the ancestral superclass; using the data processor to generate a second homoeologous estimate for the subjects independent of the ancestral superclass of the subjects; using the data processor to cluster the subjects into a first first degree family network based on one or more of the second homoeologous estimates; using the data processor to generate a third homoeologous estimate for the subjects within the first first degree family network; using the data processor to combine the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; using the data processor to construct nuclear families based on the combined homoeologous estimate; using the data processor to identify variants in the nuclear families; assigning genotype likelihood scores to variants in samples from parents and children of the constructed nuclear family trios to determine whether the variants are related to the de The method includes using a data processor to calculate the probability that a variant is a de novo mutation and, independently, to simply identify variants in the child sample that are not present in either parent sample of the trio and calculate the probability that the variant is a de novo mutation, and then combine both probabilities, thereby forming a dataset of putative de novo mutations.
[0078] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to filter the identified variants before an ancestral superclass assignment for the subject is established.
[0079] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to filter the identified variants before a second homoeologous estimate for the subject is generated.
[0080] In some exemplary embodiments, filtering variants comprises removing variants with an alternative allele frequency of greater than about 10% across samples from multiple human subjects, or a p-value of greater than about 10. -6 The present invention includes removing variants that are above 100% confidence interval (CI) or violate Hardy-Weinberg equilibrium, or removing variants that are missing calls in more than about 5% of samples from multiple human subjects, or a combination of these.
[0081] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to remove low quality samples after the identified variants have been filtered.
[0082] In some exemplary embodiments, low quality samples are those with a D-stat greater than 0.12 or less than 75% read coverage at 20x, or both.
[0083] In some exemplary embodiments, integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to augment the first homoeologous estimate.
[0084] In some exemplary embodiments, homoeologous estimation involves the calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. In some exemplary embodiments, the genotype likelihood scores are based on DNA sequence samples from multiple human subjects in multiple nuclear families.
[0085] In some exemplary embodiments, the non-transitive computer-implemented method further comprises using a data processor to filter the variants after the probability of the variant being a de novo mutation is calculated based on the genotype likelihood score.
[0086] In some exemplary embodiments, the non-transitive computer-implemented method further comprises using a data processor to filter the variants after the probability of the variant being a de novo mutation has been calculated based on simply identifying variants in the child sample that are not present in either parent sample.
[0087] In some exemplary embodiments, filtering the variants comprises removing variants with a genotype quality (GQ) annotation of less than about 35 in the child sample, or variants with an alternative allele count (AC) of 10 or more across samples, or variants with a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or variants with an allele balance (AB) of more than about 2% in either parent sample, or variants with an allele balance (AB) of less than about 15% in the child sample, or variants with an AB of more than about 90% in the child sample, or variants with alternative allele homozygosity in either parent sample, or combinations thereof.
[0088] In some exemplary embodiments, the non-transitory computer-implemented method further comprises using the data processor to annotate the variants with quality control metrics.
[0089] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to filter variants based on the sample BAM file data after putative de novo mutations are identified based simply on identifying variants in the child sample that are not present in either parent sample.
[0090] In some exemplary embodiments, the non-transitive computer-implemented method further includes using the data processor to classify the putative de novo mutation as a medium confidence de novo mutation if an allele balance of the putative de novo mutation is greater than or equal to about 0.15 in the child sample.
[0091] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to classify the putative de novo mutation as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is less than or equal to about 0.02 in each parent sample.
[0092] In some exemplary embodiments, the non-transitive computer-implemented method further includes using the data processor to classify the putative de novo mutation as a medium confidence de novo mutation if the mapping quality of the putative de novo mutation is not less than about 40.
[0093] In some exemplary embodiments, the non-transitive computer-implemented method further includes using the data processor to classify the putative de novo mutation as a medium confidence de novo mutation if a quality by depth (QD) value of the putative de novo mutation is not less than about 2.
[0094] In some exemplary embodiments, the non-transitive computer-implemented method further includes using the data processor to classify the putative de novo mutation as a medium confidence de novo mutation if the MAC of the putative de novo mutation is less than about 20 across all samples.
[0095] In some exemplary embodiments, the non-transitive computer-implemented method further includes using the data processor to classify the putative de novo mutation as a medium confidence de novo mutation if a soft clip read of the putative de novo mutation is about 3 or less at the variant site in a carrier of the putative de novo mutation.
[0096] In some exemplary embodiments, the non-transitive computer-implemented method further includes using a data processor to classify the putative de novo mutation as a medium confidence de novo mutation if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs.
[0097] In some exemplary embodiments, the non-transitory computer-implemented method further includes using a data processor to classify the putative de novo mutation as a medium-confidence de novo mutation if the allele balance (AB) of the putative de novo mutation is greater than or equal to about 0.15 in the child sample and less than or equal to about 0.02 in each parent sample, and if the mapping quality (MQ) of the putative de novo mutation is not less than about 40, and if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if the minor allele count (MAC) of the putative de novo mutation is less than about 20 across all samples, and if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in a carrier of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs.
[0098] In some exemplary embodiments, the non-transitory computer-implemented method further includes using the data processor to classify the medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality (GQ) annotation of about 90 or more in the parent samples, and if the medium-confidence de novo mutation has a read depth (DP) of about 10 or more in each parent sample, and if the medium-confidence de novo mutation has an alternative DP of about 7 or more in the child sample, and if the medium-confidence de novo mutation has a SNP with a QD of greater than about 3, and if the medium-confidence de novo mutation has an INDEL with a QD of greater than about 5.
[0099] In some exemplary embodiments, the DNA sequence sample comprises an exome sequence. In some exemplary embodiments, the plurality of human subjects comprises more than 10K subjects.
[0100] In some exemplary embodiments, the D-stat of a low quality sample is determined by comparing the distribution of the actual allele balance of the sample to the distribution of the expected allele balance using a KS test.
[0101] In another aspect, the present disclosure provides a system, which may be used, for example, to perform the methods and non-transitory computer-implemented methods. The system generally includes a data processor, a memory coupled to the data processor, and a program stored in the memory, the program identifying variants in DNA sequence samples from a plurality of human subjects, establishing an ancestral superclass designation for the subjects based on the identified variants, generating a first homoeologous estimate for the subjects within the ancestral superclasses, generating a second homoeologous estimate for the subjects independent of the ancestral superclasses for the subjects, clustering the subjects into a first first degree family network based on one or more of the second homoeologous estimates, generating a third homoeologous estimate for the subjects within the first first degree family network, merging the first homoeologous estimate and the third homoeologous estimate to obtain a merged homoeologous estimate, constructing nuclear families based on the merged homoeologous estimate, identifying variants in the nuclear families, assigning genotype likelihood scores to variants in samples from parents and children of the constructed nuclear family trios to calculate a probability that the variant is a de novo mutation, and independently identifying variants in the child sample that are not present in either parent sample of the trio to determine whether the variant is a de novo mutation. The method includes instructions for calculating the probability of a mutation being a novo mutation and then combining both probabilities, thereby forming a dataset of putative de novo mutations.
[0102] In some exemplary embodiments, the program includes instructions for filtering the identified variants before an ancestral superclass assignment for the subject is established. In some exemplary embodiments, the program includes instructions for filtering the identified variants before a second homoeologous estimate for the subject is generated.
[0103] In some exemplary embodiments, filtering variants comprises removing variants with an alternative allele frequency of greater than about 10% across samples from multiple human subjects, or a p-value of greater than about 10. -6These include removing variants that are above 100% confidence interval (CI) or that violate Hardy-Weinberg equilibrium (HWE), or that are missing calls in more than about 5% of samples from multiple human subjects, or a combination of these.
[0104] In some exemplary embodiments, the program includes instructions for removing low quality samples after the identified variants have been filtered. In some exemplary embodiments, low quality samples are those with a D-stat greater than 0.12 or less than 75% read coverage at 20x, or both.
[0105] In some exemplary embodiments, integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to augment the first homoeologous estimate.
[0106] In some exemplary embodiments, homoeologous estimation involves the calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. In some exemplary embodiments, the genotype likelihood scores are based on DNA sequence samples from multiple human subjects in multiple nuclear families.
[0107] In some exemplary embodiments, the program includes instructions for filtering the variants after the probability that the variant is a de novo mutation is calculated based on the genotype likelihood score.
[0108] In some exemplary embodiments, the program includes instructions for filtering the variants after the probability that the variant is a de novo mutation has been calculated based on simply identifying variants in the child sample that are not present in either parent sample.
[0109] In some exemplary embodiments, filtering the variants comprises removing variants with a genotype quality (GQ) annotation of less than about 35 in the child sample, or variants with an alternative allele count (AC) of 10 or more across samples, or variants with a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or variants with an allele balance (AB) of more than about 2% in either parent sample, or variants with an allele balance (AB) of less than about 15% in the child sample, or variants with an AB of more than about 90% in the child sample, or variants with alternative allele homozygosity in either parent sample, or combinations thereof.
[0110] In some exemplary embodiments, the program includes instructions for annotating the variants with quality control metrics. In some exemplary embodiments, the program includes instructions for filtering variants based on the sample BAM file data after putative de novo mutations have been identified based simply on identifying variants in the child sample that are not present in either parent sample.
[0111] In some exemplary embodiments, the program includes instructions for classifying the putative de novo mutation as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is greater than or equal to about 0.15 in the child sample.
[0112] In some exemplary embodiments, the program includes instructions for classifying the putative de novo mutation as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is less than or equal to about 0.02 in each parent sample.
[0113] In some exemplary embodiments, the program comprises instructions for classifying the putative de novo mutation as a medium confidence de novo mutation if the mapping quality of the putative de novo mutation is not less than about 40.
[0114] In some exemplary embodiments, the program includes instructions for classifying the putative de novo mutation as a medium confidence de novo mutation if the quality by depth (QD) value of the putative de novo mutation is not less than about 2.
[0115] In some exemplary embodiments, the program includes instructions for classifying the putative de novo mutation as a medium confidence de novo mutation if the MAC of the putative de novo mutation is less than about 20 across all samples.
[0116] In some exemplary embodiments, the program includes instructions for classifying a putative de novo mutation as a medium confidence de novo mutation if a soft clip read of the putative de novo mutation is about 3 or less at a variant site in a carrier of the putative de novo mutation.
[0117] In some exemplary embodiments, the program includes instructions for classifying the putative de novo mutation as a medium confidence de novo mutation if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs.
[0118] In some exemplary embodiments, the program includes instructions for classifying the putative de novo mutation as a medium confidence de novo mutation if the allele balance (AB) of the putative de novo mutation is greater than or equal to about 15% in the child sample and less than or equal to about 2% in each parent sample, and if the mapping quality (MQ) of the putative de novo mutation is not less than about 40, and if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if the minor allele count (MAC) of the putative de novo mutation is less than about 20 across all samples, and if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs.
[0119] In some exemplary embodiments, the program includes instructions for classifying a medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality (GQ) annotation of about 90 or more in the parent samples, and if the medium-confidence de novo mutation has a read depth (DP) of about 10 or more in each parent sample, and if the medium-confidence de novo mutation has an alternative DP of about 7 or more in the child sample, and if the medium-confidence de novo mutation has a SNP with a QD of greater than about 3, and if the medium-confidence de novo mutation has an INDEL with a QD of greater than about 5.
[0120] In some exemplary embodiments, the DNA sequence sample comprises an exome sequence. In some exemplary embodiments, the plurality of human subjects comprises more than 10K subjects.
[0121] In some exemplary embodiments, the D-stat of a low quality sample is determined by comparing the distribution of the actual allele balance of the sample to the distribution of the expected allele balance using a KS test.
[0122] In some exemplary embodiments, the method, non-transitory computer-implemented method, or system includes assigning genotype likelihood scores to variants in samples from parents and children of the constructed nuclear family trio to calculate the probability that the variant is a de novo mutation, selecting variants with a significantly high probability that the variant is a de novo mutation, and independently simply identifying variants called in the child sample that are not called in either parent sample of the trio, and then combining the two sets of de novo mutations, thereby forming a dataset of putative de novo mutations.
[0123] In another aspect, the present disclosure provides a predictive model of relatedness in human populations. The predictive model may be prepared by a process including establishing a first population data set, performing a 120-year burn-in phase to establish a second population data set, and modifying the second population data set by performing the following steps: (a) moving individuals in the second population data set into age pools according to the age of the individuals; (b) randomly selecting pairs of single men and women from the single men and women in the second population data set who are more distantly related than first cousins, marrying them by a marriage specified by an age parameter, and selecting pairs until a number specified by a marriage rate parameter is reached; (c) divorcing couples at a specified divorce rate. Couples are randomly selected from the second population data set and marked as single upon divorce; (d) randomly selecting pairs of single men and women or married couples from the second population data set at a specified ratio and allowing them to reproduce according to a specified birth rate until a target number of successful conceptions is reached. Parents are restricted to being more distantly related than first cousins, and all individuals in the second population dataset are restricted to having one child per year. (e) Allowing individuals in the second population dataset to die at a specific mortality rate with specific mortality and age parameters. (f) Allowing individuals to move in and out of the second population dataset, thereby maintaining the population age and sex distribution and the proportion of those of reproductive age who are married in the second population dataset. (g) Allowing individuals to move within the second population dataset, whereby individuals are randomly selected from subpopulations and randomly assigned to another population, if any, until a specified migration rate between the subpopulations is achieved. Steps (a)-(g) are repeated for a predetermined number of years with one-year intervals, where the steps are applied to the population dataset obtained from the previous iteration.
[0124] In some exemplary embodiments, establishing the dataset of the first population further includes specifying a number of subpopulations and sizes. In some exemplary embodiments, establishing the first population dataset further includes assigning ages to individuals in the first population dataset between zero and the maximum fertility age.
[0125] In some exemplary embodiments, the maximum age for conception is 49 years. In some exemplary embodiments, performing the burn-in phase further comprises keeping the number of births and deaths of individuals in the second population dataset equal and keeping the net migration rate of individuals at zero.
[0126] In some exemplary embodiments, performing the burn-in phase further includes moving the dataset of the second population of individuals from the young pool to the mating pool when the individuals age above the minimum fertility age, and moving individuals from the mating pool to the old pool when the individuals age above the maximum fertility age, as well as removing individuals from all age groups if they emigrate or die.
[0127] In some exemplary embodiments, the minimum age for birth is 15 years and the maximum age for conception is 49 years. In another aspect, the present disclosure provides a method using a predictive model, where the verification of the individual is done randomly.
[0128] In another aspect, the present disclosure provides a method of using a predictive model, wherein identifying the individual is performed in a clustered manner. In some exemplary embodiments, verifying the individuals further includes collecting kinship data and association statistics regarding the verified individuals, including first-degree or second-degree relationships, or both, between the verified individuals. Effect of the Invention
[0129] As described above, according to the present invention, a system and method for utilizing relatedness in genome data analysis are provided. [Brief description of the drawings]
[0130] [Figure 1] 1 is a flow chart of an exemplary embodiment in which phasing / identification of genetic variants in a population is performed by leveraging population relatedness. [Diagram 2] 1 is a flow chart of an exemplary embodiment in which compound heterozygous mutations (CHMs) are identified in a population. [Diagram 3] 1 is a flow chart of an exemplary embodiment in which de novo mutations (DNMs) in a population are identified. [Figure 4] 1 is a flowchart of a method for creating a predictive model of relatedness in a human population, according to an exemplary embodiment. [Figure 5A] 1 depicts a flowchart of an exemplary embodiment in which homoeologousity is determined. [Figure 5B] 1 depicts a flowchart of an exemplary embodiment in which homoeologousity is determined. [Figure 5C] 1 depicts a flowchart of an exemplary embodiment in which homoeologousity is determined. [Figure 5D] 1 depicts a flowchart of an exemplary embodiment in which homoeologousity is determined. [Figure 6A] 1 depicts a flow chart of an exemplary embodiment in which compound heterozygous mutations (CHMs) are identified / phased in a population. [Figure 6B] 1 depicts a flow chart of an exemplary embodiment in which compound heterozygous mutations (CHMs) are identified / phased in a population. [Figure 6C] 1 depicts a flow chart of an exemplary embodiment in which compound heterozygous mutations (CHMs) are identified / phased in a population. [Figure 7A]1 depicts a flow chart of an exemplary embodiment in which de novo mutations (DNMs) are identified in a population. DNM calling, filtering, and confidence ranking workflow. GQ=genotype quality. MAC is DiscovEHR minor allele count. DP=read depth of DNM site. AD=alternative allele depth. AB=alternative allele balance. MQ=mapping quality. QD=quality by depth of junction called DiscovEHR pVCF. Homopolymer INDELs are INDELs with four or more consecutive base pairs of the same nucleotide. Blacklisted genes include PDE4DIP, PRAMEF1, PABPC3, NBPF10, NBPF14, olfactory genes (OR*), MUC genes (MUC*), and HLA genes (HLA-*). [Figure 7B] 1 depicts a flow chart of an exemplary embodiment in which de novo mutations (DNMs) are identified in a population. DNM calling, filtering, and confidence ranking workflow. GQ=genotype quality. MAC is DiscovEHR minor allele count. DP=read depth of DNM site. AD=alternative allele depth. AB=alternative allele balance. MQ=mapping quality. QD=quality by depth of junction called DiscovEHR pVCF. Homopolymer INDELs are INDELs with four or more consecutive base pairs of the same nucleotide. Blacklisted genes include PDE4DIP, PRAMEF1, PABPC3, NBPF10, NBPF14, olfactory genes (OR*), MUC genes (MUC*), and HLA genes (HLA-*). [Figure 8]A, B, C and D represent a scheme providing an overview of different types of population-based genomic studies and corresponding sampling methods, showing that strict ascertainment increases family structure and impacts the statistical analysis approach that needs to be used. Panel A shows a schematic diagram of (1) conventional population-based genomic studies (gray boxes), (2) medical population-based genomic (HPG) studies (green boxes), and (3) family-based genomic studies (yellow boxes). Panel B shows a line graph of family structure in the three aforementioned ascertainment approaches. Panel C shows a scatter plot of family structure in the three aforementioned ascertainment approaches (lines indicate first- and second-degree pairwise relationships ascertained from the three aforementioned ascertainment approaches). Panel D shows the statistical analysis approaches classified into four categories based on the level of family structure. [Figure 9] FIG. 2 is a flow chart of an exemplary embodiment outlining the cascade analysis performed to determine the phase of potential compound heterozygous mutations (pCHMs) in the analyzed dataset (DiscovEHR dataset; see Examples). [Figure 10] 1 is an exemplary operating environment. [Figure 11] 1 illustrates several system components configured to perform the disclosed methods. [Figure 12A] A, B, C, and D show relatedness seen in the first 61K sequenced individuals from the DiscovEHR cohort according to an exemplary embodiment. Panel A shows a plot of IBD0 vs. IBD1. [Figure 12B] Figure 1 shows relatedness found in the first 61K sequenced individuals from the DiscovEHR cohort according to an exemplary embodiment. Panel B shows a histogram plotting the size distribution of first-degree family networks within the analyzed cohort. [Figure 12C] Figure 1 shows relatedness seen in the first 61K sequenced individuals from the DiscovEHR cohort according to an exemplary embodiment. Panel C shows a pedigree of a first-degree family network including the 25 sequenced individuals, reconstructed from pairwise IBD inference. [Figure 12D] FIG. 1 shows relatedness found in the first 61K sequenced individuals from the DiscovEHR cohort according to an exemplary embodiment. Panel D shows a scheme representing the largest second-degree family network of 7,084 individuals. [Figure 13A] FIG. 1 shows the cumulative relatedness within a DiscovEHR cohort as a function of the number of ascertained individuals in a dataset ascertained by an exemplary embodiment. [Figure 13B] FIG. 1 shows the cumulative relatedness within a DiscovEHR cohort as a function of the number of ascertained individuals in a dataset ascertained by an exemplary embodiment. [Figure 14A] Figure 1 shows a comparison of ascertainment of first-degree relatives among 61K DiscovEHR participants with random ascertainment of a simulated population according to an exemplary embodiment. Panel A shows ascertainment of first-degree relative pairs. [Figure 14B] Figure 1 shows an example embodiment of ascertaining first-degree relatives among 61K DiscovEHR participants compared to random ascertainment of a simulated population. Panel B shows ascertainment of the number of individuals with multiple first-degree relatives. [Figure 15A] Figure 1 shows ascertainment fit to the accumulation of first-degree relatedness within a simulated population and a validated DiscoverEHR cohort according to an exemplary embodiment. Panel A shows the accumulation of pairs of first-degree relatives. [Figure 15B] Figure 1 shows ascertainment matched to the accumulation of first-degree relatedness within a simulated population and a validated DiscoverEHR cohort, according to an exemplary embodiment. Panel B shows the proportion of confirmed participants with one or more first-degree relatives. [Figure 15C] Figure 1 shows ascertainment fit to the cumulative accumulation of first-degree relatedness within a simulated population and a validated DiscoverEHR cohort, according to an exemplary embodiment. Panel C shows simulated ascertainment projections along with upper and lower bounds on the number of first-degree relationships. [Figure 15D]Figure 1 shows ascertainment fit to cumulative first-degree relatedness within the simulated population and ascertained DiscoverEHR cohorts according to an exemplary embodiment. Panel D shows simulated projections along with upper and lower bounds on the proportion of ascertained participants with one or more first-degree relatives. [Figure 16A] FIG. 1 shows the first 92K sequenced individuals from the validated and expanded DiscovEHR cohort according to an exemplary embodiment. Panel A shows a plot of IBD0 vs. IBD1. [Figure 16B] Panel A shows the first 92K sequenced individuals from the validated and expanded DiscovEHR cohort according to an exemplary embodiment. Panel B shows a histogram plotting the size distribution of first-degree family networks within the analyzed cohort. [Figure 16C] Panel C shows the first 92K sequenced individuals from the validated and expanded DiscovEHR cohort according to an exemplary embodiment. Panel D shows a pedigree of the first-degree family network including the 25 sequenced individuals, reconstructed from pairwise IBD inference. [Figure 16D] Panel D shows the first 92K sequenced individuals from the validated and expanded DiscovEHR cohort according to an exemplary embodiment. Panel E shows a scheme representing the largest second-degree family network of 7,084 individuals. [Figure 17A] Figure 1 shows a comparison between ascertainment of first-degree relatives among 92K extended DiscovEHR participants compared to random ascertainment of a simulated population according to an exemplary embodiment. Panel A shows ascertainment of first-degree relative pairs. [Figure 17B] Figure 1 shows a comparison between ascertainment of first-degree relatives among 92K extended DiscovEHR participants compared to random ascertainment of a simulated population according to an exemplary embodiment. Panel B shows ascertainment of the number of individuals with multiple first-degree relatives. [Figure 18A]Figure 1 shows population and ascertainment matching simulating the accumulation of first-degree relatedness in an extended DiscovEHR cohort according to an exemplary embodiment. Panel A shows the accumulation of pairs of first-degree relatives. [Figure 18B] Figure 1 shows population and ascertainment matching simulating the accumulation of first-degree relatedness in an expanded DiscovEHR cohort according to an exemplary embodiment. Panel B shows the proportion of confirmed participants with one or more first-degree relatives. [Figure 18C] Figure 1 shows population and ascertainment matching simulating the accumulation of first-degree relatedness in the extended DiscovEHR cohort according to an exemplary embodiment. Panel C shows simulated ascertainment projections along with upper and lower bounds on the number of first-degree relationships. [Figure 18D] Figure 1 shows population and ascertainment matching for simulated accumulation of first-degree relatedness in the expanded DiscovEHR cohort, according to an exemplary embodiment. Panel D shows simulated projections along with upper and lower bounds on the proportion of ascertained participants with one or more first-degree relatives. [Figure 19A] Figure 1 shows simulated population and ascertainment matching of the pooling of first- and second-degree relatedness within a validated DiscovEHR cohort according to an exemplary embodiment. Panel A shows the pooling of first- and second-degree related pairs. [Figure 19B] Figure 1 shows population and ascertainment matching simulating the accumulation of first- and second-degree relatedness within a validated DiscovEHR cohort, according to an exemplary embodiment. Panel B shows the proportion of confirmed participants with one or more first- and second-degree relatives. [Figure 19C] Figure 1 shows population and ascertainment matching simulating the accumulation of first- and second-degree relatedness within a validated DiscovEHR cohort, according to an exemplary embodiment. Panel C shows simulated ascertainment projections along with upper and lower bounds on the number of first- and second-degree relationships. [Figure 19D]Figure 1 shows population and ascertainment matching simulating the accumulation of first- and second-degree relatedness within a validated DiscovEHR cohort, according to an exemplary embodiment. Panel D shows simulated projections along with upper and lower bounds on the proportion of validated participants with one or more first- or second-degree relatives. [Figure 20A] Figure 1 shows simulated population and ascertainment matching of the pooling of first- and second-degree relatedness within a validated extended DiscovEHR cohort according to an exemplary embodiment. Panel A shows the pooling of first- and second-degree related pairs. [Figure 20B] Figure 1 shows population and ascertainment matching simulating the accumulation of first- and second-degree relatedness within a validated expanded DiscovEHR cohort according to an exemplary embodiment. Panel B shows the proportion of confirmed participants with one or more first- and second-degree relatives. [Figure 20C] Figure 1 shows population and ascertainment matching simulating the accumulation of first- and second-degree relatedness within a validated expanded DiscovEHR cohort according to an exemplary embodiment. Panel C shows simulated ascertainment projections along with upper and lower bounds on the number of first- and second-degree relationships. [Figure 20D] Figure 1 shows population and ascertainment matching simulating the accumulation of first- and second-degree relatedness within the ascertained expanded DiscovEHR cohort according to an exemplary embodiment. Panel D shows simulated projections along with upper and lower bounds on the proportion of ascertained participants with one or more first- or second-degree relatives. [Figure 21A] Figure 1 shows the number of compound heterozygous mutations (CHM) identified per individual in the DiscovEHR cohort according to an exemplary embodiment. Panel A shows the number of CHMs per individual in the DiscovEHR cohort. [Figure 21B] Figure 1 shows the number of compound heterozygous mutations (CHM) identified per gene in the DiscovEHR cohort according to an exemplary embodiment. Panel B shows the number of CHMs per gene in the DiscovEHR cohort. [Figure 21C]Figure 1 shows the number of de novo mutations (DNMs) identified per individual in the DiscovEHR cohort according to an exemplary embodiment. Panel C shows the distribution of the number of exonic high-confidence DNMs among trio children in the DiscovEHR cohort. [Figure 21D] FIG. 1 shows the number of de novo mutations (DNMs) identified per gene in the DiscovEHR cohort, according to an exemplary embodiment. Panel D shows several non-synonymous DNMs per gene. [Figure 22] 1 is a chart showing the range of genomic distances between phased compound heterozygous variant (CHM) variants identified for the DiscovEHR dataset, according to an exemplary embodiment. [Figure 23] Reconstructed pedigrees from the DiscovEHR cohort showing segregation of known disease-causing variants, including variants for (A) aortic aneurysm, (B) long QT syndrome, and (C) thyroid cancer. [Figure 24] Pedigrees reconstructed from sequenced DiscovEHRs including 22 / 29 carriers of the tandem duplication in LDLR and from 10 unaffected related (first- or second-degree) individuals from the sequenced cohort. [Diagram 25] 1 is a decision cascade of an exemplary embodiment for determining the phase of potential compound heterozygous mutations (pCHM) among 92KDiscover participants. [Figure 26A] Figure 1 shows the resulting expanded DiscovEHR cohort for identified compound heterozygous mutations (CHMs) and de novo mutations (DNMs) according to an exemplary embodiment. Panel A shows the distribution of the number of CHMs per individual in the DiscovEHR cohort. [Figure 26B] FIG. 1 shows the resulting expanded DiscovEHR cohort for identified compound heterozygous mutations (CHMs) according to an exemplary embodiment. Panel B shows the distribution of the number of CHMs per gene. [Figure 26C]Figure 1 shows the resulting expanded DiscovEHR cohort for de novo mutations (DNMs) according to an exemplary embodiment. Panel C shows the distribution of 3,415 exonic high- and medium-confidence DNMs among trio children in the DiscovEHR cohort. [Figure 26D] FIG. 1 shows the resulting expanded DiscovEHR cohort for de novo mutations (DNMs) according to an exemplary embodiment. Panel D shows the distribution of nonsynonymous DNMs across one or more of 2,802 genes. [Figure 27] 1 is a chart showing the range of genomic distances between phased compound heterozygous variant (CHM) variants identified for the extended DiscovEHR, according to an exemplary embodiment. [Figure 28] FIG. 13 is a cohort profile showing the number of family trios, family trios with parental ages, probands with one or more exonic DNMs, exonic DNMs, medium / high confidence DNMs, single nucleotide DNMs, medium / high confidence variants, and random variants identified in the extended DiscovEHR dataset, according to an exemplary embodiment. [Figure 29] Figure 1 shows the number of DNMs identified per confidence level and per person in the extended DiscovEHR cohort according to an exemplary embodiment. Panel A shows the distribution of DNM numbers per confidence level in the extended DiscovEHR cohort. Panel B shows the distribution of DNM numbers per individual in the extended DiscovEHR cohort in which DNMs were identified according to an exemplary embodiment. [Diagram 30] 1 is a histogram plotting the distribution of DNM numbers identified per functional effect DNM in the expanded DiscovEHR cohort, according to an exemplary embodiment. [Diagram 31] 1 is a histogram plotting the distribution of DNM numbers identified by type of DNM (transitions, transversions, and indels) in the extended DiscovEHR cohort, according to an exemplary embodiment. [Diagram 32]1 is a histogram plotting the distribution of DNM numbers identified by single nucleotide DNM type (→) in the expanded DiscovEHR cohort, according to an exemplary embodiment. [Diagram 33] 1 is a histogram plotting the distribution of DNM numbers identified per 10M exonic base pairs per chromosome in the extended DiscovEHR cohort, according to an exemplary embodiment. [Diagram 34] FIG. 1 is a bar graph plotting the distribution of the percentage of DNMs or randomly selected variants occurring in regions of the genome known to be enriched for CG dinucleotides (conventionally written as CpG, where "p" represents a phosphate between two bases) in the expanded DiscovEHR cohort, according to an exemplary embodiment. [Diagram 35] Shown is an image of pedigree prediction reconstructed from the sequenced expanded DiscovEHR cohort including 25 / 37 carriers of a de novo FH-causing tandem duplication in LDLR and 20 non-carrier, related (first- or second-degree) individuals. [Figure 36A] FIG. 1 shows the relationship between paternal and maternal age at birth and the number of exonic DNMs identified in children in the DiscovEHR cohort, according to an exemplary embodiment. [Figure 36B] FIG. 1 shows the relationship between paternal and maternal age at birth and the number of exonic DNMs identified in children in the DiscovEHR cohort, according to an exemplary embodiment. [Figure 37] 1 is a chart showing the correlation between maternal and paternal age at the child's birth and DNMs identified in the child in the DiscovEHR cohort, according to an exemplary embodiment. [Figure 38] 1 is a histogram plotting pathogenicity predictions of DNMs and random variants identified in the extended DiscovEHR cohort, according to an exemplary embodiment. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
[0131] The term "a" should be understood to mean "at least one." The terms "about" and "approximately" should be understood to allow for standard variations that would be understood by one of ordinary skill in the art. When ranges are provided, the endpoints are included.
[0132] Previous large-scale human genome studies typically collected human samples across many different geographic regions and / or health care systems and combined them to generate cohorts for analysis. In such cohorts, the total number of individuals sampled is often large, but the extent of relatedness and family structure in these cohorts has tended to be relatively small. Many statistical methods commonly used in the context of genomic analysis, including association studies and principal component analysis, require that all samples are unrelated. Otherwise, the statistical power of these tests will be biased, resulting in inflated p-values and false positive findings (Figure 8) (Kang et al. (2010), Nature Publishing Group 42, 348-354; Sun and Dimitromanolakis (2012), Methods Mol. Biol. 850, 47-57; Devlin and Roeder (1999), Biometrics 55, 997-104; and Voight and Pritchard (2005), PLoS Genet 1, e32-10).
[0133] Removing family structure from a dataset is a viable option when closely related samples are rare in the dataset (Lek, et al. (2016), Nature Publishing Group 536, 285-291; Fuchsberger et al. (2016), Nature Publishing Group 536, 41-47; Locke et al. (2015), Nature 518, 197-206; and Surendran et al. (2016) Nat Genet 48, 1151-1161). It is also a feasible option when a subset of unrelated data is sufficient for statistical analysis, such as calculating principal components (PCs) and then projecting the remaining samples onto these PCs (Dewey et al. (2016), Science 354, aaf6814-aaf6814). Several methods exist to help researchers retain the largest possible set of unrelated individuals (Staples at al. (2013), Genet. Epidemiol. 37, 136-141; Chang at al. (2015), Gigascience 4, 7). Unfortunately, removing related individuals not only reduces sample size, but also discards useful relationship information. Indeed, such loss of information is unacceptable for many analyses, even when the dataset has a moderate family structure.
[0134] The present disclosure is based, at least in part, on the recognition that information regarding family and ancestry structure and relatedness in datasets of genomic samples from multiple subjects is useful because it provides numerous analytical opportunities that allow, inter alia, to investigate links between rare genetic variations (e.g., compound heterozygous mutations and / or de novo mutations) and disease.
[0135] The present disclosure is also based, at least in part, on the recognition that genome-wide homozygosity (IBD) estimates are a superior metric for quantifying the level of relatedness within a dataset of multi-subject genomic samples, and between two sets of individuals.
[0136] Several statistical methods have been developed to model accurate pairwise relationships. For example, genome-wide association studies using mixed models are more powerful and outperform methods that do not model confounding relatedness (Kang et al. (2010), Nature Publishing Group 42, 348-354; Zhang et al. (2010), Nat Genet 42, 355-360; Yang et al. (2014), Nat Genet 46, 100-106; Kirkpatrick and Bouchard-Cote (2016), arXiv q-bio.QM), but mixed models do not fully exploit the information contained in the family structure and do not scale practically to datasets containing hundreds of thousands of samples and hundreds to thousands of phenotypes. Pairwise relationships are also sometimes used for pedigree-free QTL linkage analysis (Day-Williams et al. (2011), Genet. Epidemiol. 35, 360-370). Further software packages exist to model population and family structure, for pairwise relationship estimation (PCrelate) (Conomos et al. (2016), Am. J. Hum. Genet. 98, 127-148) and principal component analysis (PC-AiR) (Conomos et al. (2015), Genet. Epidemiol. 39, 276-293).
[0137] In contrast to traditional genome-wide association studies, recent and future large-scale genomic studies, such as those embodied in the present disclosure, sample tens to hundreds of thousands of participants from distinct geographic regions. As a result, these studies will identify a much larger proportion of people from the same geographic region, and thus family and ancestry structure within the sample dataset, identifying rare variants that segregate in families that are not appreciated in traditional population-wide association studies.
[0138] For several reasons, data from such large-scale genomic studies are rich in family structure and distant hidden relatedness. First, these studies sample extensively from a particular geographic region, for example through a health care system population, and the number of pairs of related individuals identified increases combinatorially as more samples are identified from a single population (Figure 8A). Second, families living in the same geographic region are likely to be served by the same physicians in the same health care system due to common insurance coverage and convenience. Third, common genetic and environmental factors may increase the frequency of medical collaboration targeted to certain families. In populations with low migration rates, both family structure and distant hidden relatedness are even more pronounced (Henn et al. (2012), PLoS ONE 7, e34267). The influence of family structure can be observed through the effect of sampling methods on linkage, pedigree-based analysis, IBD modeling, and analysis of non-related individuals (Figure 8, Panel D). "Linkage" refers to traditional linkage analysis using one or more informative pedigrees. "Pedigree-based analyses" refer to statistical methods beyond linkage that use pedigree structure within a larger cohort that includes unrelated individuals. "IBD modeling" refers to analyses that model pairwise relationships between individuals without using overall pedigree structure. "Unrelated-relative analyses" refer to analyses that assume all individuals in a cohort are unrelated.
[0139] In this disclosure, we focus on family structure and use both real and simulated data to reveal high levels of family structure. One of the improvements of this disclosure is to identify and / or phase compound heterozygous mutations (CHM) and / or de novo mutations (DNM) more accurately and reliably than traditional approaches (see data disclosed in the Examples section).
[0140] Thus, the present disclosure provides a method for phasing genetic variants in an ascertained population by leveraging the relatedness of the populations. A flow outlining an exemplary phasing method is provided in FIG.
[0141] This method can be applied to various genetic variants in various populations.Non-limiting examples of the types of genetic variants that can be evaluated include point mutations, insertions, deletions, inversions, duplications, and multimerization.Non-limiting examples of types of populations include single medical network populations, multiple medical network populations, racially, culturally, or socially homogeneous or heterogeneous populations, mixed age populations or homogeneous populations in terms of age, geographically concentrated or geographically dispersed populations, or combinations thereof.Non-limiting examples of the means that can obtain genetic variants include the following steps:
[0142] - Sample preparation and sequencing (Dewey et al. (2016), Science 354, aaf6814-1~aaf6814-10). Upon completion of sequencing, raw data from each sequencing run is collected in local buffer storage and uploaded to the DNAnexus platform (Reid et al. (2014); BMC Bioinformatics 15, 30) for automated analysis.
[0143] - Sample-level read files were generated with CASAVA (Illumina Inc., San Diego, CA) and aligned to GRCh38 (Li and Durbin (2009); Bioinformatics 25, 1754-176; Li (2013); arXiv q-bio.GN) with BWA-mem.
[0144] - The resulting BAM files are processed using GATK (McKenna et al. (2010); Genome Res. 20, 1297-1303) and Picard to sort, mark duplicates, and perform local alignment corrections of reads around putative indels.
[0145] - Sequenced variants are annotated with snpEFF (Cingolani et al. (2012); Fly (Austin) 6, 80-92) using Ensembl85 gene definitions to reveal functional effects on transcripts and genes.
[0146] It will be understood that the method is not limited to any of the steps described above and obtaining sequence variants may be accomplished by any suitable means. Figure 1 is a flow chart of an exemplary embodiment in which genetic variant phasing / identification in a population is performed by leveraging population relatedness. Low quality sequence variants may be removed in step 1 from a dataset of nucleic acid sequence samples obtained from multiple human subjects by any suitable means. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7) and those disclosed in the present Examples.
[0147] Ancestor superclass assignments for each of one or more of the samples may be established in step 2 by any suitable means. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7) and those disclosed herein.
[0148] Low quality samples may be removed from the dataset in step 3 by any suitable means. Non-limiting examples of such means include those disclosed in Dewey et al. (2016), Science 354, aaf6814-1 to aaf6814-10, and those disclosed in the present Examples.
[0149] An initial homozygosity estimate for subjects within ancestral superclasses may be generated in step 4 by any suitable means. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7), and those disclosed herein.
[0150] A second homoeologous estimate of the subject, independent of the subject's ancestral superclass, may be generated in steps 5 and 6, and the subject may be clustered into a first first-degree family network based on one or more of the second homoeologous estimates by any suitable means. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7), and those disclosed herein.
[0151] A third homozygosity estimate for the subject within the first-degree family network may be generated in step 7 by any suitable means. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7), and those disclosed in the present Examples.
[0152] The first homoeologous estimate and the third homoeologous estimate may be combined in step 8 by any suitable means to obtain a combined homoeologous estimate. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7), and those disclosed herein.
[0153] A second first-degree family network for the subject based on the pooled homoeologous estimates may be constructed in step 9 by any suitable means. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7), and those disclosed in the present Examples.
[0154] At step 10, the variants may be phased as being compound heterozygous mutations (CHM) or not CHM according to the combined homozygosity estimates and the second first-degree family network by any suitable means, or the variants may be identified as being de novo mutations (DNM) according to the combined homozygosity estimates and the second first-degree family network by any suitable means. Non-limiting examples of such means include those disclosed in Figures 6 and 7 and in the present Examples.
[0155] 5-7 show basic operating logic for purposes of illustration, not limitation, of a method for generating inter-ancestry (IBD) estimates and using the IBD estimates to phase genetic variants as compound heterozygous mutations (CHMs) or potential compound heterozygous mutations (pCHMs), or de novo mutations (DNMs). The programs (e.g., EAGLE, PLINK, etc.) identified in the logic are exemplary of the steps for which they are identified, but it is understood that such programs are not the only way to perform such steps.
[0156] Phasing variants as compound heterozygous mutations (CHM) may include (1) phasing variants according to population allele frequency; (2) removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both; and removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of about 5 or less, or an alternative allele balance (AB) of about 10% or less, or a combination thereof; and removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of about 5 or less, or an AB of about 10% or less, or a combination thereof; (3) selecting remaining variants as potential compound heterozygous mutations (pCHM), where one or more pairs of variants exist in the same sample and in the same gene; and (4) phasing pCHM as cis pCHM or trans pCHM, and then classifying pCHM phased as trans pCHM as CHM. Phasing variants according to population allele frequency can be facilitated by any suitable means, including but not limited to EAGLE (Loh et al. (2016), Nat Genet 48, 1443-1448). Variants that do not meet certain selection criteria can be removed, and the remaining variants are selected as potential compound heterozygous mutations, which are phased by any suitable means, including those described in this example. These exemplary embodiments are also shown in FIG. 6.
[0157] Phasing variants as compound heterozygous mutations may include removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of about 7 or less, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of about 10 or less, or an AB of about 20% or less, or a combination thereof. These steps may be performed as described elsewhere herein, except that the exclusion parameters are set to a more stringent level.
[0158] In some exemplary embodiments, the method further comprises: (1) assigning scores to the CHMs according to functional effect priority; and (2) selecting the CHM with the highest functional effect priority score per gene per sample, such that when a human has multiple CHMs within the same gene, the CHM most likely to result in protein function inhibition is identified. These steps may be performed by any suitable means, including but not limited to SIFT (Loh et al. (2016); Nat Genet 48, 1443-1448) (damaging), PolyPhen2 HDIV45 (damaging and possibly damaging), PolyPhen2 HVAR (damaging and possibly damaging), LRT46 (deleterious), and MutationTaster (Schwarz et al. (2014); Nat. Methods 11, 361-362) (disease causing automatic and disease causing).
[0159] Phasing a variant as a de novo mutation includes (1) identifying the variant in the second first-degree family network and the trio samples; (2) assigning genotype likelihood scores to the variant in the parent samples and the corresponding child samples in the trio, calculating the probability that the variant is a de novo mutation, and identifying the variant as a putative de novo mutation if the calculated probability is statistically significant; (3) identifying the variant in the child sample of the trio and identifying the variant as a putative de novo mutation if the variant is not present in any parent sample of the trio; (4) identifying a putative de novo mutation with a genotype quality (GQ) annotation of less than about 35 in the child sample, or a putative de novo mutation with an alternative allele count (AC) of 10 or more across samples, or a putative de novo mutation with a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or a putative de novo mutation with an allele balance (AB) of more than about 2% in either parent sample. (4) filtering the identified putative de novo mutations by removing putative de novo mutations that have novo mutations, or putative de novo mutations with less than about 15% allele balance (AB) in the child sample, or putative de novo mutations with more than about 90% AB in the child sample, or putative de novo mutations with alternative allele homozygosity in either parent sample, or combinations thereof; and (5) combining the filtered identified putative de novo mutations, thereby forming a putative de novo mutation dataset. These steps may be performed by any suitable means, including those described in this example. These exemplary embodiments are also shown in FIG. 7.
[0160] In some exemplary embodiments, the method further comprises classifying the putative de novo mutation in the putative de novo mutation dataset as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is greater than or equal to about 0.15 in the child sample and less than or equal to about 0.02 in each parent sample, and if the mapping quality of the putative de novo mutation is not less than about 40, and if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if the MAC of the putative de novo mutation is less than about 20 across all samples, and if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with a monopolymer run greater than about 4. In some exemplary embodiments, the method further includes classifying the medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality annotation of about 90 or more in the parent samples, and if the medium-confidence de novo mutation has a read depth of about 10 or more in each parent sample, and if the medium-confidence de novo mutation has an alternative read depth of about 7 or more in the child sample, and if the medium-confidence de novo mutation has a SNP of more than about 3 QD, and if the medium-confidence de novo mutation has an INDEL of more than about 5 QD. Either exemplary embodiment may be performed in any manner, including but not limited to those disclosed in the present examples.
[0161] The present disclosure also provides a method for identifying compound heterozygous mutations (CHM) in a population. A flow chart showing an example of a method for identifying CHM is provided in FIG. The method may be applied to any type of DNA sequence sample from any type of human subject obtained by any means. Non-limiting examples of variants include point mutations, insertions, deletions, inversions, duplications, and multimerizations. Non-limiting examples of types of human subjects include human subjects from a single medical network population, a multi-medical network population, a racially, culturally, or socially homogenous or heterogeneous population, a mixed age population or a homogenous population with respect to age, a geographically concentrated population or a geographically dispersed population, or a combination thereof. DNA sequence samples may be obtained in any of a number of ways, including but not limited to those disclosed in Dewey et al. (2016), Science 354, aaf6814-1 to aaf6814-10.
[0162] In some exemplary embodiments, the DNA sequence sample comprises an exome sequence. Exome DNA can be isolated by any of the commonly used methods or as described in Dewey et al. (2016), Science 354, aaf6814-1 to aaf6814-10.
[0163] Variants in DNA sequence samples from multiple human subjects may be identified in step 11 by any suitable means. Non-limiting examples of means by which variants may be identified include the following steps:
[0164] Upon completion of sequencing, raw data from each sequencing run is collected in local buffer storage and uploaded to the DNAnexus platform (Reid et al. (2014); BMC Bioinformatics 15, 30) for automated analysis.
[0165] - Sample-level read files were generated with CASAVA software (Illumina Inc., San Diego, CA) and aligned to GRCh38 (Li and Durbin (2009); Bioinformatics 25, 1754-176; Li (2013); arXiv q-bio.GN) with BWA-mem.
[0166] - The resulting BAM files are processed using GATK (McKenna et al. (2010); Genome Res. 20, 1297-1303) and Picard to sort, mark duplicates, and perform local alignment corrections of reads around putative indels.
[0167] - Sequenced variants are annotated with snpEFF (Cingolani et al. (2012); Fly (Austin) 6, 80-92) using Ensembl85 gene definitions to reveal functional effects on transcripts and genes.
[0168] It will be understood that the method is not limited to any of the steps described above and obtaining sequence variants may be accomplished by any suitable means. FIG. 2 is a flow chart of an exemplary embodiment in which compound heterozygous mutations (CHMs) are identified in a population. An assignment of an ancestral superclass of a subject based on the identified variants may be established at step 12 by any suitable means. A first homoeologous estimate of the subject within the ancestral superclass may be generated at step 13 by any suitable means. A second homoeologous estimate of the subject independent of the ancestral superclass of the subject may be generated at step 14 by any suitable means. At step 15, the subjects may be clustered into a first first-degree family network based on one or more of the second homoeologous estimates by any suitable means. A third homoeologous estimate of the subject within the first first-degree family network may be generated at step 16 by any suitable means. The first homoeologous estimate and the third homoeologous estimate may be combined at step 17 by any suitable means to obtain a combined homoeologous estimate. At step 18, a second first-degree family network may be constructed based on the combined homoeologous estimate by any suitable means. Non-limiting examples of such tools include PLINK (Chang et al. (2015); Gigascience 4, 7) and those disclosed in the present Examples. In some exemplary embodiments, homoeologous estimation involves the calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs.
[0169] Variants in the samples may be phased according to population allele frequency in step 19 by any suitable means, including but not limited to EAGLE (Loh et al. (2016), Nat Genet 48, 1443-1448).
[0170] A pair of phased variants may be classified as a potential CHM in step 20 based on the presence of the two variants in the same subject and gene, which was confirmed by testing all possible combinations of heterozygous pLoF and / or deleterious missense variants within the gene in the same person.
[0171] Potential CHMs may be phased as cis or trans in step 21, and potential CHMs phased as trans may be classified as CHMs. Potential CHMs may be phased by any suitable means. In a non-limiting example, a combination of population allele frequency-based phasing and ancestry / relationship-based phasing using EAGLE is used to determine whether the phase of a potential CHM is cis or trans (this exemplary process is also shown in FIG. 9).
[0172] In some examples, the method further comprises filtering the identified variants before an ancestral superclass assignment for the subject is established, and in some exemplary embodiments, the method further comprises filtering the identified variants before a second homoeologous estimate for the subject is generated. Variants may be filtered by any suitable means. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7) and those disclosed in the present Examples.
[0173] In some exemplary embodiments, filtering variants comprises removing variants with an alternative allele frequency of greater than about 10% across samples from multiple human subjects, or a p-value of greater than about 10. -6 The methods include removing variants that are above 100% of all variants in the sequence, that violate Hardy-Weinberg Equilibrium (HWE), or that are missing calls in more than about 5% of samples from multiple human subjects, or combinations thereof. Variants that do not meet certain selection criteria can be removed, and the remaining variants are selected as potential compound heterozygous mutations, and the potential compound heterozygous mutations are phased by any suitable means, including those described in this example. These exemplary embodiments are also shown in FIG. 6.
[0174] In some exemplary embodiments, the method includes removing low-quality samples after the identified variants are filtered. The low-quality samples may be removed by any suitable means. Non-limiting examples of such means include those disclosed in Dewey et al. (2016), Science 354, aaf6814-1 to aaf6814-10, which are generally known, and those disclosed in the present examples. In some exemplary embodiments, parameters are adjusted such that samples with D-stat greater than 0.12 or with 20x read coverage less than 75%, or both, are the low-quality samples to be removed.
[0175] Integrating the first homoeologous estimate with the third homoeologous estimate may include supplementing the pairwise homoeologous estimate specific to the third homoeologous estimate to enhance the first homoeologous estimate, which may be facilitated by, for example, but not limited to, PLINK (Chang et al. (2015); Gigascience 4, 7) and the means disclosed in this embodiment.
[0176] In some exemplary embodiments, the method further comprises filtering the variants after they have been phased according to population allele frequency, which in some exemplary embodiments may comprise dividing the DNA sequence sample of the human subject into genomic segments of approximately equal size and with substantial segment overlap and breakpoints in intergenic regions. Phasing variants according to population allele frequency may be facilitated by any suitable means, including but not limited to EAGLE (Loh et al. (2016), Nat Genet 48, 1443-1448). Filtering the phased variants according to population allele frequency may include removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of about 5 or less, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of about 5 or less, or an AB of about 10% or less, or a combination thereof. Filtering the phased variants according to population allele frequency may include removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of about 7 or less, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of about 10 or less, or an AB of about 20% or less, or a combination thereof. Variants that do not meet certain selection criteria may be removed, and the remaining variants are selected as potential compound heterozygous variants, which are phased by any suitable means, including those described in this example.These exemplary embodiments are also illustrated in FIG.
[0177] Potential CHMs can be phased based on trio data, or parent-child data, or full sibling data, or distant relative data, or combinations thereof, or phased based on minor allele count (MAC), or phased based on population allele frequency, or combinations thereof. Phasing can be facilitated by any suitable method commonly used in the art. In a non-limiting example, a combination of population allele frequency-based phasing and ancestry / relationship-based phasing using EAGLE is used to phase potential CHMs. This exemplary process is also shown in Figure 9.
[0178] In some exemplary embodiments, the method further includes scoring the CHMs according to functional effect priority and selecting the CHMs with the highest functional effect priority score per gene per sample, thereby obtaining a set of medically important mutations. These steps may be performed by any suitable means, including but not limited to SIFT (Loh et al. (2016); Nat Genet 48, 1443-1448) (damaging), PolyPhen2 HDIV (damaging and possibly damaging), PolyPhen2 HVAR (damaging and possibly damaging), LRT (deleterious), and MutationTaster (Schwarz et al. (2014); Nat. Methods 11, 361-362) (disease causing automatic and disease causing).
[0179] In some exemplary embodiments, the D-stat of a low quality sample is determined by comparing the distribution of the actual allele balance of the sample to the distribution of the expected allele balance using a KS test.
[0180] The present disclosure also provides a method for identifying de novo mutations (DNMs) in a population. A flow chart showing an example of a method for identifying DNMs is provided in FIG. The method may be applied to any type of DNA sequence sample from any type of human subject obtained by any means. Non-limiting examples of variants include point mutations, insertions, deletions, inversions, duplications, and multimerizations. Non-limiting examples of types of human subjects include human subjects from a single medical network population, a multi-medical network population, a racially, culturally, or socially homogenous or heterogeneous population, a mixed age population or a homogenous population with respect to age, a geographically concentrated population or a geographically dispersed population, or a combination thereof. DNA sequence samples may be obtained in any of a number of ways, including but not limited to those disclosed in Dewey et al. (2016), Science 354, aaf6814-1 to aaf6814-10.
[0181] The DNA sequence sample includes or is an exome sequence. Exome DNA can be isolated by any of the commonly used methods or as described in Dewey et al. (2016), Science 354, aaf6814-1 to aaf6814-10.
[0182] Variants in DNA sequence samples from multiple human subjects may be identified in step 22 by any suitable means. Non-limiting examples of means by which variants may be identified include the following steps:
[0183] Upon completion of sequencing, raw data from each sequencing run is collected in local buffer storage and uploaded to the DNAnexus platform (Reid et al. (2014); BMC Bioinformatics 15, 30) for automated analysis.
[0184] - Sample-level read files were generated with CASAVA (Illumina Inc., San Diego, CA) and aligned to GRCh38 (Li and Durbin (2009); Bioinformatics 25, 1754-176; Li (2013); arXiv q-bio.GN) with BWA-mem.
[0185] - The resulting BAM files are processed using GATK (McKenna et al. (2010); Genome Res. 20, 1297-1303) and Picard to sort, mark duplicates, and perform local alignment corrections of reads around putative indels.
[0186] - Sequenced variants are annotated with snpEFF (Cingolani et al. (2012); Fly (Austin) 6, 80-92) using Ensembl85 gene definitions to reveal functional effects on transcripts and genes.
[0187] It is understood that the present disclosure is not limited to any of the steps described above, and obtaining sequence variants may be accomplished by any suitable means. 3 is a flow chart of an exemplary embodiment in which de novo mutations (DNMs) in a population are identified. An assignment of an ancestral superclass of a subject based on the identified variants may be established at step 23 by any suitable means. A first homoeologousness estimate of the subject within the ancestral superclass may be generated at step 24 by any suitable means. A second homoeologousness estimate of the subject independent of the ancestral superclass of the subject may be generated at step 25 by any suitable means. At step 26, the subjects may be clustered into a first first-degree family network by any suitable means based on one or more of the second homoeologousness estimates. A third homoeologousness estimate of the subject within the first first-degree family network may be generated at step 27 by any suitable means. The first homoeologousness estimate and the third homoeologousness estimate may be combined at step 28 by any suitable means to obtain a combined homoeologousness estimate. Non-limiting examples of such tools include PLINK (Chang et al. (2015); Gigascience 4, 7), and those disclosed herein. Homologousness estimation can include calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs.
[0188] Further, based on the integrated homoeologous estimation, a nuclear family may be constructed in step 29. Variants in the nuclear family may be identified in step 30. In step 31, the variants of the samples from the parents and children of the constructed nuclear family trio may be assigned genotype likelihood scores, and the probability that the variant is a de novo mutation may be calculated, and independently, variants of the child sample that are not present in either parent sample of the trio may simply be identified, and the probability that the variant is a de novo mutation may be calculated, and then both sets of putative de novo mutations may be combined, thereby forming a data set of putative de novo mutations. Non-limiting examples of means for carrying out the above steps include those disclosed in the present embodiment.
[0189] In some exemplary embodiments, the method further comprises filtering the identified variants before an ancestral superclass assignment for the subject is established, and in some exemplary embodiments, the method further comprises filtering the identified variants before a second homoeologous estimate for the subject is generated. Variants may be filtered by any suitable means. Non-limiting examples of such means include PLINK (Chang et al. (2015); Gigascience 4, 7) and those disclosed in the present Examples.
[0190] Filtering variants includes removing variants with alternative allele frequencies greater than about 10% across samples from multiple human subjects, or a p-value of about 10 -6 The selection criteria may include removing variants that are above 100% of all variants in the sequence, that violate Hardy-Weinberg equilibrium (HWE), or that are missing calls in more than about 5% of samples from multiple human subjects, or a combination thereof. Variants that do not meet certain selection criteria may be removed, and the remaining variants are selected as potential compound heterozygous mutations, and the potential compound heterozygous mutations are phased by any suitable means, including those described in this example.
[0191] In some exemplary embodiments, the method includes removing low-quality samples after the identified variants are filtered. The low-quality samples may be removed by any suitable means. Non-limiting examples of such means include those disclosed in Dewey et al. (2016), Science 354, aaf6814-1 to aaf6814-10, which are generally known and therefore not further detailed herein, and those disclosed in the present examples.
[0192] In some exemplary embodiments, parameters are adjusted such that samples with D-stats greater than 0.12 or less than 75% read coverage at 20x, or both, are low-quality samples that are removed. In some exemplary embodiments, the D-stats of low-quality samples are determined by comparing the distribution of the actual allele balance of the sample to the distribution of the expected allele balance using a KS test.
[0193] Integrating the first homoeologous estimate with the third homoeologous estimate may include supplementing the pairwise homoeologous estimate specific to the third homoeologous estimate to augment the first homoeologous estimate, which may be facilitated by, for example, but not limited to, PLINK (Changetal. (2015); Gigascience 4, 7) and the means disclosed in this embodiment.
[0194] Filtering the variants may include removing variants with a genotype quality (GQ) annotation of less than about 35 in the child sample, or variants with an alternative allele count (AC) of 10 or more across samples, or variants with a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or variants with an allele balance (AB) of more than about 2% in either parent sample, or variants with an allele balance (AB) of less than about 15% in the child sample, or variants with an AB of more than about 90% in the child sample, or variants with alternative allele homozygosity in either parent sample, or combinations thereof. A putative de novo mutation can be classified as a medium-confidence de novo mutation if the allele balance (AB) of the putative de novo mutation is greater than or equal to about 15% in the child sample and less than or equal to about 2% in each parent sample, if the mapping quality (MQ) of the putative de novo mutation is not less than about 40, if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, if the minor allele count (MAC) of the putative de novo mutation is less than about 20 across all samples, if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs. Classifying a medium-confidence de novo mutation as a high-confidence de novo mutation may be performed if the medium-confidence de novo mutation has a genotype quality (GQ) annotation of about 90 or more in the parent samples, and if the medium-confidence de novo mutation has a read depth (DP) of about 10 or more in each parent sample, and if the medium-confidence de novo mutation has an alternative DP of about 7 or more in the child sample, and if the medium-confidence de novo mutation has a SNP with a QD of greater than about 3, and if the medium-confidence de novo mutation has an INDEL with a QD of greater than about 5. These steps may be performed by any suitable means.
[0195] These exemplary embodiments include those described in this example and are also shown in FIG. The term "D-stat" as used herein refers to a QC metric that can be generated and used to identify low quality samples. The poor quality of a sample may be due to contamination, which may cause problems in downstream analysis. The D-stat of a sample may be calculated, for example, by comparing the distribution of the actual allele balance of the sample to a reference allele balance distribution (e.g., expected allele balance distribution). This reference distribution may be calculated, for example, from multiple samples without any evidence of contamination that are captured and sequenced using the same platform used to query the sample to be analyzed. The QC metric value of the D-stat as used herein is equivalent to the D-statistic generated from the KS (Kolmogorov-Smirnov) test prior to the calculation of the p-value. The D-stat is unitless. The D-statistic value of the KS test can be a value between 0 and 1, with 1 meaning the maximum difference between the cumulative distribution of the reference distribution and the sample distribution. In some exemplary embodiments, low quality samples are identified by comparing the distribution of the actual allele balance of the sample to the expected / reference distribution of the allele balance calculated according to the KS test. In some exemplary embodiments, samples determined to have a particular D-stat value are considered to be low quality samples and are excluded from further analysis. In some exemplary embodiments, samples that are considered to be low quality and should be excluded have a D-stat value greater than 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, or 0.12. In preferred embodiments, samples that are considered to be low quality and should be excluded have a D-stat value greater than about 0.12. In even more preferred embodiments, samples that are considered to be low quality and should be excluded have a D-stat value greater than 0.12.
[0196] Any of the methods described or illustrated may be implemented as a non-transitory computer-implemented method and / or system, and any suitable computer system known to those skilled in the art may be used for this purpose.
[0197] 10 illustrates various aspects of an exemplary environment 201 in which the method and system may operate. The method may be used in a variety of networks and systems that use both digital and analog equipment. Functional descriptions are provided herein, indicating that each function may be performed by software, hardware, or a combination of software and hardware.
[0198] The environment 201 may include a local data / processing center 210. The local data / processing center 210 may include one or more networks, such as a local area network, to facilitate communication between one or more computing devices. The one or more computing devices may be used to store, process, analyze, output, and / or visualize biological data. The environment 201 may optionally include a medical data provider 220. The medical data provider 220 may include one or more sources of biological data. For example, the medical data provider 220 may include one or more medical systems from which medical information of one or more patients is available. The medical information may include, for example, medical history, medical professional opinions and annotations, laboratory reports, diagnoses, physician orders, prescriptions, vital signs, fluid balance, respiratory function, blood parameters, electrocardiograms, x-rays, CT scans, MRI data, laboratory test results, diagnoses, prognoses, evaluations, admission and discharge records, and patient registration information. The medical data provider 220 may include one or more networks, such as a local area network, to facilitate communication between one or more computing devices. The medical information may be stored, processed, analyzed, output, and / or visualized using one or more computing devices. The medical data provider 220 may de-identify the medical information and provide the de-identified medical information to the local data / processing center 210. The de-identified medical information may include a unique identifier for each patient to distinguish one patient's medical information from another patient's medical information while keeping the medical information de-identified. The de-identified medical information ensures that a patient's personal information is not associated with the patient's specific medical information. The local data / processing center 210 may analyze the de-identified medical information to assign one or more phenotypes to each patient (e.g., by assigning an International Classification of Diseases "ICD" code and / or a Physician Procedural Terminology "CPT" code).
[0199] The environment 201 may include an NGS sequencing facility 230. The NGS sequencing facility 230 may include one or more sequencers (e.g., Illumina HiSeq 2500, Pacific Biosciences PacBio RS II, etc.). The one or more sequencers may be configured for exome sequencing, whole exome sequencing, RNA-seq, whole genome sequencing, targeted sequencing, etc. In an exemplary embodiment, the medical data provider 220 may provide a biological sample from a patient associated with de-identified medical information. A unique identifier may be used to maintain an association between the biological sample and the de-identified medical information corresponding to the biological sample. The NGS sequencing facility 230 may sequence the exome of each patient based on the biological sample. The NGS sequencing facility 230 may include a biobank (e.g., manufactured by Liconic Instruments) to store the biological sample prior to sequencing. Biological samples may be contained in tubes (each tube assigned to a patient), and each tube may include a barcode (or other identifier) that may be scanned to automatically record the sample in the local data / processing center 210. The NGS sequencing facility 230 may include one or more robots for use in one or more phases of sequencing to ensure homogenous data and efficient non-disruptive operation. In this manner, the NGS sequencing facility 230 may sequence tens of thousands of exomes per year. In one embodiment, the NGS sequencing facility 230 has a functional capacity to sequence at least 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10,000, 11,000, or 12,000 whole exomes per month.
[0200] Biological data (e.g., raw sequencing data) generated by the NGS sequencing facility 230 may be transferred to the local data / processing center 210, which may then transfer the biological data to the remote data / processing center 240. The remote data / processing center 240 may comprise a cloud-based data storage and processing center including one or more computing devices. The local data / processing center 210 and the NGS sequencing facility 230 may communicate data with the remote data / processing center 240 via one or more high-capacity fiber lines, although other data communication systems (e.g., the Internet) are also contemplated. In an exemplary embodiment, the remote data / processing center 240 may include a third-party system, such as Amazon Web Services (DNAnexus). The remote data / processing center 240 may facilitate automation of analysis steps and allow data to be shared with one or more collaborators 250 in a secure manner. Upon receiving the biological data from the local data / processing center 210, the remote data / processing center 240 may use bioinformatics tools to execute an automated series of pipeline steps for primary and secondary data analysis, resulting in the generation of annotated variant files for each sample. The results of such data analysis (e.g., genotypes) may be sent back to the local data / processing center 210 and may be incorporated, for example, into a Laboratory Information Management System (LIMS) configured to maintain the status of each biological sample.
[0201] The local data / processing center 210 may then utilize the biological data (e.g., genotypes) obtained by the NGS sequencing facility 230 and the remote data / processing center 240 in combination with the de-identified medical information (including the identified phenotypes) to identify associations between genotypes and phenotypes. For example, the local data / processing center 210 may apply a phenotype-first approach, in which phenotypes that may have treatable potential in a particular disease area, e.g., extreme blood lipids for cardiovascular disease, are defined. Another example is the study of obese patients to identify individuals who may be protected from a range of typical comorbidities. Another approach starts from genotypes and a hypothesis that, for example, gene X is involved in the development of disease Y or protection from disease Y.
[0202] In an exemplary embodiment, one or more collaborators 250 may access some or all of the biological data and / or de-identified medical information via a network, such as the Internet 260.
[0203] 11 , one or more of local data / processing center 210 and / or remote data / processing center 240 may comprise one or more computing devices including one or more of genetic data component 300, phenotypic data component 310, genetic variant-phenotype association data component 320, and / or data analysis component 330. Genetic data component 300, phenotypic data component 310, and / or genetic variant-phenotype association data component 320 may be configured for one or more of sequence data quality assessment, read alignment to a reference genome, variant identification, variant annotation, phenotype identification, variant-phenotype association identification, data visualization, combinations thereof, and the like.
[0204] In an exemplary aspect, one or more of the components may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the method and system may take the form of a computer program product on a computer-readable storage medium having computer-readable program instructions (e.g., non-transitory computer software) embodied in the storage medium. More specifically, the method and system may take the form of computer software executed on the web. Any suitable computer-readable storage medium may be utilized, including hard disks, CD-ROMs, optical storage devices, or magnetic storage devices.
[0205] In an exemplary embodiment, the genetic data component 300 may be configured to functionally annotate one or more genetic variants. The genetic data component 300 may also be configured to store, analyze, receive, etc., one or more genetic variants. The one or more genetic variants may be annotated based on sequence data (e.g., raw sequence data) obtained from one or more patients (subjects). For example, the one or more genetic variants may be annotated from each of at least 100,000, 200,000, 300,000, 400,000, or 500,000 subjects. The functional annotation of the one or more genetic variants results in the generation of genetic variant data. By way of example, the genetic variant data may include one or more variant call format (VCF) files. A VCF file is a text file format for configuring calls for SNPs, indels, and / or structural variations. Variants are evaluated for functional impact on transcripts / genes, and potential loss-of-function (pLoF) candidates are identified. Variants are annotated with snpEff using Ensembl75 gene definitions, and then functional annotation is further processed for each variant (and gene).
[0206] Sequential number and / or letter labeling of method steps provided herein is not meant to limit the method or any example embodiments thereof to the particular order presented.
[0207] Various publications are cited throughout this specification, including patents, patent applications, published patent applications, accession numbers, technical papers, journal articles, etc. Each of these cited documents is incorporated by reference in its entirety into this document for all purposes.
[0208] The present disclosure will be more fully understood by reference to the following examples, which are provided to more fully illustrate the present disclosure and are intended to be illustrative and should not be construed as limiting the scope of the disclosure. EXAMPLES
[0209] Example 1.1 Estimation of associations and relatedness content in a cohort of 61K human exomes We analyzed a cohort of 61K human exomes. The cohort originated from a study by the Regeneron Genetics Center (RGC) and Geisinger Health System (GHS) that began in 2014 (Dewey et al. (2016), Science 354, aaf6814-aaf6814). The DiscovEHR study densely sampled patients from a single health care system serving a low migration population. This 61K human exome cohort is referred to herein as the DiscovEHR dataset. A significant amount of family structure was identified within the DiscovEHR dataset, and simulations disclosed herein estimated that 70%-80% of individuals in the dataset would have first- or second-degree relatives when the study confirmed its 250K target.
[0210] Idiosyncratic by descent (IBD) estimation was used to identify different types of family relationships in the dataset, and PRIMUS (Staples et al. (2014), Am. J. Hum. Genet. 95, 553-564) was used to classify pairwise relationships into different family classes and to reconstruct pedigrees (as further described in Example 8). Due to limitations in accurately estimating IBD rates in distant relatives from whole exome sequencing (WES) data, only estimated first-, second-, and high-confidence third-degree relationships were included among DiscovEHR dataset samples.
[0211] In total, we identified 20 identical twins, 8,802 parent-child relationships, 6,122 full-sib relationships, and approximately 20,000 second-degree relationships in the dataset (Figure 12A). Because the IBD sharing distributions of second- and third-degree relationships overlap with each other, in this study we set a hard cutoff between the two expected average values. Third-degree relationships (marked with an asterisk in Figure 12A) are difficult to estimate accurately due to technical limitations of the exome data, as well as the extensive and overlapping variations around the expected mean IBD ratios of more distant relationship classes (e.g., fourth- and fifth-degree). We then treated individuals as nodes and relationships as edges to generate an undirected graph. Using only first-degree relationships, we identified 7,684 connected components, referred to as first-degree family networks. Figure 12B shows the distribution of the size of the first-degree family networks, which range from 2 to 25 sequenced individuals. Similarly, 7,136 second-degree family networks were found, with the largest containing 7,123 individuals (~12% of the entire dataset, Figure 12D). In Figure 12D, first-degree family networks within a second-degree family network are represented as red boxes of a size proportional to the number of individuals in the network (including the descent of the first-degree family network shown in Figure 12C). Single individuals are represented as black nodes connected by second-degree ties drawn as blue edges.
[0212] We were also able to identify approximately 4,500 third-degree relationships within the second-degree family network. By relaxing the minimum IBD cutoff for IBD estimation within the represented ancestry groups, we identified significantly more than 50K third-degree relationships in the DiscovEHR dataset. Although individuals of European ancestry only constitute 96.5% of the DiscovEHR dataset (see Table 1a below), the vast majority (>99%) of pairwise relationships found in the dataset involve individuals of European ancestry (see Table 1b below). In any case, many relationships exist between people of European ancestry, between people of non-European ancestry, and between individuals of various ancestries. For example, we found a trio in the DiscovEHR dataset that included a European father and an East Asian mother whose child was assigned unknown ancestry because the child did not closely match the reference population.
[0213] [Table 1-1]
[0214] [Table 1-2]
[0215] [Table 1-3]
[0216] The recruitment rate of relatives far exceeded the experimental sample ascertainment rate (Figure 13A) and throughout the simulations (Figure 14A), the recruitment rate of relatives far exceeded the sample ascertainment rate. The recruitment of additional pairwise relationships resulted in more individuals being involved in these relationships. Currently, 50.4% of 61K individuals have at least one first- or second-degree relative in the DiscovEHR dataset (Figure 13B).
[0217] Example 1.2 Estimation of associations and relatedness content in a cohort of 92K human exomes We analyzed a larger clinical cohort of 92,455 human exomes that originated from an ongoing study by the Regeneron Genetics Center (RGC) and Geisinger Health System (GHS) that began in 2014 (Staples et al. (2018), Am. J. Hum. Genet. 102(5):874-889). Similarly, this expanded DiscoverEHR cohort represents a dense sample of participants from a single health care system serving a largely rural population with low migration rates in central Pennsylvania.
[0218] The set containing the first 61K samples prepared and sequenced (Example 1.1) was designated the "VCRome set". The remaining set of 31K samples was prepared by the same process, except that instead of capture with NimbleGen probes, a slightly modified version of IDT's xGen probes was used, with supplementary probes used to capture regions of the genome that are covered by NimbleGen's VCRome capture reagents but are poorly covered by the standard xGen probes. The captured fragments were bound to streptavidin-conjugated beads, and non-specific DNA fragments were removed by a series of stringent washes according to the manufacturer's (IDT) recommended protocol. This second set of samples was designated the "xGen set". Variant calls were generated with GATK, which was used to locally correct alignment of the aligned and duplicate-marked reads of each sample around the putative indels. GATK's HaplotypeCaller was used to process the INDEL alignment corrected duplicate marked reads to identify all exon positions where the sample differed from the genomic reference in a genomic variant calling format (gVCF). Genotyping was accomplished using GATK's GenotypeGVCF on each sample and a training set of 50 randomly selected samples that output single sample variant calling format (VCF) files that identify both single nucleotide variants (SNVs) and indels compared to the reference. The single sample VCF files were then used to create pseudo samples that contained all variant sites from both sets of single sample VCF files. Additionally, a separate pVCF file for the VCRome set was created by jointly calling the 200 single sample gVCF files with the pseudo sample to force a call or no call for each sample at all variant sites across the two capture sets. All 200 sample pVCF files were combined to create a VCRome pVCF file, and this process was repeated to create an xGen pVCF file.The VCRome and xGen pVCF files were combined to create a combined pVCF file. Variants were annotated by aligning sequence reads to GRCh38 and using Ensembl85 gene definitions. Gene definitions were restricted to 54,214 transcripts corresponding to 19,467 protein-coding genes with annotated starts and stops. After sample QC process, 92,455 exomes were left for analysis.
[0219] From the extended DiscovEHR dataset of 92,455 individuals, we identified 43 identical twins, 16,476 parent-child relationships, 10,479 full sib relationships, and 39,000 second-degree relationships (Figure 16, Panel A). To generate an undirected graph, we treated individuals as nodes and relationships as edges. Using only first-degree relationships, we identified 12,594 connected components, called first-degree family networks. Figure 16, Panel B shows the distribution of the size of the first-degree family networks, which range from 2 to 25 sequenced individuals. Similarly, we identified 10,173 second-degree family networks, with the largest containing 19,968 individuals (22% of the entire dataset, Figure 16, Panel C). We also identified approximately 5,300 third-degree relationships within the second-degree family networks. By using a low IBD cutoff (PIJ_HAT>0.09875) for IBD estimation within ancestry groups without considering second-degree family networks, we identified over 100,000 third-degree relationships within the extended DiscovEHR cohort. Given that 95.9% of the extended DiscovEHR individuals were of European ancestry (Table 2a), it is not surprising that the majority of pairwise relationships found (98.6%) were between two individuals of European ancestry (Table 2b). Nevertheless, we identified many relationships between people of European ancestry, between people with the same non-European ancestry, and between individuals with different ancestries. For example, there were several trios with one European parent and one East Asian parent, as well as children whose ancestry could not be assigned to a superpopulation due to mixed genome types. Importantly, experimentally (Figure 17A) and through simulations (Figure 18A), we found that the rate of kinship accumulation far exceeded the rate of ascertaining samples. This is expected given that the number of possible pairwise relationships in a dataset combinatorially increases as the dataset size increases, and also the likelihood that previously unrelated individuals in the dataset will become involved in newly identified relationships.Currently, 39% of individuals in the expanded DiscovEHR cohort are likely to have at least one first-degree relative in the dataset, and 56% of participants are likely to have one or more first- or second-degree relatives in the dataset (Figure 17, Panel B).
[0220] [Table 2-1]
[0221] [Table 2-2]
[0222] [Table 2-3]
[0223] [Table 2-4]
[0224] Example 2 Simulation and relatedness prediction using SimProgeny In an attempt to model, understand, and predict the growth of relationship networks within the DiscovEHR and extended DiscovEHR datasets, we developed a simulation framework (hereafter "SimProgeny") that is capable of simulating the lineages of millions of people over hundreds of years, distributed across multiple subpopulations. From these simulated populations, SimProgeny can model various sampling approaches and estimate the amount of relatedness a researcher would expect to find for a given population and set of sampling parameters (see Example 17).
[0225] SimProgeny was used to simulate DiscovEHR and Extended DiscovEHR cohorts and the ascertainment of the first 61K and first 92K participants from them, respectively. The simulations show that DiscovEHR and Extended DiscovEHR participants were enriched for close relatives rather than randomly sampled from the cohort. The real data was calculated in periodic "freezes", shown by punctuation marks connected by a faint line, as shown in Figures 14A and 14B. We also took the samples and relationships identified in the 61K freezes and then shuffled the ascertainment order to show that the first half of the 61K DiscovEHR participants were enriched for first-degree relationships compared to the second half. Cohorts of various sizes were simulated using parameters similar to the real cohorts ascertained by DiscovEHR. We then ran random ascertainments from each of these cohorts to see which cohort size best fitted the real data. An important note is that none of these cohort sizes fit the real data, and the random ascertainment approach is a poor fit. Alternative ascertainment approaches that enrich for first-degree relatives may yield better fits compared to random ascertainment. Figure 14A shows that ascertainment of first-degree relative pairs in an effective sampling population of size 270K fits a shuffled version of the real data well, but underestimates the number of relative pairs below the 61K participants ascertained, and significantly overestimates the number of relative pairs above the 61K participants ascertained. Figure 14B shows that the 270K population fits the shuffled real data best for the number of individuals with one or more first-degree relatives, but fits the real data poorly.
[0226] Similar results were observed using the extended DiscovEHR dataset (Figures 17A and 17B). We then shuffled the samples and relationships identified in the 92K freeze, showing that the first half of the 92K extended DiscovEHR participants were enriched for first-degree relationships compared to the second half. Random ascertainments were then performed from each of these populations to see which population size best fit the real data. Figure 17A shows that ascertainment of first-degree relative pairs in an effective sampling population of size 403K fits the shuffled version of the real data well, but underestimates the number of relative pairs below the 92K participants we ascertained, and significantly overestimates the number of relative pairs above the 92K participants we ascertained. Figure 17B shows that the 403K population fits the shuffled real data best for the number of individuals with one or more first-degree relatives, but fits the real data poorly, suggesting that the extended DiscovEHR participants were not ascertained randomly.
[0227] Enrichment of close relatives was modeled by using a clustered ascertainment approach (see Example 17), which produced simulations that better fit real data in DiscovEHR (Figures 15A and 15B) and Extended DiscovEHR (Figures 18A and 18B). In both Figures 15 and 18, real data was calculated with periodic "freezes", shown as punctuation marks connected by light lines. Most simulation parameters were set based on information about the demographics of the real population and the DiscovEHR ascertainment approach. However, two parameters were unknown and were chosen in light of the fit to the real data: 1) the effective population size from which the sample was ascertained, and 2) the increased likelihood that someone would be ascertained given previously ascertained first-degree relatives, called "clustered ascertainment". All panels show the same three simulated population sizes spanning the estimated effective population sizes. Clustered ascertainment was simulated by randomly ascertaining individuals in conjunction with a random number of Poisson-distributed first-degree relatives (the lambda values for the Poisson distribution are shown in the legend). These simulation results suggested that the effective sampling population size was approximately 475K individuals, and that a Poisson distribution with lambda 0.2 best matched enrichment for first-degree relatives. This was consistent with the understanding that the majority of current participants live in a specific geographic area, such as Danville, PA (approximately 500K individuals) in this example, rather than being evenly distributed across the entire GHS jurisdiction (>2.5M individuals).
[0228] After identifying simulation parameters that provided a reasonable fit to the real data, SimProgeny was used to obtain aggregate predictions of the first-degree relationships that should be expected when the DiscovEHR and Extended DiscovEHR studies are expanded to a goal of 250K participants. Results showed that if participant ascertainment continues in the same manner, approximately 150K first-degree relationships should be expected in DiscovEHR (Figure 15C) and Extended DiscovEHR (Figure 18C), involving approximately 60% of DiscovEHR participants (Figure 15D), and involving approximately 60% of Extended DiscovEHR participants (Figure 18D).
[0229] The simulation analysis was then extended to include second-degree relationships, and the simulation results suggested that with 250K participants, one would expect well over 200K combined first- and second-degree relationships involving over 70% of the individuals in DiscovEHR (Figure 19) and Extended DiscovEHR (Figure 20). In this analysis, real data was calculated with periodic "freezes", indicated by punctuation marks connected by faint lines in the figures. Most simulation parameters were set based on information about the demographics of the real population and the DiscovEHR ascertainment approach. All panels show the same three simulated population sizes, spanning the estimated effective population sizes. Clustered ascertainment was simulated by randomly ascertaining individuals in conjunction with a random number of Poisson-distributed first-degree relatives and a separate random number of second-degree relatives (both Poisson distributions have lambdas indicated in the figure legends).
[0230] Simulation results showed a clear enrichment of relatedness in the DiscovEHR HPG study, providing important insight into the vast amount of relatedness that is expected to be seen as we continue to ascertain additional participants.
[0231] Example 3.1 Exploiting relatedness instead of treating it as a random variable in the DicoverEHR dataset Using the family reconstruction tool PRIMUS (Staples et al. (2014), Am. J. Hum. Genet. 95, 553-564.), we reconstructed all 7,684 first-degree family networks in the DiscovEHR dataset and found that 98.9% of these pedigrees were uniquely reconstructed when considering IBD estimates and reported ages. These pedigrees included 1,081 nuclear families (925 trios, 134 quartets, 19 quintets, and 3 sextets). Table 3 below shows the breakdown of trios by ancestry. We divided the 1,081 nuclear families into their individual trio components. For example, a quartet would split into two separate trios with the same parents. As the DiscovEHR cohort was mostly European, the majority of trios included individuals of European ancestry. Individuals with unknown ancestry were generally children of parents with different lineages. For example, all three of the East Asian-European-unknown trios contain a European father and an East Asian mother, resulting in mixed children. There were no closely matched reference populations for these mixed European-East Asian individuals, resulting in unknown ancestry.
[0232] [Table 3]
[0233] Figure 12C shows the largest first-degree pedigree identified in the DiscovEHR dataset, which contains 25 sequenced individuals. These relationships and pedigrees were used in several ways, including the following:
[0234] Compound heterozygous mutations A major goal of human genetics is to gain a deeper understanding of the function of all genes in the human genome. Homozygous loss-of-function mutations (LoF) are a powerful tool to gain insight into gene function by analyzing the phenotypic effects of these "human knockouts" (KOs). In recent large-scale sequencing studies, rare (MAF<1%) homozygous LoFs have been of particular interest and have been important in identifying many gene-phenotype interactions (Lek et al.(2016), Nature Publishing Group 536, 285-291, Dewey et al.(2016), Science 354, aaf6814-aaf6814, Saleheen et al.(2017), Nature Publishing Group 544, 235-239, and Narasimhan et al.(2016), Science 352, 474-477). Rare compound heterozygous mutations (CHMs) of two heterozygous LoFs are functionally equivalent to rare homozygous KOs, but they are rarely investigated in such large-scale sequencing studies (Lek et al. (2016), Nature Publishing Group 536, 285-291; Dewey et al. (2016), Science 354, aaf6814-aaf6814; and Saleheen et al. (2017), Nature Publishing Group 544, 235-239). Accurate identification of rare CHMs in LoF is valuable because (1) rare CHMs substantially increase the number of human gene KOs and improve statistical power, (2) KOs of rare CHMs may contain extremely rare heterozygous mutations that may be absent in homozygous carriers, and (3) rare CHMs provide a more complete set of KOs for "human KO projects" (Saleheen et al. (2017), Nature Publishing Group 544, 235-239; Perdigoto (2017), Nat. Rev. Genet. 18, 328-329).
[0235] A survey of rare CHMs in the DiscovEHR dataset was conducted. First, 39,459 high-quality potential CHMs (pCHMs) were identified, consisting of pairs of rare heterozygous variants that were either putative LoFs (pLoFs, i.e., nonsense, frameshift, or splice site mutations) or missense variants with strong evidence of being deleterious (see Example 10). Second, pCHMs were phased using a combination of allele frequency-based phasing using EAGLE and ancestry-based phasing using reconstructed ancestry and relationship data (Figure 9). EAGLE phased pCHMs with 91.4% accuracy based on trio validation (see Table 4 below). However, because there was extensive ancestry and relationship data in this cohort, nearly one-third of the pCHMs could be phased with approximately 100% accuracy based on these data (see Table 4 below), allowing for a reduction in inaccurate phasing by approximately 31%. The phased pCHMs spanned the full range from singleton to 1% MAF (see Table 5 below).
[0236] [Table 4]
[0237] All pCHMs with MAF<1% and MAC>1 occurring in the children of the trio were phased using the reconstructed trio and considered to be "true". pCHMs in which one or more of the contributing variants were determined to be de novo in the children were excluded. The phased pCHMs of the trio were then used to evaluate other phasing methods. The accuracy of EAGLE was assessed by removing all first-degree relatives of one child from each reconstructed nuclear family and then phasing all variants in the remaining dataset. The phased pCHMs of EAGLE were compared to the phased pCHMs of the trio.
[0238] [Table 5]
[0239] Because the accuracy of pCHM tended to decrease for extremely rare variants, the MAF of the rarer of the two pCHM variants was used to place pCHM into its respective frequency bin. pCHMs with a MAC of 1 were phased using the association data and considered not to be de novo mutations in pCHM carriers. Unknown phases for pCHM were the result of one or both pCHM variants being filtered by EAGLE (MAC=1 or missing >10%) and there being no association data for phasing.
[0240] After processing, we phased 39% of pCHMs in trans, yielding a high-confidence set of 13,335 rare and deleterious CHMs distributed among 11,375 of the 61K individuals (mean = 0.22; max = 6; Figure 21, Panel A). The median genomic distance between pCHM variants in cis (5,308 bps) was slightly less than half the median distance between variants in trans (11,201 bps; Figure 22). Nearly one-third of CHMs contained at least one pLoF, and 9.8% of CHMs were composed of two pLoF variants (see Table 6 below). Over 3,385 of the 19,467 targeted genes contained one or more CHM carriers (see Table 7 below), of which 1,555 (46%) had multiple carriers. Eleven genes with over 85 CHM carriers were predicted to be the most LoF resistant in the genome based on the ExAC pLI scores (Lek et al. (2016), Nature Publishing Group 536, 285-291) (see Table 8 below).
[0241] [Table 6]
[0242] This table provides a breakdown of CHMs that are composed of rare (<1% MAF) pLoFs and missense variants, and also indicates how many of these CHMs comprise indel-indel, indel-SNP, and SNP-SNP pairings.
[0243] [Table 7]
[0244] [Table 8]
[0245] For the 11 genes with the most CHM, we display their pLI scores as reported by ExAC, as well as the percentile of each gene relative to the LoF tolerance, calculated by ranking all genes by their pLI score and dividing by the total number of genes with a reported pLi score.
[0246] * The pLI score for SSPO is not reported by ExAC. To obtain a more robust set of human knockout genes and demonstrate the added value of CHM, we combined CHM with 3,915 homozygous pLoFs found among 61K DiscovEHR participants. pLoF-pLoF CHM increased the number of genes in ≥1 and ≥10 individuals with putative KO by 15% and 54%, respectively (see Table 6 above). The benefit of including CHM in the KO analysis became even more significant when considering missense variants predicted to disrupt protein function: CHM increased carriers by 28% for ≥1 gene and carriers by 246% for ≥10 genes, where both copies of the gene are predicted to be fully knocked out or disrupted.
[0247] The results of the trio validation showed that familial relationship-based phasing was 100% accurate (750 / 750 pCHM) and EAGLE phasing was less accurate at 91.4% (459 / 502 pCHM; see Table 3 above). Visual validation of the Illumina read data was also performed for 190 pCHMs (115 cis and 79 trans; 126 phased with EAGLE and 74 phased with familial / relationship). Visual validation yielded overall accuracies of 95.8% and 89.9% for familial / relationship and EAGLE phasing, respectively (see Table 9 below). The Illumina read-based validation results were consistent with the trio validation results, although it should be noted that the Illumina read-based validation accuracy results were lower than the phasing accuracy determined by phasing with trio. This difference is likely due to enrichment of false-positive pCHMs in small problematic regions of exons that are prone to sequencing and variant calling errors.
[0248] [Table 9]
[0249] We randomly selected 200 pCHMs among 61K DiscovEHR participants where both variants were within 75 base pairs of each other, and then visually verified the phase by looking at the read stacks that spanned the two variants. Ten (5%) could not be confidently phased using the read stacks, either because there were no overlapping reads in both variants, or the reads gave conflicting results (i.e., some reads showed cis and some showed trans).
[0250] De novo mutations De novo mutations (DNMs) are a class of rare variations that are likely to generate extreme phenotypes in humans due to reduced purifying selection. Recent sequencing studies have shown that DNMs are a major contributor to human genetic disease (de Ligt et al. (2012), N. Engl. J. Med. 367, 1921-1929; Deciphering Developmental Disorders Study (2017). Prevalence and architecture of de novo mutations in developmental disorders. Nature Publishing Group 542, 433-438; and Fromer et al. (2014), Publishing Group 506, 179-184), indicating that DNMs are a valuable tool for gaining a deeper understanding of gene function.
[0251] Using nuclear families reconstructed from the DiscovEHR dataset, we confidently called 1,800 medium- and high-confidence exonic DNMs distributed among 887 of 1,262 valid children in trios (see Example 12). The average number of DNMs per individual was 1.43, with a maximum of 49 (Figure 21C). PolyPhen2 predicted 28.2% (N=507) of the DNMs as "probably damaging" and an additional 8.6% (N=154) as "possibly damaging". DNMs were distributed across 1,597 genes (Figure 21D), with only one gene harboring more than 5 DNMs. The most common type of DNM was nonsynonymous SNVs (57.17%), followed by synonymous SNVs (25.56%). Table 10 below provides a full breakdown of DNM types and shows that the proportions of DNMs falling into the various functional classes correspond well with those found in recent studies of DNMs in developmentally disabled children.
[0252] [Table 10]
[0253] * Deciphering Developmental Disorders Study (DDD) (Deciphering Developmental Disorders Study (2017). Prevalence and architecture of de novo mutations in developmental disorders. Nature Publishing Group 542, 433-438). The DDD paper also reports 57 DNMs from other classes that are not included in our analysis nor in this table. Percentages have been adjusted accordingly.
[0254] We attempted to perform visual validation of 23 high-confidence, 30 medium-confidence, and 47 low-confidence DNMs across all functional classes. Eight medium-confidence and two low-confidence variants could not be confidently called as true-positive or false-positive DNMs. The remaining 23 / 23 (100%) high-confidence DNMs, 19 / 22 (86%) medium-confidence DNMs, and 12 / 43 (28%) low-confidence DNMs could be confirmed as true positives. Visual validation also confirmed that the majority (40 / 49) of potential DNMs in individuals with more than 10 DNMs were likely false-positive calls.
[0255] Segregation of variants with phenotype in families We used reconstructed pedigree data from within the DiscovEHR dataset to distinguish between novel / rare population variations and familial variants, and leveraged this pedigree data to identify highly penetrant disease variants segregating in families that are not duly represented in population-wide association studies. While this is not intended to be a survey of all known Mendelian disease-causing variations transmitted through these families, we did identify several illustrative examples, including familial aortic aneurysm (Figure 23, Panel A), long QT syndrome (Figure 23, Panel B), thyroid cancer (Figure 23, Panel C), and familial hypercholesterolemia (FH; Figure 24) (Maxwell, EK, et al. (2017). Profiling copy number variation and disease associations from 50,726 DiscovEHR Study exomes). The FH example was particularly interesting as it reconstructed 27 / 29 carriers of a de novo familial hypercholesterolemia-causing tandem duplication of the LDLR. Five additional carriers (not depicted) were also included in this pedigree. Elevated LDL and total cholesterol, and increased prevalence of coronary artery disease and early-onset ischemic heart disease ('ischemic heart disease age' < 55 years in men and < 65 years in women) segregate with the duplication carriers. Their common ancestry history indicates that they all inherited this duplication event from a common ancestor approximately six generations ago. The remaining two samples were first-degree relatives of each other, but they were not successfully genotyped and could not be connected to the larger pedigree.
[0256] Sequencing studies continue to capture and sequence ever-increasing sizes of human populations, revealing the highly complex and intertwined nature of human relatedness. In the DiscovEHR dataset, we identified approximately 35K first- and second-degree relationships, reconstructed 7,684 pedigrees, and revealed second-degree family networks for over 7,000 participants. Studies of founder populations have long emphasized the complexity of the relationships (Old Order Amish (McKusick, VA, HOSTETLER, JA, and EGELAND, JA (1964). GENETIC STUDIES OF THE AMISH, BACK-LOOK AND POTENTIALITIES. Bull Johns Hopkins Hosp 115, 203-222), Hutterites (Ober et al. (2001). The American Journal of Human Genetics 69, 1068-1079), and Ashkenazi Jews (Gusev et al. (2012). Mol. Biol. Evol. 29, 473-486), and recent studies of non-founder populations have demonstrated extensive levels of relatedness (UK Biobank (Bycroft et al. (2017). Genome-wide genetic data on ~500,000 UK Biobank participants), NHAMES (Malinowski et al. (2015). Frontier Genetics 2016). 6,317) and AncestryDNA (Han et al. (2017), Nat Commun 8,14238.). Whereas once only a small number of individuals were involved in large sequencing cohorts, close ties are likely to involve a large proportion, if not the majority, of individuals in large-scale health population-based genomic (HPG) studies. Herein, we show that it is possible to obtain a large number of close family ties, nuclear families, and informative pedigrees through simulation and real data. This observation was likely to be more pronounced in datasets collected for HPG studies, as families tend to visit the same health care system and have similar genetic and environmental disease risks.It is becoming clear that it is no longer possible to simply exclude pairs of closely related individuals from association studies, knowing that they represent only a small fraction of the entire cohort. Traditional approaches that obtain a maximum size set of unrelated relationships dramatically reduce the size of HPG cohorts, making them unsuitable for many important disease phenotyping analyses performed in these types of cohorts. Instead, new methods are needed to exploit the relatedness information outlined herein.
[0257] In this study, we demonstrate several ways of leveraging relatedness information. First, we improved the accuracy of phasing rare compound heterozygous mutations (CHMs). Although accurate phasing of CHMs was obtained using EAGLE, ancestry- and relationship-based phasing was much more accurate, reducing phasing error for pCHMs by an estimated 31%. Because phasing using pairwise relationships assumes that if two variants appear together in two relatives, they are in cis and have segregated together from a common ancestor, the accuracy of relationship-based phasing for pCHMs could be slightly degraded when variants with MAFs >1% were included. It is very likely that two common variants that segregated independently will appear together in multiple people and, as a result, be erroneously phased as cis by the algorithm. For common variants, phasing using population allele frequencies may be more appropriate than relationship-based phasing.
[0258] Second, pedigree reconstruction of the relationships identified in the HPG study provided informative trios and informative pedigrees that could be used in several ways. The 1,262 reconstructed trios were used to find 1,800 DNMs, allowing for tracing of known disease-causing mutations through the extended pedigrees. The number and size of informative pedigrees will continue to increase as larger portions of the population are sequenced, providing richer pedigree data sets. Pedigrees and relationships are particularly useful for very rare variants, as transmission of the rare variant through pedigrees provides strong evidence that it is real and allows the use of more traditional Mendelian genetic approaches. Pedigrees proved particularly useful when combined with DiscovEHR's ability to recontact patients and recruit additional family members to augment small to medium-sized pedigrees in follow-up studies.
[0259] Rather than viewing relatedness as a vexing problem that needs to be addressed, it should be seen as an opportunity to harness a valuable and untapped source of genetic insight. As the era of genome-based precision medicine begins, innovative methods and tools are needed that can effectively mine family structure and distant relatedness in ever-growing sequencing cohorts.
[0260] Example 3.2 Exploiting relatedness instead of treating it as a random variable in the extended DicoverEHR dataset The pedigree structure of 12,574 first-degree family networks in the extended DiscovEHR dataset was reconstructed using the pedigree reconstruction tool PRIMUS. We found that 98.9% of these pedigrees could be unambiguously reconstructed into a single pedigree structure when considering LBD estimates and reported ages of participants. These pedigrees include 2,192 nuclear families (1,841 trios, 297 quartets, 50 quintets, three sextets, and one septet). Table 11 shows the breakdown of trios by ancestry. Figure 14, panel C shows the largest first-degree pedigree, including 34 sequenced individuals.
[0261] [Table 11]
[0262] Compound heterozygous mutations The 57,355 high-quality pCHMs, consisting of pairs of rare heterozygous mutations, were accepted to be putative LoFs (pLoFs, i.e., nonsense, frameshift, or splice site mutations) or missense variants with strong evidence of being deleterious. We then performed phasing of pCHMs by using a combination of allele frequency-based phasing by EAGLE and ancestry-based phasing with reconstructed pedigree and relationship data (Figure 25). Trio validation showed that EAGLE phased pCHMs with an average accuracy of 89.1% (Table 12 below). However, in this cohort, due to extensive ancestry and relationship data, we phased 25.2% of pCHMs and 33.8% of trans-CHMs with highly accurate trio and relationship phasing data (R 98.0%; Table 12), reducing inaccurate phasing of trans-CHMs by approximately one-third. The phased pCHMs spanned the full range from singleton to 1% MAF (see Table 13 below).
[0263] [Table 12]
[0264] All pCHMs with MAF<1% and MAC>1 occurring in the offspring of a trio were phased using the reconstructed trio and considered to be "true". pCHMs in which one or more of the contributing variants were determined to be de novo in the offspring were excluded. The phased pCHMs of the trio were then used to evaluate the other phasing methods. The accuracy of EAGLE was assessed by removing all first-degree relatives of one child from each reconstructed nuclear family and then phasing all variants in the remaining dataset. The phased pCHMs of EAGLE were compared to the phased pCHMs of the trio.
[0265] [Table 13]
[0266] After processing, we phased 40.3% of pCHMs in trans, yielding a high-confidence set of 20,947 rare and deleterious CHMs distributed among 17,533 of the 92K individuals (mean ~0.23 per individual; maximum ~10 per individual; Figure 26, Panel A). The median genomic distance between pCHM variants in cis (5,955 bp) was slightly more than half the median distance between pCHMs in trans (11,600 bp; Figure 27). Nearly one-third of CHMs contained at least one pLoF, and 8.9% of CHMs were composed of two pLoF variants (see Table 14 below). Over 4,216 of the 19,467 targeted genes contained one or more CHM carriers (see Table 15 below), of which 2,468 had multiple carriers (Figure 26, Panel B). The ExAC pLI scores indicate that the 10 genes with a value of >125 in CHM carriers are most likely to be LoF resistant in the genome (see Table 16 below).
[0267] [Table 14]
[0268] This table provides a breakdown of CHMs that are composed of rare (<1% MAF) pLoFs and missense variants, and also indicates how many of these CHMs comprise indel-indel, indel-SNP, and SNP-SNP pairings.
[0269] [Table 15]
[0270] [Table 16]
[0271] For the 10 genes with the most CHM, we display their pLI scores as reported by ExAC3, as well as the percentile of each gene relative to the LoF tolerance, calculated by ranking all genes by their pLI score and dividing by the total number of genes with a reported pLi score.
[0272] * The pLI score for SSPO is not reported by ExAC. To obtain a more robust set of genes where both copies of a gene are knocked out or disrupted in the same individual and to demonstrate the added value of CHM, CHM was combined with 6,560 rare (MAF<1%) homozygous pLoFs found among 92K DiscovEHR participants. pLoF-pLoF CHM increased the number of genes knocked out in R1 and R20 individuals by 15% and 61%, respectively (see Table 16 below). The benefit of including CHMs in the KO analysis becomes even more significant when considering missense variants predicted to disrupt protein function. Combined, 20,364 rare homozygous pLOFs and deleterious missense variants were found in 92K participants. Homozygous pLoFs or carriers of predicted deleterious missense variants provided a large number of genes predicted to be fully knocked out or disrupted. However, when carriers of CHM were included, 26% more genes were knocked out or disrupted in R1 individuals and 397% more genes were knocked out or disrupted in R20 individuals (Table 15).
[0273] De novo mutations Nuclear families reconstructed from 92K extended DiscovEHR participants were able to covertly call 3,415 medium- and high-confidence exonic DNMs distributed across 1,783 of 2,602 valid children in trios (mean ~1.31; max ~48; Figure 26, Panel C). PolyPhen2 predicts 29.1% (n~995) of DNMs as "probably damaging" and an additional 9.2% (n~316) as possibly damaging. DNMs are distributed across 2,802 genes (Figure 26, Panel D), accommodating the most (9) TTNs. The most common type of DNMs is nonsynonymous SNVs (58.5%), followed by synonymous SNVs (24.3%). Table 17 provides a full breakdown of DNM types and shows that the proportions of DNMs falling into various functional classes are roughly consistent with those found in a recent study of DNMs in developmentally disabled children. The workflow for calling, filtering, and confidence ranking DNMs was followed as described in Figure 7. From a cohort of 92,455 sequenced GHS exomes, 2,602 trios were identified (Figure 28). From trios sorted based on low-, medium-, and high-confidence DNMs, 6,645 exonic DNMs were identified. Using families reconstructed from the extended DiscovEHR dataset, 3,409 medium- and high-confidence exonic DNMs, as well as 3,045 single-nucleotide DNMs, were confidently called from 2,602 family trios (Figure 29, panels A and B). Most individuals in the cohort had less than 5 DNMs. Furthermore, from a cohort of 92,455 sequenced GHS exomes, we identified 2,602 trios and sorted these based on low-, medium-, and high-confidence variants to generate 10,000 random variants, providing 73,192 medium- / high-confidence variants.
[0274] The most common type of DNM was nonsynonymous SNV, followed by synonymous SNV. Stop-loss SNVs were the least common DNM. This result was similar to that obtained in the DiscovEHR cohort containing 61K exome sequencing data (see Table 17 below). Figure 30 provides a full breakdown of the types of medium and high confidence exonic DNMs (n=3409) found in the extended DiscovEHR cohort and shows the proportion of DNMs classified into different functional effect classes.
[0275] [Table 17]
[0276] * Deciphering Developmental Disorders Study (DDD) (Deciphering Developmental Disorders Study (2017). Prevalence and architecture of de novo mutations in developmental disorders. Nature Publishing Group 542, 433-438). The DDD paper also reports 57 DNMs from other classes that are not included in our analysis nor in this table. Percentages have been adjusted accordingly.
[0277] Figure 31 provides a breakdown of the types of medium- and high-confidence exonic DNMs (n=3409) found in the extended DiscovEHR cohort, showing the proportion of DNMs caused by transitions, transversions, and indels. Among the medium- and high-confidence exonic DNMs (n=3409) found in the extended DiscovEHR, the number of mutations due to transitions was 2038, the number of mutations due to transversions was 1007, and the number of mutations due to indels was 364. Thus, the ratio of transitions to transversions (Ti:Tv) was 2:1, which was similar to the ratio of transitions to transversions obtained from other studies. Among the single-nucleotide DNMs (n=3045), cysteine to thymine and guanine to adenine were the most common mutations (Figure 32).
[0278] Medium- and high-confidence DNMs were evenly distributed across the autosomes. One-way chi-square test (χ 2 Test) showed that the DNMs per 10M exonic base pairs did not deviate significantly from a random distribution (p=0.045) (Figure 33).
[0279] Mutations in CG dinucleotides (conventionally written CpG, where "p" stands for the phosphate between the two bases) account for one-third of disease-causing germline mutations in humans (Cooper and Krawczak (1990); Hum. Genet. 85:55-74). Of the medium / high confidence DNMs (n=3,409), approximately 13% of DNMs were explained due to DNMs in CpG islands. Of the random variants (n=10,000), approximately 10% of DNMs were explained due to DNMs in CpG islands. DNMs were more likely to occur in CpG islands than random variants (χ 2 = 32.3661, df value = 1, p value = 1.28E-08) (Figure 34). This is more expected due to the high variability of these sites.
[0280] An attempt was made to perform visual validation of 23 high-confidence, 30 medium-confidence, and 47 low-confidence DNMs across all functional classes. Eight medium-confidence and two low-confidence variants could not be confidently called as true-positive or false-positive DNMs. The remaining 23 / 23 (100%) high-confidence DNMs, 19 / 22 (86%) medium-confidence DNMs, and 12 / 43 (28%) low-confidence DNMs were confirmed as true positives. Visual validation also confirmed that the majority (40 / 49) of potential DNMs in individuals with more than 10 DNMs were most likely false-positive calls.
[0281] Segregation of variants with phenotype in families Using reconstructed pedigree data from among 92K extended DiscovEHR participants, we distinguished between rare population variants and familial variants and leveraged this to identify highly penetrant disease variants segregating in families. While this is not intended to investigate all known Mendelian disease-causing variants transmitted through these pedigrees, as with the DiscoverEHR dataset, it has identified familial aortic aneurysm, long QT syndrome, thyroid cancer, and familial hypercholesterolemia (FH[MLM:143890]; Figure 35). During the CNV call update, we found 37 carriers of the FH-causing tandem duplication among 92K exomes. Based on this, we reconstructed 30 of the 37 carriers into a single extended pedigree. The carriers' common ancestral history showed that they all inherited this duplication event from a common ancestor approximately six generations ago. Two of the remaining seven carriers are second-degree relatives to each other, but genotyping sequence data were not available to confirm that the remaining seven carriers are also distantly related to the other carriers in Figure 36. In the pedigree depicted in Figure 36, carrier and non-carrier status was determined from exome data from each individual, and elevated maximum LDL levels (values below the symbols) were found, as well as an increased prevalence of coronary artery disease (CAD, filled red) and pure hypercholesterolemia (ICD 272.0; blue) isolated in the overlapping carriers. Five additional carriers (not depicted) were also found to be distant relatives (seventh to ninth degree relatives) of individuals in this pedigree (Figure 36).
[0282] Example 4 Patients and samples We collected two sets of data by applying the predictive models to the cohorts: (A) the DicovEHR cohort with 61,720 de-identified patient exomes, and (B) the extended DicovEHR cohort with 92,455 de-identified patient exomes.
[0283] We sequenced all de-identified patient participants from both cohorts obtained from the Geisinger Health System (GHS). All participants consented to participate in the MyCode® Community Health Initiative (Carey et al. (2016), Genet. Med. 18, 906-913) and provided DNA samples for genomic analysis in the Regeneron-GHS DiscovEHR study (Dewey et al. (2016), Science 354, aaf6814-aaf6814). All patient exomes were linked to the corresponding de-identified electronic health records (EHR). Details of the first 50,726 sequenced individuals have been published previously (Dewey et al. (2016), Science 354, aaf6814-aaf6814; Abul-Husn et al. (2016), Science 354, aaf7000-aaf7000).
[0284] The study did not specifically target families participating in the study, but rather targeted adults with chronic health problems who have frequent interactions with the health care system, as well as participants from coronary catheterization laboratories and bariatric care units.
[0285] Example 5 Sample preparation, sequencing, variant calling, and sample QC Sample preparation and sequencing were previously described by Dewey et al. (Dewey et al. (2016), Science 354, aaf6814-aaf6814).
[0286] Upon completion of sequencing, raw data from each Illumina Hiseq 2500 run was collected into local buffer storage and uploaded to the DNAnexus platform (Reid et al. (2014)15,30) for automated analysis. Sample-level read files were generated with CASAVA (Illumina Inc., San Diego, CA) and aligned with BWA-mem to GRCh38 (Li and Durbin (2009); Bioinformatics 25, 1754-1760, Li, H. (2013); arXiv q-bio.GN). The resulting BAM files were processed using GATK and Picard to sort and mark duplicates, and locally correct alignments of reads around putative indels. Sequenced variants were annotated in snpEFF (Cingolani et al. (2012); Fly (Austin) 6, 80-92) using Ensembl85 gene definitions to reveal functional effects on transcripts and genes. Gene definitions were restricted to 54,214 transcripts that are protein-coding with annotated starts and stops, corresponding to 19,467 genes.
[0287] We excluded individuals with low-quality DNA sequence data, as indicated by high frequencies of homozygosity, low sequence data coverage, or genetically identified duplications that could not be confirmed as true monozygotic twins. 61,019 exomes remained for analysis. Additional information regarding sample preparation, sequencing, variant calling, and variant annotation is reported in Dewey et al. (2016), Science 354, aaf6814-1 to aaf6814-10.
[0288] Example 6 Principal components and ancestry estimation We merged the dataset with HapMap3 (International HapMap 3 Consortium, Altshuler et al. (2010), Nature Publishing Group 467, 52-58) using PLINK v1.9, retaining only SNPs that were present in both datasets. We also applied the following PLINK filters: --maf 0.1 --geno 0.05 --snps-only --hwe 0.00001. Principal component (PC) analysis of the HapMap3 samples was computed and each sample was populated into their PCs in the dataset using PLINK. The PCs of the HapMap3 samples were used to train kernel density estimates (KDEs) for each of the five ancestry superclasses, namely African (AFR), Admixed American (AMR), East Asian (EAS), European (EUR), and South Asian (SAS). We used the KDE to calculate the likelihood that each sample belonged to each superclass. For each sample, an ancestry superclass was assigned based on the likelihood. If a sample had two ancestry groups with a likelihood greater than 0.3, it was assigned AFR in preference to EUR, AMR in preference to EUR, AMR in preference to EAS, SAS in preference to EUR, and AMR in preference to AFR; otherwise it was assigned “unknown” (this was done to provide a rigorous estimate of EUR and EAS populations, as well as a comprehensive estimate of highly admixed populations in the dataset). If zero or more than two ancestry groups had a sufficiently high likelihood, the sample was assigned “unknown” ancestry. Samples with unknown ancestry were excluded from ancestry-based homosomal ancestry (IBD) calculations.
[0289] Example 7 IBD Prediction We filtered out high quality and common variants by running PLINK on the entire dataset with the following flags: --maf 0.1 --geno 0.05 --snps-only --hwe 0.00001. We then took a two-pronged approach to obtain accurate IBD estimates from the exome data. First, IBD estimates between individuals were calculated within ancestry superclasses determined from the ancestry analysis (e.g., AMR, AFR, EAS, EUR, and SAS). The following PLINK flags were used to obtain IBD estimates up to second-degree relationships: --genome --min 0.1875. This allowed for more accurate relatedness estimates, as all samples shared similar ancestral alleles. However, this approach could not predict relationships between individuals of different ancestries, for example, children with a European father and an Asian mother.
[0290] Next, IBD estimates were calculated between all individuals using the --min 0.3 PLINK option to capture first-degree relationships between individuals of different ancestry. Individuals were then grouped into first-degree family networks, where network nodes are individuals and edges are first-degree relationships. Each first-degree family network was run through the prePRIMUS pipeline (Staples et al. (2014); Am. J. Hum. Genet. 95, 553-564) to match sample ancestry to appropriate ancestral minor allele frequencies to improve IBD estimates. This process accurately estimated first- and second-degree relationships between individuals within each family network (minimum PI_HAT of 0.15).
[0291] Finally, we combined the IBD estimates from the above two approaches by adding the missing relationships from the family network-derived IBD estimates to the ancestry-based IBD estimates. This approach yielded accurate IBD estimates for second-degree relationships among all samples of similar ancestry and for first-degree relationships among all samples.
[0292] The IBD rate for third-degree relatives is difficult to accurately estimate from large exome sequencing datasets with diverse lineages, as the analysis often results in an over-number of expected third-degree relationships due to artificially inflated IBD estimates. During the ancestry-specific IBD analysis, a --min 0.09875 cutoff was used to capture the number of third-degree relationships that may be present in the DiscovEHR and extended DiscovEHR cohorts, but these were not used in any phasing or pedigree-based analyses. Rather, the relationship-based analyses disclosed herein used only high-confidence third-degree relationships identified within first- and second-degree family networks.
[0293] Example 8 Family reconstruction All first-degree family networks identified within the DiscovEHR and Extended DiscovEHR cohorts were reconstructed in PRIMUS v1.9.0. Joint IBD estimates were provided to PRIMUS within the range of genetically derived sex and EHR-reported age. A relatedness cutoff of PI_HAT>0.375 was specified to restrict reconstruction to first-degree family networks, and a cutoff of 0.1875 was specified to define second-degree networks.
[0294] Example 9 Allele Frequency-Based Phasing EAGLEv2.3 (Loh et al. (2016); Nat Genet 48, 1443-1448) was used to phase all biallelic variants from 61,019 exomes. To parallelize the analysis within DNAnexus, the genome was divided into overlapping segments of approximately 40K variants with a minimum overlap of 500 variants and 250K base pairs. Because the goal was to phase putative compound heterozygous mutations within genes, care was taken to ensure that segment breakpoints occurred in intergenic regions.
[0295] We performed a liftover of the EAGLE-provided genetic_map_hg19.txt.gz file from hg19 to GRCh38 to remove all variants that swapped chromosomes or changed the relative order within a chromosome such that both the chromosomal location and the cM location were not in ascending order when sorted. In most cases, this QC step was to remove pericentromeric inversions. We also removed SNPs that mapped to alternative chromosomes. In total, of the 3.3 million SNPs, only 2,783 were removed from the genetic map file. Data for each segment were provided to EAGLE as PLINK format files and run on DNAnexus with the following EAGLE command line parameters:
[0296] --geneticMapFile=genetic_map_hg19_withX.txt.GRCh38_liftover.txt.gz --maxMissingPerIndiv 1 --genoErrProb 0.01 --numThreads=16 Example 10 Compound heterozygosity calls The objective was to obtain high-confidence compound heterozygous mutation (CHM) calls for putative loss-of-function (pLoF) variants to identify humans with both copies of a potentially knocked-out or perturbed gene. Variants were classified as pLoF if they were a frameshift, stop codon gain, stop codon loss, start codon gain, start codon loss, or splicing acceptor or donor altering variant. A second, expanded set of potentially deleterious variants was created, which included pLOF and likely perturbing missense variants, defined by being predicted to be deleterious by all of the following five methods: SIFT (Loh et al. (2016); Nat Genet 48, 1443-1448) (damaging), PolyPhen2 HDIV (damaging and possibly damaging), PolyPhen2 HVAR (damaging and possibly damaging), LRT (deleterious), and MutationTaster (Schwarz et al. (2014); Nat. Methods 11, 361-362) (disease causing automatic and disease causing).
[0297] Rare (alternative allele counts <1%) cryptic compound heterozygous variants (pCHMs) were identified by testing all possible combinations of heterozygous pLoFs and / or deleterious missense variants within the same gene in the same individual. -15(Chang et al.(2015);Gigascience4,7.)) We filtered out all variants with >10% deletion across 61K samples or with another variant within 10 base pairs of the same individual. We also filtered out SNPs with QD<3, AB<15%, and read depth<7, and INDELs with QD<5, AB<20%, and read depth<10. After filtering, we obtained 39,459 high-quality pCHMs distributed across 25,031 individuals that could knock out or disrupt the function of both copies of a human gene when the pCHM variants are phased in trans.
[0298] The next step was to phase the pCHMs. A combination of population allele frequency-based phasing and ancestry / relationship-based phasing using EAGLE was used to determine whether the pCHMs were cis or trans. Figure 9 shows the pCHM phasing workflow used to obtain the most accurate phasing of each pCHM in the DiscovEHR dataset. Figure 2 shows the pCHM phasing workflow used to obtain the most accurate phasing of each pCHM in the extended DiscovEHR dataset. Ancestry and relationship phasing proved to be more accurate than EAGLE phasing, so ancestry and relationship data were used for phasing. Table 18 below shows the logic used to determine the phase of pCHMs for different types of familial relationships. For all remaining pCHMs, the EAGLE phase data described above was used. We excluded pCHMs phased with EAGLE where either or both variants were singletons because EAGLE's phasing accuracy on singletons was not significantly different from random guessing (see Table 19 below for the DiscovEHR dataset and Table 20 for the Extended DiscovEHR dataset). In the DiscovEHR dataset, when two variants in pCHM have the same minor allele count (MAC) below 100, they are almost always in cis (36 of 37 among our trios), exceeding the accuracy of EAGLE's pCHM phasing. In the Extended DiscovEHR dataset, when two variants in pCHM have the same minor allele count (MAC) below 100, they are almost always in cis (22 of 22 among trios), exceeding the accuracy of EAGLE's pCHM phasing.
[0299] [Table 18-1]
[0300] [Table 18-2]
[0301] These rules were applied in order, starting with the top rule for each relationship. A "?" result means that the pCHM could not be phased. A "NA" result indicates that the result should not have occurred and is likely the result of a sequencing error or the propagation of other non-Mendelian variants. PC_rel refers to non-pCHM carriers in the parent-child relationship. "Rare" refers to a MAF less than 1%, which includes all variants as used herein.
[0302] [Table 19]
[0303] All pCHMs with MAF below 1% were binned by the less frequent of the two variants that make up the pCHM. Correct calls and accuracy were determined by comparing EAGLE phasing of pCHMs to phasing determined in trios. The number of pCHMs with incorrect EAGLE phasing determined to be in cis or trans within trios is also provided. pCHMs where one or both variants were determined to be de novo in the child of the trio were excluded. All pCHMs with MAC >6 had similar accuracy in the mid to high 90s, but a decrease in accuracy was observed for MACs between 2 and 6. EAGLE singleton phasing did not perform significantly better than random guessing, so EAGLE-phased singletons were excluded from the results of phased pCHMs, as was done when measuring the overall accuracy of EAGLE-phased pCHMs.
[0304] [Table 20]
[0305] All pCHMs with MAF less than 1% were binned by the less frequent of the two variants that make up the pCHM. Correct calls and accuracy were determined by comparing EAGLE phasing of pCHMs to phasing determined in trios. We also provided the number of pCHMs with incorrect EAGLE phasing determined to be in cis or trans within trios. pCHMs where one or both variants were determined to be de novo in a child of a trio were excluded. It can be seen that all pCHMs with MAC>9 have similar accuracy90, but EAGLE pCHM phasing accuracy decreases with MACs between 2 and 9. Because phasing of pCHMs with EAGLE singletons does not perform significantly worse than random guessing, EAGLE-phased singletons were excluded from the results of phased pCHMs, as was done when measuring the overall accuracy of EAGLE-phased pCHMs. * Due to low phasing accuracy of singleton variants, 2,838 pCHMs containing singleton variants were removed, and the remaining 401 singleton variants were therefore phased with trio and relationship data only.
[0306] To measure the exact accuracy of EAGLE's pCHM phasing across the entire dataset, we ran EAGLE across the entire dataset, excluding all first-degree relatives of one child in each nuclear family before phasing. This removal was necessary because including parental haplotypes improves the accuracy of phasing children of trios compared to samples with no parents in the dataset.
[0307] Finally, when an individual had multiple pCHMs in the same gene, we retained only the pCHM with the most deleterious profile (see Table 21 below).We were able to phase >99% of all pCHMs and identify 13,335 rare compound heterozygous mutations (CHMs).
[0308] [Table 21]
[0309] When an individual had two or more trans pCHMs in the same gene, the values in this table were used to identify and retain the most deleterious pCHM. The effect score was calculated by adding the functional effect scores of the two variants and penalizing the pair if they did not affect all gene transcripts. The pCHM with the lower score was predicted to be the most deleterious and was retained.
[0310] Example 11.1 Compound heterozygous mutation validation of DiscovEHR dataset The accuracy of phasing was evaluated by comparing the phasing predictions with phasing performed using trios and Illumina reads. First, the phasing accuracy of pCHMs was evaluated using the phased pCHMs of the trios as the truth. The phasing approach for each family relationship was performed independently of the phasing of the trios, so the accuracy of phasing for each relationship class could be properly measured as long as the pCHM carrier was a child of the trio. Tables 4 and 12 above show that the accuracy of family relationship-based phasing was 100% accurate for these rare pCHMs. EAGLE phasing was less accurate at 91.4% and 89.1% for the DiscovEHR and Extended DiscovEHR datasets, respectively. For the DiscovEHR dataset, the accuracy of EAGLE in phasing pCHMs was evaluated over a range of different minor allele frequencies, consistently yielding an accuracy of 95% or higher for MACs greater than 6 and around 77% for MACs between 2 and 6 (see Table 19 above). EAGLE fading did not work with singletons as expected.
[0311] We next attempted to validate the 200 pCHMs with short Illumina reads (~75 bp) by looking at the read stacks in the Integrated Genomics Viewer (IGV) (Robinson et al. (2011); Nat. Biotechnol. 29, 24-26) to see if the two variants occur in the same read or independently. During this validation process, we found that a pCHM consisting of two deletions where the end of the first deletion is within 10 bps of the second deletion was in fact one large deletion miscalled as two separate deletions (N=1,109 out of 39,459 pCHMs). As only 15 were phased as in trans (~0.1% of the entire CHM dataset), these pCHMs were not excluded from the overall analysis, but were excluded when we selected the 200 pCHMs for validation. Of the 200 randomly selected pCHMs using short reads, we were able to use the reads to definitively phase 190. The remaining 10 cases showed lead evidence of both cis and trans phasing, most likely indicating that one or both of the variants were false positive calls.
[0312] Example 11.2 Compound heterozygous mutation validation in the extended DiscovEHR dataset For the DiscovEHR dataset, Table 12 above shows that the accuracy of family-based phasing of rare pCHMs was 99.6% (1060 / 1064 pCHMs). EAGLE phasing was less accurate at 89.1% (766 / 860 pCHMs; Table 12 above). Evaluating EAGLE's pCHM phasing accuracy at various ranges of minor allele frequencies, we found that EAGLE consistently achieved >90% accuracy for MACs >9 and approximately 77% accuracy for MACs between 2 and 9 (see Table 20 above). EAGLE phasing performed poorly for singletons.
[0313] We next attempted to validate the 200 pCHMs with short (975 bp) Illumina reads by looking at the read stacks in the Integrated Genomics Viewer (IGV) (Robinson et al. (2011); Nat. Biotechnol. 29, 24-26) to see if the two variants occurred in the same read or independently. Using the short reads, we selected 190 pCHMs (115 cis and 79 trans, 126 EAGLE phased and 74 ancestry or relationship phased). The remaining 10 cases showed read evidence of both cis and trans phasing, most likely because one or both of the variants were false positive calls. Visual validation showed an overall accuracy of 95.8% and 89.9% for ancestry and relationship phasing and EAGLE phasing, respectively (see Table 22). Although Illumina read-based validation results were consistent with the Trio validation results, Illumina read-based validation accuracy results were lower than the Trio phasing accuracy, and this difference is likely due to the enrichment of false-positive pCHMs in problematic small exonic regions that are prone to sequencing and variant calling errors.
[0314] [Table 22]
[0315] 200 pCHMs were randomly selected from 92K extended DiscovEHR participants, where both variants were within 75 base pairs of each other, and phased visually by examining the read stacks that spanned the two variants. Ten (5%) could not be confidently phased using the read stacks, either because there were no overlapping reads in both variants or because the reads gave conflicting results (i.e., some reads showed cis and some showed trans).
[0316] Example 12 De novo mutation (DNM) detection We combined the results of two different approaches to detect DNMs. The first method is TrioDeNovo (Wei et al. (2015); Bioinformatics 31, 1375-1381), which reads the likelihood of the child and parent genotypes for each variant site in the child. These likelihoods were input into a Bayesian framework to calculate the posterior likelihood that the variant in the child is a DNM. The second program is DeNovoCheck (https: / / sourceforge.net / projects / denovocheck), which is described in the Supplementary Methods of de Ligt, et al. (de Ligt et al. (2012); N. Engl. J. Med. 367, 1921-1929). DeNovoCheck takes in a set of candidate DNMs identified as not present in either parent and called in the child. It then checks whether the variant is present in the child and whether it is absent in both parents by examining the BAM file. These potential DNMs were filtered and the confidence level of each DNM in the combined set was assessed using various QC metrics. Figure 7 illustrates this DNM calling process, shows the variant filters that were applied, and provides the criteria used to classify each DNM as either low, medium, or high confidence.
[0317] Example 13.1 Testing the correlation between parental age at conception and number of DNMs in children in the DiscovEHR dataset In this analysis, samples with 10 or more DNMs were excluded as outliers (N=6 excluded samples). They are likely to represent technical artifacts, or somatic mutations. Maternal and paternal ages were highly correlated (rho=0.78, p=1.2x10^-262); when modeled jointly, neither was significant due to collinearity (0.0053 maternal DNMs / year, p=0.48; 0.0076 paternal DNMs / year, p=0.26; Poisson regression) (Figures 36A and 36B). Next, parental age difference (paternal and maternal ages) was tested along with either maternal or paternal age at birth. Both paternal and maternal ages were found to be equally predictive of DNM number (i.e., age difference was not significantly associated with number of DNMs, given maternal or paternal age).
[0318] An increase in exonic DNM numbers was also observed with respect to both maternal (0.012 DNM / year, p=0.011; Poisson regression; Figure 37) and paternal (0.011 DNM / year; p=0.007) age at birth, consistent with other reports (Deciphering Developmental Disorders Study (2017). Nature 542, 433-438; Kong et al. (2012) Nature 542, 433-438; Rahbari et al. (2016) Nat. Genet. 48, 126-133; and Wong et al. (2016) Nat. Commun. 7, 10486). In particular, maternal and paternal age at birth were highly correlated in the dataset (rho=0.78, p=1.2x10^-262; Figure 38), and thus the ratios were not additive and no significant differences were identified to distinguish them as driving factors.
[0319] Example 13.2 Testing the correlation between parental age at conception and number of DNMs in children in the extended DiscovEHR dataset The Expanded DiscoverEHR cohort tested for a correlation between parental age at conception and the number of DNMs in the child, and showed similar results to the DiscovEHR cohort. Consistent with other reports, maternal age at birth (0.011 DNM / year, p=7.3x10-4 ; Poisson regression; Figure 37) and paternity (0.010 DNM / year; p=5.6x10 -4 An increase in the number of exonic DNMs was observed for both maternal and paternal age at birth. Notably, maternal and paternal age at birth were highly correlated in the data set (r=0.79; Figure 39), and thus the ratios are not additive, nor are there significant differences that could be distinguished as driving factors.
[0320] Paternal age was correlated with the number of DNMs per capita using a Poisson distribution (n=2587, coefficient=0.010, p=5.67E-4). Similarly, maternal age was correlated with the number of DNMs per capita using a Poisson distribution (n=2587, coefficient=0.011, p=7.35E-4). Furthermore, paternal and maternal ages were also correlated with each other (R 2 =0.79;p<10E-308).
[0321] Pathogenicity of DNMs was predicted using functional prediction algorithms, namely SIFT (damaging), PolyPhen2HDIV (damaging and possibly damaging), PolyPhen2HVAR (damaging and possibly damaging), LRT (deleterious), and MutationTaster (Schwarz et al. (2014); Nat. Methods 11, 361-362) (disease causing automatic and disease causing). Pathogenicity predictions of DNMs differ significantly from those of random variant distribution (Figure 38). A higher percentage of DNMs also showed concordant predictions of non-pathogenicity. DNMs are 1.8 times more likely to be predicted as pathogenic by the 5 / 5 algorithm. Random variants are 1.5 times more likely to have discordant pathogenicity predictions.
[0322] Example 14 LDLR tandem duplication inference for distant relatives Although it is not possible to know the true family history of de-identified individuals in this cohort, using PRIMUS-reconstructed pedigrees (Staples et al. (2014); Am. J. Hum. Genet. 95, 553-564.), ERSA distant relationship inference (Huff et al. (2011); GenomeRes. 21, 768-774.), and the ability of PADRE to connect pedigrees (Staples et al. (2016); The American Journal of Human Genetics 99, 154-162), we identified the best pedigree representation of mutation carriers of a novel tandem duplication in LDLR (Maxwell et al. (2017). Profiling copy number variation and disease associations from 50,726 DiscovEHR Study exomes). HumanOmniExpress sequence data have previously been used to infer more distant relationships.
[0323] Example 15 SimProgeny SimProgeny can simulate populations of millions of people distributed across one or more subpopulations and track their descendants over hundreds of years. To find a good balance between simplicity and realism, we selected several key population-level parameters that users can adjust (see Table 23 below). These parameters were chosen to provide a good approximation of real populations and family descent structures while keeping the simulation tool relatively simple. Default values are based on US demographics. The default values are set to work for both cohorts, and these parameters can be easily customized to model various populations by modifying the configuration files included in the SimProgeny code (web resource). For more details on the population simulation process, see Example 17.
[0324] [Table 23]
[0325] For the framework developed for the DiscovEHR cohort, the end of reproduction was age 49 years, and for the framework developed for the extended DiscovEHR cohort, the end of reproduction was age 50 years.
[0326] In addition to modeling populations, SimProgeny simulates two ascertainment approaches that model the selection of individuals from a population for genetic studies: random ascertainment and clustered sampling. Random ascertainment gives each individual in a population an equal chance of being ascertained without replacement. Clustered sampling is an approach that enriches for close relatives, and is done by randomly selecting individuals with many first- and second-degree relatives. The number of first-degree relatives is determined by sampling values from a Poisson distributed with a user-specified first-degree ascertainment lambda (default 0.2). The number of second-degree relatives is determined in the same way, with a default second-degree ascertainment lambda of 0.03. For more information on SimProgeny's ascertainment options, see Example 17.
[0327] Example 16 Simulation and validation of the underlying DiscovEHR population To avoid overcomplicating the simulation model, the simulation included populations with starting sizes of 200K, 300K, 400K, 450K, 500K, 550K, 600K, and 1000K individuals. The SimProgeny parameters (see Table 23 above) were adjusted to match publicly available national, state, and county-level data and our own understanding of how individuals were ascertained through the GHS. Sources for selected parameters can be found in the supplemental file Simulation_parameters.xls. We reduced the average migration rate for the state of Pennsylvania (PA) because the GHS primarily targets rural areas, which tend to have lower migration rates than urban areas. The simulation was run with a burn-in period of 120 years and continued for 101 years afterward. The simulated population grew by approximately 15%, which is similar to the growth of PA since the mid-20th century.
[0328] We performed both random and clustered ascertainment. For both ascertainment approaches, the ascertainment order (specified by the ordered_sampling_proportance parameter) of the first 5% of the population was shuffled to model the random sequencing order of individuals in the GHS biobank at the start of the collaboration. This parameter choice had no effect on random ascertainment and negligible effect on the accumulation of pairwise relationships in clustered ascertainment, but it did affect the proportion of individuals with one or more relatives in the dataset ascertained by clustered sampling by producing an inflection point (more pronounced at higher lambda values). This inflection point would be less pronounced when modeling the freezing process of real data and a smoother transition between the sequencing samples from the biobank and the newly ascertained individuals.
[0329] Example 17 The SimProgeny population and validation simulation process The simulation began by initializing a user-specified number of subpopulations and sizes. Ages were initially assigned between 0 and the maximum reproductive age (default 49). Individuals in the population existed in one of three age-specific pools: adolescents, fertile individuals, and elderly individuals. Individuals were assigned to the subpopulation's adolescent pool if they were under reproductive age (default 15), and to the subpopulation's mating pool if they were within the reproductive age range (default 15-49). Individuals moved from the juvenile pool to the mating pool when they reached or exceeded the minimum reproductive age. Similarly, they moved from the mating pool to the elderly population once they exceeded reproductive age. Individuals who emigrated or died were removed from all age groups. After establishing the initial population, the simulation underwent a 120-year burn-in phase to establish family relationships and age distributions that more closely matched the input parameters, while requiring equal births and deaths and zero net migration rates. After burn-in, the simulation was run for a specific number of years with a given population growth and migration rate. Simulations were run in 1-year increments, with the following steps performed within each subpopulation for each year unless otherwise stated.
[0330] 1. Age - Moving individuals who have aged out of an age pool to the next age pool. 2. Mating - Simulates single men and women participating in monogamous marriages. This process is important to obtain a realistic number of full-sib relationships. Male and female pairs are randomly selected from a single reproductive-age male-female pool, and they successfully marry based on the likelihood of marrying at the age specified in the male and female "marriage by age" parameters. Pairs are drawn until the number of successful marriages defined by the marriage rate is reached. Couples are limited to relatives more distant than first cousins. During the burn-in phase, the marriage rate is doubled until it reaches the initial marriage rate specified by the user (default is 66% of the mating pool being married).
[0331] 3. Splitting - Simulates a man and woman breaking up a marriage with a specified divorce rate. Couples are randomly selected and both individuals are marked as single. 4. Mating - Simulates all possible reproductions that can occur in a population over a one year period. Mother and father pairs are randomly selected from either a single birth age pool or a married pool with a ratio defined by the full sibship rate (default is 88% of all births are to married couples). Pairs are drawn and attempts to reproduce are made until a target number of successful conceptions is reached (default birth rate is 0.0219 births per person). The chance of successful conception is based on the mother's age and corresponding birth rate. Parents are restricted to being related no more than first cousins, and all individuals are restricted to having one child per year.
[0332] 5. Selection - Simulate individuals who die. Mortality rates (default 0.0095 deaths per person) are used to determine the expected number of deaths in the population in a given year. Age-specific male and female mortality parameters are used to weight the likelihood of a randomly selected individual dying. If a random number between 0 and 1 exceeds the probability of dying at that person's age, the individual is kept and another individual is selected. Unlucky individuals are added to the pool of dead and removed from the other pool of living. All individuals over the age of 120 are automatically added to the pool of dead and counted towards the annual target number of deaths.
[0333] 6. Emigration - Simulate migration between populations. Emigration is done by randomly selecting one individual from the mating pool and, if the individual is married and of reproductive age, removing him from the population along with his spouse. The proportion of young and old people leaving is recorded along with the number of reproductive-age couples. Emigration is done in a way that maintains the age distribution and the number of reproductive-age couples. First, a juvenile is randomly selected from the existing population and a new individual of the same sex and age is added to the youth pool, and the process is repeated until the appropriate proportion of juveniles has been added. The same process is repeated for the old. Next, two reproductive-age individuals are selected from the existing population and two new individuals are added with the corresponding ages. One is assigned to male and the other to female, and the two migrants marry. This step is repeated until the number of couples is replenished. Finally, reproductive-age individuals are added in the same process as adding new young people, and this is repeated until the target number of migrants is achieved. This process helps to maintain the age and sex distribution of the population, as well as the proportion of reproductive-age individuals who are married.
[0334] 7. Transplantation - Simulates people moving within a subpopulation. To simulate the lack of genetic isolation between subpopulations, individuals are allowed to move between subpopulations within the whole population. A single migration rate is used across the whole population. Individuals from a subpopulation are randomly selected and randomly assigned to one of the other subpopulations until the desired number of transplants is achieved. This procedure does not occur if there is only one subpopulation or if the transplant rate is 0 (default is 1% of the whole population transplants).
[0335] The simulation proceeds for a specified period of time, tracking each founder and its descendants. Both random and clustered ascertainment were performed. For both ascertainment approaches, the ascertainment order (specified by the ordered_sampling_proportance parameter) of the first 5% of the population was shuffled to model the random sequencing order of individuals in the GHS biobank at the start of the collaboration. This parameter choice had no effect on random ascertainment and negligible effect on the accumulation of pairwise relationships in clustered ascertainment, but it did affect the proportion of individuals with one or more relatives in the dataset ascertained with clustered sampling by producing an inflection point (more pronounced at higher lambda values). This inflection point would be less pronounced if modeling a real data freezing process and a smoother transition between sequencing samples from the biobank and newly ascertained individuals. Users can specify a subpopulation ascertainment order if they wish to simulate ascertaining from one or more subpopulations before moving on to the next set of subpopulations. By default, we first group all subpopulations together and ascertain them as if they were a single set of individuals. The user can also specify the initial percentage of the population that is identified before moving on to other subpopulations or the entire population. The program produced output for the entire population in the form of a ped file, a list of samples that were identified in the order in which they were identified, and several result files summarizing useful population and ascertainment statistics.
[0336] Example 18 Methods that use pedigree structure to aid in the identification of the genetic cause of a given phenotype typically involve innovative variations on association mapping, linkage analysis, or both. Such methods include MORGAN31, pVAAST15, FBAT (www.hsph.harvard.edu / fbat / fbat.htm), QTDT (csg.sph.umich.edu / abecasis / qtdt / ), ROADTRIPS, rareIBD, and RV-GDT. The appropriate method to use depends on the phenotype, mode of inheritance, ancestry, pedigree structure / size, number of pedigrees, and size of the unrelated dataset. In addition to using associations and pedigrees to directly examine gene-phenotype associations, they can also be used in many other methods to generate additional or improved data, namely pedigree-aware imputation, pedigree-aware phasing, Mendelian error checking, compound heterozygous knockout detection and de novo variant calling, and variant call validation.
[0337] The present disclosure is not limited to the exemplary embodiments described and illustrated above, but variations and modifications are possible within the scope of the appended claims. The technical ideas that can be understood from the above-described embodiment will be described below as supplementary notes. [Appendix 1] 1. A method for phasing genetic variants in a population by leveraging relatedness of the population, comprising: establishing an ancestor superclass designation for each of one or more of said samples; generating a first homozygosity estimate for the subject within an ancestral superclass; generating a second homoeologous estimate for the subject independent of the subject's ancestral superclass; clustering the subjects into a first first-degree family network based on one or more of said second homozygosity estimates; generating a third homozygosity estimate for the subject within the first-degree family network; combining the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; constructing a second first-degree family network for the subject based on pooled homoeologous inference; and The variants are phased as either compound heterozygous mutations (CHM) or not CHM according to the combined homoeologous estimate and the second first-degree family network, or the variants are phased as either compound heterozygous mutations (CHM) or not CHM according to the combined homoeologous estimate and the second first-degree family network, The method further comprises identifying the mutation as a de novo mutation (DNM). [Appendix 2] The method of claim 1, wherein integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to enhance the first homoeologous estimate. [Appendix 3] Phasing variants as compound heterozygous mutations (CHMs) phasing variants according to population allele frequency; removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both; removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of about 5 or less, or an alternative allele balance (AB) of about 10% or less, or a combination thereof; removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of about 5 or less, or an AB of about 10% or less, or a combination thereof; Selecting the remaining variants as potential compound heterozygous mutations (pCHMs), where more than one pair of variants exist within the same sample and within the same gene; and 3. The method of claim 1 or 2, comprising phasing pCHM as cis pCHM or trans pCHM, and then classifying said pCHM phased as trans pCHM as CHM. [Appendix 4] Phasing variants as compound heterozygous mutations The method of claim 3, comprising removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of less than about 10, or an AB of about 20% or less, or a combination thereof. [Appendix 5] To assign a score to CHM according to the severity of the phenotype, and The method of Claim 3 or Claim 4, further comprising selecting the CHM with the highest phenotypic severity score per gene per sample, such that if the person has multiple CHMs within the same gene, the CHM most likely to result in protein function inhibition is identified. [Appendix 6] Phasing the variant as a de novo mutation includes identifying the variant in a sample of the second first-degree family network and its trios; assigning genotype likelihood scores to variants in parent samples and corresponding child samples in a trio, calculating the probability that the variant is a de novo mutation, and identifying the variant as a putative de novo mutation if the calculated probability is statistically significant; identifying a variant in a child sample of the trio and, if the variant is not present in either parent sample of the trio, identifying the variant as a putative de novo mutation; filtering the identified putative de novo mutations by removing putative de novo mutations having a genotype quality (GQ) annotation of less than about 35 in the child sample, or a putative de novo mutation having an alternative allele count (AC) of 10 or more across the samples, or a putative de novo mutation having a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or a putative de novo mutation having an allele balance (AB) of more than about 2% in either of the parent samples, or a putative de novo mutation having an allele balance (AB) of less than about 15% in the child sample, or a putative de novo mutation having an AB of more than about 90% in the child sample, or a putative de novo mutation having alternative allele homozygosity in either of the parent samples, or a combination thereof; 6. The method of any one of claims 1 to 5, comprising combining the filtered identified putative de novo mutations, thereby forming a putative de novo mutation dataset. [Appendix 7] 7. The method of claim 6, further comprising classifying a putative de novo mutation in the putative de novo mutation dataset as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is greater than or equal to about 0.15 in the child samples and less than or equal to about 0.02 in each of the parent samples, and if the mapping quality of the putative de novo mutation is not less than about 40, and if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if the MAC of the putative de novo mutation is less than about 20 across the samples, and if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with a monopolymer run greater than about 4. [Appendix 8] 8. The method of claim 7, further comprising classifying a medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality annotation of about 90 or more in the parent samples, and if the medium-confidence de novo mutation has a read depth of about 10 or more in each of the parent samples, and if the medium-confidence de novo mutation has an alternative read depth of about 7 or more in the child sample, and if the medium-confidence de novo mutation has a SNP with a QD of greater than about 3, and if the medium-confidence de novo mutation has an INDEL with a QD of greater than about 5. [Appendix 9] The method of any one of claims 1 to 8, further comprising removing low-quality samples from the dataset, the low-quality samples being samples with a D-stat greater than 0.12 or a 20x read coverage less than 75%, or both. [Appendix 10] 1. A method for identifying compound heterozygous mutations (CHM) in a population, comprising: identifying variants in DNA sequence samples from a plurality of human subjects; establishing an ancestral superclass designation for the subject based on the identified variants; generating a first homozygosity estimate for the subject within an ancestral superclass; generating a second homoeologous estimate for the subject independent of the subject's ancestral superclass; clustering the subjects into a first first-degree family network based on one or more of said second homozygosity estimates; generating a third homozygosity estimate for the subject within the first-degree family network; combining the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; constructing a second first-degree family network based on the integrated homoeologous estimates; phasing variants in the samples according to population allele frequency; classifying the phased variants as potential CHMs based on the presence of two or more variants in the same subject and gene; and The method includes phasing potential CHMs as being in cis or in trans with another variant in the same subject and gene, and then classifying the potential CHMs phased as being in trans as CHMs. [Appendix 11] 11. The method of claim 10, further comprising filtering the identified variants before an ancestral superclass assignment for the subject is established. [Appendix 12] 11. The method of claim 10, further comprising filtering the identified variants before a second homoeologous estimate for the subject is generated. [Appendix 13] Filtering variants may include removing variants having an alternative allele frequency of greater than about 10% across the samples from a plurality of the human subjects, or a p-value of greater than about 10. -6 13. The method of claim 11 or 12, comprising removing variants that are greater than 100% of the total number of samples from said human subjects, that violate Hardy-Weinberg Equilibrium (HWE), or removing variants that are missing calls in more than about 5% of said samples from a plurality of said human subjects, or a combination thereof. [Appendix 14] 14. The method of any one of claims 11 to 13, further comprising removing low quality samples after the identified variants have been filtered. [Appendix 15] The method of claim 14, wherein a low quality sample is a sample having a D-stat greater than 0.12 or a read coverage of less than 75% at 20x, or both. [Appendix 16] The method according to any one of appendix 10 to 15, wherein integrating the first homoeologous estimation and the third homoeologous estimation comprises supplementing the third homoeologous estimation with a pairwise homoeologous estimation specific to the first homoeologous estimation to enhance the first homoeologous estimation. [Appendix 17] 17. The method of any one of claims 10 to 16, wherein homoeologous estimation includes calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. [Appendix 18] 18. The method of any one of claims 10 to 17, further comprising filtering the variants after they have been phased according to population allele frequency. [Appendix 19] 19. The method of claim 18, wherein filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of less than about 5, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of less than about 5, or an AB of about 10% or less, or a combination thereof. [Appendix 20] 20. The method of any one of claims 10 to 19, wherein phasing variants according to population allele frequency comprises dividing a DNA sequence sample of a human subject into genomic segments of approximately equal size and having substantial segment overlap and breakpoints in intergenic regions. [Appendix 21] 21. The method of any one of claims 10 to 20, wherein the potential CHM is phased based on trio data, or parent-child data, or full sibling data, or distant relative data, or a combination thereof, or is phased based on minor allele counts (MAC), or is phased based on population allele frequencies, or a combination thereof. [Appendix 22] 21. The method of any one of claims 10 to 20, further comprising assigning a score to the CHMs according to phenotypic severity, and selecting the CHMs with the highest phenotypic severity score per gene per sample, thereby obtaining a set of medically important mutations. [Appendix 23] 23. The method of any one of appendix 10 to 22, wherein the DNA sequence sample comprises an exome sequence. [Appendix 24] 24. The method of any one of claims 10-23, wherein the plurality of human subjects comprises more than 10K subjects. [Appendix 25] The method of claim 15, wherein the D-stat of a low quality sample is determined by comparing the distribution of actual allele balance of the sample to the distribution of expected allele balance using the Kolmogorov-Smirnov test. [Appendix 26] 20. The method of claim 19, wherein filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of less than about 10, or an AB of about 20% or less, or a combination thereof. [Appendix 27] 1. A non-transitive computer-implemented method for identifying compound heterozygous mutations (CHMs) in a population, comprising: using a data processor of a computing device to identify variants in DNA sequence samples from a plurality of human subjects; using said data processor to establish an ancestry superclass assignment for the subject based on the identified variants; using the data processor to generate a first homozygosity estimate for the subject within an ancestral superclass; using the data processor to generate a second homoeologous estimate for the subject independent of the subject's ancestral superclass; using said data processor for clustering subjects into a first first-degree family network based on one or more of said second homozygosity estimates; using said data processor to generate a third homozygosity estimate for the subject within a first first-degree family network; using the data processor to combine the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; using said data processor to construct a second first-degree family network based on the pooled homoeologous estimates; using said data processor for phasing variants in the sample according to population allele frequencies; using said data processor to classify the phased variants as potential CHMs based on the presence of two or more variants in the same subject and gene; and The non-transitory computer-implemented method includes using the data processor to phase potential CHMs as being in cis or in trans with another variant in the same subject and gene, and then classifying the potential CHMs phased as being in trans as CHMs. [Appendix 28] 28. The non-transitory computer-implemented method of claim 27, further comprising using the data processor to filter identified variants before an ancestral superclass assignment for the subject is established. [Appendix 29] 29. The non-transitive computer-implemented method of claim 27 or claim 28, further comprising using the data processor to filter identified variants before a second homoeologous estimate for the subject is generated. [Appendix 30] Filtering variants may include removing variants having an alternative allele frequency of greater than about 10% across the samples from a plurality of the human subjects, or a p-value of greater than about 10. -6 29. The non-transitory computer-implemented method of claim 28, further comprising removing variants that violate Hardy-Weinberg Equilibrium (HWE), that are greater than 100% or that are missing calls in more than about 5% of said samples from a plurality of said human subjects, or a combination thereof. [Appendix 31] 31. The non-transitory computer-implemented method of any one of claims 28 to 30, further comprising using the data processor to remove low quality samples after the identified variants have been filtered. [Appendix 32] 32. The non-transitive computer-implemented method of claim 31, wherein a low quality sample is a sample having a D-stat greater than 0.12 or a read coverage of 20x less than 75%, or both. [Appendix 33] 33. The non-transitory computer-implemented method of any one of appendices 27 to 32, wherein integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the third homoeologous estimate with a pairwise homoeologous estimate specific to the first homoeologous estimate to augment the first homoeologous estimate. [Appendix 34] 32. The non-transitive computer-implemented method of any one of appendices 27-31, wherein homoeologous estimation comprises calculation of genome-wide IBD values of 0, 1, and 2 between the sample pairs. [Appendix 35] 33. The non-transitory computer-implemented method of any one of claims 27-32, further comprising using the data processor to filter variants after they have been phased according to population allele frequency. [Appendix 36] 36. The non-transitory computer-implemented method of claim 35, wherein filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of less than about 5, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of less than about 5, or an AB of about 10% or less, or a combination thereof. [Appendix 37] 36. The non-transitory computer-implemented method of any one of claims 27 to 35, wherein phasing variants according to population allele frequency comprises dividing the human subject's DNA sequence sample into genomic segments of approximately equal size and having substantial segment overlap and breakpoints in intergenic regions. [Appendix 38] 38. The non-transitory computer-implemented method of any one of claims 27 to 37, wherein the latent CHM is phased based on trio data, or parent-offspring data, or full sibling data, or distant relative data, or a combination thereof, or is phased based on minor allele counts (MACs), or is phased based on population allele frequencies, or a combination thereof. [Appendix 39] 39. The non-transitory computer-implemented method of any one of claims 27 to 38, further comprising using the data processor to score the CHMs according to phenotypic severity and select the CHMs with the highest phenotypic severity score per gene per sample, thereby obtaining a set of medically important mutations. [Appendix 40] 40. The non-transitory computer-implemented method of any one of claims 27 to 39, wherein the DNA sequence sample comprises an exome sequence. [Appendix 41] 40. The non-transitory computer-implemented method of any one of Clauses 27-39, wherein the plurality of human subjects comprises more than 10K subjects. [Appendix 42] 33. The non-transitory computer-implemented method of claim 32, wherein the D-stat of a low quality sample is determined by comparing the distribution of actual allele balance of the sample to the distribution of expected allele balance using a Kolmogorov-Smirnov test. [Appendix 43] 37. The non-transitory computer-implemented method of claim 36, wherein filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of less than about 10, or an AB of about 20% or less, or a combination thereof. [Appendix 44] a data processor; a memory coupled to the data processor; and a program stored in the memory, the program comprising: identifying variants in DNA sequence samples from a plurality of human subjects; establishing an ancestral superclass designation for the subject based on the identified variants; generating a first homozygosity estimate for the subject within an ancestral superclass; generating a second homoeologous estimate for the subject independent of the subject's ancestral superclass; clustering the subjects into a first first-degree family network based on one or more of said second homozygosity estimates; generating a third homozygosity estimate for the subject within the first-degree family network; combining the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; constructing a second first-degree family network based on the integrated homoeologous estimates; phasing variants in the samples according to population allele frequency; classifying the phased variants as potential CHMs based on the presence of two or more variants in the same subject and gene; and The system includes instructions for phasing potential CHMs as being in cis or in trans with another variant in the same subject and gene, and then classifying the potential CHMs phased as being in trans as CHMs. [Appendix 45] 45. The system of claim 44, wherein the program includes instructions for filtering the identified variants before an ancestral superclass assignment for the subject is established. [Note 46] 46. The system of claim 44 or claim 45, wherein the program comprises instructions for filtering the identified variants before a second homoeologous estimate for the subject is generated. [Appendix 47] Filtering variants may include removing variants having an alternative allele frequency of greater than about 10% across the samples from a plurality of the human subjects, or a p-value of greater than about 10. -6 47. The system of any one of claims 44 to 46, comprising removing variants that are greater than 100% of the total number of samples from the human subjects, that are against Hardy-Weinberg Equilibrium (HWE), or that are missing calls in more than about 5% of the samples from the human subjects, or a combination thereof. [Appendix 48] 48. The system of any one of claims 45 to 47, wherein the program includes instructions for removing low quality samples after the identified variants have been filtered. [Appendix 49] The system of claim 48, wherein a low quality sample is a sample having a D-stat greater than 0.12, or a 20x read coverage less than 75%, or both. [Appendix 50] The system of any one of appendix 44 to 49, wherein integrating the first homoeologous estimate with the third homoeologous estimate includes supplementing the third homoeologous estimate with a pairwise homoeologous estimate specific to the first homoeologous estimate. [Appendix 51] 51. The system of any one of claims 44 to 50, wherein homoeologous estimation includes calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. [Appendix 52] 52. The system of any one of appendix 44 to 51, wherein the program includes instructions for filtering variants after the variants have been phased according to population allele frequency. [Appendix 53] 53. The system of claim 52, wherein filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of less than about 5, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 2 or less, or a DP of less than about 5, or an AB of about 10% or less, or a combination thereof. [Appendix 54] 54. The system of any one of claims 44 to 53, wherein phasing variants according to population allele frequency comprises dividing the DNA sequence sample of the human subject into genomic segments of approximately equal size and having substantial segment overlap and breakpoints in intergenic regions. [Appendix 55] The system of any one of claims 44 to 53, wherein the potential CHM is phased based on trio data, or parent-child data, or full sibling data, or distant relative data, or a combination thereof, or is phased based on minor allele counts (MACs), or is phased based on population allele frequencies, or a combination thereof. [Appendix 56] 56. The system of any one of claims 44 to 55, wherein the program includes instructions for scoring the CHMs according to phenotypic severity and selecting the CHMs with the highest phenotypic severity scores per gene per sample, thereby obtaining a set of medically important mutations. [Appendix 57] 57. The system of any one of claims 44 to 56, wherein the DNA sequence sample comprises an exome sequence. [Appendix 58] 58. The system of any one of claims 44 to 57, wherein the plurality of human subjects comprises more than 10K subjects. [Appendix 59] 49. The system of claim 49, wherein the D-stat of a low quality sample is determined by comparing the distribution of actual allele balance of the sample to the distribution of expected allele balance using a Kolmogorov-Smirnov test. [Appendix 60] 54. The system of claim 53, wherein filtering the phased variants according to population allele frequency includes removing variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both, and removing SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or a combination thereof, and removing insertions or deletions (INDELs) with a QD of about 5 or less, or a DP of less than about 10, or an AB of about 20% or less, or a combination thereof. [Note 61] 1. A method for identifying de novo mutations (DNMs) in a population, comprising: identifying variants in DNA sequence samples from a plurality of human subjects; establishing an ancestral superclass designation for the subject based on the identified variants; generating a first homozygosity estimate for the subject within an ancestral superclass; generating a second homoeologous estimate for the subject independent of the subject's ancestral superclass; clustering the subjects into a first first-degree family network based on one or more of said second homozygosity estimates; generating a third homozygosity estimate for the subject within the first-degree family network; combining the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; Constructing nuclear families based on pooled homoeologous estimates; Identifying variants in nuclear families; the method comprising: assigning genotype likelihood scores to variants in samples derived from parents and children of the constructed nuclear family trio to calculate the probability that the variants are de novo mutations; selecting the variants that have a significantly high probability of being de novo mutations; and independently, simply identifying variants called in the child sample that are not called in either parent sample of the trio; and then combining the two sets of de novo mutations, thereby forming a dataset of putative de novo mutations. [Appendix 62] 62. The method of claim 61, further comprising filtering the identified variants before an ancestral superclass assignment for the subject is established. [Appendix 63] 63. The method of claim 61 or claim 62, further comprising filtering the identified variants before a second homoeologous estimate for the subject is generated. [Appendix 64] Filtering variants may include removing variants having an alternative allele frequency of greater than about 10% across the samples from a plurality of the human subjects, or a p-value of greater than about 10. -6 64. The method of claim 62 or 63, comprising removing variants that are greater than 100% of the total number of samples from said human subjects, that are against Hardy-Weinberg Equilibrium (HWE), or that are missing calls in more than about 5% of said samples from a plurality of said human subjects, or a combination thereof. [Appendix 65] 65. The method of any one of claims 62 to 64, further comprising removing low quality samples after the identified variants have been filtered. [Note 66] The method of claim 65, wherein a low quality sample is a sample having a D-stat greater than 0.12, or a read coverage at 20x less than 75%, or both. [Appendix 67] 66. The method of any one of appendix 61 to 65, wherein integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the third homoeologous estimate with a pairwise homoeologous estimate specific to the first homoeologous estimate to enhance the first homoeologous estimate. [Appendix 68] 67. The method of any one of claims 61 to 66, wherein homoeologous estimation includes calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. [Appendix 69] 69. The method of any one of appendix 61-68, wherein the genotype likelihood scores are based on DNA sequence samples from multiple human subjects in multiple nuclear families. [Appendix 70] 69. The method of any one of claims 61 to 68, further comprising filtering the variants after the probability that the variant is a de novo mutation is calculated based on the genotype likelihood score. [Appendix 71] 71. The method of any one of claims 61 to 70, further comprising filtering the variants after the probability that the variant is a de novo mutation is calculated based on simply identifying variants in the child sample that are not present in either parent sample. [Appendix 72] 72. The method of claim 70 or 71, wherein filtering variants comprises removing variants with a genotype quality (GQ) annotation of less than about 35 in the child sample, or variants with an alternative allele count (AC) of 10 or more across the samples, or variants with a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or variants with an allele balance (AB) of more than about 2% in either parent sample, or variants with an allele balance (AB) of less than about 15% in the child sample, or variants with an AB of more than about 90% in the child sample, or variants with alternative allele homozygosity in either parent sample, or combinations thereof. [Appendix 73] 72. The method of any one of claims 61-71, further comprising annotating the variants with quality control metrics. [Appendix 74] 74. The method of any one of claims 61 to 73, further comprising filtering variants based on sample BAM file data after putative de novo mutations are identified based on simply identifying variants in the child sample that are not present in either parent sample. [Appendix 75] 75. The method of any one of appendix 61 to 74, further comprising classifying the putative de novo mutation as a medium-confidence de novo mutation if the allele balance of the putative de novo mutation is about 0.15 or greater in the child sample. [Appendix 76] 76. The method of any one of appendix 61 to 75, further comprising classifying the putative de novo mutation as a medium confidence de novo mutation if the allele balance of the putative de novo mutation is less than or equal to about 0.02 in each of the parent samples. [Appendix 77] 77. The method of any one of appendix 61 to 76, further comprising classifying the putative de novo mutation as a medium confidence de novo mutation if the mapping quality of the putative de novo mutation is not less than about 40. [Appendix 78] 78. The method of any one of appendix 61 to 77, further comprising classifying the putative de novo mutation as a medium confidence de novo mutation if the putative de novo mutation has a quality by depth (QD) value of not less than about 2. [Appendix 79] 79. The method of any one of appendix 61-78, further comprising classifying the putative de novo mutation as a medium confidence de novo mutation if the MAC of the putative de novo mutation is less than about 20 across the samples. [Appendix 80] 80. The method of any one of appendix 61 to 79, further comprising classifying the putative de novo mutation as a medium confidence de novo mutation if the soft clipped read of the putative de novo mutation is about 3 or less at the variant site in a carrier of the putative de novo mutation. [Appendix 81] 81. The method of any one of appendix 61-80, further comprising classifying the putative de novo mutation as a medium confidence de novo mutation if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs. [Appendix 82] 82. The method of any one of appendix 61-81, further comprising classifying the putative de novo mutation as a medium-confidence de novo mutation if the allele balance (AB) of the putative de novo mutation is greater than or equal to about 0.15 in the child sample and less than or equal to about 0.02 in each of the parent samples, and if the mapping quality (MQ) of the putative de novo mutation is not less than about 40, and if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if the minor allele count (MAC) of the putative de novo mutation is less than about 20 across the samples, and if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs. [Appendix 83] 83. The method of claim 82, further comprising classifying a medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality (GQ) annotation of about 90 or greater in the parent samples, and if the medium-confidence de novo mutation has a read depth (DP) of about 10 or greater in each of the parent samples, and if the medium-confidence de novo mutation has an alternative DP of about 7 or greater in the child sample, and if the medium-confidence de novo mutation has a SNP with a QD of greater than about 3, and if the medium-confidence de novo mutation has an INDEL with a QD of greater than about 5. [Appendix 84] 84. The method of any one of claims 61 to 83, wherein the DNA sequence sample comprises an exome sequence. [Appendix 85] 85. The method of any one of claims 61-84, wherein the plurality of human subjects comprises more than 10K subjects. [Appendix 86] 67. The method of claim 66, wherein the D-stat of a low quality sample is determined by comparing the distribution of actual allele balance of the sample to the distribution of expected allele balance using the Kolmogorov-Smirnov test. [Appendix 87] A non-transitive computer-implemented method for identifying de novo mutations (DNMs) in a population, comprising: using a data processor of a computing device to identify variants in DNA sequence samples from a plurality of human subjects; using the data processor to establish ancestral superclass designations for the subjects based on the identified variants; using a data processor to generate a first homozygosity estimate for the subject within an ancestral superclass; using the data processor to generate a second homoeologous estimate for the subject independent of the subject's ancestral superclass; using a data processor to cluster subjects into a first first-degree family network based on one or more of said second homozygosity estimates; using the data processor to generate a third homozygosity estimate for the subject within the first-degree family network; using a data processor to combine the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; using a data processor to construct nuclear families based on pooled homoeologous estimates; using a data processor to identify variants in nuclear families; the non-transitive computer-implemented method comprising: assigning genotype likelihood scores to variants of samples from parents and children of the constructed nuclear family trio to calculate the probability that the variants are de novo mutations; selecting the variants that have a significantly high probability that they are de novo mutations; and independently, simply identifying variants called in the child sample that are not called in either parent sample of the trio; and then combining the two sets of de novo mutations, thereby forming a dataset of putative de novo mutations. [Appendix 88] 88. The non-transitory computer-implemented method of claim 87, further comprising using the data processor to filter identified variants before an ancestral superclass assignment for the subject is established. [Appendix 89] 89. The non-transitive computer-implemented method of claim 87 or claim 88, further comprising using the data processor to filter identified variants before a second homoeologous estimate for the subject is generated. [Appendix 90] Filtering variants may include removing variants having an alternative allele frequency of greater than about 10% across the samples from a plurality of the human subjects, or a p-value of greater than about 10. -6 89. The non-transitory computer-implemented method of any one of claims 88-89, comprising removing variants that violate Hardy-Weinberg Equilibrium (HWE), that are greater than 100% or that are missing calls in more than about 5% of the samples from a plurality of said human subjects, or a combination thereof. [Appendix 91] 91. The non-transitory computer-implemented method of any one of claims 88-90, further comprising using the data processor to remove low quality samples after the identified variants have been filtered. [Appendix 92] 92. The non-transitive computer-implemented method of claim 91, wherein a low quality sample is a sample having a D-stat greater than 0.12, or a 20x read coverage less than 75%, or both. [Appendix 93] 93. The non-transitory computer-implemented method of any one of appendices 87-92, wherein integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the third homoeologous estimate with a pairwise homoeologous estimate specific to the first homoeologous estimate to augment the first homoeologous estimate. [Appendix 94] 94. The non-transitive computer-implemented method of any one of claims 87-93, wherein homoeologous estimation comprises calculation of genome-wide IBD values of 0, 1, and 2 between the sample pairs. [Appendix 95] 95. The non-transitory computer-implemented method of any one of Clauses 87-94, wherein said genotype likelihood scores are based on DNA sequence samples from multiple human subjects in multiple nuclear families. [Appendix 96] 96. The non-transitory computer-implemented method of any one of clauses 87-95, further comprising using the data processor to filter variants after a probability that a variant is a de novo mutation has been calculated based on the genotype likelihood score. [Appendix 97] 97. The non-transitory computer-implemented method of any one of clauses 87-96, further comprising using the data processor to filter variants after a probability that a variant is a de novo mutation has been calculated based on simply identifying variants in the child sample that are not present in either parent sample. [Appendix 98] 98. The non-transitory computer-implemented method of claim 96 or 97, wherein filtering variants comprises removing variants having a genotype quality (GQ) annotation of less than about 35 in the child sample, or variants having an alternative allele count (AC) of 10 or greater across the samples, or variants having a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or variants having an allele balance (AB) of greater than about 2% in either parent sample, or variants having an allele balance (AB) of less than about 15% in the child sample, or variants having an AB of greater than about 90% in the child sample, or variants having alternative allele homozygosity in either parent sample, or combinations thereof. [Appendix 99] 99. The non-transitory computer-implemented method of any one of clauses 87-98, further comprising using the data processor to annotate variants with quality control metrics. [Appendix 100] 99. The non-transitory computer-implemented method of any one of clauses 87-99, further comprising using the data processor to filter variants based on sample BAM file data after putative de novo mutations have been identified based on simply identifying variants in the child sample that are not present in either parent sample. [Appendix 101] 101. The non-transitory computer-implemented method of any one of appendix 87 to 100, further comprising using the data processor to classify a putative de novo mutation as a medium confidence de novo mutation if an allele balance of the putative de novo mutation is greater than or equal to about 0.15 in the child sample. [Appendix 102] 102. The non-transitory computer-implemented method of any one of appendix 87-101, further comprising using the data processor to classify a putative de novo mutation as a medium confidence de novo mutation if an allele balance of the putative de novo mutation is less than or equal to about 0.02 in each of the parent samples. [Appendix 103] 103. The non-transitory computer-implemented method of any one of appendixes 87-102, further comprising using the data processor to classify a putative de novo mutation as a medium confidence de novo mutation if a mapping quality of the putative de novo mutation is not less than about 40. [Appendix 104] 104. The non-transitory computer-implemented method of any one of appendixes 87-103, further comprising using the data processor to classify a putative de novo mutation as a medium confidence de novo mutation if a quality by depth (QD) value of the putative de novo mutation is not less than about 2. [Appendix 105] The non-transitory computer-implemented method of any one of Appendices 87-104, further comprising using the data processor to classify a putative de novo mutation as a medium confidence de novo mutation if the MAC of the putative de novo mutation is less than about 20 across all of the samples. [Appendix 106] 106. The non-transitory computer-implemented method of any one of appendixes 87-105, further comprising using the data processor to classify a putative de novo mutation as a medium confidence de novo mutation if a soft clip read of the putative de novo mutation is about 3 or less at the variant site in a carrier of the putative de novo mutation. [Appendix 107] 107. The non-transitory computer-implemented method of any one of clauses 87-106, further comprising using the data processor to classify a putative de novo mutation as a medium confidence de novo mutation if the putative de novo mutation is not an INDEL with more than about 4 monopolymer runs. [Appendix 108] 108. The non-transitory computer-implemented method of any one of clauses 87-107, further comprising using the data processor to classify a putative de novo mutation as a medium confidence de novo mutation if an allele balance (AB) of the putative de novo mutation is greater than or equal to about 0.15 in the child sample and less than or equal to about 0.02 in each of the parent samples, and if a mapping quality (MQ) of the putative de novo mutation is not less than about 40, and if a quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if a minor allele count (MAC) of the putative de novo mutation is less than about 20 across the samples, and if soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with a monopolymer run greater than about 4. [Appendix 109] 109. The non-transitory computer-implemented method of claim 108, further comprising using the data processor to classify a medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality (GQ) annotation of about 90 or greater in the parent samples, and if the medium-confidence de novo mutation has a read depth (DP) of about 10 or greater in each of the parent samples, and if the medium-confidence de novo mutation has an alternative DP of about 7 or greater in the child sample, and if the medium-confidence de novo mutation has a SNP with a QD of greater than about 3, and if the medium-confidence de novo mutation has an INDEL with a QD of greater than about 5. [Appendix 110] 109. The non-transitory computer-implemented method of any one of claims 87 to 109, wherein the DNA sequence sample comprises an exome sequence. [Appendix 111] 111. The non-transitory computer-implemented method of any one of Clauses 87-110, wherein the plurality of human subjects comprises more than 10K subjects. [Appendix 112] 93. The non-transitory computer-implemented method of claim 92, wherein the D-stat of a low quality sample is determined by comparing the distribution of actual allele balance of the sample to the distribution of expected allele balance using a Kolmogorov-Smirnov test. [Appendix 113] a data processor; a memory coupled to the data processor; and a program stored in the memory, the program comprising: Identifying variants in DNA sequence samples from multiple human subjects; Establishing an ancestral superclass designation for the subject based on the identified variants; generating a first homozygosity estimate for the subject within an ancestral superclass; generating a second homozygosity estimate for the subject independent of the subject's ancestral superclass; clustering the subjects into a first first-degree family network based on one or more of the second homoeologous estimates. Generate a third homozygosity estimate for the subject within the first-degree family network; combining the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; Constructing nuclear families based on integrated homoeologous estimates Identifying variants in nuclear families, The system includes instructions for assigning genotype likelihood scores to variants in samples from parents and children of the constructed nuclear family trio to calculate the probability that the variants are de novo mutations, selecting the variants that have a significantly high probability of being a de novo mutation, and independently simply identifying variants called in the child sample that are not called in either parent sample of the trio, and then combining the two sets of de novo mutations, thereby forming a dataset of putative de novo mutations. [Appendix 114] 114. The system of claim 113, wherein the program includes instructions for filtering the identified variants before an ancestral superclass assignment for the subject is established. [Appendix 115] 115. The system of claim 113 or 114, wherein the program includes instructions for filtering the identified variants before a second homoeologous estimate for the subject is generated. [Appendix 116] Filtering variants may include removing variants having an alternative allele frequency of greater than about 10% across the samples from a plurality of the human subjects, or a p-value of greater than about 10. -6 16. The system of claim 114 or 115, comprising removing variants that are greater than 100% of the total number of variants in the human subject population, that are against Hardy-Weinberg Equilibrium (HWE), or that are missing calls in more than about 5% of the samples from a plurality of the human subjects, or a combination thereof. [Appendix 117] The system of any one of claims 114 to 116, wherein the program includes instructions for removing low quality samples after the identified variants have been filtered. [Appendix 118] The system of claim 117, wherein a low quality sample is a sample having a D-stat greater than 0.12, or a 20x read coverage less than 75%, or both. [Appendix 119] The system of any one of appendix 113 to 118, wherein integrating the first homoeologous estimate with the third homoeologous estimate includes supplementing the third homoeologous estimate with a pairwise homoeologous estimate specific to the first homoeologous estimate. [Appendix 120] The system of any one of appendix 113 to 119, wherein homoeologous estimation includes calculation of genome-wide IBD values of 0, 1, and 2 between sample pairs. [Appendix 121] 121. The system of any one of appendix 113-120, wherein the genotype likelihood scores are based on DNA sequence samples from multiple human subjects in multiple nuclear families. [Appendix 122] 122. The system of any one of appendix 113-121, wherein the program includes instructions for filtering variants after a probability that the variant is a de novo mutation is calculated based on the genotype likelihood score. [Appendix 123] The system of any one of claims 113 to 122, wherein the program includes instructions for filtering variants after a probability that the variant is a de novo mutation has been calculated based on simply identifying variants in the child sample that are not present in either parent sample. [Appendix 124] 124. The system of claim 122 or 123, wherein filtering variants comprises removing variants having a genotype quality (GQ) annotation of less than about 35 in the child sample, or variants having an alternative allele count (AC) of 10 or greater across the samples, or variants having a read depth (DP) of less than about 7 and an alternative DP of less than about 4 in the child sample, or variants having an allele balance (AB) of greater than about 2% in either parent sample, or variants having an allele balance (AB) of less than about 15% in the child sample, or variants having an AB of greater than about 90% in the child sample, or variants having alternative allele homozygosity in either parent sample, or combinations thereof. [Appendix 125] The system of any one of Appendices 113 to 124, wherein the program includes instructions for annotating variants with quality control metrics. [Appendix 126] The system of any one of claims 113 to 125, wherein the program includes instructions for filtering variants based on sample BAM file data after putative de novo mutations have been identified based on simply identifying variants in the child sample that are not present in either parent sample. [Appendix 127] 127. The system of any one of appendices 113 to 126, wherein the program includes instructions for classifying a putative de novo mutation as a medium-confidence de novo mutation if an allele balance of the putative de novo mutation is greater than or equal to about 0.15 in the child sample. [Appendix 128] 128. The system of any one of appendices 113 to 127, wherein the program includes instructions for classifying a putative de novo mutation as a medium-confidence de novo mutation if the allele balance of the putative de novo mutation is less than or equal to about 0.02 in each of the parent samples. [Appendix 129] 129. The system of any one of appendix 113 to 128, wherein the program includes instructions for classifying a putative de novo mutation as a medium-confidence de novo mutation if the mapping quality of the putative de novo mutation is not less than about 40. [Appendix 130] 130. The system of any one of appendices 113 to 129, wherein the program comprises instructions for classifying a putative de novo mutation as a medium-confidence de novo mutation if the putative de novo mutation has a quality by depth (QD) value that is not less than about 2. [Appendix 131] 131. The system of any one of appendix 113-130, wherein the program includes instructions for classifying a putative de novo mutation as a medium confidence de novo mutation if the MAC of the putative de novo mutation is less than about 20 across the samples. [Appendix 132] 132. The system of any one of appendices 113 to 131, wherein the program includes instructions for classifying a putative de novo mutation as a medium confidence de novo mutation if a soft clip read of the putative de novo mutation is about 3 or less at the variant site in a carrier of the putative de novo mutation. [Appendix 133] 133. The system of any one of appendix 113 to 132, wherein the program includes instructions for classifying the putative de novo mutation as a medium confidence de novo mutation if the putative de novo mutation is not an INDEL having more than about 4 monopolymer runs. [Appendix 134] 134. The system of any one of appendix 113 to 133, wherein the program comprises instructions for classifying a putative de novo mutation as a medium confidence de novo mutation if the allele balance (AB) of the putative de novo mutation is greater than or equal to about 15% in the child sample and less than or equal to about 2% in each of the parent samples, and if the mapping quality (MQ) of the putative de novo mutation is not less than about 40, and if the quality by depth (QD) value of the putative de novo mutation is not less than about 2, and if the minor allele count (MAC) of the putative de novo mutation is less than about 20 across the samples, and if the soft clipped reads of the putative de novo mutation are less than or equal to about 3 at the variant site in carriers of the putative de novo mutation, and if the putative de novo mutation is not an INDEL with a monopolymer run greater than about 4. [Appendix 135] 135. The system of claim 134, wherein the program comprises instructions for classifying a medium-confidence de novo mutation as a high-confidence de novo mutation if the medium-confidence de novo mutation has a genotype quality (GQ) annotation of about 90 or more in the parent samples, and if the medium-confidence de novo mutation has a read depth (DP) of about 10 or more in each of the parent samples, and if the medium-confidence de novo mutation has an alternative DP of about 7 or more in the child sample, and if the medium-confidence de novo mutation has a SNP with a QD of greater than about 3, and if the medium-confidence de novo mutation has an INDEL with a QD of greater than about 5. [Appendix 136] 136. The system of any one of claims 113 to 135, wherein the DNA sequence sample comprises an exome sequence. [Appendix 137] 17. The system of any one of claims 113 to 136, wherein the plurality of human subjects comprises more than 10K subjects. [Appendix 138] The system of claim 118, wherein the D-stat of a low quality sample is determined by comparing the distribution of actual allele balance of the sample to the distribution of expected allele balance using a Kolmogorov-Smirnov test.
Claims
1. A method for manufacturing a computer-implemented computer-implemented system, comprising: a processor; a memory coupled to the processor; and a program stored in the memory, the program comprising: identifying variants in DNA sequence samples from a plurality of human subjects; generating a first homozygosity estimate for the subject within an ancestral superclass; generating a second homozygosity estimate for the subject independent of the subject's ancestral superclass; clustering subjects into first-degree family networks based on the first homozygosity estimate and the second homozygosity estimate; generating third homozygosity estimates within a set of first-degree family networks; combining the first homoeologous estimate and the third homoeologous estimate to generate a combined homoeologous estimate; constructing a second first-degree family network for the subject based on the pooled homoeologous estimates; phasing the identified variants as hLoF based on the combined homoeologous estimate and the second first-degree family network; The system further comprises instructions executable by the processor.
2. 2. The system of claim 1, wherein the program comprises instructions for filtering identified variants before the first homoeologous estimate and the second homoeologous estimate are generated.
3. Filtering variants may include removing variants having an alternative allele frequency of greater than about 10% across the samples from a plurality of the human subjects, or a p-value of greater than about 10. -6 3. The system of claim 1 or claim 2, further comprising removing variants that violate Hardy-Weinberg Equilibrium (HWE), that are greater than 100%, or that are missing calls in more than about 5% of the samples from a plurality of said human subjects, or a combination thereof.
4. 4. The system of claim 2 or claim 3, wherein the program comprises instructions for removing samples with D-stat greater than 0.12 or read coverage at 20x less than 75%, or both, after the identified variants have been filtered.
5. The system of any one of claims 1 to 4, wherein integrating the first homoeologous estimate with the third homoeologous estimate comprises supplementing the first homoeologous estimate with a pairwise homoeologous estimate specific to the third homoeologous estimate to enhance the first homoeologous estimate.
6. 1. A computer-implemented method for phasing heterozygous loss-of-function mutations (hLoF), said computer-implemented method comprising: generating, by a processor, a first homoeologous estimate for the subject within an ancestral superclass; generating, by the processor, a second homoeologous estimate for the subject independent of the subject's ancestral superclass; clustering, by the processor, subjects into a first first-degree family network based on the first homozygosity estimate and the second homozygosity estimate; generating, by the processor, a third homoeologous estimate within the set of first first-degree family networks; combining, by the processor, the first homoeologous estimate and the third homoeologous estimate to generate a combined homoeologous estimate; constructing, by the processor, a second first-degree family network for the subject based on the pooled homoeologous estimates; phasing, by the processor, variants as hLoF based on the combined homoeologous estimate and the second first-degree family network; 4. A computer-implemented method comprising:
7. The first homoeologous estimate and the second homoeologous estimate are: The computer-implemented method of claim 6 , wherein the samples are generated by the processor by establishing, for each of one or more of the samples, a designation of an ancestor superclass.
8. Phasing the variant as hLoF further comprises: identifying, by the processor, variants as compound heterozygous mutations (CHMs) based on the combined homoeologous estimates and the second first-degree family network; phasing variants identified by said processor as CHM as putative, missense, or non-putative; and phasing variants identified by said processor as putative or missense as hLoF; 8. The computer-implemented method of claim 6 or claim 7, further comprising:
9. Phasing the identified variant as hLoF further comprises: classifying the variant as hLoF by said processor if the variant is a frameshift, stop codon gain, stop codon loss, start codon gain, start codon loss, splicing acceptor alteration, or splicing donor alteration variant; The computer-implemented method of claim 8 , further comprising:
10. Phasing the variant as hLoF further comprises: classifying compound heterozygous mutations (CHMs) in which a missense mutation occurs as disruptive or non-disruptive by the processor; and classifying, by said processor, the CHM in which the perturbing missense mutation occurred as hLoF; The computer-implemented method of claim 6 , further comprising:
11. Phasing a variant as being CHM or not CHM involves: phasing, by the processor, variants according to population allele frequency; removing, by said processor, variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both; removing single nucleotide polymorphisms (SNPs) with a quality by depth (QD) of about 2 or less, or a read depth (DP) of less than about 5, or an alternative allele balance (AB) of about 10% or less, or a combination thereof, by the processor; removing insertions or deletions (INDELs) with QD of about 2 or less, or DP of less than about 5, or AB of about 10% or less, or a combination thereof, by said processor; selecting, by said processor, the remaining variants as potential compound heterozygous mutations (pCHMs), where one or more pairs of variants exist within the same sample and within the same gene; phasing, by said processor, pCHM as cis pCHM or trans pCHM, and then classifying said pCHM phased as trans pCHM as CHM; The computer-implemented method of claim 10, comprising:
12. 1. A non-transitory computer-implemented method for fading hLoF, the non-transitory computer-implemented method comprising: generating a first homozygosity estimate for the subject within an ancestral superclass; generating a second homozygosity estimate for the subject independent of the subject's ancestral superclass; clustering each of the subjects into a first-degree family network based on the first homozygosity estimate and the second homozygosity estimate; generating a third homoeologous estimate within the first-degree family network; combining the first homoeologous estimate with the third homoeologous estimate; constructing a second first-degree family network of the samples based on the pooled homoeologous estimates; phasing the identified variants as hLoF based on the combined homoeologous estimate and the second first-degree family network; A non-transitory computer-implemented method comprising instructions executable by a processor.
13. 1. A computer-implemented method for fading hLoF, the computer-implemented method comprising: establishing, by the processor, a designation of an ancestor superclass for each of one or more of the samples; generating, by the processor, a first homoeologous estimate for the subject within an ancestral superclass; generating, by the processor, a second homoeologous estimate for the subject independent of the subject's ancestral superclass; clustering, by the processor, subjects into a first first-degree family network based on one or more of the second homozygosity estimates; generating, by the processor, a third homozygosity estimate for the subject within the first first-degree family network; combining, by the processor, the first homoeologous estimate and the third homoeologous estimate to obtain a combined homoeologous estimate; constructing, by the processor, a second first-degree family network for the subject based on the pooled homoeologous estimates; phasing, by the processor, variants as being either compound heterozygous mutations (CHM) or not CHM according to the combined homoeologous estimate and the second first-degree family network; phasing variants identified by said processor as CHM as putative, missense, or non-putative; and phasing variants identified by said processor as CHM and as putative or missense as hLoF; 4. A computer-implemented method comprising:
14. Phasing a putative variant as hLoF: classifying a variant identified by said processor as CHM as hLoF if said variant is a frameshift, stop codon gain, stop codon loss, start codon gain, start codon loss, splicing acceptor alteration, or splicing donor alteration variant; The computer-implemented method of claim 13, further comprising:
15. Phasing a missense variant as hLoF: classifying, by the processor, the CHM in which the missense mutation occurs as perturbed or non-perturbed; and classifying, by said processor, the CHM in which the perturbing missense mutation occurred as hLoF; 14. The computer-implemented method of claim 13, comprising:
16. Phasing a variant as being CHM or not CHM involves: phasing, by the processor, variants according to population allele frequency; removing, by said processor, variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both; removing, by the processor, SNPs with a quality by depth (QD) of about 2 or less, or a read depth (DP) of less than about 5, or an alternative allele balance (AB) of about 10% or less, or any combination thereof; removing insertions or deletions (INDELs) with QD of about 2 or less, or DP of less than about 5, or AB of about 10% or less, or a combination thereof, by said processor; selecting, by said processor, the remaining variants as potential compound heterozygous mutations (pCHMs) where one or more pairs of variants exist within the same sample and within the same gene, phasing the pCHMs as cis pCHMs or trans pCHMs, and then classifying the pCHMs phased as trans pCHMs as CHMs; 14. The computer-implemented method of claim 13, comprising:
17. Phasing a variant as being CHM or not CHM involves: removing, by said processor, variants that are out of Hardy-Weinberg equilibrium (HWE) or within 10 base pairs of another variant in the same sample, or both; removing, by the processor, SNPs with a quality by depth (QD) of about 3 or less, or a read depth (DP) of less than about 7, or an alternative allele balance (AB) of about 15% or less, or any combination thereof; and removing insertions or deletions (INDELs) with QD of about 5 or less, or DP of less than about 10, or AB of about 20% or less, or any combination thereof, with said processor; 20. The computer-implemented method of claim 16, comprising:
Citation Information
Patent Citations
Systems and methods for genomic annotation and distributed variant interpretation
WO2013067001A1