Inferring compound heterozygosity from large-scale exome sequencing data

Abstract

Recessive diseases arise when both the maternal and the paternal copies of a gene are impacted by a damaging genetic variant in the affected individual. When a patient carries two different potentially causal variants in a gene for a given disorder, accurate diagnosis requires determining that these two variants occur on different copies of the chromosome (i.e., are in trans) rather than on the same copy (i.e. in cis). However, current approaches for determining phase, beyond parental testing, are limited in clinical settings. We developed a strategy for inferring phase for rare variant pairs within genes, leveraging genotypes observed in exome sequencing data from the Genome Aggregation Database (gnomAD v2, n=125,748). When applied to trio data where phase can be determined by transmission, our approach estimates phase with 95.7% accuracy and remains accurate even for very rare variants (allele frequency < 1×10^-4). We also correctly phase 95.9% of variant pairs in a set of 293 patients with Mendelian conditions carrying presumed causal compound heterozygous variants. We provide a public resource of phasing estimates from gnomAD, including phasing estimates for coding variants across the genome and counts per gene of rare variants in trans, that can aid interpretation of rare co-occurring variants in the context of recessive disease.

Fig. 1:. Overview of phasing approach using Expectation-Maximization method in gnomAD.

a, Schematic of phasing approach. b, Histogram of ${\text{P}}_{\text{trans}}$ scores for variant pairs in cis (top, blue) and in trans (bottom, red). c, Proportion of variant pairs in each ${\text{P}}_{\text{trans}}$ bin that are in trans. Each point represents variant pairs with ${\text{P}}_{\text{trans}}$ bin size of 0.01. Blue dashed line at 10% indicates the ${\text{P}}_{\text{trans}}$ threshold at which ≥ 90% of variant pairs in bin are on the same haplotype $\left({\text{P}}_{\text{trans}}\le 0.02\right)$. Red dashed line at 90% indicates the ${\text{P}}_{\text{trans}}$ threshold at which ≥ 90% of variant pairs in bin are on opposite haplotypes $\left({\text{P}}_{\text{trans}}\ge 0.55\right)$. Calculations are performed using variant pairs with population AF ≥ 1×10⁻⁴. d, Performance of ${\text{P}}_{\text{trans}}$ for distinguishing variant pairs in cis and trans. Accuracy is calculated as the proportion of variant pairs correctly phased (green bars) divided by the proportion of variant pairs phased using ${\text{P}}_{\text{trans}}$ (orange plus green bars). b-d, ${\text{P}}_{\text{trans}}$ scores are population-specific.

Fig. 2:. Phasing accuracy as a function of variant allele frequency (AF).

Phasing accuracy at different AF bins for all variant pairs (a), variant pairs in trans (b), and variant pairs in cis (c). Shading of squares and numbers in each square represent the phasing accuracy. Y-axis labels refer to the more frequent variant in each variant pair and X-axis labels refer to the rarer variant in each variant pair. Accuracy is the proportion of correct classifications (i.e., correct classifications / all classifications) and is calculated for all unique variant pairs seen in the trio data across all populations using population-specific ${\text{P}}_{\text{trans}}$ calculations.

Fig. 3:. Phasing accuracy using population-specific versus cosmopolitan Ptrans estimates.

Population-specific ${\text{P}}_{\text{trans}}$ estimates are shown in light blue and cosmopolitan ${\text{P}}_{\text{trans}}$ estimates are shown in medium blue. Accuracies are shown separately for variants in trans (a, left) and variants in cis (b, right)

Fig. 4:. Phasing accuracy as a function of distance between variant pairs.

a, Phasing accuracy (y-axis) as a function of physical distance (in base pairs on log10 scale) between variants (x-axis). Blue represents variants on the same haplotype (in cis), and red represents variants on opposite haplotypes (in trans). b, Same as a, except the x-axis shows genetic distance (in centiMorgans). Accuracies for a and b are based on unique variant pairs observed across all genetic ancestry groups and are calculated using population-specific ${\text{P}}_{\text{trans}}$ estimates.

Fig. 5:. Counts of genes with variants in trans in gnomAD.

a, Proportion of genes with one or more individuals in gnomAD carrying predicted compound heterozygous (in trans) variants or a homozygous variant at ≤ 1% and ≤ 5% AF stratified by predicted functional consequence. b, Number of genes with ≥ 1 individual in gnomAD carrying compound heterozygous (in trans) or homozygous predicted damaging variants at ≤ 1% AF, stratified by predicted functional consequence and Mendelian disease-association in the Online Mendelian Inheritance in Man database (OMIM). In total, 28 genes (25 non-disease, 2 AD, and 1 AR) carried predicted compound heterozygous loss-of-function variants at ≤ 1% AF, only seven of which were high confidence “human knock-out” events following manual curation. For predicted compound heterozygous variants, both variants in the variant pair must be annotated with a consequence at least as severe as the consequence listed (i.e., a compound heterozygous loss-of-function variant would be counted under the pLoF category but also included with a less deleterious variant under the other categories). All homozygous pLoF variants previously underwent manual curation as part of Karczewski et al. AF, allele frequency; comp het, compound heterozygous; hom, homozygous; AD, autosomal dominant; AR, autosomal recessive.

Fig. 6:. Publicly-available browser for sharing of phasing data.

a, Sample gnomAD browser output for two variants (1-55505647-G-T and 1-55523855-G-A) in the gene PCSK9. On the top, a table subdivided by genetic ancestry group displays how many individuals in gnomAD from that genetic ancestry are consistent with the two variants occurring on different haplotypes (trans), and how many individuals are consistent with their occurring on the same haplotype (cis). Below that, there is a 3×3 table that contains the 9 possible combinations of genotypes for the two variants of interest. The number of individuals in gnomAD that fall in each of these combinations are shown and are colored by whether they are consistent with variants falling on different haplotypes (red) or the same haplotype (blue), or whether they are indeterminate (purple). The estimated haplotype counts for the four possible haplotypes for the two variants as calculated by the EM algorithm is displayed on the bottom right. The probability of being in trans for this particular pair of variants is > 99%. b, Variant co-occurrence tables on the gene landing page. For each gene (GBA1 shown), the top table lists the number of individuals carrying pairs of rare heterozygous variants by inferred phase, AF, and predicted functional consequence. The number of individuals with homozygous variants are tabulated in the same manner and presented as a comparison below. AF thresholds of ≤ 5%, ≤ 1%, and ≤ 0.5% are displayed across six predicted functional consequences (combinations of pLoF, various evidence strengths of predicted pathogenicity for missense variants, and synonymous variants). Both variants in the variant pair must be annotated with a consequence at least as severe as the consequence listed (i.e., pLoF + strong missense also includes pLoF + pLoF).