A Combinatorial Partitioning Method to Identify Multilocus Genotypic Partitions That Predict Quantitative Trait Variation

Step 1: Search and Evaluation of the State Space

Searching the state space of all possible ways to partition m-locus genotypes into sets of genotypic partitions for all subsets of l loci can be separated into two nested combinatorial operations, illustrated in Figure 2 for m = 2. The first operation, described in Figure 1, Step 1, and illustrated in Figure 2A, consists of systematically selecting all possible subsets of variable loci for the desired values of m, of which there are ways. For each subset M, the m-locus genotypes are identified, illustrated in Figure 2B by the 3 × 3 grid for two diallelic loci A and B. The second operation, depicted in Figure 2C, is to evaluate the possible sets of genotypic partitions over the desired range of k. The number of ways to partition g_M genotypes into a set of k genotypic partitions is known as a Stirling number of the second kind (Comtet 1974), computed from the sum

For diallelic loci, the maximum number of m-locus genotypes that could be observed is g_M =3^m. For two diallelic loci with m = 2 and g_M = 9, there are 21,146 ways to partition G_M into k = 2, 3, …, 9 partitions. Combining these two combinatorial operations (choosing M from L and enumerating all K possible from G_M) in the example application below where the number of variable loci l = 18, if m = 2 and g_M = 9 for all M in L, results in 3,235,338 possible sets of genotypic partitions. At this size it is practical to search the two-locus state space exhaustively. However, by adding even one more variable locus such that g_M = 27, and we partition G_M into k = 2, 3, …, 27, the number of sets of genotypic partitions for a single combination of three loci is ∼10²⁵, clearly out of the range of what can be exhaustively enumerated. The computational feasibility of the CPM is further addressed in the discussion below. In this article, to illustrate the CPM, we restrict our application to the two-locus case, and evaluate directly all possible sets of genotypic partitions. This complete enumeration of the state space provides additional advantages. As we have no a priori knowledge of how the state space of sets of genotypic partitions, as defined above, relates to the measured phenotypic characteristics, enumeration in low dimensions allows us to examine this relationship directly. Also, observing the nature of this relationship over the entire state space can aid in the development of algorithms for searching higher-dimensional spaces.

To evaluate the state space, we must select a statistical function that provides a measure of the value of each set of genotypic partitions. This statistical function is often referred to as the objective function. The sum of the squared deviations of the trait means of the partitions from the overall sample mean (SS_K) is a natural choice for an objective function of a set of genotypic partitions. The value of this measure increases as the individuals within genotypic partitions become more phenotypically similar and as the phenotypic differences between partitions increase. A disadvantage of the partition sum of squares is that it will tend to increase as k increases, favoring the division of genotypes into a greater number of partitions. The bias-corrected estimate of genotypic variance (Boerwinkle and Sing 1986) is used to compensate for this bias. It can be written as

where n is the total sample size, Y is the sample grand mean, n_i is the sample size and Y_i is the mean of partition i, Y_ij is the phenotype of the jth individual in the ith partition, and MS_W is the mean squared estimate of the phenotypic variability among individuals within genotypic partitions. This statistic was derived using the expected mean squares from a standard one-way analysis of variance model to obtain an unbiased estimator of the measured genotypic variance. This correction has the effect of penalizing the scaled partition sum of squares by a quantity that increases with k if the estimate of MS_W does not decrease as additional partitions are considered. For comparative purposes, the proportion of variability explained by a set of genotypic partitions is preferred for general usage over s²_K. The proportion, denoted pv_K, is computed as

Note that the pv_K statistic, which is a function of unbiased estimators, is not an unbiased estimator of the proportion of variability explained by set K (Wijsman and Nur 2001). However, it is appropriate for our purpose of comparing sets of genotypic partitions.

For most traits and combinations of variable loci, most of the sets evaluated will explain very little phenotypic variability and will not be useful for predictive purposes. For this reason, some filter is needed to decide which sets of genotypic partitions should be retained for further consideration in Step 2 of the CPM approach (Fig. 1). Many criteria could be used in setting this filter, including criteria based on the significance level of an F-test (e.g., MS_K/MS_W > F[α = 0.005; k − 1, n − k]), biological significance (e.g., pv_K > 0.05), or some proportion of all of the sets considered (e.g., the top 1%, the top 100, or simply the best). In the application described below, we use a cutoff based on the distribution of the F statistic.

