A unified test of linkage analysis and rare-variant association for analysis of pedigree sequence data

pVAAST

pVAAST searches through the personal genomic data from disease pedigrees, sporadic cases and unaffected controls to identify genes associated with disease. To do so, it combines logarithm of odds (lod) scores with association signals to generate a unified test statistic that offers a higher power compared to either method alone. Unlike lod scores in traditional parametric linkage analysis, the lod score in pVAAST is designed for sequence data. Specifically, the statistical model assumes that the dysfunctional variants influencing disease-susceptibility can be directly detected. As a result, the pVAAST lod score is in general more powerful than traditional linkage analysis with sequencing data, as we show below. Moreover, this assumption allows us to calculate lod scores for de novo mutations, which is not possible with traditional linkage analysis, given that de novo mutations are not in linkage with other markers. pVAAST is built upon the CLRT used in VAAST, but in addition integrates the linkage information (quantified by a lod score) as a separate log likelihood ratio in the pVAAST CLRT (CLRT_p) (Fig. 1). pVAAST evaluates the significance of the CLRT_P score using a combination of a randomization test and a gene-drop simulation¹⁴ (Online Methods).

Simulated family data

We first evaluated the performance of pVAAST to identify variants causing rare Mendelian diseases using simulated family data and unaffected control genomes (we recorded the parameterization of all pVAAST experiments in this manuscript in Supplementary Note 1). We investigated three disease models using both association- and pedigree-based approaches: dominant, recessive and dominant resulting from de novo mutations. In all models, we compared pVAAST with two rare-variant association tests, VAAST¹² and SKAT-O³ (version 0.91; using the ‘linear.weighted’ kernel and ‘optimal.adj’ method). For comparison, we also included a nonparametric linkage method based on an idealized scenario with perfect knowledge of identity-by-descent (IBD) states in all families and a two-point parametric linkage analysis using Superlink¹⁵ for dominant and recessive models (Supplementary Note 2). For the de novo model, we included a Poisson-based test, which detects excess inheritance error in cases (Supplementary Note 2). pVAAST correctly controls for type I error in all three scenarios (Supplementary Fig. 1).

We used each method to analyze the required sample size at four different levels of population-attributable risk (PAR)¹⁶ (Fig. 2a–c). Under all disease models, pVAAST was consistently the most powerful approach. The required sample size of pVAAST was usually an order of magnitude lower than for nonparametric linkage analysis, demonstrating the value of case-control sequencing data in the identification of genes associated with rare Mendelian diseases. Under dominant and de novo models, pVAAST typically required half the sample size of VAAST, and one-fifth the sample size of SKAT-O. Under the de novo model, the Poisson-based test was more powerful than rare-variant association tests alone (VAAST and SKAT-O), but substantially less powerful than pVAAST. In general, parametric linkage analysis performed worse than nonparametric (Fig. 2a–b and Supplementary Table 1), which is expected given that our nonparametric test was based on perfect knowledge of IBD states.

We also benchmarked the performance of pVAAST in common, complex diseases by simulating four-generation families (Fig. 2d). We compared the relative performance of four different choices of sequenced pedigree members: affected parent-offspring pairs, affected first-cousin pairs, affected second-cousin pairs and the entire pedigree. We simulated mildly deleterious risk alleles with a selection coefficient of 0.001, which resulted in an average MAF of 1.9 × 10⁻³ (Fig. 2e). With all pedigree members shown in Figure 2e, pVAAST required only 66% of the sample size of VAAST, and with affected first- or second-cousin pairs, pVAAST required 79% the sample size of VAAST (Fig. 2e). We observed no performance improvement with affected parent-offspring pairs in pVAAST compared to VAAST. With a selection coefficient of 0.01 (average MAF = 2.2 × 10⁻⁴) (Fig. 2f), we observed a similar trend but with slightly better pVAAST performance in all scenarios. pVAAST correctly controlled for type I error in all scenarios (Supplementary Fig. 2). For both the rare and common disease simulations, we also compared the performance of pVAAST to ASKAT⁹ (version 1.2d, build 2013-09-05), an extension of SKAT that accommodates family-based studies. However, ASKAT controls for familial relationships through asymptotic assumptions, and for the relatively small sample sizes that we evaluated, the type I error of ASKAT was inflated (Supplementary Fig. 3a–f).

