The FBAT (Laird et al., 2000) statistic Z_c is based on the covariance U_c between the offspring's traits and genotypes: (1)Uc=∑iTi(Xi-E[Xi])     (2)Variance(Uc)=∑iVariance(TiXi)     (3)Ti=Yi- μ     (4)Zc=Uc∕Variance (Uc)     Here, X_i denotes the offspring genotype in trio i at the genomic marker being tested. For a nuclear family with multiple offspring, there will be as many father-mother-offspring trios contributing to the test independently. The subscript “c” in the above formula denotes that FBAT is based on only “complete” trios in the data. X_i is defined by the genetic model (additive, dominant, recessive) under consideration. For example, for additive model, X_i counts the number of non-reference alleles observed in the offspring, and can take a value of 0, 1, or 2 for a bi-allelic genomic marker (Laird et al., 2000).

T_i is the coded trait defined as Y_i − μ, where Y_i denotes the observed trait of the offspring in trio i. Although Y can take several types of values, in this paper we focus on dichotomous traits where the observed trait Y_i is “1” for affected offspring and “0” for unaffected offspring.

μ∈[0, 1] is an offset value that can be chosen to maximize the power of the test (Laird et al., 2000). When μ = 0, T_i = Y_i, implying that only affected trios are used in the test (since Y_i is 0 for unaffected offspring). When μ > 0, affected trait T_i > 0 and unaffected trait T_i < 0, so both affected and unaffected trios are used in the test. For the analyses presented in this paper, we used μ = 0.5 in order to assign equal but opposite weights to affected and unaffected trios.

Figure 1A shows an example of an informative complete trio for autosomal chromosomes (Sebastiani et al., 2004). Figures 1B,C show examples of informative trio types with female and male offspring respectively for X chromosome. A comprehensive list of informative complete trios for autosomal chromosomes, as well as the X chromosome, is shown shown in Figure S1. The corresponding statistics X − E[X] and Variance(X) shown in Figure S1 are for the additive genetic model. Statistics for dominant and recessive models are in Table S1 (autosomal chromosomes) and Table S2 (X chromosome).

Here we describe an example to explain computation of the statistics X − E[X] and Variance(X). For the trio type shown in Figure 1A where one parent is reference homozygous, the other parent is heterozygous and the offspring is heterozygous, under additive genetic model, Xi=1E[Xi]=0* 1∕2+1* 1∕2= 1∕2Xi-E[Xi]=1∕2Variance(Xi)=E[Xi2]-(E[Xi])2Variance(Xi)=(0*1∕2+1*1∕2)-(1∕2)2Variance(Xi)=1∕4

Trio types 6 and 7 for autosomal chromosomes (Figure S1) are not informative under the dominant model (Sebastiani et al., 2004) because both the possible offspring genotypes—heterozygous and alternate allele homozygous—have equal penetrance under dominant genetic model. Similarly, trio types 1 and 2 are not informative under the recessive model.

Once X_i − E[X_i], Variance(X_i), and T_i are computed for each trio, U_c and Variance(U_c) are computed by summation over all the trios and, finally, the FBAT statistic Z_c is computed as ratio of U_c and standard deviation of U_c.

Z_c is essentially a z-score measuring deviation from the null hypothesis of no linkage and no association. When evaluating bi-allelic markers, a positive Z_c indicates that the allele being tested was over-transmitted to the affected offspring, whereas a negative Z_c indicates under-transmission to affected offspring. P-values are computed considering this as a two-sided test.

The FBAT test described above does not account for incomplete trios in the study cohort. This might lead to undetected bias in the test statistic. We have implemented CIFBAT to detect potential bias in the FBAT statistic in presence of incomplete trios and to identify significant genomic markers with higher precision. CIFBAT considers all the valid completions of each incomplete trio equally likely and computes QIs of the FBAT statistic over many iterations of randomized completions.

