The efficiency of mapping of quantitative trait loci using cofactor analysis in half-sib design

This simulation study was designed to study the power and type I error rate in QTL mapping using cofactor analysis in half-sib designs. A number of scenarios were simulated with different power to identify QTL by varying family size, heritability, QTL effect and map density, and three threshold levels for cofactor were considered. Generally cofactor analysis did not increase the power of QTL mapping in a half-sib design, but increased the type I error rate. The exception was with small family size where the number of correctly identified QTL increased by 13% when heritability was high and 21% when heritability was low. However, in the same scenarios the number of false positives increased by 49% and 45% respectively. With a liberal threshold level of 10% for cofactor combined with a low heritability, the number of correctly identified QTL increased by 14% but there was a 41% increase in the number of false positives. Also, the power of QTL mapping did not increase with cofactor analysis in scenarios with unequal QTL effect, sparse marker density and large QTL effect (25% of the genetic variance), but the type I error rate tended to increase. A priori, cofactor analysis was expected to have higher power than individual chromosome analysis especially in experiments with lower power to detect QTL. Our study shows that cofactor analysis increased the number of false positives in all scenarios with low heritability and the increase was up to 50% in low power experiments and with lower thresholds for cofactors.


INTRODUCTION
The purpose of mapping quantitative trait loci (QTL) in livestock is to identify chromosomal regions affecting quantitative traits and ultimately to use the existing variation in those chromosomal segments to select superior individuals from a population. Detection of QTL in outbred half-sib family structure has mainly been based on interval mapping of single QTL on individual chromosomes [6,7,14]. These methods do not take possible QTL on other chromosomes into account. Jansen [10,11] and Zeng [20] proposed methods to account for linked and unlinked QTL by fitting markers as cofactors. These methods were developed for inbred line cross experiments, and Jansen et al. [12] and Kao et al. [13] described methods for multiple QTL mapping in outcrossing species. Cofactor analysis is expected to have higher power than individual chromosome analysis and especially in experiments with lower power to detect QTL. The argument put forward for increased power is the decrease in the residual variance by correcting for variance explained by the cofactors [3,20]. Consequently, the test statistic, which is a function of residual variance, increases and the power of detecting additional QTL increases. The increase can be substantial, since putative QTL become part of the complete model in which all cofactor effects are estimated jointly to give the best fit of the data. De Koning et al. [3] identified initially five QTL for milk yield, significant at 5% chromosome-wise threshold with individual chromosome analysis. Using cofactors, they observed eight QTL exceeding the 5% genome-wise significance threshold. Bennewitz et al. [1] also observed a 39% increase (from 18 to 25) in the number of QTL detected when cofactors were included in the model in comparison to single chromosome analyses. Schulman et al. [16] reported 67% increase (from 6 to 10) in the number of QTL detected with cofactor analysis in comparison with single chromosome analysis. Viitala et al. [18] used the cofactor analysis to identify QTL affecting milk production traits. They found only two genome-wise significant QTL when chromosomes were analysed individually. The cofactor analyses detected a total of 31 genome-wise significant QTL from the same data set. Holmberg and Andersson-Eklund [9] also used the cofactor analysis in half-sib design to identify QTL affecting health traits in dairy cattle. The increase in QTL number may be due to an increase in power or it may be due to an increase in the type 1 error rate. Churchill and Doerge [2] suggested that taking account of a priori known major gene(s) could significantly increase the power for detecting unlinked QTL effects secondary to the major gene. However, they cautioned that if the major gene(s) are determined in the light of the data from the present experiment, the type I error level of the procedure can not be guaranteed. In cofactor analysis, the cofactors are determined from the same experimental data. It is also not clear how adjusting for a nonexistent QTL effect in non-segregating families affects the analyses of other chromosomes for these families. The cofactor effects are overestimated in the first round of analysis i.e. individual chromosome analysis, which may result in some non-significant regions of the chromosome crossing the significance threshold. Fitting non-significant effects as cofactors could add noise to the analyses and might result in loss of power rather than in gain [3]. Since the number of QTL detected in half-sib designs is substantially higher in cofactor analysis, it is crucial to investigate if the expected type I error rate is preserved in cofactor analysis.
Power to detect segregating QTL is a function of the number of individuals genotyped for the genetic markers and phenotyped for the quantitative traits, the effect of segregating QTL, heritability of the trait, type I error rate allowed and marker density as well as residual variance [17,19]. The threshold level for inclusion of a cofactor will also influence the power as well as type I error rate. Therefore, this simulation study was designed to study the power and type I error rate in QTL mapping including the cofactor with varying family size, heritability, QTL effect, map density and threshold level of the cofactor in order to test the following two hypotheses. Hypothesis I: using cofactors, as implemented in the half-sib design by De Koning et al. [3] increases the power of QTL mapping, especially when the power for detection is low and with liberal threshold for cofactor. Hypothesis II: using cofactors increases type I error rate in QTL mapping in half-sib designs.