De novo inheritance in an enteropathy pedigree

We performed whole-genome sequencing on a family quartet and used pVAAST to identify the potential causal mutation for a child with undiagnosed enteropathy (Fig. 3a). The proband was a 12-year-old male with severe diarrhea, total villous atrophy and hypothyroidism. Both parents and the sibling of the proband were unaffected. The phenotype was most consistent with the IPEX syndrome (OMIM 304790), but clinical sequencing of the FOXP3 and IL2RA genes revealed no pathogenic mutations.

We analyzed this pedigree using both the dominant and recessive models in pVAAST. Under the dominant model, the highest-ranking gene, STAT1, had a P value of 3.97 × 10⁻⁶. The only variant in this gene is a de novo mutation in the affected child, with a lod score of 0.70 and a CLRT_p score of 11.724. The second ranking gene was PAX3 (P = 3.33 × 10⁻³; lod score = 0 and CLRT_p score = 11.047). STAT1 was the only gene in the genome with a lod score >0.1; genes with lod scores between 0.1 and 0 fit an inheritance pattern of dominance with incomplete penetrance. Under the recessive model, no gene has a P value <1.18 × 10⁻³ (Fig. 3b). We validated the de novo inheritance pattern by genotyping the offspring and parental genotypes with Sanger sequencing. Other than this mutation, we did not identify any exonic variation in STAT1 in the family. This heterozygous mutation is observed only in the proband but not in the parents or unaffected sibling.

The de novo mutation found in the affected child is a single-nucleotide guanine-to-adenine mutation, causing the amino acid change T385M in the DNA-binding motif of STAT1; the reference allele–encoded threonine is conserved among almost all sequenced vertebrate genomes¹⁷. STAT1 encodes a transcription factor belonging to the signal transducers and activator of transcription family; both gain- and loss-of-function mutations in STAT1 cause human disease¹⁸. Gain-of-function mutations in STAT1 cause autosomal dominant chronic mucocutaneous candidiasis (CMC)^19–21 and an IPEX-like phenotype²². The T385M mutation was reported as a cause of CMC in a Japanese patient²³ and a Ukrainian patient²⁴. These data support T385M as the causative mutation for this patient’s phenotype, and demonstrate pVAAST s ability to identify a causal de novo mutation from a family quartet with a single affected proband.

Dominant inheritance in a cardiac septal defect pedigree

We analyzed whole-genome sequencing data from a previous study²⁵ on a single pedigree affected with cardiac septal defects and having an autosomal dominant pattern of inheritance (Fig. 4a). Previously²⁵, the G296S mutation in GATA-binding protein 4 (encoded by GATA4) was identified as the cause of cardiac septal defects in this pedigree using genome-wide linkage mapping followed by sequencing of the GATA4 coding region and functional studies. pVAAST successfully identified GATA4 with genome-wide significance (P = 2.0 × 10⁻⁹; Fig. 4b). The mutation encoding G296S had a CLRT_p score of 38.4 (CLRT_V score = 13.2; lod score = 5.47), and no other variants received a positive CLRT_V or lod score in GATA4. The second-ranking gene was ITIH2, with a P value of 2.3 × 10⁻⁵ and a lod score of 1.51. Because the prevalence parameter (disease prevalence in general population) was set to 0.01 to match that of cardiac septal defects, no other gene received a positive lod score in the pedigree. ASKAT was not applicable to this example owing to the small sample size (Supplementary Fig. 3g). When VAAST analyzed the genomic sequence of a single affected individual in the cardiac septal defect pedigree (the affected individual in the second generation), GATA4 was ranked forty-first genome-wide, with a P value of 2.0 × 10⁻³ (Supplementary Fig. 4).

We also analyzed the cardiac septal defect pedigree using a two-point parametric linkage test implemented in Superlink¹⁵. The mutation encoding G296S in GATA4 has a lod score of 5.13 and was the highest-scoring variant genome-wide. Assuming 2ln(10^lod) is χ² distributed with two degrees of freedom (penetrance and recombination frequency), the P value of the mutation encoding G296S from two-point linkage analysis was 7.32 × 10⁻⁶.