In addition to the proportion of variability that a set is predicting, we include an additional criterion in the evaluation of each set of genotypic partitions. Variable loci that have alleles of low relative frequency are expected to result in some multilocus genotypes being represented by only a few individuals in any given sample. It follows that some of the sets will also have few individuals in a genotypic partition. Because our method of evaluating sets is based on both within- and between-partition sums of squares, sparse partitions with only one or a few individuals do not have sufficient degrees of freedom to reliably estimate both partition means and partition sums of squares. To avoid these cases, we set a lower bound on the number of individuals observed for a partition to be included in the CPM evaluation. The value of this lower bound is dependent on the judgement of the investigator. Setting this bound too low could lead to spurious results, as could be the case if a partition contained one or a few outliers. However, setting this bound too high could result in failing to consider sets that could be useful for predicting phenotypic variability and providing insights into the trait biology. In the application described in this article, we set a lower bound of five individuals in a valid partition.

Step 2: Validating the Selected Sets of Genotypic Partitions

The next step in the CPM is to validate the retained sets of genotypic partitions (Fig. 1, Step 2). A very large number of sets are considered in the process of searching and evaluating the state space of all sets of genotypic partitions. This problem is common to the methods usually applied in the fields of data mining, pattern recognition, machine learning, and artificial intelligence. These loosely related fields have long struggled to find appropriate ways to validate the model or collection of models that were constructed as a consequence of considering the large number of possible relationships that are present within a particular data set (Ripley 1998, section 2.7). One of the most common methods of model validation is multifold cross-validation (Stone 1978), where the number of folds is typically 10 (Kohavi 1995). This method simulates the process of going back into the population of inference and collecting an independent sample with which to validate the constructed models. Briefly, 10-fold cross validation is carried out by randomly dividing the sample into 10 approximately equal-sized groups. The first group is removed from the sample, and the remaining nine are used to estimate the parameters of the model, which in the case of the CPM are the means of the k genotypic partitions of a set. Then the one group that was withheld is used to compute a portion of the predicted within partition sum of squares. This process is repeated 10 times, with each of the 10 groups being withheld from the parameter estimation and used to compute a portion of the predicted within partition sum of squares exactly once for each of the 10 groups. The 10 portions of the predicted within partition sum of squares are then summed to provide an estimate (SS_W,CV) of the cross-validated within partition sum of squares. To reduce the possibility that a particular random assignment of the sample into 10 cross validation groups might favor one set over another, the random 10-fold random division of the sample and cross validation is repeated 10 times and the resulting cross-validated within-partition sums of squares are averaged for each set of genotypic partitions.

This cross-validated estimate (SS_W,CV) of the trait variability within each genotypic partition can then be used to judge the predictive ability of a given set of genotypic partitions. For consistency and comparability, we use SS_W,CV and SS_K,CV = SS_Total −  SS_W,CV in place of SS_W and SS_K in the equations above to calculate a cross-validated proportion of variability explained, denoted as pv_K,CV. The larger pv_K,CV becomes, the more predictive the set K is said to be. Note that SS_W,CV must be ≥SS_W, which implies that pv_K,CV will be ≤pv_K.

Step 3: Select the 'Best' Sets of Genotypic Partitions and Make Inferences about the Relationship between Phenotypic and Genotypic Variation

