Block network mapping approach to quantitative trait locus analysis

We demonstrate our method by applying it to simulated data with 742 genotypic sequences of Diversity Outbred (DO) Mice and 1000 phenotypes. We simulated phenotypes on the 19 autosomes and did not add a sex effect. We selected 19 genomic locations, one on each autosome, and generated 19 QTL effect sizes from an exponential distribution. Using the genotypes at each location, we created the QTL effects and scaled the variance to be 1. Then we added N(0,1) noise and the QTL effects together.

We compare our results to that of simulations using Interval Mapping (IM) with expectation maximization from the R/qtl package. For the R/qtl simulations, we use scanone with the default settings and calc.genoprob with step = 0 and error.prob = 0. For the P-values calculations we run scanone with 1000 permutations (n.perm = 1000). The simulated data are obtained by choosing only one SNP that influences a particular phenotype on each of the 19 autosomes with varying effect sizes. These effect sizes range between 1.65×10^?5 and 10.03 (see Additional file 1: Figure S1). We compare the power [5, 27, 28] of our method with that of IM for different effect sizes. With 1000 phenotypes, 19 autosomes, and 1 “true signal” on each autosome, we have 19,000 effect size data points. We arrange them in order of increasing effect size and then divide them into 76 groups of 4000 data points with 200 points offset. For example, the first group is composed of the first 4000 data points with the lowest effect sizes, and then the second group is composed of data points 200 to 4200, and so on. Then the power and false discovery rate (FDR) are calculated within each group separately (Fig. 4). While the R/qtl scanone implementation of IM assigns a P-value for each SNP, BNM assigns a empirical likelihood r
_s
(as described in the “Methods” Section) and thus an R-value = 1?r
_s
. To compare the power of the two methods at matching thresholds, we choose a P-value for IM and look for the BNM R-value with comparable average FDR over the 76 groups (Additional file 1: Figure S2). We compare the power of our method with that of IM at 3 different P-value thresholds (P-value = 0.001, 0.03, 0.05) and their FDR-matching BNM R-values (R-value = 0.146, 0.322, 0.383). In all cases BNM has a higher power (Fig. 4
a–c). This is more prominent at higher effect sizes even though BNM has a monotonically decreasing FDR with increasing effect size (Fig. 4
d–f).

We also compare the ROC curves in 6 different effect size ranges (Fig. 5). In this case the data points are divided into 6 groups of 3800 data points each with an offset of 3100 points. BNM outperforms IM at moderate and high effect sizes. At very low effect sizes, both IM and BNM do not have much predictive power.

In all of the above results, the power and FDR calculations are based on the precise location of the true SNP. In other words, even if the method predicts a signal (at a specified threshold) in the neighborhood of the “true signal”, it is considered a false positive. This is more stringent than the usual expectation of QTL mapping [29] so we repeated the above analysis after uniformly dividing each autosome into blocks of SNPs within d Mb from each other. For example, if one or more signals are obtained in a particular block, this counts as 1 true positive (if the true signal is in this block) or 1 false positive (if the true signal is not in this block). We did the analysis at three different block sizes, d=2 (Additional file 1: Figures S3, S4 and S5), d=3 (Additional file 1: Figures S6, S7 and S8), and d=4 (Additional file 1: Figures S9, S10 and S11). BNM still shows a higher power than IM at all block sizes and (P-value, R-value) thresholds, despite the fact that adding this freedom improves the IM R/qtl results much more than it does the results from BNM. Even with the use of blocks of varying sizes, IM shows a decreasing and then increasing FDR as effect sizes are increased while the BNM FDR continues to be a monotonically decreasing function of effect size (Additional file 1: Figures S3, S6 and S9).

Next we examine how the power and false discovery rates change with the choice of different samples of phenotypes and with decreasing number of mice (Fig. 6). To examine the variation with choice of phenotype sets, we use three samples of 500 phenotypes. In two samples we randomly select 500 of the 1000 phenotypes and in the third we select the 500 phenotypes with the highest average effect sizer over the 19 autosomes. As above, we arrange the data points in order of increasing effect size and then divide them into 76 groups of 2000 data points with 100 points offset. Then the power and FDR are calculated within each group separately (Fig. 6
a, d). At low threshold values we see high variation in the average FDR between the samples (Fig. 7
a) which is due to the low number of predicted signals, making for larger statistical uncertainties. Except for the lowest threshold, this variation decreases when the analysis is repeated after setting blocks of size 2 Mb (Additional file 1: Figure S13a), 3 Mb (Additional file 1: Figure S15a), and 4 Mb (Additional file 1: Figure S17a).

The effect of population size on QTL detection has been demonstrated [30, 31], so we investigated the performance of BNM with a change in the number of mice. To examine the change and variation in FDR as we decrease the number of mice, we randomly select three samples of 600 mice and three samples of 400 mice out of the total number of 742 mice. For all of the 6 samples we used the 500 phenotypes with the highest average effect size over the 19 autosomes. Choosing the phenotypes in this manner slightly increases the fraction of high effect signals which will allow us to go to slightly higher average effect sizes in the 76 groups of data points.

As the number of mice decreases, we see more variation in the FDR between the three samples (Figs. 6
d–f), particularly for IM at low P-values (Fig. 7
e). As is to be expected, the FDR increases as the number of mice decreases for both IM (Fig. 7
e) and BNM (Fig. 7
d). For each of the 742 (all mice), 600 and 400 mice samples, we match the IM P-values to R-values of comparable average FDRs and compare the powers at P-value = 0.05 (Fig. 6). In all cases, BNM shows higher power and less variation in FDR than IM. Applying the block analysis with d=2,3,4 improves the IM FDR and removes the effect of the lower number of mice but this is not as much the case for the BNM FDR (Additional file 1: Figures S13, S15 and S17). In fact, for the block analysis, BNM’s FDR increases with smaller numbers of mice while IM’s FDR is relatively insensitive, making the matching BNM R-value much lower at the chosen IM P-value. Now when we lower the number of mice to 400, BNM shows less power (Additional file 1: Figures S12, S14 and S16). Overall, however, BNM shows better ROC curves in all cases with and without the block analysis (see Additional file 1: Figures S18, S19, S20 and S21).

Finally, we use the same 742 DO mice to map neutrophil counts in whole blood obtained from [21]. We set our R-value threshold to 0.383 since this is the FDR-matching R-value to P-value = 0.05 in our simulated data. We found signals on loci in chromosomes 1, 11, 12, 15, 16, 17, and 19 (Fig. 8). The loci we found on chromosome 1 are between 123.301336 Mb and 132.515233 Mb. This interval included Cxcr4 which is involved in neutrophil trafficking [32].