BACKGROUND
the availability of high density genotypes of single nucleotide polymorphism  markers for plants and livestock species, in conjunction with phenotypic data for complex traits, allows the calculation of: 1) estimates of genomic breeding values   <cit>  for genomic evaluation  <cit> , and 2) estimates of the effects of genomic regions associated with the genetic variability in genome wide association studies   <cit> .
there is an increasing number of gwas data sets analyzed by mixed models and multiple testing procedures  <cit> , after fitting all individual effects of genomic regions into the model  <cit> .
the model may be difficult to fit when both, the number of individuals and snp effects, are large.
we propose to use a linear transformation of genomic breeding values to estimate the marker effects from a simpler equivalent mixed model, and then testing those effects using a standardized test statistic that employs the variance  of the same effects.
the method of genomic selection proposed by meuwissen et al.
<cit>  to estimate gebvs starts by fitting the snp effects to a given data set.
next is to estimate gebv of any individual using its genotype , by adding across the entire genome those solutions corresponding to the individual's snp.
the mixed model employed conveys vectors of fixed effects, and random effects of markers or snps  assumed to be normally distributed with null mean and a covariance matrix proportional to the identity matrix times the variance of snp effects .
errors are assumed to be gaussian, independent and identically distributed with null mean and covariance matrix .
an equivalent mixed model discussed by garrick  <cit>  and stranden  <cit>  is fitted after the linear transformation a = z g where a is a random vector of breeding values, and z an incidence matrix that relates elements in a to those in g. each column of z is associated with a given snp and the elements are standardized by functions of snp allele frequencies and by the total number of snp.
it is worth noting that the same z is used in our implementation of the model of meuwissen et al.
<cit>  to relate the vector of marker effects in g to the data phenotypes.
moreover, gebvs in the equivalent model have variance-covariance matrix .
the procedure requires that the variances are equal, i.e.
.
once the equivalent model is fit, snp effects are calculated by the transformation g = z'g− 1a, and individual snp effects in g are divided by the square root of its variance ) to get the so called snpej test statistics.
we also provide a formula to calculate var without having to fit the model with snp effects.
the next step is to select genome segments that may be highly associated with the genetic variability of the trait for each chromosome.
in doing so, we look for the snp having the highest value of minus the logarithm of the p-value throughout the chromosome.
once the snp is located, a segment of one mb to the left and one to the right is defined, and a relationship matrix is calculated using only the information from those markers.
the relationship matrix is used as the proportional variance-covariance of the segment effects in a model that also includes fixed effects and random gebvs.
in a final step, the likelihood ratio is calculated to test the significance of the largest effect segment of each chromosome by comparing against a reduced model with fixed effects and gebvs.
the critical value  is adjusted by the bonferroni correction.
the algorithm not only delivers genome wide associations and genomic predictions efficiently, but it also minimizes computing time and memory requirements.
moreover, the specific variance of the snp effects is used in calculating the test, thus taking into account the amount of information of any given marker.
instead, other testing approaches rely on a prior variance or a constant estimate of the additive variance.
methods
dataset
the experimental population was raised at the michigan state university swine teaching and research farm, east lansing, mi  <cit> .
parents from the initial generation  were four duroc boars mated to  <dig> pietrain sows by artificial insemination.
from all resulting f <dig> animals,  <dig> females and  <dig> males  were selected as parents for the f <dig> generation, by avoiding full or half sib matings.
a total of  <dig>  f <dig> piglets were born alive from  <dig> litters out of  <dig> farrowing groups.
phenotypic data for growth, carcass merit and meat quality traits were collected for approximately 950 f <dig> pigs .
data used for the study were measures of the growth trait 13 week tenth rib backfat  .
the trait was chosen as it displays a sizable heritability  and a normal distribution.
animal protocols were approved by the michigan state university all university committee on animal use and care .
genotyping and data editing
dna was isolated from white blood cells using standard procedures as previously described for this population  <cit> .
quantity and quality of dna samples were determined using a qubit fluorometer .
the experimental population was genotyped with two marker snp panels.
1)  <dig> animals were genotyped  with a commercial panel, the illumina porcinesnp <dig> beadchip   <cit>  and 2) 612 f <dig> animals were genotyped with a second panel composed of a 9 k tagsnp set referred to as the geneseek genomic profiler for porcine ld   <cit>  .
a set of  <dig>  snp out of m =  <dig> , were eliminated from all analyses as their physical positions were unknown.
mendelian inconsistencies  were taken as missing genotypes, and  <dig> animals  with more than 10% of snp missing were not used for any analysis.
by similar considerations,  <dig>  snp were removed from the analyses as they had more than 10% missing data.
additionally,  <dig>  snp were excluded as their minor allele frequency  was below  <dig> .
this editing procedure followed that of badke et al.
<cit>  and gualdrón et al.
<cit> , and the program plinkv <dig>   <cit>  was used for the task.
f <dig> animals genotyped with the 9 k panel were imputed to 60 k following procedures discussed by gualdrón et.
al  <cit> , by means of the software alphaimpute  <cit> , resulting in imputation accuracy of around  <dig>   <cit> .
genotypes imputed in the f <dig> had a second editing procedure by maf <  <dig> , which excluded  <dig> virtually monomorphic snp.
the editing policies and genotype imputation resulted in a data set with records from  <dig> pigs  having  <dig>  snp per animal.
estimation of genomic relationship matrix
the genomic relationship matrix was estimated from observed and imputed high density  snp genotypes.
genotypes were expressed as allelic dosage  <cit> , such that genotypes were entered into a marker matrix m of dimension , where n is the number of animals and m the number of snp, having elements in the interval  <cit> , i.e.
the count of the allele used as reference.
in the sequel, we will use the sub index i to refer to the individual.
matrix m was standardized to matrix z that has generic elements equal to
elements of z are then calculated by subtracting twice the frequency of the reference allele at the jth marker , to the corresponding element of m
<cit> , and then dividing the resulting difference by the square root of the expected variance 2pj of each element in the column multiplied by the number of columns  in m. the allele frequency pj was calculated from the f <dig> generation .
the genomic relationship matrix was finally calculated as:
<dig>
prediction model
using the genomic relationship matrix from equation , the centered animal model for genomic evaluation can be written as:
2a
where y is the phenotypic vector containing the data on 13-week tenth rib backfat , x is the incidence matrix that relates records to the fixed effects of sex in β, vector a contains the random breeding values such that , e is the random error vector such that , and i is the identity matrix.
variance components were estimated with reml using the regress version  <dig> - <dig> r package  <cit> .
following stranden et al.
<cit>  an equivalent model to  is
2b
every element in  is defined as before except for the vector g of snp effects.
to show that  and  are equivalent models, we employ the fact that a = z g. then, the variances of a and g are related in the following manner:
necessary conditions for models  and  to be equivalent  are that g = z z ' and .
variance of snp effects
in this section, we describe the algorithm to calculate the variance of the estimated snp effects g. the snp effects were obtained from a linear transformation of breeding values in
<cit> , as follows:
<dig>
the last step results from the fact that model equivalence involves .
now, from equation   is obtained as follows:
<dig>
now, we know that the predictor error variance  of  from model  is equal to:
so that
matrix caa results from inverting the coefficient matrix of the mixed model equations  <cit>  such that:
then, on replacing with the latter expression into  ), we have:
<dig>
expression  results in a large matrix of dimension  with m the number of snp.
however, we only need its diagonal elements.
also notice that the first term in , z ' g− 1 z, can be computed and stored to be reused for the different traits, whereas caa has to be computed for each trait.
standardization of snp effects
the estimated snp effects in  were standardized by dividing with their corresponding  obtained from  as follows:
<dig>
p-values and genome screening
the p-values were assessed as  <dig> minus the cumulative probability density of the absolute value of snpe j, a number that was then multiplied by  <dig> so as to obtain:
where Φ is the cumulative density function of the normal distribution for the random variable x. when analyzing the trait 13 week tenth rib backfat , the p-values for each snp were plotted across the genome as –log <dig>  using the physical position of the snp in mega-bases .
standardization of snp effects using the pev of the marker
a second standardization of the jth snp effect  was performed using the pev  as follows:
<dig>
as discussed above, .
the p-values and genome screening for snpep j were assessed and plotted in the same fashion as for snpe j.
simulation
a plasmode simulation was performed to compare how the standardized values snpe j and snpep j affected the nominal size of the test for the effect to be equal to zero.
data on  <dig> animals with  <dig>  snp each were used for the study, and the  <dig> snp on chromosome  <dig> were reshuffled.
two scenarios were considered: 1) dependency: rows of the genotype matrix were permuted for columns corresponding to snp on chromosome  <dig>  thus keeping linkage disequilibrium  within chromosomes but breaking the relationship between genotypes and phenotypes for the  <dig> snp on the chromosome.
2) independency: the genotype of any animal was permuted independently by marker  for those snp on chromosome  <dig>  and the relationship with the phenotype was broken too.
for both scenarios model  was fitted to the data, and two tests were calculated for each scenario: test1 = snpej and test2 = snpepj.
permutations were repeated  <dig> times per scenario, and in each permutation the g matrix was calculated while fitting model .
as a result, the heritability of the trait was similar to the original heritability due to relationships in the other  <dig> chromosomes being kept intact, and p-values for those snp  on chromosome  <dig> were obtained for the different tests.
under the null hypothesis and assuming independence , an approach that controls for type i error appropriately  <cit> , the  <dig> test p-values follow a uniform distribution.
consequently, to estimate the empirical quantiles of the distribution for the null hypothesis, we used a uniform density u ∼  to generate  <dig> replicated sets for the  <dig> p-values.
snp effects and tests obtained by a single marker model
the snp effects were tested on a one by one basis.
the model approach used for testing purposes is better known as “efficient mixed-model association”   <cit> .
the model included fixed effects of sex and one-marker-at-a-time; random variable was the animal effect with variance-covariance equal to the genomic relationship matrix using all markers, which was calculated as described before.
the r package rrblup  <cit>  was used for fitting the different models and for calculating the tests and p-values.
proportion of variance explained by segments with large effect
after the genome screen using model 2a, the snp with the smallest p-values were selected to form snp segments.
these segments were defined by taking all snp within one mb upstream and one mb downstream of the snp with smallest p-value on each chromosome.
the size of the segment was chosen using a criterion similar to the one employed by hayes et al.
<cit> .
the point of change in the rate of decay in linkage disequilibrium in this population was about r2 =  <dig>  at 1 mb , which essentially would imply a minimal contribution to the additive variance from markers located beyond such distance.
moreover, segment sizes about two mb have been reported to be significant in association studies .
the proportion of variance associated with each segment was estimated by building a genomic relationship matrix g <dig> ) using all snps that belonged to the segment, whereas genomic relationship matrix g <dig> was built using all remaining snps.
the model fitted can be represented as:
<dig>
where a <dig> is the vector of additive random effects associated with those snp located in the segment, such that , and a <dig> is the vector of additive random effects associated with all snps except those involved with a <dig>  such that .
model  assesses the proportion of variance explained by the segment of interest  from the genome variance explained by all snps .
the variances estimated in  were compared with those estimates from model .
hayes et al.
<cit>  used a similar model to assess the segment variance.
applying either model , or the approach of hayes et al.
<cit>  gave similar estimated variance components.
in practice, the advantage of fitting model  is that g <dig> is computed by subtracting from g the columns of z related to the segment being tested.
let zs be a matrix having as columns those related to the segment being tested, then .
on the contrary, in the model of hayes et al.
<cit>
gis different from segment to segment.
additionally, the calculation of g <dig> and  is fast and involves only those snps located in the segment.
to adjust the level of significance for multiple comparisons, a bonferroni correction  was performed.
in this context, if the pig genome is ~2800 mb long and the average size of the segment is 2 mb, there are  <dig> segments along the genome with corresponding multiple tests.
thus, for α =  <dig> , the bc was equal to  <dig> /1400 =  <dig> e− <dig> .
hence, in order to evaluate the significance of the segments, a second p-value for the likelihood ratio test  was calculated to compare against bc.
this p − valuelrt was assessed as  <dig> minus the distribution function of a chi-square  random variable with  <dig>  degrees of freedom  <cit>  as follows:
where Ω is the distribution function of a random variable having the χ <dig> as density, and lrt is the likelihood ratio test obtained by contrasting appropriate models.
RESULTS
genome screening
the p-values of the  <dig> snp were obtained as described in the methods section.
first, the p-values for snpej, i.e.
using , were plotted along the genome  to identify genomic positions that are associated with variation in 13-week tenth rib backfat .
large peaks  above  <dig> can be seen at chromosomes  <dig> and  <dig>  suggesting noticeable genetic variation for the trait.
on the other hand, p-values for snpepj  were very large, with a maximum − log <dig> of  <dig> .
in essence, the pattern observed in figure  <dig> is the result of dividing the non-standardized snp effects by a constant.
specifically, the normalizing value was , with var  =  <dig> .
the use of the square root of the difference between those two values resulted in a practically constant denominator for the test-statistic that was equal to  <dig> .
also, a look at figure  <dig> suggests signals at chromosomes  <dig>   <dig>   <dig>  and  <dig>  a fact that is not observed in figure  <dig>  however, this might be an artefact of the constant denominator that tends to overestimate the true variability for some snp, thus resulting in corresponding false positives across the genome.figure  <dig>
manhattan plot for trait 13-week tenth rib backfat  by standardization snp
ej
.
genome screening for  <dig> snp using standardization .
−log10   versus the absolute snp position in mb .
the red line represents a genome-wide significance threshold .
numbers from  <dig> to  <dig> represent the chromosome id.
manhattan plot for trait 13-week tenth rib backfat  by standardization snp
epj
.
genome screening for  <dig> snp using standardization .
−log10   versus the absolute snp position in mb .
numbers from  <dig> to  <dig> represent the chromosome id.
in order to study the type i error rate of the two proposed tests we performed a plasmode simulation  <cit> .
a plasmode is a dataset created from real data where some of the truth is known.
in brief, our plasmode is a simulation that uses reshuffling in a portion of the data as explained in the methods section.
we performed a simulation assuming independent snp, and another one keeping the dependency between snp  intact.
simulation results were plotted into a quantil-quantil plot graph  using the number –log for each case of standardization.
first, the p-values for test <dig>  obtained in the scenario under independent snps  displayed an identical distribution of p-values when obtained by the reference distribution u ∼ .
in contrast, under dependency  less extreme p-values were observed, a fact that was not reflected in a uniform distribution.
this is a well known fact in human genetic epidemiology  <cit> , where the implementation of the bonferroni correction of p-values from associated snp under the assumption of independence results in tests that are too conservative.
on the other hand, for test  <dig>  even p-values obtained for independent snp  displayed a distribution that was too conservative.
furthermore, the results from the dependent scenario  were even more conservative than those from the independent scenario , thus indicating that the use of the square root of  as the denominator of the test-statistic results in a more powerful and not too liberal choice when compared to the use of the square root of pev = .figure  <dig>
quantil-quantil plot of the observed and expected –log obtained by simulation.
reference distribution was an independent and uniform distribution u ∼  for  <dig> p-values simulated .
test1 = under dependent  and standardization by  .
test1 = under independent  and standardization by  .
test2 = under independent  and standardization by pev .
each scenario has  <dig> p-values permuted  <dig> times.
bands represent confidence intervals of 95% , green band = test <dig>  pink band = test <dig>
snp effects and tests obtained by the marker model
the analyses of one snp tested at a time using the emma procedure  <cit>  resulted in p-values that were almost identical  to those of snpej .
the time taken to compute  <dig> snp tests one at a time was 84 minutes.
in comparison, the algorithm used to fit model  and to perform the tests of standardized effects took a total time of 29 minutes .
this time includes the computation of the g matrix, the fit of the animal model, the back transformation to calculate the snp effects, and the calculation of the standard errors that are needed to compute the test-statistics.
tests of segment effects
we also compared the results from our proposed method to those obtained with a segment-based likelihood ratio test that has been used by animal breeders  <cit> .
due to computational demand, we only performed the lrt to test for segment effects.
thus, the snp with the smallest p-values ) on each chromosome were chosen, whereas no segments were tested using lrt for regions with snpepj resulting in exceedingly low p-values.
the three segments from chromosomes with the smallest p-values are displayed in table  <dig>  and the remaining segments from the  <dig> other chromosomes are shown in the additional files .
all segments measured 2 mb .
the estimates of the variance components and the loglikelihood obtained from model equation  were compared with those from model equation .
these results are displayed in table  <dig> table  <dig>
snp selected by smallest p-value per chromosome
snp name = snp marker name, position mb = marker physical position in mega-bases, −log10 = −logarithm in base  <dig> of the smallest p-value,  = absolute value of the snp effect estimated for the trait 13 week tenth rib backfat .
variance components and loglikelihood for models with or without the segment
snp − log
10
lk_m1
lk_m2
lrt
p-value
lrt
vare_m1
vara_m1
vare_m2
vara_m2
segmva
%segmva
seg-chromosome = number of chromosome where segment is located, m1 = model without the segment: y = xβ + a + e, m2 = model  with the segment y = x
β + a
1 + a
2 + e, snp − log
10
= −logarithm in base  <dig> of the snp p-value selected to create a segment, lk_m1 = −loglikelihood for m <dig>  lk_m2 = −loglikelihood for m <dig>  lrt = likelihood ratio test for m <dig> and m <dig>  p-value
lrt = p-value for lrt, vare_m1 = error variance  of m <dig>  vara_m1 = additive variance  of m <dig>  vare_m2 = error variance  of m <dig>  vara_m2 = additive variance  of m <dig>  segmva = additive variance segment  of m <dig>  %segmva = proportion in% of the total variance explained by the segment.
results from the lrt indicated that the segment on chromosome  <dig> was significant: p − valuelrt ‒ 6 =  <dig> e− <dig>  a number smaller than the critical  <dig>  bonferroni threshold for  <dig> segments .
on the contrary, the segments located on all other chromosomes were not significant.
the proportion of variance explained by the segment from chromosome  <dig>  =  <dig> ) was 11% of the total variance, a fact that was reflected in a similar reduction of the estimated additive variance  in model :  <dig>  +  <dig>  =  <dig> .
this latter value is close to  <dig> , i.e.
the estimated value of  from model  .
for all other chromosomal segments, the estimated value of  did not decrease to a significant amount.
discussion
the main goal of this research was to develop a novel procedure to perform a rapid genome scan, or gwas analysis, from a genomic evaluation.
moreover, the sufficient statistics of our methodology are: the best linear unbiased prediction  of the breeding values from an animal model, g as the covariance matrix , z as the standardized marker effects matrix, variance components, and caa.
this setting makes the implementation extremely feasible after the genomic evaluation has been performed as discussed by legarra et al.
<cit> .
variance of the snp effect
first, the snp effects  were calculated by a linear transformation of  using expression .
then, we calculated  using an expression derived from mixed model theory ).
next, we divided  by the square root of  to standardize the effect, and referred the statistics as snpej.
the p-values for the tests of specific genome regions were calculated with a level of significance − log10 =  <dig>  additionally, prediction error variance () was employed for a second standardization, and it was called the snpepj statistic.
after the analyses, we obtained higher p-values  =  <dig> ) and detected stronger signals  for snpepj than with snpej.
furthermore, a simulation was carried out with the same structure of snps markers and animal data as in the current study, in order to compare the performance of empirical p-values of both standardized tests.
the snps markers of chromosome  <dig> were reshuffled, and two scenarios were simulated: 1) dependent genotypes , and 2) independent genotypes .
neither scenario displayed a relationship with the phenotype, whereas both standardized tests were calculated at each scenario.
the reference distribution for the p-values considered was the uniform.
in the independent scenario , standardization with  gave an empirical distribution of p-values that resembled the uniform density, but in the dependent scenario  the snpej performed conservatively.
instead, the standardization with  produced conservative results in the independent scenario , and very conservative tests in the dependent scenario .
in this context, standardizing snp effects with  resulted in p-values that were closer to the simulated ones.
moreover, the performance of snpej under ld was not too conservative, a scenario that could be extrapolated to the genotypes in the current study.
in addition, the p-values calculated using the emma procedure  <cit>  were similar to those obtained with snpej.
these results suggest that snpej behaves reasonably to control type i error rate or false positives.
also, the computing time for fitting model  and then calculating  using expressions - was  <dig>  to  <dig> times less than the computing time for the emma model.
in order to identify snp with important phenotypic associations  <cit> , the calculation of snp effects  from genomic breeding values
<cit>  has been used in several studies  <cit> .
in this context, the variance of snp effects has been estimated using different approaches.
wang et al.
<cit>  employed the classical definition of the variance of additive effects from quantitative genetics  <cit> , so that the variance for each jth marker was obtained as follows: .
whereas, mcclure et al.
<cit>  proposed equating the variance of snp effects to , and then normalizing the snp effects with the square root of this estimated and constant variance.
this test performed similar to snpep j , when the estimated snp effects  was divided by a constant denominator, a value almost equal to the prior variance  <dig> , and resulted in a very conservative test.
in contrast, the advantage of the standardized test  presented here was that each snp effect was scaled by its own  standard deviation rather than the use of a prior variance  <cit>  or by the square of each specific snp effect
<cit>  as variance.
furthermore, the computation of snpej, involves the same variance for the same snps markers and animals, i.e.
, and the use of the standardized incidence matrix z, a function of 2pj, takes into account this latter quantity into snpej.
additionally, the matrix z uses the allele frequencies from the f <dig> generation calculated with unrelated individuals, and a proper expected variance by marker .
in addition, the test statistics snpej that standardizes snp effects produces a p-value, a result that is appealing to many researchers that are more familiar with the method of testing one snp at the time rather than with the proportion of additive variance that is explained by a genomic region.
a further advantage of the method is that detection of many false positives are avoided, and genome positions with sizeable effects are highlighted.
candidate segment approach
later in the research, genome segments that expressed higher signals were located.
to this purpose, snps with the smallest p-values from snpej  were selected, and for each of these snp a segment of 2 mb long  was created.
the next step was to estimate the variance components and the log-likelihood from the centered animal models  and .
the latter model includes the random vector of snp segments a <dig>  lastly, we compare the performance of both models.
hayes et al.
<cit>  used a similar model to , although the random snp effect was taken from the breeding value and fitted as a separate segment effect.
we observed similar results from the use of either approach.
the advantage of fitting model  is that matrix g is the same for all segments, so that it was calculated only once, and stored in memory for the calculations, whereas in the model of hayes et al.
<cit>  a different g has to be calculated for each segment.
this implies an extended computing time and higher requirements of cpu memory to obtain similar results to those from model .
to evaluate the significance of the segments, the effects of each chromosome segment were tested by the likelihood ratio test.
the size of the test was adjusted by the bonferroni correction.
as a result, the segment located on chromosome  <dig>  was significant, and explained 11% of the trait total variance.
previous studies by edwards et al.
<cit>  and choi et al.
<cit> , using microsatellites and a small number of snp, found significant regions  on chromosome  <dig> for 13 week tenth rib backfat in the current population under study.
additionally, forty eight markers between the physical position between 128 mb and 139 mb on chromosome  <dig> , have been reported to be associated with the trait.
furthermore, recent studies showed the importance of chromosome  <dig>  <cit>  in the expression of the trait.
therefore, our results confirm the presence of genetic variability in the trait from chromosome  <dig>
CONCLUSIONS
fast genome screening of snp effects linearly transformed from genomic breeding values is advantageous, as a by-product of genomic evaluations for different species of farm animals.
moreover, the standardized tests of snp effects using their own variance  developed in this study helps in detecting specific genomic regions involved in the additive variation of the trait and reducing false positive locations using less computing time.
additionally, genome segments of about 2 mb formed by surrounding the snp with the smallest p-values on each chromosome, and tested with a standardized test involving  and with the bonferroni correction, could detect genome regions responsible for sizeable fractions of the trait genetic variance.
this methodology involving genome scan and candidate segment approach is a useful method for meta-analyses of genome-wide association studies, as it enables the detection of specific genome regions that affect an economically relevant trait when using multiple populations.
code and data to obtain and reproduce the results presented is publicly available at https://www.msu.edu/~steibelj/jp_files/gblup.html.
electronic supplementary material
additional file 1:
highest − log
10
on each chromosome for trait 13-week tenth rib backfat  by standardization snp
ej
and emma.
the blue and red circle represents highest − log <dig> on each chromosome by the standardization snp
ej and efficient mixed-model association  using rrblup.
respectively.
additional file 2:
dispersion plot of − log
10
for trait 13-week tenth rib backfat  by emma and standardization snp
ej
.
dispersion plot for  <dig> –log <dig>  by efficient mixed-model association  using the rrblup r package in the x axis, and by the standardization snp
ej in the y axis.
red straight line is the reference line 0– <dig>
additional file 3:
variance components and loglikehood for models with or without the segment for all chromosomes.
.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
jps, rjcc, jlgd: performed and supervised statistical and simulation analyses and wrote the manuscript.
rob, cwe: designed the resource population and led collection of phenotypic data.
cwe, ner: performed dna extraction and coordinated genotyping with commercial laboratory.
jps, rob, cwe: designed high density genotyping scheme.
all authors read and approved the paper.

 the availability of high density genotypes of single nucleotide polymorphism markers for plants and livestock species, in conjunction with phenotypic data for complex traits, allows the calculation of: 1) estimates of genomic breeding values for genomic evaluation, and 2) estimates of the effects of genomic regions associated with the genetic variability in genome wide association studies. there is an increasing number of gwas data sets analyzed by mixed models and multiple testing procedures, after fitting all individual effects of genomic regions into the model. the model may be difficult to fit when both, the number of individuals and snp effects, are large. we propose to use a linear transformation of genomic breeding values to estimate the marker effects from a simpler equivalent mixed model, and then testing those effects using a standardized test statistic that employs the variance of the same effects. the method of genomic selection proposed by meuwissen et al. to estimate gebvs starts by fitting the snp effects to a given data set. next is to estimate gebv of any individual using its genotype, by adding across the entire genome those solutions corresponding to the individual's snp. the mixed model employed conveys vectors of fixed effects, and random effects of markers or snps assumed to be normally distributed with null mean and a covariance matrix proportional to the identity matrix times the variance of snp effects. errors are assumed to be gaussian, independent and identically distributed with null mean and covariance matrix. an equivalent mixed model discussed by garrick and stranden is fitted after the linear transformation a = z g where a is a random vector of breeding values, and z an incidence matrix that relates elements in a to those in g. each column of z is associated with a given snp and the elements are standardized by functions of snp allele frequencies and by the total number of snp. it is worth noting that the same z is used in our implementation of the model of meuwissen et al. to relate the vector of marker effects in g to the data phenotypes. moreover, gebvs in the equivalent model have variance-covariance matrix. the procedure requires that the variances are equal, i.e.. once the equivalent model is fit, snp effects are calculated by the transformation g = z'g− 1a, and individual snp effects in g are divided by the square root of its variance ) to get the so called snpej test statistics. we also provide a formula to calculate var without having to fit the model with snp effects. the next step is to select genome segments that may be highly associated with the genetic variability of the trait for each chromosome. in doing so, we look for the snp having the highest value of minus the logarithm of the p-value throughout the chromosome. once the snp is located, a segment of one mb to the left and one to the right is defined, and a relationship matrix is calculated using only the information from those markers. the relationship matrix is used as the proportional variance-covariance of the segment effects in a model that also includes fixed effects and random gebvs. in a final step, the likelihood ratio is calculated to test the significance of the largest effect segment of each chromosome by comparing against a reduced model with fixed effects and gebvs. the critical value is adjusted by the bonferroni correction. the algorithm not only delivers genome wide associations and genomic predictions efficiently, but it also minimizes computing time and memory requirements. moreover, the specific variance of the snp effects is used in calculating the test, thus taking into account the amount of information of any given marker. instead, other testing approaches rely on a prior variance or a constant estimate of the additive variance. methods dataset the experimental population was raised at the michigan state university swine teaching and research farm, east lansing, mi. parents from the initial generation were four duroc boars mated to pietrain sows by artificial insemination. from all resulting f animals, females and males were selected as parents for the f generation, by avoiding full or half sib matings. a total of  f piglets were born alive from litters out of farrowing groups. phenotypic data for growth, carcass merit and meat quality traits were collected for approximately 950 f pigs. data used for the study were measures of the growth trait 13 week tenth rib backfat. the trait was chosen as it displays a sizable heritability and a normal distribution. animal protocols were approved by the michigan state university all university committee on animal use and care. genotyping and data editing dna was isolated from white blood cells using standard procedures as previously described for this population. quantity and quality of dna samples were determined using a qubit fluorometer. the experimental population was genotyped with two marker snp panels. 1) animals were genotyped with a commercial panel, the illumina porcinesnp beadchip and 2) 612 f animals were genotyped with a second panel composed of a 9 k tagsnp set referred to as the geneseek genomic profiler for porcine ld. a set of snp out of m = , were eliminated from all analyses as their physical positions were unknown. mendelian inconsistencies were taken as missing genotypes, and animals with more than 10% of snp missing were not used for any analysis. by similar considerations, snp were removed from the analyses as they had more than 10% missing data. additionally, snp were excluded as their minor allele frequency was below. this editing procedure followed that of badke et al. and gualdrón et al., and the program plinkv was used for the task. f animals genotyped with the 9 k panel were imputed to 60 k following procedures discussed by gualdrón et. al, by means of the software alphaimpute, resulting in imputation accuracy of around. genotypes imputed in the f had a second editing procedure by maf < , which excluded virtually monomorphic snp. the editing policies and genotype imputation resulted in a data set with records from pigs having snp per animal. estimation of genomic relationship matrix the genomic relationship matrix was estimated from observed and imputed high density snp genotypes. genotypes were expressed as allelic dosage, such that genotypes were entered into a marker matrix m of dimension, where n is the number of animals and m the number of snp, having elements in the interval, i.e. the count of the allele used as reference. in the sequel, we will use the sub index i to refer to the individual. matrix m was standardized to matrix z that has generic elements equal to elements of z are then calculated by subtracting twice the frequency of the reference allele at the jth marker, to the corresponding element of m, and then dividing the resulting difference by the square root of the expected variance 2pj of each element in the column multiplied by the number of columns in m. the allele frequency pj was calculated from the f generation. the genomic relationship matrix was finally calculated as: prediction model using the genomic relationship matrix from equation, the centered animal model for genomic evaluation can be written as: 2a where y is the phenotypic vector containing the data on 13-week tenth rib backfat, x is the incidence matrix that relates records to the fixed effects of sex in β, vector a contains the random breeding values such that, e is the random error vector such that, and i is the identity matrix. variance components were estimated with reml using the regress version - r package. following stranden et al. an equivalent model to is 2b every element in is defined as before except for the vector g of snp effects. to show that and are equivalent models, we employ the fact that a = z g. then, the variances of a and g are related in the following manner: necessary conditions for models and to be equivalent are that g = z z' and. variance of snp effects in this section, we describe the algorithm to calculate the variance of the estimated snp effects g. the snp effects were obtained from a linear transformation of breeding values in, as follows: the last step results from the fact that model equivalence involves. now, from equation is obtained as follows: now, we know that the predictor error variance of from model is equal to: so that matrix caa results from inverting the coefficient matrix of the mixed model equations such that: then, on replacing with the latter expression into ), we have: expression results in a large matrix of dimension with m the number of snp. however, we only need its diagonal elements. also notice that the first term in, z' g− 1 z, can be computed and stored to be reused for the different traits, whereas caa has to be computed for each trait. standardization of snp effects the estimated snp effects in were standardized by dividing with their corresponding obtained from as follows: p-values and genome screening the p-values were assessed as minus the cumulative probability density of the absolute value of snpe j, a number that was then multiplied by so as to obtain: where Φ is the cumulative density function of the normal distribution for the random variable x. when analyzing the trait 13 week tenth rib backfat, the p-values for each snp were plotted across the genome as –log using the physical position of the snp in mega-bases. standardization of snp effects using the pev of the marker a second standardization of the jth snp effect was performed using the pev as follows: as discussed above,. the p-values and genome screening for snpep j were assessed and plotted in the same fashion as for snpe j. simulation a plasmode simulation was performed to compare how the standardized values snpe j and snpep j affected the nominal size of the test for the effect to be equal to zero. data on animals with snp each were used for the study, and the snp on chromosome were reshuffled. two scenarios were considered: 1) dependency: rows of the genotype matrix were permuted for columns corresponding to snp on chromosome thus keeping linkage disequilibrium within chromosomes but breaking the relationship between genotypes and phenotypes for the snp on the chromosome. 2) independency: the genotype of any animal was permuted independently by marker for those snp on chromosome and the relationship with the phenotype was broken too. for both scenarios model was fitted to the data, and two tests were calculated for each scenario: test1 = snpej and test2 = snpepj. permutations were repeated times per scenario, and in each permutation the g matrix was calculated while fitting model. as a result, the heritability of the trait was similar to the original heritability due to relationships in the other chromosomes being kept intact, and p-values for those snp on chromosome were obtained for the different tests. under the null hypothesis and assuming independence, an approach that controls for type i error appropriately, the test p-values follow a uniform distribution. consequently, to estimate the empirical quantiles of the distribution for the null hypothesis, we used a uniform density u ∼  to generate replicated sets for the p-values. snp effects and tests obtained by a single marker model the snp effects were tested on a one by one basis. the model approach used for testing purposes is better known as “efficient mixed-model association”. the model included fixed effects of sex and one-marker-at-a-time; random variable was the animal effect with variance-covariance equal to the genomic relationship matrix using all markers, which was calculated as described before. the r package rrblup was used for fitting the different models and for calculating the tests and p-values. proportion of variance explained by segments with large effect after the genome screen using model 2a, the snp with the smallest p-values were selected to form snp segments. these segments were defined by taking all snp within one mb upstream and one mb downstream of the snp with smallest p-value on each chromosome. the size of the segment was chosen using a criterion similar to the one employed by hayes et al.. the point of change in the rate of decay in linkage disequilibrium in this population was about r2 =  at 1 mb, which essentially would imply a minimal contribution to the additive variance from markers located beyond such distance. moreover, segment sizes about two mb have been reported to be significant in association studies. the proportion of variance associated with each segment was estimated by building a genomic relationship matrix g ) using all snps that belonged to the segment, whereas genomic relationship matrix g was built using all remaining snps. the model fitted can be represented as: where a is the vector of additive random effects associated with those snp located in the segment, such that, and a is the vector of additive random effects associated with all snps except those involved with a such that. model assesses the proportion of variance explained by the segment of interest from the genome variance explained by all snps. the variances estimated in were compared with those estimates from model. hayes et al. used a similar model to assess the segment variance. applying either model, or the approach of hayes et al. gave similar estimated variance components. in practice, the advantage of fitting model is that g is computed by subtracting from g the columns of z related to the segment being tested. let zs be a matrix having as columns those related to the segment being tested, then. on the contrary, in the model of hayes et al. gis different from segment to segment. additionally, the calculation of g and is fast and involves only those snps located in the segment. to adjust the level of significance for multiple comparisons, a bonferroni correction was performed. in this context, if the pig genome is ~2800 mb long and the average size of the segment is 2 mb, there are segments along the genome with corresponding multiple tests. thus, for α = , the bc was equal to /1400 =  e−. hence, in order to evaluate the significance of the segments, a second p-value for the likelihood ratio test was calculated to compare against bc. this p − valuelrt was assessed as minus the distribution function of a chi-square random variable with degrees of freedom as follows: where Ω is the distribution function of a random variable having the χ as density, and lrt is the likelihood ratio test obtained by contrasting appropriate models.  genome screening the p-values of the snp were obtained as described in the methods section. first, the p-values for snpej, i.e. using, were plotted along the genome to identify genomic positions that are associated with variation in 13-week tenth rib backfat. large peaks above can be seen at chromosomes and suggesting noticeable genetic variation for the trait. on the other hand, p-values for snpepj were very large, with a maximum − log of. in essence, the pattern observed in figure  is the result of dividing the non-standardized snp effects by a constant. specifically, the normalizing value was, with var =. the use of the square root of the difference between those two values resulted in a practically constant denominator for the test-statistic that was equal to. also, a look at figure  suggests signals at chromosomes and a fact that is not observed in figure  however, this might be an artefact of the constant denominator that tends to overestimate the true variability for some snp, thus resulting in corresponding false positives across the genome.figure manhattan plot for trait 13-week tenth rib backfat by standardization snp ej. genome screening for snp using standardization. −log10  versus the absolute snp position in mb. the red line represents a genome-wide significance threshold. numbers from to represent the chromosome id. manhattan plot for trait 13-week tenth rib backfat by standardization snp epj. genome screening for snp using standardization. −log10  versus the absolute snp position in mb. numbers from to represent the chromosome id. in order to study the type i error rate of the two proposed tests we performed a plasmode simulation. a plasmode is a dataset created from real data where some of the truth is known. in brief, our plasmode is a simulation that uses reshuffling in a portion of the data as explained in the methods section. we performed a simulation assuming independent snp, and another one keeping the dependency between snp intact. simulation results were plotted into a quantil-quantil plot graph using the number –log for each case of standardization. first, the p-values for test obtained in the scenario under independent snps displayed an identical distribution of p-values when obtained by the reference distribution u ∼. in contrast, under dependency less extreme p-values were observed, a fact that was not reflected in a uniform distribution. this is a well known fact in human genetic epidemiology, where the implementation of the bonferroni correction of p-values from associated snp under the assumption of independence results in tests that are too conservative. on the other hand, for test even p-values obtained for independent snp displayed a distribution that was too conservative. furthermore, the results from the dependent scenario were even more conservative than those from the independent scenario, thus indicating that the use of the square root of as the denominator of the test-statistic results in a more powerful and not too liberal choice when compared to the use of the square root of pev =.figure quantil-quantil plot of the observed and expected –log obtained by simulation. reference distribution was an independent and uniform distribution u ∼  for p-values simulated. test1 = under dependent and standardization by. test1 = under independent and standardization by. test2 = under independent and standardization by pev. each scenario has p-values permuted times. bands represent confidence intervals of 95%, green band = test pink band = test snp effects and tests obtained by the marker model the analyses of one snp tested at a time using the emma procedure resulted in p-values that were almost identical to those of snpej. the time taken to compute snp tests one at a time was 84 minutes. in comparison, the algorithm used to fit model and to perform the tests of standardized effects took a total time of 29 minutes. this time includes the computation of the g matrix, the fit of the animal model, the back transformation to calculate the snp effects, and the calculation of the standard errors that are needed to compute the test-statistics. tests of segment effects we also compared the results from our proposed method to those obtained with a segment-based likelihood ratio test that has been used by animal breeders. due to computational demand, we only performed the lrt to test for segment effects. thus, the snp with the smallest p-values ) on each chromosome were chosen, whereas no segments were tested using lrt for regions with snpepj resulting in exceedingly low p-values. the three segments from chromosomes with the smallest p-values are displayed in table  and the remaining segments from the other chromosomes are shown in the additional files. all segments measured 2 mb. the estimates of the variance components and the loglikelihood obtained from model equation were compared with those from model equation. these results are displayed in table  table snp selected by smallest p-value per chromosome snp name = snp marker name, position mb = marker physical position in mega-bases, −log10 = −logarithm in base of the smallest p-value, = absolute value of the snp effect estimated for the trait 13 week tenth rib backfat. variance components and loglikelihood for models with or without the segment snp − log 10 lk_m1 lk_m2 lrt p-value lrt vare_m1 vara_m1 vare_m2 vara_m2 segmva %segmva seg-chromosome = number of chromosome where segment is located, m1 = model without the segment: y = xβ + a + e, m2 = model with the segment y = x  β + a 1 + a 2 + e, snp − log 10  = −logarithm in base of the snp p-value selected to create a segment, lk_m1 = −loglikelihood for m lk_m2 = −loglikelihood for m lrt = likelihood ratio test for m and m p-value lrt = p-value for lrt, vare_m1 = error variance of m vara_m1 = additive variance of m vare_m2 = error variance of m vara_m2 = additive variance of m segmva = additive variance segment of m %segmva = proportion in% of the total variance explained by the segment. results from the lrt indicated that the segment on chromosome was significant: p − valuelrt ‒ 6 =  e− a number smaller than the critical bonferroni threshold for segments. on the contrary, the segments located on all other chromosomes were not significant. the proportion of variance explained by the segment from chromosome  =  ) was 11% of the total variance, a fact that was reflected in a similar reduction of the estimated additive variance in model:  +   = . this latter value is close to, i.e. the estimated value of from model. for all other chromosomal segments, the estimated value of did not decrease to a significant amount. discussion the main goal of this research was to develop a novel procedure to perform a rapid genome scan, or gwas analysis, from a genomic evaluation. moreover, the sufficient statistics of our methodology are: the best linear unbiased prediction of the breeding values from an animal model, g as the covariance matrix, z as the standardized marker effects matrix, variance components, and caa. this setting makes the implementation extremely feasible after the genomic evaluation has been performed as discussed by legarra et al.. variance of the snp effect first, the snp effects were calculated by a linear transformation of using expression. then, we calculated using an expression derived from mixed model theory ). next, we divided by the square root of to standardize the effect, and referred the statistics as snpej. the p-values for the tests of specific genome regions were calculated with a level of significance − log10 =  additionally, prediction error variance () was employed for a second standardization, and it was called the snpepj statistic. after the analyses, we obtained higher p-values  =  ) and detected stronger signals for snpepj than with snpej. furthermore, a simulation was carried out with the same structure of snps markers and animal data as in the current study, in order to compare the performance of empirical p-values of both standardized tests. the snps markers of chromosome were reshuffled, and two scenarios were simulated: 1) dependent genotypes, and 2) independent genotypes. neither scenario displayed a relationship with the phenotype, whereas both standardized tests were calculated at each scenario. the reference distribution for the p-values considered was the uniform. in the independent scenario, standardization with gave an empirical distribution of p-values that resembled the uniform density, but in the dependent scenario the snpej performed conservatively. instead, the standardization with produced conservative results in the independent scenario, and very conservative tests in the dependent scenario. in this context, standardizing snp effects with resulted in p-values that were closer to the simulated ones. moreover, the performance of snpej under ld was not too conservative, a scenario that could be extrapolated to the genotypes in the current study. in addition, the p-values calculated using the emma procedure were similar to those obtained with snpej. these results suggest that snpej behaves reasonably to control type i error rate or false positives. also, the computing time for fitting model and then calculating using expressions - was to times less than the computing time for the emma model. in order to identify snp with important phenotypic associations, the calculation of snp effects from genomic breeding values has been used in several studies. in this context, the variance of snp effects has been estimated using different approaches. wang et al. employed the classical definition of the variance of additive effects from quantitative genetics, so that the variance for each jth marker was obtained as follows:. whereas, mcclure et al. proposed equating the variance of snp effects to, and then normalizing the snp effects with the square root of this estimated and constant variance. this test performed similar to snpep j, when the estimated snp effects was divided by a constant denominator, a value almost equal to the prior variance, and resulted in a very conservative test. in contrast, the advantage of the standardized test presented here was that each snp effect was scaled by its own standard deviation rather than the use of a prior variance or by the square of each specific snp effect as variance. furthermore, the computation of snpej, involves the same variance for the same snps markers and animals, i.e., and the use of the standardized incidence matrix z, a function of 2pj, takes into account this latter quantity into snpej. additionally, the matrix z uses the allele frequencies from the f generation calculated with unrelated individuals, and a proper expected variance by marker. in addition, the test statistics snpej that standardizes snp effects produces a p-value, a result that is appealing to many researchers that are more familiar with the method of testing one snp at the time rather than with the proportion of additive variance that is explained by a genomic region. a further advantage of the method is that detection of many false positives are avoided, and genome positions with sizeable effects are highlighted. candidate segment approach later in the research, genome segments that expressed higher signals were located. to this purpose, snps with the smallest p-values from snpej were selected, and for each of these snp a segment of 2 mb long was created. the next step was to estimate the variance components and the log-likelihood from the centered animal models and. the latter model includes the random vector of snp segments a lastly, we compare the performance of both models. hayes et al. used a similar model to, although the random snp effect was taken from the breeding value and fitted as a separate segment effect. we observed similar results from the use of either approach. the advantage of fitting model is that matrix g is the same for all segments, so that it was calculated only once, and stored in memory for the calculations, whereas in the model of hayes et al. a different g has to be calculated for each segment. this implies an extended computing time and higher requirements of cpu memory to obtain similar results to those from model. to evaluate the significance of the segments, the effects of each chromosome segment were tested by the likelihood ratio test. the size of the test was adjusted by the bonferroni correction. as a result, the segment located on chromosome was significant, and explained 11% of the trait total variance. previous studies by edwards et al. and choi et al., using microsatellites and a small number of snp, found significant regions on chromosome for 13 week tenth rib backfat in the current population under study. additionally, forty eight markers between the physical position between 128 mb and 139 mb on chromosome, have been reported to be associated with the trait. furthermore, recent studies showed the importance of chromosome in the expression of the trait. therefore, our results confirm the presence of genetic variability in the trait from chromosome  fast genome screening of snp effects linearly transformed from genomic breeding values is advantageous, as a by-product of genomic evaluations for different species of farm animals. moreover, the standardized tests of snp effects using their own variance developed in this study helps in detecting specific genomic regions involved in the additive variation of the trait and reducing false positive locations using less computing time. additionally, genome segments of about 2 mb formed by surrounding the snp with the smallest p-values on each chromosome, and tested with a standardized test involving and with the bonferroni correction, could detect genome regions responsible for sizeable fractions of the trait genetic variance. this methodology involving genome scan and candidate segment approach is a useful method for meta-analyses of genome-wide association studies, as it enables the detection of specific genome regions that affect an economically relevant trait when using multiple populations. code and data to obtain and reproduce the results presented is publicly available at https://www.msu.edu/~steibelj/jp_files/gblup.html. electronic supplementary material additional file 1: highest − log 10 on each chromosome for trait 13-week tenth rib backfat by standardization snp ej and emma. the blue and red circle represents highest − log on each chromosome by the standardization snp ej and efficient mixed-model association using rrblup. respectively. additional file 2: dispersion plot of − log 10 for trait 13-week tenth rib backfat by emma and standardization snp ej. dispersion plot for –log by efficient mixed-model association using the rrblup r package in the x axis, and by the standardization snp ej in the y axis. red straight line is the reference line 0– additional file 3: variance components and loglikehood for models with or without the segment for all chromosomes.. competing interests the authors declare that they have no competing interests. authors’ contributions jps, rjcc, jlgd: performed and supervised statistical and simulation analyses and wrote the manuscript. rob, cwe: designed the resource population and led collection of phenotypic data. cwe, ner: performed dna extraction and coordinated genotyping with commercial laboratory. jps, rob, cwe: designed high density genotyping scheme. all authors read and approved the paper. 