Recessive inheritance in a Miller’s syndrome pedigree

We investigated the performance of pVAAST on a recessive disease, Miller’s syndrome, using previously generated⁷ whole-genome sequencing data from a two-generation pedigree (Fig. 5a). The two offspring are affected with Miller’s syndrome and primary ciliary dyskinesia, both of which are rare recessive Mendelian diseases. The two diseases are caused by compound heterozygous mutations in the DHODH and DNAH5 genes, respectively⁷. All four individuals in the family quartet were sequenced. pVAAST identified only five genes with positive lod scores, and the two disease-causal genes (DHODH and DNAH5) were ranked first and second genome-wide (Fig. 5b), with P values of 3.3 × 10⁻⁵ and 1.3 × 10⁻⁴, and CLRT_V scores of 27.9 and 30.8, respectively. The lod scores were 1.204 in both cases. In both genes, only the two causal mutations received positive scores; all other variants had scores of 0.

We also explored the performance of pVAAST after removing one affected child (B01) from the pedigree. That is, we converted the original Miller’s syndrome pedigree to a trio family with two unaffected parents and one affected child. In this scenario, DHODH and DNAH5 were ranked first and thirteenth genome-wide, respectively (Supplementary Fig. 5a), both with lod scores of 0.602. We also ran VAAST over the genome-sequencing data of only one affected child (i.e., not using the data from the parents and the affected sibling). DHODH and DNAH5 were ranked tenth and twenty-seventh, respectively (Supplementary Fig. 5b). In our previous work, by enforcing a strict filtering method based on inheritance patterns and minor allele frequencies, VAAST was also able to identify the correct causal genes in this pedigree but was unable to produce an accurate P value that accounted for the familial relationships¹².

Challenging situations in pedigree studies

In linkage analysis, factors such as incomplete penetrance, locus heterogeneity and missing phenotypes negatively affect linkage signals and thus reduce disease-gene identification power. The cardiac septal defect pedigree data presented above (Fig. 4) is a large pedigree with no locus heterogeneity and very high penetrance (93.3%) for the mutation encoding G296S. We modified the genotype and phenotype data from this pedigree (Supplementary Note 3) to benchmark pVAAST in four scenarios: (i) missing phenotypes, (ii) reduced penetrance, (iii) locus heterogeneity and (iv) reduced number of informative meioses in the family. For each test case, we evaluated the lod score and the genome-wide ranking of GATA4 (ranked by P values). The lod score reported by pVAAST was approximately a monotonic function of each of the four parameters and was highly correlated with the classic two-point parametric lod score (Fig. 6). pVAAST was robust to pedigrees with missing phenotype data. For example, when 82% of pedigree members had unknown phenotypes, the lod score of GATA4 was 1.5 and genome-wide ranking was first (Fig. 6a,b). Reduced penetrance generally decreased the lod score without significantly compromising the genome-wide ranking (Fig. 6c,d). Specifically, the genome-wide ranking of GATA4 was consistently first until the penetrance dropped below 40%; even with penetrance of 20%, GATA4 was ranked eighth genome-wide. In comparison, locus heterogeneity had a greater impact on power (Fig. 6e,f). When locus heterogeneity was modest, GATA4 always ranked first or second. However, when the proportion of affected individuals carrying G296S fell to 50%, the lod score dropped below 0.2, and the genome-wide ranking was beyond fiftieth. The original family has 20 informative meiosis events, and our results show pVAAST ranked GATA4 first genome-wide even when there are only 11 informative meioses in the family (Fig. 6g). Furthermore, with only six meioses, pVAAST still ranked GATA4 second genome-wide. This suggests that for a rare Mendelian disease with high penetrance and low locus heterogeneity within the family, the risk gene can often be identified among the top hits genome-wide using a typical three-generation pedigree.

For comparison, we evaluated the genome-wide ranking of GATA4 with three alternative approaches. In the first approach, we calculated a two-point parametric lod score at each polymorphism site with Superlink¹⁵ and designated the lod score from the best-scoring site overlapping a protein-encoding gene as the gene lod score. We then ranked all genes by the gene lod scores. We also attempted to perform multipoint linkage analysis with Merlin²⁶, but this proved computationally infeasible. In the second approach, we applied the same procedure to the pVAAST lod score to calculate the ranking (Fig. 6). Finally, we evaluated a hard-filtering approach that only considered variants that perfectly fit the expected inheritance pattern with minor allele frequencies below 0.5% (Supplementary Note 3).