Simulation of data
Phenotypes and marker data were simulated for 15 half-sib families. The parameters for nine different scenarios are presented in Table I. Marker alleles were sampled on 15 chromosomes, each with 12 tetra-allelic markers with equal allele frequencies placed with a distance of 5 cM between loci, except in the sparse-marker density scenario where four markers were placed at 20 cM intervals. Bi-allelic QTL alleles were simulated on five chromosomes situated halfway between markers 7 and 8, except in the sparse marker density scenario where the QTL location was halfway between markers 3 and 4. The QTL alleles were assumed to have equal frequency. In each scenario the QTL explained 50% of the genetic variance and the other 50% of the genetic variance was due to polygenes unlinked to the QTL. In scenario 1-6, each of the five QTL explained 10% and in scenario 7 each of the two QTL explained 25% of the genetic variance. In scenario 8 (unequal QTL effect I), five QTL with diminishing effects were simulated, explaining 13.3%, 11.7%, 10.0%, 8.3% and 6.7%  Table. such that the smallest QTL explained half the variance of the biggest QTL and collectively explained half of the genetic variance. Similarly, in scenario 9 (unequal QTL effect II), the five QTL explained 12, 11, 10, 9 and 8% of genetic variance, such that the smallest QTL explained two-thirds of the variance of the biggest QTL. Each of the scenarios was replicated 100 times. Therefore, there were 500 chromosomes with a QTL and 1000 chromosomes without a QTL simulated in each scenario except in scenario 7 (large QTL effect) where it was 200 and 1300 respectively.
The total phenotypic variance of the trait was assumed to be 100. Each of the nine scenarios was simulated with two different "heritabilities", 0.90 and 0.30. The phenotypes were simulated on each offspring for 15 half-sib families by adding the sire effect, effect of QTL allele inherited from the sire and error component. Therefore, the heritability of 0.90 resembles a granddaughter design with a heritability normally observed for yield traits in granddaughter designs in dairy cattle. Similarly, the heritability 0.30 resembles the low heritability traits like disease resistance traits in a granddaughter design. This also resembles the heritability observed in milk yield traits in a daughter-design in dairy cattle. Three family sizes were considered with 100, 50 and 25 offspring per half-sib family. There were three threshold levels for considering cofactors in the analysis: 0.05, 0.01 and 0.10.