the marker with the highest − log in each chromosome is selected, and the segment is expanded one mb upstream and one mb downstream of the marker. currently, association studies are analysed using statistical mixed models, with marker effects estimated by a linear transformation of genomic breeding values. however, approaches used to estimate the parameters rely on a prior variance or on a constant estimate of the additive variance. moreover, genome scanning of candidate segments can be used in meta-analyses of genome-wide association studies, as it enables the detection of specific genome regions that affect an economically relevant trait when using multiple populations. the standardized test of marker effects using their own variance helps in detecting specific genomic regions involved in the additive variance, and in reducing false positives. furthermore, genomic predictions are solved efficiently by the procedure, and the p-values are virtually identical to those calculated from tests for one marker effect at a time. moreover, the standardized test reduces computing time and memory requirements. the likelihood ratio is then calculated to test for the effect in every chromosome against a reduced model with fixed effects and genomic breeding values. the following steps are used to locate genome segments displaying strong association. the standardized test was neither conservative nor liberal with respect to type i error rate, compared to a similar test using predictor error variance, a method that was too conservative.
 currently, association studies are analysed using statistical mixed models, with marker effects estimated by a linear transformation of genomic breeding values. the variances of marker effects are needed when performing the tests of association. however, approaches used to estimate the parameters rely on a prior variance or on a constant estimate of the additive variance. alternatively, we propose a standardized test of association using the variance of each marker effect, which generally differ among each other. random breeding values from a mixed model including fixed effects and a genomic covariance matrix are linearly transformed to estimate the marker effects.  the standardized test was neither conservative nor liberal with respect to type i error rate, compared to a similar test using predictor error variance, a method that was too conservative. furthermore, genomic predictions are solved efficiently by the procedure, and the p-values are virtually identical to those calculated from tests for one marker effect at a time. moreover, the standardized test reduces computing time and memory requirements. the following steps are used to locate genome segments displaying strong association. the marker with the highest − log in each chromosome is selected, and the segment is expanded one mb upstream and one mb downstream of the marker. a genomic matrix is calculated using the information from those markers only, which is used as the variance-covariance of the segment effects in a model that also includes fixed effects and random genomic breeding values. the likelihood ratio is then calculated to test for the effect in every chromosome against a reduced model with fixed effects and genomic breeding values. in a case study with pigs, a significant segment from chromosome explained 11% of total genetic variance.  the standardized test of marker effects using their own variance helps in detecting specific genomic regions involved in the additive variance, and in reducing false positives. moreover, genome scanning of candidate segments can be used in meta-analyses of genome-wide association studies, as it enables the detection of specific genome regions that affect an economically relevant trait when using multiple populations. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords genome wide associationmarker variancepig genotypeissue-copyright-statement© the author 2014 