Figure 2A shows an example of an admissible incomplete trio type for autosomal chromosomes (Sebastiani et al., 2004). Figures 2B,C show examples of admissible incomplete trio types with female and male offspring respectively for X chromosome. A complete list of all admissible incomplete trio types for autosomal chromosomes (Sebastiani et al., 2004) as well as X chromosome is shown in Figure S2. We computed X − E[X] and Variance(X) for all valid completions of these incomplete trios under additive, dominant, and recessive models respectively. Table S3 lists these statistics for autosomal chromosomes, and Tables S4, S5 list these statistics for the X chromosome for trios with male and female offspring respectively. For non-informative completions (both homozygous parents), both X − E[X] and Variance(X) are equal to 0.

We will now describe how CIFBAT computes QIs of the FBAT statistic. In the following explanation, subscript “c” denotes a complete trio or a statistic related to complete trios, subscript “m” denotes missing (incomplete) trio or a statistic related to incomplete trios, and subscript “r” denotes a random variable.

Suppose for a genomic marker under evaluation, our data set consists of k complete trios and d incomplete trios. The total U_c and Variance(U_c) for all the complete trios are computed as: (5)Uc=∑kTk (Xk-E[Xk]) (6)Variance (Uc)=∑kTk2*Variance (Xk) Next, for each incomplete trio, CIFBAT chooses a completion assuming uniform distribution of all valid completions as listed in Tables S3–S5. Let X' denote the offspring's genotype for a randomly chosen completion. The total contribution of incomplete trios is a random variable U_mr computed as summation of contributions based on their random completions by CIFBAT: (7)Umr=∑dTd(Xd′-E[Xd′])     (8)Variance(Umr)=∑dTd2*Variance(Xd′)     For a single iteration of random completion of all incomplete trios, the corresponding total U statistic and the variance are computed as the sum of the statistics from complete and incomplete trios. (9)Utotal,r=Uc+Umr     (10)Variance(Utotal,r)=Variance(Uc)+Variance(Umr)     The corresponding Z statistic is computed as: (11)Zr=Utotal,r∕Variance(Utotal,r)   For each genomic marker, CIFBAT executes a pre-defined number (1000 by default) of iterations of computing U_mr and the corresponding Z_r, each time with a randomly selected set of completions for all incomplete trios. QIs of Z_r (α/2 and 100- α/2) and the corresponding p-values are then computed for a pre-defined confidence level (α = 0.05 by default).

In order to explore the statistical response of CIFBAT under various scenarios of missing data, we simulated family genotype data under no association and a genetic model of disease explained in detail in Data Sheet 1. The following three standard scenarios of missing data, as defined in the statistics literature, were simulated.

When data is “missing completely at random” (MCAR), the probability of an observation being missing does not depend on observed or unobserved measurements. This scenario was simulated by introducing missing genotypes for randomly chosen samples in the data. This type of missing data was use in the study of markers with no association.

Data can also be “missing at random” (MAR), which means that the missingness mechanism depends only on an observed measurement, and not on any unobserved measurements. This scenario was simulated with respect to the following observed variables separately: gender, subpopulation, and disease status.

Lastly, data can be “missing not at random” (MNAR) when the missingness mechanism depends on unobserved measurements. This can be difficult to detect and can lead to invalid inference if the missing data is ignored or incorrectly modeled. We simulated this scenario by concentrating missing data on heterozygous or homozygous samples, i.e., samples that were heterozygous or homozygous for the non-reference allele were more likely to be missing.

The simulated binary phenotype is either uniformly random or based on a bi-allelic additive genetic model closely following the work of Yang et al. (2003), where each individual starts with a low probability for disease and each causative SNP additively increases the probability. The model for log probability of disease, π_ij, is defined as: (12)log(πij)=αij+∑mβmX(Gijm)   Here, α_ij is the initial log probability of disease for offspring j in family i. α represents environmental or other hidden factors that contribute to disease. For each causative SNP m, β is a log ratio of penetrance values comparing having one or two alleles to none. X(G) is a vector of dummy variables indicating a count of alternative alleles. Disease status was probabilistically determined using π_ij; offspring within a family have correlated probabilities for disease, but are independent across families.

In each simulated scenario, 4000 pedigrees were generated and the mean population prevalence was 12.5%. An equal number of cases and controls were used, and missing data was generated at rates of 0, 1, 5, and 10% at the sample level (Table 2).