Steps 1 and 2 provide a strategy for identifying sets of genotypic partitions that predict variation in quantitative trait levels. The third and final step in this approach (Fig. 1) is to select some subset of the validated sets of genotypic partitions on which inferences can be made. The utility of information about these sets depends on the goals of the investigator and the questions that are being asked. One use could be to identify the overall “winner,” for eaxmple, the set of partitions that best predicts phenotypic variation. However, many other goals cannot be addressed by considering only one “best” set. For instance, if we are interested in enhancing our understanding of the relationship between phenotypic variability and genotypic variability at multiple loci, there may be information about the penetrance function or norm of reaction to be gained by scrutinizing a subgroup of multiple sets of genotypic partitions that may be almost as predictive as the overall best predictive set. When making inferences about such a group of multiple sets, we may ask the following questions: How much trait variability does each of the selected sets explain? Which combinations of loci are involved in the selected subgroup of all possible sets? If there are multiple combinations of loci represented in these selected sets, how do the genotype–phenotype relationships in one combination compare with another? Are the selected loci acting additively on the variance of the trait? How do the genotype–phenotype relationships observed in the selected sets compare with the relationships expected using traditional multilocus models? What criteria should be used in making the subjective decision of how many sets should be selected for drawing inferences about genotype–phenotype relationships? An application of the CPM to illustrate it's utility in providing analytical results to deal with these questions is presented next.

Sample

To illustrate the combinatorial partitioning method, we present an application to identify sets of genotypic partitions of two-locus combinations of 18 diallelic loci (16 single nucleotide polymorphism [SNPs] and two insertion-deletions [InDels]) located in six candidate CHD susceptibility gene regions (Table 1) that predict interindividual variability in plasma triglycerides (Trig) levels, a known CHD risk factor. The sample used here comes from the Rochester Family Heart Study (RFHS), a population-based study of multigeneration pedigrees sampled without regard to health status from the Rochester, Minnesota, population (Turner et al. 1989). Because the genetic influence on variation in most quantitative trait levels is influenced by age (Reilly et al. 1992; Zerba et al. 1996; Jarvik et al. 1997; Nelson et al. 1999), we restricted our analyses to adults ages 20–60 yr. From 281 pedigrees, we selected 188 males for this sample application. All individuals were within the prescribed age window and had Trig measurements and complete genotypes for all 18 variable loci. In this sample, the distribution of Trig was right skewed (skewness = 2.7) and leptokurtic (kurtosis = 9.6). In accordance with common practice in the study of lipid metabolism, we applied a natural logarithm transformation (lnTrig). This transformation greatly reduced both the skewness (0.83) and the kurtosis (1.04). Finally, to reduce the influence of age and body size on lnTrig variability, we adjusted lnTrig using a linear regression model that included age, height, and weight (each up to third-order polynomial terms) and waist-to-hip ratio and body mass index.

Analysis

The single-locus contribution to adjusted lnTrig (henceforth referred to as lnTrig) variability (pv) is shown in Table 1 for each of the 18 diallelic loci. The bias-corrected estimate of the proportion of lnTrig variability explained ranges from 0.000 for nine of the loci to 0.021 for the StuI locus in the LDLR gene. The variability among single-locus genotypes was not statistically significant at the level of α = 0.01 for any of the 18 variable loci.

Step 1: Evaluation of the State Space of Sets of Genotypic Partitions

As discussed above, we limited the state space in this application to all possible sets of genotypic partitions (K: k = 2, 3, …, 9) defined by all possible pairs (m = 2) of the l = 18 diallelic loci. We set five as the lower bound for the number of individuals that must be present to constitute a valid partition. Sets containing a partition with fewer than five individuals were excluded from further evaluation. To filter the potentially millions of valid sets of genotypic partitions for validation in Step 2, we selected a criterion based on the test-wise significance of each set. All sets with an F statistic that exceeded the 0.995 quantile of the F distribution corresponding to each k considered, that is, F(α = 0.005; k − 1, n −k), were retained for further validation. This combination of criteria for k = 2, 3, …, 9 resulted in the consideration of 794,699 sets and the retention of 7710 sets (1.0%).

The number of sets considered is much lower than the 3,235,338 possible sets associated with all two-locus combinations of 18 diallelic loci. There are two reasons for this. The first is that for most of the pairwise combinations of loci, fewer than nine two-locus genotypes were observed, such that the total number of enumerable sets of genotypic partitions for this sample is 1,001,270. The second reason is that 206,571 sets were not evaluated because they contained partitions with fewer than five individuals. The proportion of variability explained by each of the retained sets is displayed in Figure 3. Each line in the plot corresponds to the retained sets for a different number of partitions (k = 2, 3, …, 8). The retained partitions for each k are sorted by the proportion of variability explained for each set, with the sets that explain the least at the left of the plot and the sets that explain the most at the right.

