segregation distortion is a phenomenon that the genotypic frequency array of a locus does not follow a typical mendelian ratio . depending on the population under investigation , mendelian ratio of a locus varies from 1 : 1 for a backcross to 1 : 2 : 1 for an f2 and to 1 : 1 : 1 : 1 for a four - way cross . , the effects of distorted markers on the result of qtl mapping were not known . for the reason of precaution , people simply discarded all the distorted markers in qtl mapping . recently , we found that distorted markers can be safely used for qtl mapping with no detrimental effect on the result of qtl mapping . this finding can help qtl mappers save tremendous resources by using all available markers , regardless whether they are mendelian or not . we also found that if distorted markers are handled properly , they can be beneficial to qtl mapping . the reason behind the distortion is due to one or more segregation distortion loci ( sdl ) . these loci are subject to gametic selection , zygotic selection , or both and their ( unobservable ) distorted segregation causes the observed markers to deviate from the mendelian ratio . several investigators [ 411 ] have attempted to map these segregation distortion loci using molecular markers . agricultural scientists are interested in mapping qtl for economically important traits while evolutionary biologists are interested in mapping sdl that respond to natural selection . since the theory of segregation distortion has been introduced and discussed in previous studies [ 7 , 8 ] and our own research , this study only presents the em ( expectation - maximization ) implementation of the statistical method . the variance - covariance matrix of estimated parameters under the em algorithm is also derived and presented in appendix a for interested readers . we only investigate interval mapping where a model contains a single qtl at a time and the entire genome is scanned through repeated calling of the same program for different locations of the genome . the technical difference between the joint mapping and qtl mapping occurs only in one place . in the traditional interval mapping of qtl , the conditional probabilities of genotypes for a qtl are calculated using flanking marker genotypes with the prior probabilities of qtl genotypes being substituted by the mendelian ratio . for the joint mapping , the genotypic frequencies ( segregation ratios ) are treated as unknown parameters that are subject to estimation . let m and n be the left and right flanking markers bracketing the qtl ( denoted by g for short ) . the interval of the genome carrying the three loci is labeled by a segment mgn . the three genotypes of the qtl are denoted by g1g1 , g1g2 , and g2g2 , respectively . let r1 and r2 be the recombination fractions for segments mg and gn , respectively . the joint distribution of the marker genotypes conditional on the qtl genotype can be derived using the markov chain property under the assumption of no interference between consecutive loci in segregation . let us denote the three ordered genotypes , g1g1 , g1g2 , and g2g2 , by genotypes 1 , 2 , and 3 , respectively . if individual j takes the   th genotype for the qtl , we denote the event by gj =  , for all  = 1 , 2 , 3 . the joint probability of the two markers conditional on the genotype of the qtl is 
								 ( 1)pr(mj=,nj=    gj= )  = pr(mj=    gj=)pr(nj=    gj= ) 
							 for all  ,  ,  = 1 , 2 , 3 , where pr(mj =   gj =  ) = 1( ,  ) and pr(nj = we use i( ,  ) to denote the th row and the th column of the following transition matrix : 
								 ( 2)i=[(1ri)22ri(1ri)ri2ri(1ri)(1ri)2+ri2ri(1ri)ri22ri(1ri)(1ri)2 ] ,  i=1,2 . for example , 
								 ( 3)pr(mj=1,nj=2    gj=3 )  = pr(mj=1    gj=3)pr(nj=2    gj=3 )  = 1(3,1)2(3,2)=2r12r2(1r2 ) . let  = pr(g =  ) , for all  = 1 , 2 , 3 , be the probability that a randomly sampled individual from the f2 family has a genotype . we use a generic notation p for probability so that p(gj =  ) represents pr(gj =  ) and p(mj , nj  gj =  ) stands for pr(mj , nj  gj =  ) . the log likelihood function of the flanking marker genotypes in the f2 population is 
 ( 4)l(    m)=j=1nln  [ =13p(gj=)p(mj , nj    gj=)]=j=1nln  [ =131(,mj)2(,nj ) ] , 
							 where  = [ 1  2  3 ] is a vector of parameters with constraint =1 = 1 , where m in l(  m ) stands for marker information . note that without any prior information , p(gj =  ) =  , for all j = 1 ,  , n. under the assumption of mendelian segregation ,  =  where  = [ 1  2  3 ] = [ 1/4  1/2  1/4 ] . because we are dealing with the genotypic frequencies , the segregation distortion is called zygotic distortion we postulate that deviation of  from  causes a marker linked to locus g to show distorted segregation . let yj be the phenotypic value of a quantitative trait measured from individual j. the probability density of yj conditional on the genotype of individual j is normal with mean 
 ( 5)=xj+h 
							 and variance  , that is , 
 ( 6)p(yj    gj=)=122exp  [ 122(yjxjh)2 ] , 
							 where h is the th row of matrix h and  
 ( 7)h=[+110111 ] . this h matrix can be defined in a different scale , for example ,  
 ( 8)h=[+100110 ] , 
							 which does not affect the significance test . the advantage of choosing the scale in ( 7 ) is that the expectation of the dominance indicator is zero . the design matrix xj and the regression coefficients  capture non - qtl effects , for example , field location effects , year effects , and so on . the likelihood function of the phenotypic values in the f2 population is 
								 ( 9)l(,,2,    y )  = j=1nln  [ =13p(gj=)p(yj    gj= ) ]  = n2ln  ( 2)+j=1nln  { =13exp  [ 122(yj)2 ] } , 
							 where letter y in l( ,  ,  ,   y ) stands for the phenotype . this likelihood function is the one used in segregation analysis of quantitative traits   because no marker information is required . let  = [        ] be a vector of all parameters in the joint analysis . the likelihood function is obtained by combining ( 4 ) and ( 9 ) :  
 ( 10)l(    m , y)=j=1nln  [ =13p(gj=)p(yj    gj=)p(mj , nj    gj=)]=j=1nln  { =13exp  [ 122(yj)2]1(,mj)2(,nj ) }  n2ln  ( 2 ) . 
							 for qtl mapping under segregation distortion the previous two likelihood functions ( for markers and for phenotypes ) were presented as background information to introduce this joint log likelihood function . the mle ( maximum likelihood estimate ) of the parameters is solved via an em algorithm . let us define a delta function as  
 ( 11)(gj,)={1if  gj=,0if  gj , 
							 if the genotypes of all individuals are known , that is , given (gj ,  ) for all j = 1 ,  , n and  = 1 , 2 , 3 , the complete - data log likelihood is 
 ( 12)l(,)=j=1nln  [ p(yj    gj)p(mj , nj    gj)p(gj ) ] , 
							 where 
 ( 13)p(yj    gj)=122exp  [ 122=13(gj,)(yj)2],p(mj , nj    gj)=k=13p(mj , nj    gj=)(gj,)=k=13[1(,mj)2(,nj)](gj,),p(gj)==13(gj, ) . therefore , we need to take expectation of the log likelihood with respect to . the expected complete - data log likelihood function is  
 ( 14)l(    (t))=e[l(, )    (t)]=0+1()+2( ) . note that e[l( ,  )   ] stands for the expectation of l( ,  ) with respect to  conditional on parameters at the current state  and the data ( the symbol for data is suppressed for simplicity ) . the three components of ( 14 ) are  
 ( 15)0=j=1n  =13e[(gj, )    (t ) ]  [ln  1(,mj)+ln  2(,nj)],1()=n2ln  ( 2 )  122j=1n  =13e[(gj, )    (t)](yj)2,2()=j=1n  =13e[(gj, )    (t)]ln  . 
							 the first component 0 is a function of  but not a function of . therefore , it is considered as a constant . the expectation step of the em algorithm requires computing the expectation of  conditional on the data and . because  is a bernoulli variable , the expectation is simply the probability of  = 1 , that is ,  
 ( 16)e[(gj, )    (t)]=pr[(gj,)=1    (t),m , y]=p(gj=)p(yj    gj=)p(mj , nj    gj=)=13p(gj=)p(yj    gj=)p(mj , nj    gj=)=exp  [ (1/22)(yj)2]1(,mj)2(,nj)=13exp  [ (1/22)(yj)2]1(,mj)2(,nj ) . the maximization step of the em algorithm requires taking the partial derivatives of l(   ) with respect to  , setting the partial derivatives equal to zero , and solving for the parameters : 
 ( 17)l(    (t))=1()+2()=0 . the solutions of the parameters are 
 ( 18)=[j=1nxjtxj]1[j=1n=13e[(gj,)]xjt(yjh)],=[j=1n  =13e[(gj,)](hth)]1  [j=1n  =13e[(gj,)]ht(yjxj)],2=1nj=1n  =13e[(gj,)](yjxjh)2,=1nj=1ne[(gj, ) ] ,  =1,2,3 . although the wald test can be performed for testing the presence of qtl , such a test is not justified for testing the null hypothesis of mendelian segregation . therefore , the likelihood ratio tests are more justifiable . regardless what hypothesis is tested , the full model joint log likelihood function given in ( 10 ) is required . let us reintroduce this joint log likelihood function using a different notation so that different likelihood ratio tests are easily interpreted . the joint likelihood is rewritten as 
 ( 19)lqs(,)=n2log  ( 2)+j=1nln  { =13exp  [ 122(yj)2]1(,mj)2(,nj ) } , 
							 where  represents qtl effects and  stands for non - mendelian segregation . the null hypothesis for qtl detection is hqtl :  = 0 while the null hypothesis for detecting segregation distortion is hsdl :  = . the null hypothesis is hqtl :  = 0 . the likelihood ratio test statistic is 
 ( 20)qtl=2[ls(0,^)lqs(^,^ ) ] , 
									 where ls(0,^ ) is the log likelihood value under the null model  = 0 , which is 
 ( 21)ls(0,^)=l(^,^2    y)+l(^    m ) , 
									 where  
 ( 22)l(^,^2    y)=n2ln  ( ^2)12^2j=1n(yjxj^)2,l(^    m)=j=1nln  [ =13^1(,mj)2(,nj ) ] . 
									 the estimated parameters in ( 22 ) are obtained by maximizing the corresponding likelihood functions ( see appendix b for the estimation ) . the null hypothesis is hsdl :  = . the likelihood ratio test statistic is 
 ( 23)sdl=2[lq(^,)lqs(^,^ ) ] , 
									 where 
 ( 24)lq(^,)=n2ln  ( ^2)+j=1nln  { =13exp  [ 12^2(yj^)2]1(,mj)2(,nj)}. 
									 again , the mles of the parameters in ( 24 ) are obtained by maximizing this likelihood function ( see appendix b for the estimation of parameters under the null model ) . the null hypothesis is h0 :  = 0 and  = . the likelihood ratio test statistic is 
 ( 25)qs=2[l(0,)lqs(^,^ ) ] , 
									 where  
 ( 26)l(0,)=l(^,^2    y)+l( m ) . 
									 the two components of ( 26 ) are  
 ( 27)l(^,^2    y)=n2ln  ( ^2)12^2j=1n(yjxj^)2,l(    m)=j=1nln  [ =131(,mj)2(,nj ) ] . the parameters involved in these log likelihood functions are estimated using formulas given in appendix b. this hypothesis is rejected if either   0 or    or both inequalities hold . suppose that , for some reason , we know for sure that the population from which the sample is drawn is a mendelian population . the sample drawn from this population is selected based on extreme phenotypes ( selective genotyping ) . the sample is then non - mendelian regarding the qtl that control the trait subject to phenotypic selection . the reason is that segregation distortion in the sample is solely caused by selective genotyping . this example demonstrates the application of the method to joint mapping of qtl and sdl in a wheat qtl mapping experiment . the experiment was conducted by dou et al .   who made the data available to us for this analysis . a female sterile line xnd126 and an elite cultivar gaocheng 8901 with normal fertility were crossed for genetic analysis of female sterility measured as a ratio of the number of seeded spikelets to the total number of spikelets per plant . the parents and their f1 and f2 progeny were planted in the huaian experimental station in china for the 2006 - 2007 growing season under the normal autumn sowing condition . the sterility trait was transformed using the angular sine transformation , y = arcsin(x ) , where x is the phenotypic value expressed as ratio . these markers covered five chromosomes of the wheat genome with an average genome marker density of 15.5  cm per marker interval . the five chromosomes are only part of the wheat genome . the model for the female sterility is 
 ( 28)yj = xj+zj11+zj22+j , 
					 where xj = 1 for all j ,  is the intercept , zj1 = { + 1,0 ,    1 } is the genotype indicator variable for the additive effect 1 = a , and zj2 = { -1,1 ,    1 } is the genotype indicator variable for the dominance effect 2 = d. we used an interval mapping approach to scanning the entire genome . therefore , the model contains one qtl at a time . with the interval mapping , zj = [ zj1  zj2 ] is missing and can take one of three values denoted by h for  = 1 , 2 , 3 ( see definition of h ) in section 2 . the likelihood ratio test statistics were divided by 4.61 to obtain their corresponding lod scores . the top panel ( a ) shows the lod profile for qtl detection regardless whether there is segregation distortion or not . figure 3(b ) shows the lod profile for testing segregation distortion , regardless whether the qtl is present or absent . we used the quick method of piepho   to calculate the genome wide critical value of the lod for significance test . we found that the genome wide critical value was slightly less than lod = 3 criterion ( data not shown ) . in addition to the two major loci , several other regions of the genome also showed significant peaks of the lod score profile . overall we detected eight loci , four qtls and four sdls . among the detected qtl , each explains from 10% to 47% of the phenotypic variance ( listed as heritability denoted by h in table 1 ) . among the four sdl detected , all showed bias in favor of the xnd126 parent , that is , the homozygote of xnd126 allele was over represented at the cost of low representation of the other parent . the frequency of the heterozygote was close to the mendelian frequency of 0.5 , but the homozygote of gaocheng 8901 allele was almost wiped out . the estimated genotypic frequencies are plotted against the genome location as shown in figure 2 . one of the major theoretical contributions of this study is the development of the variance - covariance matrix of the estimated qtl - sdl parameters . the covariances between pairs of estimated parameters are not of interest , but the variances of the estimated parameters are important . we reported the standard errors for two selected loci , locus 4 ( qtl ) and locus 7 ( sdl ) . the standard errors are the square roots of the variances obtained from the em algorithm . the variance - covariance matrix of the estimated parameters takes the inverse of the information matrix . as a result we also performed a bootstrap analysis ( 1000 bootstrap samples ) to provide more accurate estimation of the variance . the results of the bootstrap estimates of the standard errors for the two loci are also listed in table 2 for comparison . the approximate standard errors from the em algorithm are indeed biased downward , especially for locus 4 ( qtl ) . the bootstraps method is recommended for obtaining more accurate estimates of the standard errors if the sample size is small . methods also are available for mapping viability loci or segregation distortion loci when markers do not segregate in a typical mendelian ratio [ 46 , 911 ] . however , qtl mapping and sdl mapping have never been combined in a single analysis . this study is the first attempt to combine the two seemingly different analyses into a joint one . when qtl and sdl are loosely linked or not linked , the joint analysis does not offer much advantage over the separate analyses . when they do overlap , a phenomenon called pleiotropy , joint mapping does offer some advantage . unfortunately , the wheat experiment introduced here is not a good example to demonstrate the advantage of joint analysis because the qtl and sdl detected do not overlap . an obvious situation where the joint analysis can be more powerful is qtl mapping with selective genotyping . in most designed selective genotyping experiments , directional ( one - tailed ) selection is rarely used for selective genotyping because it artificially reduces the variation of the trait and thus reduces the statistical power of qtl detection . however , with the joint analysis , the power can increase under directional selection . we simulated 500  f2 individuals for a single chromosome of 300  cm long . we placed 16 markers evenly over the chromosome with 20  cm per marker interval . we selected 300 smallest individuals out of the 500 ( one - tailed selection ) for mapping . figure 3(a ) shows the test statistic for qtl regardless whether segregation is distorted or not . the panel in the middle ( b ) shows the lod test statistic for segregation distortion , regardless whether a qtl is present or not . the panel at the bottom ( c ) shows the lod score profile of the joint analysis where the null model is no qtl and no sdl . if lod = 3 is the threshold value , neither of the separate analyses is significant . however , the joint analysis has a lod score as high as 4.5 , indicating a significant qtl in the middle of the chromosome . our model was developed under zygotic selection because we are dealing with the genotypic frequencies . however , if the true cause of segregation is gametic selection , we can still detect segregation distortion as long as the gametic selection leads to the genotypic frequencies deviating from the expected mendelian ratio . a model particularly handling gametic selection has not been developed yet , but it is not difficult . let us take the f2 population as an example to show the gametic selection model . denote the frequencies of g1 and g2 alleles from the female parent by 1 and 2 , respectively , for 1 + 2 = 1 , and the corresponding allele frequencies from the male parent . when gametic selection occurs , we treat 1 and 1 as unknown parameters for estimation . the genotypic frequencies are simply functions of the two unknown parameters , as given by 1 = 11 , 2 = 12 + 21 and 3 = 22 . qtl parameters are estimated using the same algorithm as described in the zygotic selection model because they depend only on the genotypic frequencies . estimating parameters 1 and 1 requires a modified algorithm . under the gametic selection model , we can test whether the segregation distortion is caused by the distortion of female gametes , male gametes , or both . the joint analysis developed in this study only applies to line crossing data where the marker linkage phases are known . application of the method to pedigrees warrants further investigation and it is not obvious to us at this moment . we have already extended the method to bc ( backcrosses ) , ril ( recombinant inbred lines ) , dh ( double haploids ) , and fw ( four way crosses ) and incorporated them into our qtl mapping program that is described in the next paragraph . however , interval mapping is still the quickest method of qtl mapping , even though multiple qtl mapping programs are available . compared with traditional qtl interval mapping , this joint analysis involves one additional step of updating the genotypic frequencies . this additional step presents a complication where the conditional genotypic frequencies given flanking marker genotypes can not be calculated prior to qtl mapping . fortunately , we have incorporated the joint qtl / sdl mapping into our qtl mapping program . the method =  ml  option in the proc qtl statement must be turned on with an additional suboption /distortion to invoke the joint mapping procedure . without the /distortion option proc qtl is available on our website ( http://www.statgen.ucr.edu/software.html ) and users can download the program with no charge . the program is also accompanied with a detailed user manual . let us define the individual - wise complete - data log likelihood for plant j as 
 ( a.1)lj(,)=[12ln  ( 2)122=13(gj,)(yjxjh)2]+=13(gj,)[ln  ( t1(,mj))+ln  ( t2(,nj))]+=13(gj,)ln   , 
							 where 3 = 1  1  2 so that 3 is excluded from the parameter vector . to derive the variance covariance matrix of the estimated parameters , we need to define the score vector sj( ,  ) and the hessian matrix hj( ,  ) for the individual - wise complete - data log likelihood function . the louis   information matrix of the parameters under the em algorithm is  
 ( a.2)i(^)=j=1ne[hj(^,)]j=1nvar[sj(^, ) ] , 
							 where the expectation and variance are taken with respect to the missing values . once the information matrix is define , the variance matrix of the estimated parameters simply takes  
 ( a.3)var(^)i1(^ ) . 
							 the standard error of each parameter takes the square root of each diagonal element of the above matrix . the approximation performs better for large sample sizes . the score vector is denoted by sj( ,  ) = lj( , )/ , which consists of five blocks , as follows : 
 ( a.4)sj(,)=lj(,)=12k=13(gj,)xjt(yjxjh),sj(,)=lj(,)=12=13(gj,)ht(yjxjh),sj(2,)=lj(,)2=122  + 124=13(gj,)(yjxjh)2,sj(1,)=lj(,)1=(gj,1)11(gj,3)1112,sj(2,)=lj(,)2=(gj,2)12(gj,3)1112 . 
							 concatenating the above five scores vertically the hessian matrix is denoted by hj( ,  ) = lj( , )/ , which is block diagonal with non - zero blocks given as follows : 
 ( a.6)hj()=2lj(,)t=12xjtxj , hj()=2lj(,)t=12=13(gj,)xjth,hj(2)=2lj(,)2=14=13(gj,)xjt(yjxjh),hj()=2lj(,)t=12=13(gj,)hth,hj(2)=2lj(,)2=14=13(gj,)ht(yjxjh),hj(22)=2lj(,)22=12416=13(gj,)(yjxjh)2,hj(11)=2lj(,)11=(gj,1)112(gj,3)132,hj(12)=2lj(,)12=(gj,3)132,hj(22)=2lj(,)22=(gj,2)122(gj,3)132 . the expectation of the hessian matrix e[hj( , )]and the variance matrix of the score vector var[sj( ,  ) ] can be expressed explicitly because both the hessian matrix and the score vector are linear functions of the missing value  
 ( a.8)j=[(gj,1)(gj,2)(gj,3)]t . therefore , e[hj( ,  ) ] and var[sj( ,  ) ] can eventually be expressed as functions of the expectation and variance of j , which have simple expressions because j is a multinomial variable . the hessian matrix hj( ,  ) is already expressed in linear function of j and thus the expectation can be obtained straightforwardly by replacing j by e(j ) . let us define the following matrices : 
 ( a.9)cj()=[12xjt(yjxjh1)12xjt(yjxjh2)12xjt(yjxjh3)],cj()=[12h1t(yjxjh1)12h2t(yjxjh2)12h3t(yjxjh3)],cj(2)=[124(yjxjh1)2124(yjxjh2)2124(yjxjh3)2],cj(1)=[11013 ] ,  cj(2)=[01213 ] . the score vector in matrix notation is 
 ( a.10)sj(,)=lj(,)=[cjt()cjt()cjt(2)cjt(1)cjt(2)]j+[0012200]=cjj+d . 
							 as a result ,  
 ( a.11)var[sj(,)]=cjvar(j)cjt , 
							 where  
 ( a.12)var(j)=diag  [ e(j)]e(j)e(jt ) 
							 is the variance - covariance matrix of the multinomial variable j . this appendix provides methods for estimating parameters under various null models that are required for constructing likelihood ratio test statistics .   the log likelihood for the null model is 
 ( b.1)ls(0,^)=l(^,^2    y)+l(^m ) , 
									 where 
 ( b.2)l(^,^2    y)=n2ln  ( ^2)12^2j=1n(yjxj^)2,l(^    m)=j=1nln  [ =13^1(,mj)2(,nj ) ] . the mles of the population parameters are obtained explicitly without iteration : 
 ( b.3)^=[j=1nxjtxj]1[j=1nxjtyj],^2=1nj=1n(yjxj^)2 . the estimated genotypic frequencies , however , require the following iterations :  
 ( b.4)(t+1)=1nj=1ne[(gj, ) ] ,  =1,2,3 , 
									 where 
 ( b.5)e[(gj,)]=(t)1(,mj)2(,nj)=13(t)1(,mj)2(,nj ) . the log likelihood function for the null model is 
 ( b.6)lq(^,)=n2ln  ( ^2)+j=1nln  { =13exp  [ 12^2(yj^)2]1(,mj)2(,nj)}. 
									 finding the solution of the parameters requires iterations as given as foolows : 
 ( b.7)(t+1)=[j=1nxjtxj]1  [j=1n  =13e[(gj,)]xjt(yjh(t))],(t+1)=[j=1n  =13e[(gj,)](hth)]1  [j=1n  =13e[(gj,)]ht(yjxj(t))],2(t+1)=1nj=1n  =13e[(gj,)](yjxj(t)h(t))2 , 
									 where 
 ( b.8)e[(gj,)]=exp  [ (1/22)(yj)2]1(,mj)2(,nj)=13exp  [ (1/22)(yj)2]1(,mj)2(,nj ) . the mendelian frequencies  's are constants , not parameters , as 1 = 3 = ( 1/2)2 = 1/4 . the log likelihood function for the null model is  
 ( b.9)l(0,)=l(^,^2    y)+l(    m ) , 
									 where 
 ( b.10)l(^,^2    y)=n2ln  ( ^2)12^2j=1n(yjxj^)2,l(m)=j=1nln  [ =131(,mj)2(,nj ) ] .