In the MAR scenario, for each level of missingness, we arbitrarily split the missing data 80/20 given a binary variable, such as gender or subpopulation, so that one population contains 80% of the missing data. The analysis was performed 1000 times for each level of missingness within each scenario, which results in 9 scenarios ^* 4 missingness levels ^* 1000 runs. The 9 scenarios are listed in Table 2. The parameters used and their range of values are listed in Table 3.

A false discovery rate (FDR) cut-off of 10% was used to identify significant results. The equivalent p-value threshold (from the FDR) was also used to identify significant CIFBAT QIs. An interval was considered significant if the p-value spread was below the p-value threshold.

In each simulation, the causative SNPs are known. This allows computation of performance metrics useful for comparing CIFBAT and FBAT. Among the metrics computed are sensitivity (recall), the proportion of true positives called, specificity, the proportion of true negatives called, precision, proportion of predicted trues that are actually true, and the F-measure which is the harmonic mean of recall and precision.

As more and more association studies employ whole genome sequences, it becomes important that novel association methods are evaluated for applicability and performance in genome wide association analyses. With this purpose, we used CIFBAT to analyze whole genome sequences of 784 nuclear families for uterine anomalies in mothers, represented as a binary phenotype. Samples of peripheral blood were collected from fathers, mothers and newborns at the Inova Fairfax Medical Center in Falls Church, Virginia, and sequencing was done at >40X depth using Complete Genomics' whole genome sequencing platform. Fifty-two (52) out of the 784 mothers were diagnosed with various uterine anomalies including endometriosis, bicornuate uterus, didelphic uterus, etc. (complete list in Text S1). These were used as cases in our study. Various filters were applied to the genomic markers to reduce noise in the data and to ensure applicability of family based association tests (Text S1). Altogether, we tested 3,808,482 autosomal markers under both FBAT as well as CIFBAT. A false discovery rate (FDR) cut-off of 10%, which translated to a p-value cut-off of 5.71e-05, was used to identify significant hits from FBAT. The p-value cut-off from FBAT results was also used to identify significant 95% QIs computed with CIFBAT. A QI was considered significant if the entire spread was below the p-value cut-off. A significant FBAT p-value together with a significant CIFBAT quantile interval indicated robustness against missing data. On the other hand, a significant FBAT statistic and a non-significant CIFBAT quantile interval hinted at potential bias in the FBAT statistic due to missing data.

To demonstrate the use of CIFBAT in refining candidate sets of markers, we analyzed Type 1 diabetes data that consisted of a dichotomous phenotype variable and 351 markers within 22 candidate genes that have been previously identified as contributing to the risk of Type 1 diabetes. The data were collected by the Type 1 Diabetes Genetics Consortium (T1DGC). There were 2313 pedigrees consisting of 2345 nuclear families with one or more offspring. To ensure independence of nuclear families, we excluded one randomly selected family from every pair of families that had at least one common parent. We also excluded offspring whose type 1 diabetes affectation status was unknown as well as founders that had no sequenced offspring in the data. Finally, we were left with 2314 nuclear families with one or more offspring which we then analyzed using CIFBAT. An FDR cut-off of 10% was used to indicate significant results from FBAT. The equivalent p-value cut-off of 1.10e-02 was used to indicate significant results based on CIFBAT QIs. Again, a significant FBAT statistic with a significant CIFBAT quantile interval implied robustness against missing data, whereas a significant FBAT statistic with a non-significant CIFBAT quantile interval implied potential bias due to missing data.

In addition to developing the CIFBAT method, we have extended capabilities of the original TDT and robustTDT implementations as follows. PLINK (Purcell et al., 2007) is widely used for applying Transmission Disequilibrium Test (TDT) to family based studies. PLINK's implementation of TDT handles both autosomal chromosomes as well as the X chromosome, but only under an additive model. We extended our TDT implementation to handle dominant and recessive models in addition to the additive model. Table S6 lists the informative trios for autosomal chromosomes and transmission counts for the two alleles b and c under all three models. Tables S7, S8 list the informative trios with male and female offspring respectively for the X chromosome. Similar extensions were added to robustTDT (Sebastiani et al., 2004) to handle missing genotypes for the X chromosome in addition to autosomal chromosomes under all the three genetic models (Tables S9–S11).