Cofactor analysis
Each replicate was analysed using cofactor analysis followed by analyses under a multiple QTL model using a modification of the methodology proposed by De Koning et al. [3]. In the first step of the multiple QTL analyses, the candidate regions were analysed individually fitting a single QTL within every family: where Y i j is the phenotype of j, offspring of sire i, μ i is the effect of being in sire family i, b i is the allele substitution effect of the QTL within family i, X i j is the probability that animal j inherited the (arbitrarily assigned) first haplotype of sire i, and e i j is the residual effect. In the second step, the best positions on every chromosome that exceeded a point-wise 5% threshold were identified and used as cofactors to re-analyse all the regions using the following model: where the variables are identical to (1), except for the term n k=1 b ik X i jk , which describes the multiple regression of the n cofactors that were on chromosomes other than the one under study. If this analysis revealed additional putative QTL, or the best positions of the QTL changed, the selection of cofactors was modified and the regions were re-analysed. This step is repeated until no new QTL are identified or dropped from the model, and the positions of the QTL are stable. In the present study, this step was repeated five times in each simulated data set. The difference between this analysis and that of De Koning et al. [3] is that in the present study the cofactors were fitted in the model simultaneously with the putative QTL, while De Koning et al. [3] adjusted the trait scores for cofactor effects prior to re-analysing the chromosomes. The number of significant QTL was counted in the first round (single chromosome analyses) and the fifth round and the locations of true QTL were recorded. An identified QTL was considered to be a true positive when it was identified to reside on a chromosome where a QTL was simulated. The QTL location was considered accurate when it was observed in the maker interval where it was simulated and the accuracy goes down with the QTL location drifting away from the correct marker interval. The chromosome-wise threshold level for QTL and cofactors was uniform at 0.05 except for scenario 4 (stringent threshold) and scenario 5 (liberal threshold) where the chromosome-wise threshold to cofactors were 0.01 and 0.10 respectively. When applying an empirical 5% chromosome-wide threshold, on average 5% of chromosomes without QTL in the study will be called significant. At this significance level we expect to get an average of 50 false positives (FP) in each scenario, because in 100 replications there was a total of 1000 chromosomes without QTL. The expected number of false positives in the "large QTL" scenario was 65 (5% of 1300 non-QTL chromosomes). Chromosome-wide significance thresholds were determined empirically by 1000 permutations [2], where phenotypes of half-sib offspring were randomised within half-sib families while retaining genotypes. For the cofactor analyses, permutations were performed for the chromosome under study while maintaining the link between the genotypes and phenotypes for the cofactors on other linkage groups. For a type I error rate of 0.05, a sample of 1000 permutations is usually regarded as sufficient [5,15]. A putative QTL was included as the cofactor when it exceeded the chromosome-wide threshold that was set for a given scenario. Van der Beek et al. [17] described a method for estimation of power for a half-sib design for a single interval analysis. This was used to estimate the theoretical power for each scenario for a single test taking the parameters used for simulation of data. Table II shows that cofactor analysis did not increase the power to locate QTL in the large half-sib family size scenario but there was a small increase in  Proportion (%) of QTL identified in the correct (0), ± 1, ± 2, ± 3, and beyond 3 marker interval.

RESULTS
power with medium and small family size scenarios. However, the false positives (FP) generally increased when cofactor analysis was used. The exception was with the large family in the high heritability situation, where FP decreased by 5%. In cofactor analysis with the small family size scenario, there were increases of 13% and 21% in the number of true QTL identified in high-and low heritability conditions, but also 49% and 45% increases in FP. Out of a total of 47 new significant positions found in cofactor analysis in small family size with high heritability, 24 were false positives. The accuracy of QTL location estimates remains very similar in both individual and cofactor analysis (Tab. III). With small family size (25 offspring/family) the power of detection of true QTL increased in cofactor analysis under both high-and low heritability conditions.  Table IV shows that the level of threshold (0.01, 0.05 and 0.10) for a chromosomal position to qualify as a cofactor generally did not change the power to detect QTL but a liberal threshold increased the level of FP. With a stringent threshold and high heritability there was a 13% decrease in FP in cofactor analysis, though the power to identify true QTL was the same as that observed in single chromosome analyses. However, in the case of small family size and high heritability, there was a 34% increase in FP. Out of 21 new locations, that became significant in this scenario with cofactor analysis, only six were true QTL and the rest were false positives. With low heritability and liberal threshold levels for cofactors at 0.10, there was a 13% increase in power. However, FP increased for both stringent and liberal threshold scenarios (Tab. V). The biggest increase in FP (41%) was in the case of a liberal threshold and lowheritability scenario. The theoretical power calculated for the single interval test (Tab. VI) based on the population design and the parameters used to simulate the data were very close to the power observed for single chromosome analyses. Out of the three levels of threshold used, the power was only increased in cofactor analysis in comparison to individual QTL mapping in the low heritability and liberal threshold scenario. The accuracy of the estimates of QTL location was very similar in both individual and cofactor analysis (data not shown).
The results of QTL analysis presented so far had five QTL each explaining the same amount of variance. Table VII shows that 'unequal QTL effects' does not change the power of detection of QTL, while the FP was increased or decreased in different scenarios. These two 'unequal QTL effect' scenarios with high heritability had similar power for both individual and cofactor analysis. Similar to the default scenario in the high heritability scenario, there was a 7% decrease in FP in scenario 9 (unequal QTL effect II) with the cofactor analysis. However, a 15% increase in FP was observed with scenario 8 (unequal QTL effect I) with cofactor analysis. In low heritability condition, there was an increase in FP in both scenarios 8 and 9 (Tab. V). The observed powers in both cases were lower than in the default situation (equal QTL effect). No substantial changes in QTL location estimates were observed when cofactors were included in the analysis in comparison to when they were not included (data not shown). Table VIII shows that there was no increase in the power of QTL detection with cofactor analysis in the "sparse marker" and "big QTL" scenarios. Though a decrease in FP was observed in high heritability condition in these scenarios, FP increased in the low heritability condition. With high heritability, the power   of the sparse marker scenario was 0.86 for single chromosome analyses compared to 0.97 observed in the default situation (marker interval of 5 cM). A 7% decrease in FP was observed with cofactor analysis in high heritability condition, while FP increased by 16% in the low heritability condition. In the "big QTL" scenario, power to identify QTL was 100% in high heritability condition and 0.91 with low heritability for both individual and cofactor analysis (Tab. VI) and both these two estimates were similar to their theoretical expectations. However, using cofactor analysis changed the FP for the two levels of heritability (Tab. V). With high heritability there was a 23% decrease in FP with cofactor analysis, while with low heritability, there was a 34% increase in FP.

