SNP selection and detection
SNPs for use in genotyping were selected on a weighted basis from the Celera Mouse SNP Database containing data from the DBA/2J, A/J, C57BL/6J, 129S1/SvImJ, and 129X1/SvJ strains. Sufficient SNPs were selected for coverage of at least one SNP per 300 kb on average. The 129S1/SvImJ and 129X1/SvJ strains were considered as the same strain when their alleles agreed; preference was given first to SNPs where each allele of the SNP was represented by two strains. This was done to favor selection of SNPs representing ancestral inheritance, not recent strain-specific mutations, and to favor real SNPs as opposed to errors in sequence annotation. Additional selective value was incorporated based on whether the SNP was in a gene, how many sequencing runs supported the presence of the SNP, and the proximity of the SNP to previously selected SNPs. Additional SNPs used to characterize the Tas1r3
locus were gathered from sequence from multiple strains available in GenBank (http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=nucleotide&cmd=search&term=tas1r3
). All physical positions presented in the paper are from the Celera Mouse Genome Assembly R13.
Primers for PCR and single-base extension were designed by using the SpectroDESIGNER software package (Sequenom, San Diego, California, United States). Assay designs are available as Supporting Information. All SNP assays were named for their position in the genome in the following format: the chromosomal location, the Mb position on the chromosome, and the kb position with a period separating each number.
For SNP genotyping, genomic DNA from pedigreed mice (Mouse DNA Resources, The Jackson Laboratory, Bar Harbor, Maine, United States) was diluted to 10 ng/μl, and 1 μl of DNA was combined with 2.45 μl of water, 0.1μl of 25 mM dNTPs (Invitrogen, Carlsbad, California, United States), 0.03μl of 5 units/μl HotStar Taq
(Qiagen, Valencia, California, United States), 0.625 μl of 10X HotStar PCR buffer containing 15 mM MgCl2
, 0.5μl PCR primers mixed together at a concentration of 1.25 μM for multiplexed reactions, and 0.325 μl of 25 mM MgCl2
. Reactions were heated at 95 °C for 15 min followed by 45 cycles at 95 °C for 20 s, 56 °C for 30 s, and 72 °C for 1 min and a final incubation at 72 °C for 3 min. After PCR amplification, remaining dNTPs were dephosphorylated by adding 1.5 μl of water, 0.17 μl of homogeneous mass extend reaction buffer (Sequenom), 0.3 units of shrimp alkaline phosphatase (Sequenom), and 0.03 μl of 10 units/μl exonuclease (USB Corporation, Cleveland, Ohio, United States). The reaction was placed at 37 °C for 20 min, and the enzyme was deactivated by incubating at 85 °C for 15 min. After shrimp alkaline phosphatase treatment, the genotyping reaction was combined with 0.76 μl of water, 0.2 μl of 10X Termination mix (Sequenom), 0.04 μl of 0.063 units/μl Thermosequenase (Sequenom), and 1μl of 10 mM extension primer. The MassEXTEND reaction was carried out at 94 °C for 2 min and then 99 cycles of 94 °C for 5 s, 52 °C for 5 s, and 72 °C for 5 s The reaction mix was desalted by adding 3 mg of a cationic resin, SpectroCLEAN (Sequenom), and resuspended in 30 μl of water. Completed genotyping reactions were spotted in nanoliter volumes onto a matrix arrayed into 384 elements on a silicon chip (Sequenom SpectroCHIP), and the allele-specific mass of the extension product was determined by matrix-assisted laser desorption ionization time-of-flight MS. Analysis of data was by automated allele calling from the SpectroTYPER software. All SNP data are available at NCBI dbSNP (http://www.ncbi.nih.gov/entrez/query.fcgi?db=snp
) and The Jackson Laboratory Mouse Phenome Database (http://www.jax.org/phenome
). Placement of the SNP data across the genome and major and minor allele distributions can be visualized using SNPview (http://snp.gnf.org
Statistical modeling for in silico mapping
The use of a single marker is restrictive in the sense that it only allows a representation of the genome as diallelic. The use of windows of multiple markers enables the visualization of more complex genomic relationships between multiple strains. This more accurately models actual haplotype patterns than does a binary approach. In determining the size of the SNP window to use as a definition of inferred haplotype for purposes of the algorithm, sizes of two, three, four, and five SNPs were examined. A window of only two SNPs was still found to be too limiting. Windows of three, four, and five SNPs produced similar results, but as window size is increased biologically meaningful patterns become fragmented, creating more single-strain inferred haplotypes, resulting in an increase in noise. Singly represented haplotypes can never be informative in this analysis because the commonality of haplotypes is required to achieve significant association with a phenotype. Three SNP windows were also analyzed across the whole genome to identify mirror loci. This would be a locus that has exactly the same strain distribution pattern across all 25 strains used in an in silico run. There were no mirror loci, or 1-off, or 2-off mirror loci (with one or two strains not grouped identically) that occurred outside of a 5-Mb interval.
Defining the genetic measure as a categorical unit necessitated the use of an ANOVA-based model. The type of ANOVA to use was determined by the characteristics of the phenotypic values.
The phenotypes studied here fell into two categories: binary or continuous. The coat color phenotype was considered as binary, where phenotypic values were set to 1 and 0. The HDL phenotype is an example of a continuous phenotype since the phenotypic values are measured on a continuous scale. Two different statistical methods were employed based on this distinction.
When phenotypic values are binary, the appropriate statistical approach involves first fitting a binomial generalized linear model to the data. An analysis of deviance table is then computed for the fitted model. The R language function glm with the parameter family set to binomial was used. This was followed by an application of anova.glm with the parameter test set to Chisq.
For continuous phenotypic values, a log transformation was applied to reduce the effects of outliers in the phenotypic data. Next, an F-statistic weighted for the genotypic diversity of the strains within an inferred haplotype group was used. The weighted F-statistic had the following form:
where ng is the number of strains in a given inferred haplotype, μg is the mean of phenotypic values in a given inferred haplotype , μT is the mean of all phenotypic values, k is the number of inferred haplotypes, N is the total number of data values, and wg is the weight representing the genetic diversity of the inferred haplotype. The genetic diversity ratio (wg) between two strains is the number of SNPs genome-wide in which both strains have genetic information and they disagree, divided by the total number of SNPs in which they both have genetic information. The genetic diversity coefficient for an inferred haplotype in the weighted F-statistic is the average wg between all strain pairs contained in the inferred haplotype.
The weighted F-statistic calculated at each SNP window determines if at least one of the inferred haplotypes has an average phenotypic value significantly different from the other inferred haplotypes. To assess the significance of the computed value, the null distribution of the weighted F-statistic was simulated at each SNP window by taking a million bootstrap samples of the phenotypic values. As in the algorithm used for binary phenotypes, inferred haplotype patterns present in only one strain were not included in the calculation because they are not informative in elucidating shared ancestral blocks. From this distribution of a million random F-statistics, 200 bootstrap samples of size 1 million were computed. For each bootstrap sample, a p
-value was computed by dividing the number of random F-statistics larger than the true F-statistic by the total number of random F-statistics (million). In this way 200 p
-values were collected. The vertical heights reported in the bar graphs (see Figure S2
) are the −log(p
) transform of the median of these 200 p
-values. A 95% confidence interval (CI) for the p
-value at this window was also calculated from this bootstrap distribution.
To estimate the overall false positive rate for this type of calculation, calculating a significance threshold based on the family-wise error rate (FWER) has been proposed (Churchill and Doerge 1994
). Others have noted that the traditional FWER calculation is too strict in the context of multiple testing and leads to a significant loss of power (Lander and Kruglyak 1995
). Therefore, we employed a recently developed method of bootstrap estimation of common cutoffs based on the gFWER (Dudoit et al. 2004
). Whereas the FWER method reports significance, using the most conservative criterion of only one false positive, the gFWER method controls for multiple testing while allowing for an acceptable false positive rate (in our case, α < 0.005).
The gFWER method to control for false positives as applied to in silico mapping is briefly described as follows. A null reference distribution was constructed using random bootstrap tests to determine a significance cutoff. Ten thousand bootstrap samples of phenotype values were randomly assigned to the true haplotype structure. For each random bootstrap sample, the nonparametric ANOVA approach outlined above was performed at each three-SNP window, with one difference. Whereas the initial true calculation reports the median of 200 bootstrap p
-values, the gFWER method requires an estimate of the “supremum” (least upper bound) of expected values reported at each locus, so the most significant value is reported from the 200 bootstrap p
-values (following Procedure 3 in Dudoit et al. 2004
), ensuring a conservative false positive estimate. For each bootstrap sample, the genome-wide −log(p
-value) corresponding to the (1 − α) percentile was added to the null distribution (as described in Procedure 5, Dudoit et al. 2004
). Finally, after the 10,000 bootstraps are complete, the significance threshold is set as the (1 −α) percentile in the entire null reference distribution (computed from our 10,000 randomly bootstrapped iterations). While this threshold still represents a conservative estimate of the desired false positive rate, the gFWER has significantly more power than the traditional FWER calculation.
Using this method for calculation of false positives, it is not necessary to specify the marginal distribution of the test statistic at each window of SNPs. Estimations of false positives or power that assume some parametric form of test statistic's distribution are not reliable in this context. This distribution can alter radically at each SNP window. In this context, the statistical problem of calculating quantities like discovery power (that is, ultimately the type I and type II error) is further complicated. Nearly 11,000 hypothesis tests (one at each three-SNP window) are conducted in a single run of the algorithm. Therefore, equations that currently exist for the estimation of power for QTL mapping by traditional methods cannot be applied here because they assume that the test statistic has some previously defined parametric form. Code for the above described algorithms is available upon request.
For calculation of the significance of the number of in silico QTL that overlapped with previously identified QTL for the HDL phenotype, a binomial distribution was used
given p = probability of success of 0.42 (overlap with HDL QTL in previous literature).
Therefore q = probability of failure; the 0.0025 result is the probability of at least nine successes in ten trials. Only ten loci could be assessed for this result as no information is available for traditional HDL QTL present on the X chromosome.
For the mapping of the retinal degeneration traits, 37 strains were used. This represented all of the strains for which information existed in The Jackson Laboratory database minus the most divergent wild-derived strains for which inference of haplotype would be expected to be most inaccurate. These strains were A/J, AKR/J, BALB/cByJ, BUB/BnJ, C3H/HeJ, C57BL/10J, C57BL/6J, C57BLKS/J, C57BR/cdJ, C57 l/J, C58/J, CBA/J, CE/J, DBA/1J, DBA/2J, FVB/NJ, I/LnJ, KK/HlJ, LG/J, LP/J, MA/MyJ, NOD/LtJ, NON/LtJ, NZB/BlNJ, NZW/LacJ, PERA/EiJ, PL/J, RIIIS/J, SEA/GnJ, SJL/J, SM/J, ST/bJ, SWR/J, WSB/EiJ, ZALENDE/EiJ, 129S1/SvImJ, and 129X1/SvJ. Because of the added complexity of the coat color traits, mapping was restricted to the 25 most related strains for which coat color phenotype could clearly be determined. For the albino analysis, 129S1/SvImJ, A/J, AKR/J, BALB/cByJ, C3H/HeJ, C57BL/10J, C57BL/6J, C57BLKS/J, C57BR/cdJ, C57 l/J, C58/J, CBA/J, DBA/1J, DBA/2J, I/LnJ, LP/J, MA/MyJ, NZB/BlNJ, NZW/LacJ, PERA/EiJ, PL/J, SEA/GnJ, SM/J, WSB/EiJ, and ZALENDE/EiJ strains were used. For the nonagouti mapping, the same strain set as the albino mapping was used except for the mice presenting the albino phenotype. The strains were 129S1/SvImJ, C3H/HeJ, C57BL/10J, C57BL/6J, C57BLKS/J, C57BR/cdJ, C57 l/J, C58/J, CBA/J, DBA/1J, DBA/2J, I/LnJ, LP/J, NZB/BlNJ, PERA/EiJ, SEA/GnJ, SM/J, WSB/EiJ, and ZALENDE/EiJ. Any mouse showing an agouti coat color was considered to be agouti for this analysis regardless of genotype at the agouti locus. Only limited phenotype data were available for saccharin preference, so again all strains with available data except the most divergent wild-derived strains for which inference of haplotype would be expected to be most inaccurate were used. These strains were A/J, AKR/J, BALB/cByJ, BUB/BnJ, C3H/HeJ, C57BL/6J, C57 L/J, CBA/J, CE/J, DBA/2J, FVB/NJ, I/LnJ, KK/HlJ, LP/J, NOD/LtJ, NZB/BlNJ, PL/J, RIIIS/J, SEA/GnJ, SJL/J, SM/J, ST/bJ, and SWR/J. For the mapping of the other complex traits, only the 25 strains with the closest ancestral relationship were used. These strains were 129S1/SvImJ, A/J, AKR/J, BALB/cByJ, BTBR T+ tf/J, C3H/HeJ, C57BL/10J, C57BL/6J, C57BLKS/J, C57BR/cdJ, C57 l/J, C58/J, CBA/J, DBA/1J, DBA/2J, I/LnJ, LP/J, MA/MyJ, NZB/BlNJ, NZW/LacJ, PERA/EiJ, PL/J, SEA/GnJ, SM/J, and WSB/EiJ.