There are a few notable features of the plot in Figure 3. First, the number of sets retained varies greatly by k, ranging from 11 (k = 8) to 2911 (k = 4). Second, the sets that explain the greatest proportion of variability for each k range from pv_K = 0.073 (k = 8) to pv_K = 0.093 for (k = 3). Recall that the retained sets that explain the least proportion of variability in lnTrig correspond to the minimum established by the F distribution with the appropriate degrees of freedom for each k. The third notable feature in Figure 3 is the shape of the lines. There are many sets of genotypic partitions near the cutoffs that explain roughly the same proportion of variability, and relatively few that explain substantially more. For example, for k = 3, 1357 (94%) of the retained sets explain between 0.045 and 0.069 (lower 50% of the range in pv_K), while 81 (6%) explained between 0.070 and 0.093 (upper 50% of the range in pv_K), and only 12 (0.8%) explained within the upper 25% of the range.

While Figure 3 provides a useful description of the proportion of variability being explained by the retained sets of genotypic partitions, it contains no information about which pairs of loci are involved. This information is presented in Figure 4 for k = 2, 3, 4, 5. The remaining k are not presented because of limitations in space and, as it will be shown by cross validation in the next section, all of the sets for k = 6, 7, 8 lack the ability to predict lnTrig variability. In Figure 4, there is one panel corresponding to each k. The Y-axis corresponds to the proportion of variability as in Figure 3. The X-axis corresponds to each pair of loci with a retained set of genotypic partitions, sorted on the basis of the maximum value of pv_K for each pair. The greatest proportion of variability is explained by a set with k = 3 defined by the pair {InDel (APOA1-C3-A4), HincII (LDLR)}. The set of genotypic partitions that explains the greatest proportion of ln Trig variability not involving the InDel and HincII pair of loci also corresponds to k = 3 and is defined by {InDel (APOA1-C3-A4), AvaII (LDLR)} with pv_K = 0.076. Note that HincII and AvaII are both located within LDLR, are separated by a single intron, and have previously been shown to be in very strong linkage disequilibrium in this population (Zerba et al. 1998). As a result, the most explanatory sets defined by these variable loci have very similar assignments of genotypes (and individuals) to those of partitions. The next most explanatory set also corresponds to k = 3 defined by the pair {InDel (APOA1-C3-A4), R192N (PON)}, where pv_K = 0.068. In addition to the most explanatory set overall, the pair {InDel (APOA1-C3-A4), HincII (LDLR)} identified the most explanatory set of genotypic partitions for all k. The most explanatory set in k = 4 explains nearly as much as the best set overall (pv_K = 0.090 versus pv_K = 0.093). The most explanatory sets with two and five partitions had pv_K of 0.082 and 0.087, respectively.

Step 2: Validate Retained Sets of Genotypic Partitions

The 7710 sets of genotypic partitions retained in Step 1 were validated using the repeated 10-fold cross-validation method discussed above. The cross-validation was carried out on each of the retained sets 10 times. Using the predicted within-partition sum of squares (SS_W,CV) obtained from averaging the 10 cross-validation results, pv_K,CV, a measure of the predictive ability of each retained set was obtained. The sets were then sorted by pv_K,CV and plotted in Figure 5, which, analogous to Figure 3, plots a nondecreasing line connecting the ranked sets for each number of partitions with pv_K,CV on the Y-axis. The values of pv_K,CV range from −0.023 (k = 8) to 0.067 (k = 3). All of the sets for k = 6, 7, 8 fall below 0.020, and on the basis of our inspection of Figure 5, we chose to consider these sets as not predictive. The greatest difference between the lines displaying pv_K,CV in Figure 5 and those in Figure 3 is that instead of slowly sloping toward the lower bound, there is a small proportion for each k that is noticeably less predictive than the rest, shown by the sharp drop-off at the left side of each line. As in Figure 3, a majority of the sets fall within a very narrow range of the line. This is particularly pronounced for k = 3, 4, 5. The right-hand portion of the plot showing the best sets for each k again shows that they are outliers from the rest. This is a useful feature for the selection of a subset of sets of genotypic partitions on which to make inferences, focusing attention on just a few sets from among the many possibilities.