We found that the pVAAST lod score was consistently more robust than the classic two-point parametric lod score in challenging scenarios such as low penetrance, high locus heterogeneity, small sample size and large fraction of unknown phenotypes. The ranking of GATA4 with pVAAST lod scores was usually one order of magnitude higher than with Superlink. This performance difference is perhaps not surprising given that traditional linkage analysis tests the hypothesis of disease linkage rather than disease causation and was developed for sparse marker data rather than complete sequence data. Ranking using pVAAST P values instead of lod scores further improved the accuracy of disease-gene identification, and the improvement was pronounced when the penetrance was low or the phenotypes were missing for a large fraction of the pedigree. Hard-filtering makes strict assumptions about the expected inheritance pattern and minor allele frequency of the causal mutation. When these assumptions hold, hard-filtering has comparable performance to traditional two-point linkage analysis but is less robust compared to pVAAST and pVAAST lod scores. However, hard filtering performed very poorly when any of these assumptions were violated (Fig. 6e).

We also investigated the impact of incomplete penetrance, locus heterogeneity and unknown phenotypes in conjunction with smaller family sizes. To do so, we used only a subset of the individuals in the original cardiac septal defect pedigree to reduce the number of informative meiosis. We evaluated the genome-wide ranking of GATA4 using pVAAST, pVAAST lod scores, two-point linkage analysis in Superlink, multipoint linkage analysis in Merlin²⁶ and hard filtering (Supplementary Figs. 6–8). The ranking of GATA4 was highest when using pVAAST in almost all scenarios, which is consistent with the results involving the entire family (Fig. 6).

Basic lod score calculation in pVAAST

In classic two-point parametric linkage analysis, the marker under investigation is usually assumed not to be causal but rather linked with the actual causal variant with a certain recombination probability (r). Under the null hypothesis, r = 0.5, which indicates that the marker and causal mutation are unlinked. Under the alternative hypothesis, r is a free parameter. Given the disease prevalence, allele frequency of the marker and causal allele and the penetrance of the causal allele, the likelihood of alternative and null model can be calculated for given values of r using the Elston-Stewart algorithm³⁷. The log₁₀ ratio of the maximum likelihood of the alternative and null model is the lod score.

For simplicity, we use the term causal to refer to any variant that directly increases disease risk, regardless of penetrance. Our model assumes that the disease is caused by either the locus under investigation (current locus) or some other unlinked locus in the genome (latent locus). In both models, the current and latent loci are unlinked, and there is no epistatic interaction between the alleles. The null model states that variant(s) in the latent locus cause the disease with some probability, and the current locus is not causal. The alternative model states that variants in both the current and latent loci can independently cause the disease, with different probabilities. In other words, the null model attributes the disease phenotype solely to the latent locus, and the alternative model allows variants in both the current and latent loci to be independently causal. We then maximize the likelihoods of the alternative and null models over ρ_c (genotype disease probability vector for the current variant), ρ_l (genotype disease probability vector for the latent locus) and f_l (minor allele frequency of the latent locus) and calculate the log10 likelihood ratio as the lod score. Formally,

and the likelihood for both null model and alternative model has the form

Here g_c and g_l are the genotype vectors (with values of 0, 1 and 2 corresponding to homozygous-reference, heterozygous and homozygous-nonreference genotypes) of the current and latent variant sites; p is the phenotype vector of the pedigree; f_c and f_l are the allele frequencies of the current and latent alleles. Under the null model, the expression can be further decomposed into

because only the latent allele is causal for the disease under the null model, and g_c is thus independent from p, g_l and ρ_l.