Analysis of simulated family genotype data allowed us to compare the performance of FBAT and CIFBAT under no phenotype association and the various patterns of missing data (MCAR, MAR, MNAR). The simulated scenarios of missing data are not meant to be comprehensive, but rather illustrative of the effects on algorithms that use missing data in making inferences. Comparisons were primarily based on recall and precision.

When phenotypes were uniformly random, and data was missing completely at random, the FBAT p-values remained uniformly distributed, resulting in 10% false positives at p-value threshold of 0.1 (with no FDR correction). CIFBAT Z_r score QIs, based on the interval median, were normally distributed around zero (see Figure S12). As the amount of missing data increases, the intervals widen, and fewer intervals lie completely below a given threshold. At 1% missing data 5.1% of p-value intervals were completely below 0.1 while at 10% missing data 0.3% of p-value intervals were below 0.1. However, when the CIFBAT method is used as described (computing FDR based on FBAT and comparing intervals to the p-value threshold), zero simulated markers are found to be significant at a 10% FDR level, as desired.

Across all scenarios, for both FBAT and CIFBAT, the smallest effect on the recall occurred when missing data was concentrated on the controls, and the largest effect occurred when missing data was concentrated on the cases. In general, the cases will have more occurrences of a causative variant, so when missing data is concentrated on the cases, it weakens the power of the tests to detect the overabundance of causative alleles in cases, leading to false negative results. On the other hand, when missing data is concentrated on controls, it does not affect the signal to the same extent. Figure 3A shows the difference in recall for a sliding FDR threshold when missing data is concentrated on the cases or controls.

The simulations showed that overall, in comparison to FBAT, CIFBAT tended to trade lower recall for higher precision; meaning that while fewer variants were called significant, they were more likely to be true positives (Figure 3B). For example, in the scenario where missing data was concentrated in cases, the CIFBAT recall fell from 0.48 at 1% missing data, to 0.26 at 10% missing data, compared to 0.52 at 1% to 0.40 at 10% for FBAT (Table 4). However, the CIFBAT precision rose from 0.96 to 1.0 over the same range of missing data, compared to 0.93 for all ranges of missing data with FBAT (Table 4). CIFBAT specificity and negative predictive value remained above 99% over all missingness scenarios and levels (Table S12). This suggests that CIFBAT can be useful in detecting potential bias in the FBAT statistic due to missing data and in refining the set of significant markers to be validated for downstream analyses.

Ultimately, the levels of recall and precision are influenced by the disease model parameters. For example, changes in recall and precision can be seen in Figures S3, S4 respectively when the mean penetrance is varied in the model. A complete report of changes in various performance metrics of FBAT and CIFBAT with changes in model parameters is given in Table S12.

Overall, the F-measure, which is the harmonic mean of precision and recall, was comparable between FBAT and CIFBAT, with FBAT having a median F-measure of 0.67 and CIFBAT having a median F-measure of 0.64 over various simulation scenarios (Table 4). However, the variance of F-measure was greater for CIFBAT (Figure 3C and Figure S5).

The simulation also generated results where the CIFBAT QIs were significant but the FBAT statistic was not, although these were few in number (5/300) and no clear pattern emerged to explain their occurrence. Figures S6, S7 show how precision and recall of FBAT and CIFBAT are affected by percentage of missing data. While precision of CIFBAT is higher than FBAT in presence of missing data, its recall decreases with missing data. Figures S8, S9 show effect of the FDR threshold chosen to indicate significant results.

We analyzed familial whole genome sequencing data to identify potential genomic markers for maternal uterine anomalies. Further, we utilize this section to explain interpretation and usage of CIFBAT QIs in validating robustness of FBAT results in presence of missing genotype data.