To display the relationship between pv_K and pv_K,CV, the values of pv_K for each retained set are also included in Figure 5 (jagged lines), where the color of the line is used to differentiate among k. Though the two measures are correlated, the relationship between them is not one to one. There are many sets with high values of pv_K that do not have correspondingly high values of pv_K,CV. In this particular case, the set with the highest overall value of pv_K also had the highest overall value of pv_K,CV.

The predictive ability of each of the retained sets of genotypic partitions sorted according to combinations of variable loci, analogous to Figure 4, is shown in Figure 6 for k = 2, 3, 4, 5. In comparison to Figure 4, the ranking of the pairs of variable loci is unchanged among the top three combinations. As noted previously, for each k, the top three combinations of loci are {InDel (APOA1-C3-A4), HincII (LDLR)}, {InDel (APOA1-C3-A4), AvaII (LDLR)}, and {InDel (APOA1-C3-A4), R192N (PON)}.

Step 3: Select the Best Sets of Genotypic Partitions and Make Inferences about the Relationship between Phenotypic and Genotypic Variation

Here we are faced with the challenge of restricting our attention to some subset of the 7710 retained sets with which we can make inferences. Our focus for this task is on the results plotted in Figures 5 and 6, based on the cross-validated proportion of variability explained (predictive ability). Clearly, sets that take on negative values for pv_K,CV, which occurs when SS_W,CV is greater than SS_Total, have no predictive ability and would not be useful for making inferences. Rather than select some predetermined level of pv_K,CV as being adequately predictive for inferential purposes, we instead can benefit from the information contained in the distribution of the cross-validated results. Combining the information contained in Figures 5 and 6, we see that there are three sets of genotypic partitions that stand out above the remaining 7707 sets. The values of pv_K,CV for these three best sets are 0.062 (k = 2), 0.065 (k = 3), and 0.067 (k = 3), with corresponding values of pv_K of 0.081, 0.092, and 0.093, respectively. Each of these sets are partitions of nine two-locus genotypes defined by {InDel (APOA1-C3-A4), HincII (LDLR)}. In fact, there are 57 sets defined by these two loci that predict lnTrig variability better than the best set defined by the next best pair of loci, {InDel (ApoA1-C3-A4), AvaII (LDLR)}.

The partitioning of the multilocus genotypes into the three best sets is shown in Figure 7. The genotypes are represented by a 3 × 3 grid, and the assignment of each genotype to a partition is indicated by shade, with the lightest shade for the partition with the lowest mean lnTrig and the darkest shade for the partition with the greatest mean lnTrig. In all three sets, the partition corresponding to the lightest shade (lowest mean) is unchanged. Also unchanged is the presence of (DD,−−) and (ID,+−) in the same partition as well as(ID,++) and (II,+−). The only difference between the two sets with k = 3 is the change of(II,−−) from the darkest partition to the intermediate partition. This change has only a minor effect on both pv_K and pv_K,CV. It is not surprising that this two-locus genotype is the least frequent of the nine, containing seven individuals. For the set with k = 2, the multilocus genotypes within the intermediate and high mean ln Trig in the k = 3 sets are combined into a single partition, which results in a moderate decrease in both pv_K and pv_K,CV.

To obtain an estimate of the statistical significance of the most predictive set, we performed a permutation test (Edgington 1995). The objective of a permutation test is to produce the distribution of the test statistic of interest under the null hypothesis, which in this case is proportion of variability explained by the most predictive set of genotypic partitions identified by the CPM. The test was carried out by disassociating the genotypes and phenotypes in the data set, randomly reassigning the phenotypes to genotypes via sampling without replacement, and then performing the CPM on this randomized data. This process was repeated 1000 times, and for each permutation, the set that explained the greatest proportion of variability was obtained. The proportion of permutations with results greater than the observed estimate of 0.093 for the most predictive set was 0.14.