Given ρ_c, ρ_l, f_c and f_l, all of the aforementioned probabilities can be calculated with the Elston-Stewart algorithm¹⁵ in linear computational time relative to the family size. We estimate f_c from the allele frequency in a control population and perform a grid search over ρ_c, ρ_l and f_l in the specified order to maximize the likelihood. By default, we explore ρ_c and ρ_l values ranging from 0 to 1 with increment of 0.1, and in addition the following values: 0.001, 0.01 and 0.999. We explored the following f_l values: 5 × 10⁻⁷,5 × 10⁻⁴,5 × 10⁻³, 0.01, 0.02, 0.05, 0.5 and 0.999. The resolution of the grids is tunable. In all our experiments, the aforementioned parameters offer a good balance between algorithm efficiency and statistical power, although using a finer grid may increase the power of pVAAST at the cost of longer computation time. For the dominant model, a heterozygous genotype is sufficient to be considered as a risk genotype; for a simple recessive model, a homozygous genotype is required, with the exception of sex chromosomes. Compound heterozygous scenarios are discussed below.

If more than one family is present, for each variant, we maximize the likelihood under the assumptions that ρ_c is consistent across families but ρ_l and f_l varies between families using a nested grid search. Then, within each family, the lod of one variant is chosen to be the gene lod score of this family. By default, the variant with the highest CLRT_V score is chosen, but the user can opt to use CLRT_p score or lod score alone as well. In practice, we found that in large pedigrees, using the CLRT_p score as a selection criterion may yield more favorable results. Finally, we sum the gene lod score from multiple families to generate the overall pVAAST lod score.

Extending the dominant model to de novo mutations

We accommodate de novo mutations in our model by allowing Mendelian inheritance errors to occur in the pedigree likelihood calculation. Specifically, in the Elston-Stewart algorithm³⁷, if the offspring carries a mutation absent from both parents, then this transmission has a probability of m (mutation rate per site per generation in human genome; default 1.2 × 10⁻⁸ (ref. 7)). Accordingly, we also randomly introduce Mendelian inheritance error in our gene-drop simulations¹⁴ with probability equal to the genotyping error rate.

Extending the recessive model to compound heterozygotes

Compound heterozygotes require special attention because the genotype vectors (g_l and g_c) now involve more than one variant site. Under the recessive model we are specifically interested in the situation where two deleterious mutations occur at two different chromosomes of the same gene, so that both copies are defective. To illustrate, consider a gene with three polymorphism sites, i, j and k. A straightforward approach to calculate the gene lod score would be to calculate the lod for all pairwise combination of heterozygous variant sites within the gene (i.e., i + j; i + k; and j + k) separately and then select the highest lod score. This requires the evaluation of n(n − 1)/2 combinations, where n is the number of variant sites in the gene. However, this approach is flawed because it assumes the genotype disease probabilities for all pairs of sites are independent, which is incorrect. Instead, we assume that any variant in the gene is either causal (D-variants) or neutral (N-variants)³⁸ with the same relative risk. For example, if a gene has four heterozygous sites i, j, k and l, within which i, j, and k are causal, then an individual with at least two mutations at i, j and k sites on two different chromosomes would be at risk; otherwise she will not be at risk.

Under this model, we can construct a Boolean risk vector for a gene to denote whether each variant within the gene is a D-variant or N-variant. If we know the underlying risk vector for some gene, then we can easily determine the genotype of an individual by evaluating whether he or she carries at least one D-variant on each chromosome. Then the calculation of lod score is reduced to the simple recessive case described above. However, finding the optimal risk vector is not trivial, as a brute-force approach to find the risk vector maximizing the lod score has a complexity of 0(2ⁿ), where n is the number of sites in the gene. To make this more efficient, we use an MCMC method³⁹ to approximate the optimal risk vector. Briefly, given a particular risk vector and the phenotype probability for each genotype, the joint likelihood for all sequenced and phenotyped individuals can be calculated as

where ρ_r is the probability that an individual with a risky genotype is affected; ρ_n is the probability that an individual with a neutral genotype is affected; n_a and n_b are the total numbers of affected and unaffected individuals with a risky genotype, respectively; n_c and n_d are the total number of affected/unaffected individuals with neutral genotypes, respectively. Both ρ_r and ρ_n are configurable parameters, although we found that the performance of our MCMC method was usually insensitive to these parameters.