DISCUSSION
Cofactor analysis as implemented in half-sib designs by De Koning et al. [3] does not change the power of QTL mapping but generally increases the number of false positives (FP). Generally, there was a small or no difference in power to detect QTL between single chromosome analyses and cofactor analyses. There was two exceptions in low heritability conditions where the cofactor analysis increased the power from 11% to 14% in scenarios with small family sizes and from 45% to 51% in scenarios with liberal thresholds to include a cofactor in the model. However, in those scenarios, the FP increased between 40% and 50%. This means that the expected type I error rate is not maintained in the iterative scheme of fitting cofactors in the model. Consequently, the two power measures were not comparable and the increase in both power and FP might be due to the increased probability level of type I error (α). For example, in the small family size and low heritability condition, the power increased from 0.11 to 0.14 with cofactor analysis. However, if single chromosome analyses are performed with a type I error rate of 7.4%, which is the empirical type I error rate observed with cofactor analysis in this scenario, the theoretical power would increase from 0.14 to 0.18. Consequently, the increase in power observed can be explained by the increase in type I error rate. Therefore, this study did not support the first hypothesis that cofactor analysis improves power in the half-sib design. However the FP was generally higher for cofactor analysis than for single chromosome analysis and the difference was particularly pronounced (up to 50% increase) in low power scenarios and when liberal thresholds for cofactor inclusion were used (Tab. V). Therefore, the second hypothesis supported that cofactor analysis [3] gives a higher type I error rate in half-sib design than single chromosome analysis in QTL mapping does. These results indicate that if cofactor analysis is used, there is a high risk that additionally identified QTL are false positives. We covered a wide range of scenarios with varying power. The "high heritability" in the present study, simulates the situation for most yield traits that are routinely measured in progeny testing of AI bulls in dairy cattle and "low heritability" resembles the heritability observed for disease traits. However, we did not observe an increase in power with cofactor analysis as reported earlier [1,3,16,18,20]. De Koning et al. [3] reported, using cofactor analysis, that the initial number of five suggestive QTL had increased to eight significant QTL in a study with data with an average of 41 sons per sire and ranging from 21 to 82. The current simulations would suggest that with such small family size there might be a high false positive rate among the additional QTL.
The FP was generally higher for cofactor analysis than for single chromosome analysis. The difference was particularly pronounced in low power scenarios and when a liberal threshold for cofactor inclusion was used. The main factors causing the difference in power were the family size and the genetic variance explained by the QTL (heritability). The factors influence the FP in a different manner. The powers of QTL detection were very similar in scenario 1 "low heritability and large family" (53%) and scenario 3 "high heritability and small family" (47%). Nonetheless, the increase in FP when using cofactor analysis was much more marked with the small family size, which was also the scenario with high heritability, with the largest increase in FP. This is most likely because in a small family there is a much higher probability of observing a random association between the sampled marker alleles and the sampled phenotypes just by chance. Consequently, when a QTL is included as a cofactor the markers will tend to fit some of the residual variance in families where no QTL segregates. This explanation agrees with the observation that with liberal thresholds for cofactor inclusion the FP increased both in high-and low-heritability scenarios. In this case there is an increased risk that the fitted cofactor is a false positive and even though the detection of a QTL uses a higher threshold, the expected type I error rate is not controlled in the QTL test because non genetic parts of the residual variance are partly fitted by the cofactor. Viitala et al. [18] identified two genome-wise significant QTL for milk production traits in dairy cattle when chromosomes were analysed individually. Cofactor analysis with the same data, detected a total of 31 genome-wise significant QTL. The increase in QTL number from 2 to 31 in that study shows the increase in false positives when cofactors similar with the present simulation study are fitted. The increase in FP supports the apprehension of Churchill and Doerge [2] that if gene(s) that are detected in the same experiment are fitted in the model for detecting unlinked QTL effects secondary to the major gene, the type I error level of the procedure cannot be guaranteed. Doerge and Churchill [5] carried out a sequential search for multiple QTL and reported a lower type I error (lower than expected at a 5% significant level). However, in two scenarios with 2 QTL and 3 QTL with sample size 200, they observed high "overcall" (non-simulated QTL were identified), 38 and 45 times in 500 Monte Carlo simulations, respectively. The lower type I error in some scenarios may be due to a very high QTL effect simulated on the phenotype in that study (50% of the total variance was explained by a QTL). Whereas, in the present study, a maximum 15% and minimum 3% of the total variance was attributed to a QTL in different scenarios. We also observed lower FP than expected (55 observed vs. 65 expected) in the large QTL effect scenario with high h 2 , where each QTL explained 7.5% of the total variance. The present simulation of QTL effects reflects practical situations for quantitative traits where many genes with small effects have been postulated more closely [8]. Here one specific implementation of cofactor analysis in half-sib design was studied, but in the light of the comments of Churchill and Doerge [2] and our results, it seems that an increase in type I error is a general phenomenon in cofactor analysis.
A liberal threshold for inclusion of cofactor increases FP. One of the options in cofactor analysis is to use a different significance threshold for cofactors compared to the identification of QTL. Generally, a liberal threshold is fixed for cofactors and a stringent threshold for final selection of QTL. For instance, De Koning et al. [3,4] used a 5% chromosome-wide threshold to include the cofactor, but a 5% genome-wide threshold to test significance of the QTL. They suggested that in practice, different thresholds should be used for picking cofactors and declaring a significant QTL and they recommended that a cofactor should be selected at a 5% chromosome-wide level, while for declaring that cofactor as a QTL, genome-wide thresholds should be applied [4]. In one scenario in the present study, we considered two levels of threshold, 10% for cofactors and 5% for QTL. The number of true QTL identified increase by 1.3% and 13.8% in single chromosome analysis and cofactor analysis with high-and low-heritability condition respectively. However, the results showed that using a liberal threshold for inclusion of cofactors also increases FP. With a liberal threshold and high heritability, 71% of the new QTL identified were false positives. This means that a new significant region identified using a liberal threshold has a 71% probability of being a false positive. Similarly, in a low heritability situation with a liberal threshold, the probability of a newly identified region being a true QTL is only 60%. Therefore, a liberal threshold for cofactor inflates the Type I error rate of the QTL test. This is most probably a general property and therefore, a liberal threshold for cofactors should not be used. Schulman et al. [16] used a lower threshold (P < 0.10) for the selecting cofactor. The average number of sons per grandsire in their study was 40 (21 to 82). With such small family size and liberal threshold, the chance of four additional QTL identified with cofactor analyses being false positive is quite high. However, cofactor analysis can help to decrease the type I error rate, if the experiment has very high power to detect QTL and a stringent threshold for cofactors is used. With stringent threshold FP, decreased by 13% in the high heritability scenario. A maximum decrease in FP was 22% with a cofactor analysis in a "large QTL" scenario where the power of identifying QTL was 100%. Therefore, cofactor analysis in its present form should only be used in experiments with very high power to reduce FP rather than to increase power. This recommendation is in contradiction with De Koning et al. [4] who suggested applying a liberal threshold for cofactors and a stringent threshold for declaring a significant QTL. The results indicate that one of the main reasons of increase in FP in cofactor analysis may be fitting cofactor effects in families where the QTL is not segregating. This would explain the finding that the effect is more pronounced in small families. Therefore, cofactor analysis may be improved if cofactors are only fitted within families with a significant association. With such a change, cofactor analysis may control the nominal false positive rate and thereby be applicable in more general situations.