We compared the counts of significant markers exclusively and jointly identified by FBAT and CIFBAT. Figure 4A shows that out of the 551 markers that were significant under FBAT, 242 (~44%) were also significant under CIFBAT when incomplete trios were included in the test. Validation under CIFBAT indicates the robustness of these results against missing data and serves as a way to select features with higher precision for downstream analyses. All but one of the 242 CIFBAT markers were either intergenic or were positioned in non-coding regions of their respective genes. Figure 4B compares the FBAT statistic and CIFBAT quantile interval for chr7:142008644 (GRCh37), which lies in the second exon of gene TCRBV9S1A1T. This marker was significant under FBAT with p-value less than machine epsilon (2.22e-16), as well as under CIFBAT with a 95% QI [<2.22e-16, 3.23e-10]. Figure 5 illustrates distribution of complete and incomplete trio types between case and control trios for this marker. A complete list of all the markers significant under both FBAT and CIFBAT can be found in Table S13.

The remaining 309 (~56%) markers, significant under FBAT and not under CIFBAT, indicates that these FBAT results might have been affected by missing data. In such cases, CIFBAT is useful in identifying potential false positives and refining the feature set for further downstream analyses. Example of one such marker is chr15:29443416 (GRCh37; rs7171494) which is an intronic variant in gene FAM189A1. Figure 4C compares the FBAT statistic and the CIFBAT QI for this marker, and Figure S10 illustrates the underlying distribution of complete and incomplete trio types between case and control trios. A complete list of all the markers that were significant under FBAT, but could not be validated by CIFBAT can be found in Table S14.

Additionally, there were 38 markers that were exclusively significant under CIFBAT (Table S15). While it is not possible to make conclusive remarks about this set of markers without further independent validation, these could either be categorized as markers where CIFBAT had improved statistical power due to increased sample size or as false positive results. An example marker from this set is chr11:55861880 (rs117149792) which is an upstream gene variant for gene OR8I2. Figure 4D compares the FBAT statistic and the CIFBAT QI for this marker, and Figure S11 illustrates the underlying distribution of complete and incomplete trio types between case and control trios.

Overall, CIFBAT can be a useful test to evaluate robustness of FBAT statistic against missing data. Markers that were significant under FBAT and further validated by CIFBAT after inclusion of incomplete families provide a refined set of candidate markers for downstream analysis. In contrast to refining the set of candidate markers by ~50%, CIFBAT exclusively identified only a small fraction (6.8%) of significant markers collectively identified (589) by FBAT and CIFBAT. These can be potentially interesting markers that warrant further investigation. This feature of CIFBAT is particularly useful in genome wide association studies where the initial set of markers to be evaluated for candidacy is invariably very large. CIFBAT can be used in these studies to refine the set of markers for downstream analyses. Moreover, often the underlying genetic model or the pattern of missing data is unknown; in such cases, CIFBAT can be a great advantage in evaluating the effects of missing data without any assumptions about the underlying model.

Table 5 lists results from the analysis of Type 1 Diabetes data for 351 markers within 22 candidate genes. The p-values from FBAT were corrected for multiple testing using Benjamini-Hochberg method (Benjamini and Hochberg, 1995). A cut-off of 10% false discovery rate (FDR) was used to define significance. For FBAT, the p-value cut-off at 10% FDR was 1.10e-02, which was also used to indicate significant QIs from CIFBAT.

Twelve (12) markers in the insulin (INS) gene were identified as significant by FBAT based on complete trios. Of these, 11 markers also showed significant QIs under CIFBAT when incomplete trios were included in the test (example shown in Figure 6A). For rs6356, inclusion of incomplete trios caused the distribution of z-scores to overlap with the null distribution making the lower bound of the 95% CI insignificant as shown in Figure 6B.

Seventeen (17) markers in the PTPN22 gene were identified as significant by FBAT. Of these, 4 were further validated by CIFBAT after inclusion of incomplete trios, but the remaining 13 were not significant anymore (example shown in Figure 6C).

Two (2) markers in the CTLA4 gene were significant under FBAT, but inclusion of incomplete trios under CIFBAT produced insignificant QIs for both the markers.

Four (4) markers in the interleukin (IL4R, IL2RA, IL12B) family of genes, and 1 marker in the IFIH1 gene were significant under FBAT and were further validated by CIFBAT after inclusion of incomplete trios in the test.

Overall, CIFBAT provided further validation of 20 (55.55%) out of 36 markers that were significant based on FBAT, indicating that these results were not biased by missing data. The remaining 16 markers were not significant anymore when incomplete trios were included under CIFBAT, suggesting that the FBAT test might have been biased due to missing data.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.