We start with a random risk vector, and randomly select a variant site to switch to the opposite value (neutral to risky and risky to neutral). The likelihoods for both risk vectors are calculated, and we selectively accept the new risk vector according to the Metropolis-Hastings method³⁹. This process is repeated until convergence or the maximal number of iterations is achieved. Lastly, we select the most likely risk vector from the Markov chain and calculate the lod score as described in the previous section.

Optionally, the joint likelihood can incorporate an empirical functional score. Let I_D(k) be an indicator function for whether the kth site is a D-variant. The empirical functional score (F score) is a function of VAAST CLRT_V scores across all sites in the current gene

and the updated likelihood is calculated as L* = L × e^F.

The calculation of CLRT_V score is detailed in (ref. 12). Briefly, it is twice the log-scale composite likelihood ratio of disease model versus null model, incorporating the mutation frequency in the control genome and the functional impact of the mutation on the protein sequence. This option (mcmc_use_functional_score) can be switched on or off. We used the updated likelihood function throughout the present study, although in our recessive model simulations, these two likelihood functions generated similar results.

Integrating lod scores into the CLR test

pVAAST is built on the framework of VAAST¹², which uses an extended CLRT to determine a severity score for genomic variants. The null model of the CLRT states that the frequency of a variant or variant group is the same in the control population (background genomes) and the case population (target genomes), whereas the alternative model allows these two frequencies to differ. Under a binomial distribution, the likelihood for both models can be calculated on the basis of observed allele frequencies in the control and case data sets. This likelihood is further updated by calibrated amino acid substitution and insertion and deletion (indel) severity weights.

To integrate genetic linkage information into the CLRT, we select only one sequenced and affected individual from each pedigree (pedigree representative) to establish a group of cases. The identifiers of the selected individuals can be provided, but if such information is absent, pVAAST will randomly choose one individual carrying the highest-scoring variant in the current gene. Additional affected individuals not related to any other individuals in the study can also be included among the cases. λ represents the natural log of the composite likelihood ratio calculated as previously described¹². We calculate the pVAAST CLRT (CLRT_p) score as

where LOD_j is the lod score for the ith family and

To avoid confusion, we denote the original CLRT score in VAAST (without the linkage component) as CLRT_V in this manuscript. Figure 1 provides a schematic diagram for the calculation of the CLRT_p scores in pVAAST.

Evaluating the significance of the test statistic

c represents the two parental haplotypes at the current gene locus in a particular individual. Let subscript p,pf,b and sc represent a vector of cs among all pedigree members, pedigree founders, background (control) individuals or sporadic cases, respectively, and a superscript r or s represent real data and simulated data, respectively. For example, c^r_pf represents the vector of haplotypes in all pedigree founders in the real data. T represents the unordered set of chromosomes among pedigree founders, background genomes, and sporadic cases in the real situation. Our null hypothesis is that pedigree founders, controls and sporadic cases are derived from the same population and that haplotypes in pedigree offspring randomly segregate according to Mendel’s law. When the two haplotypes in each sequenced individual are known and all pedigree founders are sequenced, a combination of a randomization test and gene-drop simulation can be used to evaluate any statistic that is a function of the genotype and phenotype data in the pedigree and controls.

We first sample (without replacement) N_pf (the cardinality of the set c^r_pf, i.e., |c^r_pf|) individuals from T as the pedigree founder (denoted by C_pf); N_ct (< = |c^r_b|) individuals as the control set for CLRT_V calculation (denoted by c_ct); and N_sc (|c^r_sc|, which can be 0) individuals as the sporadic cases (denoted by c_sc). We then generate the c_p from c_pf via gene-drop simulation¹⁴. Briefly, we simulate the two haplotypes of each offspring by randomly sampling one of each parent’s two haplotypes with equal probability. The gene dropping starts from the first generation of the pedigree and is repeated until all pedigree members are simulated. g (c_p, c_sc, c_ct) represents the desired test statistic. In pVAAST, this test statistic is CLRTp. The real data in this procedure are represented as c^r_p, c^r_ct and c^r_sc, where c^r_ct is a random subset of c^r_b, with size N_ct. If we calculate

within the described sampling space, we will have a valid P value with specified type I error under the null model. This holds because the real data are one realization of the described sampling scheme with probability equal to any other realization under the null hypothesis.

In reality, because enumerating all values of c_p, c_sc and c_ct is computationally intractable, we use a Monte Carlo method to sample n realizations of the described procedure and calculate

