Identification of rare causal variants in sequence-based studies: methods and applications to VPS13B, a gene involved in Cohen syndrome and autism

Abstract

Pinpointing the small number of causal variants among the abundant naturally occurring genetic variation is a difficult challenge, but a crucial one for understanding precise molecular mechanisms of disease and follow-up functional studies. We propose and investigate two complementary statistical approaches for identification of rare causal variants in sequencing studies: a backward elimination procedure based on groupwise association tests, and a hierarchical approach that can integrate sequencing data with diverse functional and evolutionary conservation annotations for individual variants. Using simulations, we show that incorporation of multiple bioinformatic predictors of deleteriousness, such as PolyPhen-2, SIFT and GERP++ scores, can improve the power to discover truly causal variants. As proof of principle, we apply the proposed methods to VPS13B, a gene mutated in the rare neurodevelopmental disorder called Cohen syndrome, and recently reported with recessive variants in autism. We identify a small set of promising candidates for causal variants, including two loss-of-function variants and a rare, homozygous probably-damaging variant that could contribute to autism risk.

Figure 1. (a) Median rank of causal variants among the non-synonymous variants for two disease models (M1 and M2) and three values for the NS:S ratio (R = {0.6, 1.0, 1.4}).

The proportion of causal variants in the region is 20%. HM refers to the original hierarchical model with ranking of the causal variants among the non-synonymous variants, based on their estimated effects; BE refers to the backward elimination procedure for non-synonymous variants; and HM_S refers to the ranking of causal variants only among those non-synonymous variants selected by the backward elimination procedure, with ranks based on the estimated effects from the hierarchical model. (b) The number of causal variants in Top 10 for non-synonymous variants.

Figure 2. The effect of multiple bioinformatic predictors for non-synonymous variants.

Ranking is done only within the set of variants selected by the backward elimination procedure. (a) Median rank of causal variants for two disease models (M1 and M2) and three values for the NS:S ratio (R = {0.6, 1.0, 1.4}). The proportion of causal variants in the region is 20%. HM_S refers to the hierarchical model with ranking of the causal variants among the selected non-synonymous variants, based on their estimated effects, B₁ refers to the hierarchical model with one bioinformatic predictor (B₁, Table 2), B₂ refers to the hierarchical model with one bioinformatic predictor (B₂), mB₁ refers to the hierarchical model with three bioinformatic predictors (B₁, B₁, and B₂), and mB₂ refers to the hierarchical model with four bioinformatic predictors (four B₁s). (b) The number of causal variants in Top 10 for non-synonymous variants.

Figure 3. Predicted deleteriousness scores are shown for 71 rare functional variants (non-synonymous, nonsense and splice-sites).

From the top, the first plot depicts the PolyPhen-2 score for each variant, the second depicts the GERP_RS score, and the third depicts variant counts for cases (up) and controls (down). Green tick marks indicate a variant contained in an exon, and red ticks indicate that a variant is selected by the backward elimination procedure. LoF variants are marked by a black asterisk; the homozygous probably damaging variant is marked by a red asterisk. The location of five protein domains (ChoreinN, TM2, TM4, DUF1162, Golgi targeting element, and ATG C) are depicted by boxes at the top of the plot (see Figure S10 for a complete view of VPS13B protein domains). Variants are plotted equidistantly on the x-axis.

Figure 4. Results from the backward elimination procedure for non-synonymous and splice site variants in VPS13B.

(a) The change in p value is shown as variants are being removed one by one (when the backward elimination procedure is run once on all non-synonymous variants). (b) Distribution of return counts for non-synonymous and splice site variants in VPS13B; overlaid is a fitted mixture with two components.