(I is an indicator function) and report this as the gene-level P value. Alternatively, P value can be calculated using the lod score instead of CLRT_p score as the test statistic.

We emphasize two points: (i) the number of sporadic cases can be 0 and (ii) the choice of N_ct is free and does not affect the validity of the P value.

To sample from T, the above procedure requires that the haplotypes of all pedigree founders are known. In reality, c^r_pf can be unknown or partially known because pedigree founders may not have been sequenced, thus we may not be able to directly sample from T. To accommodate this situation, we define a new set T* to be the unordered set of haplotypes among pedigree representatives (one affected sequencing individual in each pedigree, as denoted in the pVAAST parameter file), background genomes and sporadic cases in the real situation. Obviously we have

We propose sampling our test-statistics CLRT_p from the cumulative distribution function

during the simulation to approximate the distribution

The approximation becomes more accurate when the |T – T*| ≪ |T|, or in other words, when the number of unsequenced founders is small compared to the total number of sequenced background individuals, pedigree representatives and sporadic cases. Our implementation also approximates the idealized procedure owing to haplotype phase uncertainty. Despite these approximations, we observed no inflation in type I error rate in any of the experiments we evaluated (Supplementary Figs. 1 and 2).

We also documented the implementation of our simulation procedure in pVAAST in Supplementary Note 4.

Genomic data

For the enteropathy pedigree, whole-genome sequencing was performed on all four pedigree members using the Illumina HiSeq platform. We followed the Genome Analysis Toolkit (GATK) best practice to perform variant-calling steps⁴⁰. Briefly, we used Burrows-Wheeler Aligner to align reads⁴¹, GATK⁴⁰ to remove PCR duplicates and perform indel realignment and UnifiedGenotyper in GATK⁴⁰ to jointly call the genotypes in the sequenced pedigree members and 136 exomes from the 1000 Genomes Project⁴². The 136 exomes used as controls include individuals with western European ancestry (CEU) and British in England and Scotland (GBR). Potential disease-causing mutations were validated with Sanger sequencing at the University of Utah sequencing core. For the cardiac septal defect pedigree, Complete Genomics performed whole-genome sequencing and variant calling on selected pedigree members.

For the results presented the sections on cardiac septal defects, Miller Syndrome and challenging situations in pedigree studies, we used the control genome set consisting of 1,057 exomes from 1000 Genomes Project phase I data⁴³, 54 genomes from the Complete Genomics Diversity Panel⁴⁴, 184 Danish exomes⁴⁵ and nine nonduplicative genomes from the 10Gen data set⁴⁶, representing a wide variety of ethnicities and sequencing platforms. To include a wider set of variants that are unlikely to be causal for rare Mendelian diseases, we further collected high-quality variants (defined as polymorphism sites with sample sizes no smaller than 100 chromosomes) from dbSNP build 137 (http://www.ncbi.nlm.nih.gov/SNP/) and NHLBI exome sequencing data (http://evs.gs.washington.edu/EVS). We then randomly inserted these variants into the control genome set, setting the minor allele frequency equal to the reported value.

Secondary analysis studies were approved by the Western Institutional Review Board for the cardiac septal defects and Miller’s syndrome (dbGAP phs000244.v1.p1) pedigrees after initial studies were approved by local institutional review boards at sites interacting with participants. Procedures followed were in accordance with institutional and national ethical standards of human experimentation. Proper informed consent was obtained.

pVAAST runtime

pVAAST supports multithreading parallelization. The computational time is proportional to the size of pedigree and to the rounds of randomization tests being performed. On our Linux server with Intel Xeon 2.00 GHz CPUs, the enteropathy pedigree took 0.6 h (clock time) using 40 threads (1 × 10⁸ maximum randomizations). The cardiac septal defect pedigree took 181 h (clock time) using 40 threads (maximum randomizations: 1 × 10⁹). The Miller’s syndrome pedigree took 0.3 h (clock time) using 70 threads (maximum randomizations: 1 × 10⁶).

Software access

pVAAST is available for download at http://wwwyandell-lab.org/software/vaast.html with an academic user license. The source code for pVAAST is included as Supplementary Software.