# A Statistical Variance Components Framework for Mapping Imprinted Quantitative Trait Locus in Experimental Crosses

Authors: Gengxin Li; Yuehua Cui
Journal: Journal of Probability and Statistics
 (2009)

---

## Abstract

Current methods for mapping imprinted quantitative trait locus (iQTL) with inbred line
crosses assume fixed QTL effects. When an iQTL segregates in experimental line crosses,
combining different line crosses with similar genetic background can improve the accuracy
of iQTLs inference. In this article, we develop a general interval-based statistical variance
components framework to map iQTLs underlying complex traits by combining different
backcross line crosses. We propose a new iQTL variance partition method based on the
nature of marker alleles shared identical-by-decent (IBD) in inbred lines. Maternal effect is
adjusted when testing imprinting. Efficient estimation methods with the maximum likelihood
and the restricted maximum likelihood are derived and compared. Statistical properties of
the proposed mapping strategy are evaluated through extensive simulations under different
sampling designs. An extension to multiple QTL analysis is given. The proposed method
will greatly facilitate genetic dissection of imprinted complex traits in inbred line crosses.


---

## Body

## 1. Introduction

The genetic
architecture of complex phenotypes in agriculture, evolution, and biomedicine
is generally complex involving a network of
multiple genetic and environmental factors that interact with one another in
complicated ways [1]. The development of molecular markers makes it possible to
identify genetic locus (i.e., quantitative trait locus (QTL)) underlie various
traits of interest. Genetic designs with controlled crosses are generally
pursued to generate mapping populations aimed to identify QTL underlying the variation of phenotypes. Statistical method for
QTL mapping with experimental crosses dates back to the seminal work of Lander
and Botstein [2]. Various extensions have been developed since then (e.g., [3, 4]).For a diploid organism, the expression products of most
functional regions from each one of a chromosome pair are equal. A broken of
this equivalence, that is, nonequivalent genetic contribution of each parental
genome to offspring phenotype, can result ingenomic imprinting, a
phenomenon also called parent-of-origin effect [5]. Since its discovery, imprinting-like phenomena have been commonly observed in mammals and seed
plants (reviewed by Burt and Trivers [6]). However, statistical methods for
identifying imprinted genes have not been extensively studied and well
developed.The imprinted inheritance violates the Mendelian
theory and brings challenges in statistical modeling. Currently, there are two
frameworks in mapping imprinted genes. One is based on the random effect model
with pedigree-based natural population such as humans. Hanson et al. [7] first
proposed a variance components framework by partitioning the additive variance
component as two parts, a component due to maternal gene and a component due to
paternal gene. The variance component method is developed based on the
identical-by-decent (IBD) idea in which the expression of the gene for a pair
of individuals is expected to be similar if they share alleles IBD. Liu et al. [8] recently applied the model to map iQTL underlying canine hip dysplasia in a
structured canine population. However, the current IBD-based variance
components method for mapping imprinted genes assumes noninbreeding population. Their applications are immediately limited with fully or partially inbreeding
population such as the controlled inbreeding design in plants and animals. With
inbred mapping population in humans, Abney et al. [9] proposed a method to
estimate variance components of quantitative traits. However, the extension of
the method to map imprinted gene is not straightforward. No variance components
method has been proposed to map imprinted genes with inbred population in the
literature.Another general framework for mapping imprinted genes
is based on the fixed-effect model in which the effects of genetic factors are considered
as fixed. A number of studies were proposed under this framework for mapping
imprinted QTL (iQTL) with controlled crosses of outbred parents [10–12]. One
potential limitation of these methods is that allelic heterozygosity at a locus
between two outbred parents could cause confounding effects for genomic
imprinting. The genetic differences detected by such a fixed-effect model could
be caused by allelic heterozygosity of the parents rather than the imprinted
effect of iQTL [13]. A natural alternative for the mapping population is the
inbred lines. Fixed-effect models based on backcross (BC) and F2 population were
recently proposed under the maximum likelihood framework [14–17]. When inbred
lines are used, Xie et al. [18] pointed out that it is more meaningful to
inference QTL effect by its variance rather than by the allele substitution
effect. The QTL variance is generally calculated conditional on the cross, and
it, as a variable, is different from one cross to another [18]. In a
single-line cross, the estimated QTL variance cannot be simply extended to a
statistical inference space beyond that [18]. Multiple parental lines are
needed for QTL variance inference. A solution to this is to combine data from
multiple line crosses [18]. An IBD-based variance component method was proposed
by Xie et al. [18] with multiple line crosses. Extension of the IBD-based
variance component method with multiple line crosses to iQTL mapping has not
been studied.Motivated by the limitations of current methods
aforementioned and by the pressing need for efficient iQTL mapping procedure, in this article, we propose a statistical variance components framework for
iQTL mapping by combining data from multiple inbred line crosses. The proposed
model is robust in iQTL variance inference by extending the iQTL inference
space from single-line cross to multiple line crosses. A parent-specific IBD
sharing partition method is proposed by considering the inbreeding structure in
line crosses. As discussed by Cui in [14], the phenotype of an offspring is not
only controlled by its own genetic profiles, but also controlled by maternal
genotype. The effect of maternal genotype on the phenotype of her offspring, termed
maternal effect, is one potential source of confounding effect in the inference
of genomic imprinting. The existence of such parental effect may lead to
incorrect interpretations of imprinting when they are not properly accounted
for in the analysis. Parameters that model the maternal effect are also
included and adjusted when testing imprinting.With the developed model, we propose an interval-based
method for genome-wide scan and testing of iQTL. Both maximum likelihood (ML)
and restricted maximum likelihood (REML) methods are proposed and compared for
parameter estimation and power analysis. An extension to multiple QTL is also
proposed in which the multiple QTL model provides improved resolution for QTL
inference. Extensive simulations are conducted to compare the performance of
the proposed model under different sampling designs with different combinations
of family and offspring size. Coparisons of the ML and REML methods, single QTL
and multiple QTL methods are discussed. The proposed method provides a general
framework in iQTL mapping with multiple line crosses and has significant
implications in real application.


## 2. Statistical Methods

### 2.1. Genetic Design

The dissection
of imprinting effects in line crosses depends on appropriate mating designs, where the allele parental origin can be traced and distinguished. Most commonly
used inbred line crosses are the backcross,F2, and recombinant
inbred line (RIL). Reciprocal backcross design has been proposed in iQTL
mapping [14, 16]. Considering parental origin of an allele, we use the
subscripts m and f to refer to an
allele inherited from the maternal and paternal parents, respectively. The
merit of a backcross design is that two reciprocal heterozygotes in offsprings, Amaf and amAf, can be distinguished and their mean effects can be
estimated and tested to assess imprinting [14, 16]. While all individuals in an
F2 segregation
population share the same parental information, theoretically it is impossible
to distinguish the phenotypic distribution of Amaf and amAf without extra
information. Considering sex-specific recombination rates, Cui et al. [15]
recently developed an imprinting model by incorporating this information into
an interval mapping framework. No study has been reported to use RILs for iQTL
mapping.The methods proposed in Cui [14] and Cui et al. [16]
are fixed-effects QTL models where the effects of an iQTL are considered as
fixed. While only four backcross families are considered, when extending to
multiple backcross families, the inference of iQTL variance calculation is less
efficient. The variance components method, initially proposed in human linkage
analysis [19], offers a powerful alternative in assessing genomic imprinting
[7]. In this paper, we will extend the variance components method to inbred
line populations by combining different backcross lines to map iQTL.A typical backcross design often starts with the cross
between one of the parental lines and theirF1 progeny to
create a segregation population. Then large number of offsprings are collected
for QTL mapping. When imprinting effect is considered, reciprocal backcrosses
are needed. A basic design framework is illustrated in Table 1 in Cui [14]. The
two reciprocal backcrosses are treated as the base mapping units. Multiple
backcross families are sampled based on these four crosses. For simplicity, we
sample equal number of families for each backcross category. For example, a
sample of 8 families would require two of each of the four backcrosses. Noted
that the variance components method assesses the degree of allele sharing among
siblings. When it is applied to inbred line crosses, each backcross population
is considered as one family and different families are considered as
independent. For fixed total sample size, one issue is to assess whether we
should sample large number of families each with small offspring size or small
number of families each with large offspring size. For example, to sample 400
individuals, shall we sample 4 backcross families each with 100 offsprings or
100 families each with 4 progenies or other sampling strategies? The choice of
optimal designs is intensively evaluated through simulations.Table 1
The IBD sharing coefficients for full-sib
pairs in a reciprocal backcross design considering allelic parental origin.
OffspringParent-specific
IBD sharingTotal
IBDBackcrossgenotypeπmmπffπm/fπQmQfQmqfQmQfQmqfQmQfQmqfQmQfQmqfQQ×QqQmQf0.50.50.5010.521Qmqf0.50.500.50.5011QmQfqmQfQmQfqmQfQmQfqmQfQmQfqmQfQq×QQQmQf0.500.50.510.521qmQf00.50.50.50.5011qmQfqmqfqmQfqmqfqmQfqmqfqmQfqmqfqq×QqqmQf0.50.50.5000.511qmqf0.50.500.50.5112QmqfqmqfQmqfqmqfQmqfqmqfQmqfqmqfQq×qqQmqf0.500.50.500.511qmqf00.50.50.50.5112


### 2.2. The Mixed-Effect Variance Components Model

Suppose there
is a putative QTL with two segregating allelesQ and q, located in an interval responsible for the variation
of a quantitative trait. The phenotype, yik, for individual i measured in
backcross family k(=1,…,K) can be written
as a linear function of QTL, polygene, and environmental
effects:(2.1)yik=μ+aik+Gik+eik,k=1,…,K;i=1,…,nk,where nk is the number
of offsprings in the kth backcross
family; μ denotes the
overall mean; aik is the random
additive effect of the major monogenic QTL assuming normal distribution with
mean zero; Gik is the
polygenic effect that reflects the effects of unlinked genes and is assumed to
be normally distributed with mean zero; eik~N(0,σe2) is the random
environmental error uncorrelated to other effects. The phenotypic variance-covariance
for the kth family can
be expressed as(2.2)Σk=Πkσa2+Φgσg2+Iσe2,where σa2 and σg2 are the
additive and polygene variances; Πk is a matrix
containing the proportion of marker alleles shared IBD for individuals in the kth backcross
family; Φg is a matrix of
the expected proportion of alleles shared IBD; I is the identity matrix. The calculation of the IBD sharing matrix with
inbred lines can be found in Xie et al. [18] which is based on the Malécot
coefficient of coancestry [20].Noted that a backcross offspring with genotypeQmqf may be obtained
by the QQ×Qq or the Qq×QQ cross. When
there is a significant maternal effect, the mean expression for genotype Qmqf may be
different depending on whether its maternal parents carrying QQ or Qq genotype. As
described in Cui [14], maternal effect is one source of potential confounding
factor for genomic imprinting. It should be appropriately modeled and adjusted
when testing imprinting. Here, we model the cytoplasmic maternal effects as
fixed effects, and the overall mean μ is replaced by μk which models
the maternal effect of the kth distinct
backcross family.To accommodate parent-of-origin effects, the QTL
additive effect (a) can be partitioned as two terms:(1)
a component
that reflects the influence of the QTL carried on the maternally derived
chromosome (am);(2)
a component that reflects the influence of the
QTL carried on the paternally derived chromosome (af).The model that accommodates the parent-specific
effects can be expressed as(2.3)yik=μk+aikm+aikf+Gik+eik,k=1,…,K;i=1,…,nk.
Note that the
proposed design contains three distinct maternal parent genotypes. Thus the k
maternal effects indexed by μk can be compressed to three distinct maternal
effects, instead of k terms. For data vector y in family k, the above model can be reexpressed
as(2.4)yk=Xkβ+akm+akf+Gk+ek,k=1,…,K,where Xk is an indicator
matrix corresponding to the kth backcross
family and β contains
parameters associated with the three maternal
effects; akm~N(0,Πm∣kσm2), akf~N(0,Πf∣kσf2), Gk~N(0,Φgσg2), ek~N(0,Iσe2), where Πm∣k and Πf∣k are matrices
containing the proportion of marker alleles shared IBD that are derived from
the mother and father, respectively; Φg is a matrix of
the expected proportion of alleles shared IBD; I is the identity matrix; σm2 and σf2 are the
variance of alleles inherited from the maternal and paternal parents, respectively.With noninbreeding mapping population, Hanson et al. [7] expressed the phenotypic variance-covariance for the kth family
as(2.5)Σk=Πm∣kσm2+Πf∣kσf2+Φgσg2+Iσe2.However, for an inbred mapping
population, this IBD-based variance partition method cannot be directly
applied. New method considering the inbreeding structure is needed.


### 2.3. Parent-Specific Allele Sharing and Covariances between Two Inbreeding Full-Sibs

Before we get
the phenotypic variance-covariance of a pair of individualsi and j, let us first consider the parent-specific allele
sharing status. Within each BC family, there are two alleles segregating at
each locus. Because of inbreeding, the IBD values between two backcross
individuals are different from those calculated from outbred full-sibs. Consider two sibs i and j in the kth backcross
family. Without considering allelic parental origin, Xie et al. [18] proposed
to calculate the IBD value at a QTL as(2.6)πij=2θij={2forQQ−QQ,1forQQ−QqorQq−Qqwith θij being the
Malécot coefficient of coancestry [20]. Thus, for an inbred population, πij is not the
actual IBD value between individuals i and j, rather interpreted as twice the coefficient of
coancestry [18, 21]. For individuals with itself, (2.7)πii=1+Fi={2forQQ−QQ,1forQq−Qq,where Fi is the
inbreeding coefficient for individual i at the QTL. The
elements in Φg matrix are just
the expected values of πij and πii which are ϕij=5/4 and ϕii=3/2 [18].When allelic parental origin is considered, the IBD
sharing matrix can also be calculated based on the coefficient of coancestry. By definition, the coefficient of coancestry is defined as the probability that
two randomly drawn alleles from individualsi and j are identical
by descent. Figure 1 displays possible alleles shared IBD for sibs drawn in
backcross families. Consider two backcross individuals i (with two
alleles Aim and Aif) and j (with two
alleles Ajm and Ajf). Define θij as the
coefficient of coancestry between individuals i and j. By definition, θij can be
calculated as(2.8)θij=14{Pr(Aim=Ajm)+Pr(Aim=Ajf)+Pr(Aif=Ajm)+Pr(Aif=Ajf)}=14(θimjm+θimjf+θifjm+θifjf),where
θi⋅j⋅
can be
interpreted as the allelic kinship coefficient, that is, the probability that a
randomly chosen allele from individual i is IBD to a
randomly chosen allele from individual j. Note that the two terms θimjf and θifjm are not
distinguishable. However, their sum is unique and therefore the two terms can
be combined as one single term, denoted as θim/jf(=θimjf+θifjm). After the manipulation, the coefficient of coancestry
for individuals i and j can be
expressed as θij=(1/4)(θimjm+θim/jf+θifjf) which is composed of three components.Figure 1
Possible alleles-shared IBD for
individualsi and j in inbreeding
backcross families. Lines indicate alleles-shared IBD.Following Xie et al. [18], the alleles shared IBD
between individuals i and j can be
expressed as(2.9)πij=2θij=12(θimjm+θim/jf+θifjf)=πimjm+πim/jf+πifjf,where πimjm=(1/2)θimjm and πifjf=(1/2)θifjf are the alleles
shared IBD derived from the mother and father, respectively; πim/jf=(1/2)θim/jf is the alleles
shared IBD due to alleles cross sharing, a special case for inbreeding sibs. Without inbreeding, πim/jf takes value of
zero.For completely inbreeding population, the inbreeding
coefficientFi is 1 if alleles
inherited from both parents are the same since these alleles can be traced back
to the same grandparent. For example, for an individual with genotype QmQf, Pr(Qm=Qf)=1 since both
alleles Qm and Qf are inherited
from the same grandparent. Therefore, for individuals with itself, πii=1+Fi would be the
same as πij(i≠j) when i and j carry the same
genotypes. The expected proportion of alleles shared IBD ϕij can also be
calculated.Thus, the proportion of alleles-shared IBD can be
partitioned as three components for inbreeding sibs, rather than two components
considering parent-of-origin effects proposed by Hanson et al. [7]. To further
illustrate the idea, we use one backcross family to demonstrate the derivation. A full list of possible IBD sharing values for the two reciprocal backcrosses is
given in Table 1. Considering a backcross family initiated with the Qq×QQ cross. Randomly
selecting two individuals i and j with genotype QmQf and QmQf, the Malécot coefficient of coancestry can be
calculated as(2.10)πij=2θij=12{Pr(Qim=Qjm)+Pr(Qim=Qjf)+Pr(Qif=Qjm)+Pr(Qif=Qjf)}=12[1+1+1+1]=2.Thus, πimjm=πifjf=0.5 and πim/jf=1. For sib pairs i (with genotype QmQf) and j (with genotype Qmqf), πimjm=0.5, πifjf=0 and πim/jf=0.5, and πij=1 which is the
same as given in (2.6) without considering parent-of-origin partition.Considering the allelic sharing status in a complete
inbreeding population, the relationship between the maternal and paternal
alleles is no longer independent if the two alleles are in identical form. There exists a covariance term (denoted asσmf2) due to
alleles cross sharing for two inbreeding full-sibs when calculating the phenotypic
variance. Corresponding to the partition of the IBD-sharing considering allelic
parental origin, the major QTL additive variance component can be partitioned
into three components, that is, σf2, σm2, and σmf2, in which σmf2 can be
interpreted as the covariance due to alleles cross sharing in inbreeding
families. Thus, the trait covariance between two individuals i and j can be
expressed as(2.11)Cov(yi,yj)=πimjmσm2+πifjfσf2+πim/jfσmf2+ϕijσg2+Iijσe2,where Iij is an indicator
variable taking value 1 if i=j and 0 if i≠j. The variance-covariance matrix for a phenotypic
vector in the kth backcross
family can then be expressed as(2.12)Σk=Πm∣kσm2+Πm/f∣kσmf2+Πf∣kσf2+Φgσg2+Iσe2,where the elements of Πm∣k, Πf∣k, and Πm/f∣k can be found in
Table 1.For noninbreeding sib pairs with random mating,πim/jf=0 and hence Cov(am,af)=0. Model (2.12) reduces to ∑k=Πm∣kσm2+Πf∣kσf2+Φgσg2+Iσe2, the same as the variance components partition model
considering parent-of-origin effects given by Hanson et al. [7].


### 2.4. Likelihood Function and Parameter Estimation

Assuming
multivariate normality, the density function of observing a particular vector
of datay for family k is given
by(2.13)f(yk;μk,Σk)=1(2π)nk/2|Σk|1/2exp[−12(yk−μk)TΣk−1(yk−μk)],where yk=(y1k,…,ynkk)T is an nk×1 vector of
phenotypes for the kth backcross
family, and nk is the kth backcross
family size. The overall log likelihood function for K independent
backcross families is give by(2.14)ℓ=∑k=1Klog[f(yk;μk,Σk)].Note that the maternal effectμk is the same for
families with the same maternal genotype. Thus, only three maternal effects
need to be estimated. Two commonly used methods can be applied to estimate
parameters in a mixed effects model, the ML method and the REML method. Both
methods have been applied in genetic linkage analysis in a variance components
model framework [19, 22]. In general, ML estimators tend to be downwardly
biased given that it does not account for the loss in degrees of freedom
resulted from estimation of the fixed effects [23]. The REML is based on a
linear transformation of the data such that the fixed effects are eliminated
from the model, hence it provides less-biased estimators. Even though standard
softwares such as SAS have standard procedures to estimate parameters for a
mixed effects model, the estimation for the proposed model cannot be directly
fitted into a standard software. The estimation procedures for the two methods
are detailed here.


#### 2.4.1. The Ml Estimation

The phenotype
vector in thekth backcross
family follows a multivariate normal distribution, that is, yk~MVN(Xkβ,Σk). Parameters that need to be estimated are Ω=(β,σm2,σf2,σmf2,σg2,σe2) with β=(μ1,μ2,μ3)T.Defineσ2=σm2+σf2+σmf2+σg2+σe2, hm2=(σm2/σ2), hf2=(σf2/σ2), hmf2=(σmf2/σ2), hg2=(σg2/σ2), and he2=1−hm2−hf2−hmf2−hg2. σ2 is the total
phenotypic variance and hence hm2 and hf2 can be
considered as the heritability of maternal and paternal alleles, hm2+hf2+hmf2 is the total
genetic heritability due to the major QTL, hg2 is the polygene
heritability, and h2=hm2+hf2+hmf2+hg2 is the overall
heritability. The phenotypic variance-covariance between any two individuals i and j in the kth backcross
family can then be reexpressed as(2.15)Var(yikyjk)=σ2Hij∣k,  where(2.16)Hij∣k=(δiδijδijδj)with δi=πimimhm2+πififhf2+πim/ifhmf2+ϕiihg2+he2; δj is defined
similarly; δij=πimjmhm2+πifjfhf2+πim/jfhmf2+ϕijhg2.If there arenk sibs in each
backcross family, Hk={Hij∣k}nk×nk is simply an nk×nk matrix. Instead
of estimating Ω=(β,σm2,σf2,σmf2,σg2,σe2), we can estimate Ω=(β,σ2,hm2,hf2,hmf2,hg2) and solve
above equations to get the original variance estimates. Now, the log-likelihood
can be expressed as(2.17)ℓ(Ω)=∑k=1Klog[f(yk∣Ω)]∝−∑k=1K{nk2logσ2−12log|Hk|−12σ2(yk−Xkβ)′Hk−1(yk−Xkβ)}.Maximizing likelihood (2.14) is equivalent to maximize
(2.17). Here, we take an iterated estimation procedure to estimate the parameters
contained in Ω. For given values of hm2,hf2,hmf2,hg2, we can get the maximum likelihood estimates (MLEs)
of parameters (β,σ2) by setting the partial derivative of the
log-likelihood function (2.17) to zero, that is, (2.18)β^=∑k=1K(XkTHk−1Xk)−1(XkTHk−1yk),σ^2=1∑k=1Knk∑k=1K(yk−Xkβ^)THk−1(yk−Xkβ^).It can be seen thatβ^ and σ^2 are functions
of hm2, hf2, hmf2, and hg2. Plug the updated parameter values for β and σ2 into likelihood
equation (2.17), the log-likelihood function can be simplified as(2.19)ℓ(Ω)=∑k=1Klog[f(yk∣Ω)]∝−∑k=1Knk2logσ^2−12∑k=1Klog|Hk|.The simplex algorithm [24] can be applied to maximize
the function (2.19) with respect to parameters hm2, hf2, hmf2 and hg2 subject to the
the constraints that 0≤hm2,hf2,hmf2,hg2≤1 and 0≤h2≤1.To guarantee a positive definite covariance matrix
when searching for these heritability values over the constraint parameter
space, a reparameterization technique is adopted [25]. Taking δij=h2(πimjmγm2+πifjfγf2+πim/jfγmf2+ϕijγg2), where γm2=hm2/h2,γf2=hf2/h2,γmf2=hmf2/h2, γg2=hg2/h2, and h2=hm2+hf2+hmf2+hg2. We now have four new unknowns with the constraints:0≤h2≤1, γm2+γf2+γmf2+γg2=1, and γm2,γf2,γmf2,γg2≥0.The new constraints can be easily satisfied by a
reparameterization technique. Letu, vm, vf, vmf, and vg be any real
numbers. Estimating h2, γm2, γf2, γmf2, and γg2 can be done by
maximizing the likelihood function (2.19) via searching through the real domain
space with respect to u, vm, vf, vmf, and vg with the
reparameterization(2.20)h2=eu1+eu,γm2=evmevm+evf+evmf+evg,γf2=evfevm+evf+evmf+evg,γmf2=evmfevm+evf+evmf+evg,γg2=evgevm+evf+evmf+evg.MLEs ofh2, γm2, γf2, γmf2, and γg2 can be obtained
through the estimated values for u, vm, vf, vmf, and vg according to
the invariance property of MLEs. These estimated MLEs are used to update h2, hm2, hf2, hmf2, and hg2, and hence σ2 and β. The iteration steps continue until converge.


#### 2.4.2. The Reml Estimation

The REML method
was first proposed by Patterson and Thompson [26]. This method has been broadly
applied to estimate variance components in a mixed-effect model framework. Taking Ω=(β,Θ), where Θ=(σm2,σf2,σmf2,σg2,σe2). The REML
method starts with maximizing the following likelihood
function:(2.21)ℓ*(Θ)=∑k=1Klog[f(yk∣Θ)]=−12∑k=1K{log|Σk|+log(|Xk′Σk−1Xk|)+yk′Pkyk},where Pk=Σk−1−Σk−1Xk(Xk′Σk−1Xk)−1Xk′Σk−1. We can combine all family data together as one N×1 vector denoted
as y, where N=∑k=1Knk. All the Xk and the
variance-covariance matrix ∑k corresponding
to each family can be combined. The log-likelihood function for the combined
data is expressed as(2.22)ℓ*(Θ)=log[f(y∣Θ)]=−12{log|Σ|+log(|X′Σ−1X|)+y'Py},where ∑ is a block
diagonal matrix with the kth diagonal
block ∑k corresponding
to the kth family and
off-diagonal blocks being zeros; P is also a block
diagonal matrix with block elements given by Pk. The dimension of ∑ is N×N. With this combination, we develop the following REML
estimation procedure.We apply the Fisher scoring algorithm to estimate the
unknowns, which has the form(2.23)Θ(t+1)=Θ(t)+I(Θ(t))−1∂ℓ*(Θ)∂Θ∣Θ(t),where I(Θ(t)) is the Fisher
information matrix evaluated at Θ(t) which can be
expressed as (2.24)I(σm2σf2σmf2σg2σe2)=12(tr(PΠmPΠm)tr(PΠmPΠf)tr(PΠmPΠm/f)tr(PΠmPΦg)tr(PΠmP)tr(PΠfPΠm)tr(PΠfPΠf)tr(PΠfPΠm/f)tr(PΠfPΦg)tr(PΠfP)tr(PΠm/fPΠm)tr(PΠm/fPΠf)tr(PΠm/fPΠm/f)tr(PΠm/fPΦg)tr(PΠm/fP)tr(PΦgPΠm)tr(PΦgPΠf)tr(PΦgPΠm/f)tr(PΦgPΦg)tr(PΦgP)tr(PΠmP)tr(PΠfP)tr(PΠm/fP)tr(PΦgP)tr(PP)).The first derivative of the log-likelihood functionℓ* with respective
to each variance components is given by(2.25)∂ℓ*∂σm2=−12(tr(PΠm)−yTPΠmPy),∂ℓ*∂σf2=−12(tr(PΠf)−yTPΠfPy),∂ℓ*∂σmf2=−12(tr(PΠm/f)−yTPΠm/fPy),∂ℓ*∂σg2=−12(tr(PΦg)−yTPΦgPy),∂ℓ*∂σe2=−12(tr(PIN)−yTPPy).The REML estimator of β is the
generalized least squares estimator, that is, (2.26)β^=(XTΣ^−1X)−1XTΣ^−1Y.Under the current mapping framework, the likelihood
function is very complex and no global maxima is guaranteed. Thus, initial
values are very important. In both ML and REML estimation procedures, the
estimated values under the null are set as the initial parameters for the
alternative model. The additional variance component(s) to be estimated under
the alternative model is(are) set to a small positive number. In this way, we
can guarantee that the alternative model always produces larger likelihood
values than the null model does, and hence positive likelihood ratio value.


### 2.5. QTL IBD Sharing and Genome-Wide Linkage Scan

The above IBD
computation procedure assumes that a putative QTL is located right on a marker. When a QTL is located within an interval, a more efficient approach would be to
do an interval scan and to test the imprinting property of QTLs at positions
across the entire linkage group. Under the proposed framework, essentially, we
need to estimate the proportion of putative QTL alleles shared IBD at every
genome position. Here, we propose a method to calculate QTL alleles-shared IBD
inside an interval conditional on the flanking markers. The so-called expected
conditional IBD values can be derived at each test position as a function of
recombination fraction between the two flanking markers, and the one between
one flanking marker and the QTL.We use one backcross initiated with the crossQQ×Qq as an example
to illustrate the idea. For a putative QTL with two alleles Q and q, four QTL genotype pairs QQ−QQ, QQ−Qq, Qq−QQ, and Qq−Qq can be formed. If the QTL genotype is observed, the corresponding QTL alleles-shared IBD can
be calculated (see Table 1). In general, the QTL genotype is unobservable, but
its conditional distribution can be calculated from the two flanking markers. For individuals i and j with flanking
marker genotypes gi and gj, let πv∣GiGj be the IBD
values calculated at the QTL position between individual i carrying QTL
genotype Gi(=1  or 2
corresponding to QQ or Qq, resp.) and individual j carrying
genotype Gj (similarly 1 or
2), where v=imjm,ifjf or im/jf. For example, πimjm∣GiGj is the
proportion of IBD sharing between individual i carrying QTL
genotype Gi and individual j carrying
genotype Gj for alleles
derived from the mother.LetφGi∣gi and φGj∣gj be the
conditional distribution of QTL genotype Gi and Gj for individuals i and j given on the
flanking markers gi and gj, respectively. This conditional probabilities can be
easily calculated and can be found at standard QTL mapping literature (see
[27]). The probability to observe πv∣GiGj is just φGi∣giφGj∣gj. Thus, the expected IBD values between
individuals i and j at the tested
QTL position conditioning on the flanking markers gi and gj can be
calculated as π^v = E(πv∣GiGj) = ∑Gi=12∑Gj=12πv∣GiGjφGi∣giφGj∣gj. For the above example, the IBD values derived from
the maternal and paternal parents can be calculated as π^imjm = E(πimjm∣GiGj) = 0.5 φ1∣giφ1∣gj + 0.5 φ1∣giφ2∣gj +0.5 φ2∣giφ1∣gj + 0.5 φ2∣giφ2∣gj and π^ifjf = E(πifjf∣GiGj) = 0.5 φ1∣giφ1∣gj + 0.5 φ2∣giφ2∣gj. Similarly, we can calculate the conditional
expectation of IBD sharing for other backcross families.SinceφGi∣gi and φGj∣gj are functions
of recombinations, the conditional QTL IBD values vary at different testing
positions. Once the estimated IBD matrix is calculated at every 1 or
2 cm on an interval bracketed by two markers
throughout the entire genome, a grid search can be done at all testing
positions. The amount of support for a QTL at a particular map position can be
displayed graphically through the use of likelihood ratio profiles, which plot
the likelihood ratio test statistic as a function of testing positions of
putative QTL (see details in Section 2.6). The peaks of the profile plot that
pass certain significant threshold correspond
to the positions of significant QTL.


### 2.6. Hypothesis Testing

With the
estimated parameters using either the ML or REML method, we are interested in
testing the existence of QTL across the genome and assess their imprinting
mechanism. The first hypothesis is to test the existence of major QTL, termed
overall QTL test, which can be formulated as(2.27)H0:σm2=σf2=σmf2=0,H1:atleastoneparameterisnotzero.Likelihood ratio (LR) test is applied which is computed between the full (there
is a QTL) and the reduced model (there is no QTL) corresponding toH1 and H0, respectively. Let Ω˜ and Ω^ be the
estimates of the unknown parameters under H0 and H1, respectively. The log-likelihood ratio can be
calculated as(2.28)LR1=−2[logL(Ω˜∣y)−logL(Ω^∣y)].When testing the hypothesis, the polygene and the
residual variances are nuisance parameters which are constrained to be nonnegative. The three tested genetic variance components under the null are lied on the
boundaries of their alternative parameter spaces. Following Self and Liang
[28], when the null is true, LR1 may
asymptotically follows a mixture of χ2 distribution on 0,…, 3 degrees of
freedom (df) with the mixture proportion for the χk2 components
being (3k)2−3. The theoretical distribution can be used to assess
significance in linkage scan. However, since there are many point tests across
the genome, the pointwise significance value may not guarantee an appropriate
genome-wide error rate. Another approach to assess significance is to use
nonparametric permutation tests in which the critical threshold value can be
empirically calculated on the basis of repeatedly shuffling the relationships
between marker genotypes and phenotypes [29]. In simulation studies, we also
simulate the null distribution and compare it with the theoretical
distribution.For those detected QTL, the next step is to assess
their imprinting property. An identified QTL can be imprinted, completely
imprinted, partially imprinted, or not imprinted at all. These can be tested
through the following sequential tests. The first imprinting test is to assess
whether a QTL shows imprinting effect, which can be done by formulating the
following hypotheses:(2.29)H0:σm2=σf2=σ2,H1:σm2≠σf2.Rejection of H0 provides
evidence of genomic imprinting and the QTL is called iQTL. Again, likelihood
ratio test can be applied in which the log-likelihood ratio test statistics
asymptotically follows an χ2 with one df
[7]. We denote the log-likelihood ratio test statistic as LRimp. If the null is rejected, one would be interested to
test if the detected iQTL is completely maternally or paternally imprinted. The
corresponding hypotheses can be formulated as(2.30)H0:σm2=0,H1:σm2≠0for testing completely maternal
imprinting and(2.31)H0:σf2=0,H1:σf2≠0for testing completely paternal
imprinting. The likelihood ratio test statistics for the above two tests
asymptotically follow a 50 : 50 mixture of χ02 and χ12
distributions
[28]. Rejection of complete imprinting indicates partial imprinting.


### 2.7. Multiple QTL Model

In reality, more than one QTL may contribute to the phenotypic variation located in one
chromosome region or across the whole genome. The polygenic effect in model (2.4)
absorbs the effects of multiple QTL located on other chromosomes. However, when
there are multiple QTL located on the same linkage group as the tested QTL, if
their effects are not properly adjusted, the estimation could be biased due to
interference caused by these QTL outside of the
testing interval [3, 30–32]. A multiple QTL model that can test the putative
QTL effect while adjusting the effects of interference QTL deserves more
attention.Zeng [32] previously showed that IBD variables share
the same property as the indicator variables in which the shared proportion of
alleles IBD for a QTL conditional on the IBD of one flanking marker is
independent of that of a QTL on the other side of that flanking marker. Thus, conditional on one flanking marker, the interference of QTL located on the other
side of the marker can be eliminated. By conditional on the IBD of the flanking
markers, the IBD sharing of a QTL is uncorrelated with that outside this
interval. Xu and Atchley [25] showed that one marker is enough to block the
interference caused by other QTL located on the same linkage group. The authors
derived the next-to-flanking markers structure to block additional QTL effects
from both sides of testing region in one chromosome. We derive a multiple QTL
model adopting a similar idea as Xu and Atchley [25]. Assume there are total S QTL located on
a linkage group. Considering parent-specific allelic effects, the multiple QTL
model can be expressed in general as(2.32)yik=μk+∑s=1Saikms+∑s=1Saikfs+Gik+eik,k=1,…,K;i=1,…,nk.In an interval-based linkage scan, only one putative
QTL is considered at each testing position conditioning on the effects of all
other QTL. Assuming there are totalL and R QTL located on
the left and right sides of the putative QTL on a
linkage group, model (2.32) can be modified as(2.33)yik=μk+∑l=1Laikl+(aikm+aikf)+∑r=1Raikr+Gik+eik,k=1,…,K;i=1,…,nk,where aikl and aikr are the lth and rth QTL random
effects on the left and right sides of the putative
QTL, respectively. When testing the putative QTL effect, we are only interested
in blocking the total effects of QTL outside of the tested interval. Therefore, in the modified model, the effects of QTL outside of the tested interval are
not partitioned. This, however, does not affect the inference of the tested
QTL.As shown by Zeng [32] and Jansen [31, 33], one marker
is enough to block the correlation between a locus on its left and a locus on
its right. Therefore, only two additional markers flanking the current interval
are needed to block interference caused by outside QTL [25]. Let ℳl and ℳr denote two
flanking markers for the tested interval, and ℒ and ℛ denote the two
markers next to ℳl and ℳl+1 with the marker
order ℒ - ℳl - ℳl+1 - ℛ. With the modified model given in (2.33), the
covariance of phenotypes between individuals i and j in the kth backcross
family can be expressed as(2.34)Cov(yik,yjk)=∑l=1LCov(aikl,ajkl)+Cov(aikm,ajkm)+Cov(aikf,ajkf)+Cov(aikm,ajkf)+∑r=1RCov(aikr,ajkr)+ϕijσg2+Iijσe2=∑l=1Lπl∣kσl2+πimjmσm2+πim/jfσmf2+πifjfσf2+∑r=1Rπr∣kσr2+ϕijσg2+Iijσe2,where πl∣k and πr∣k are the IBD
values for QTL located on the left and right sides
of the putative QTL in the kth backcross
family, and can be calculated following (2.6) and (2.7) if their genotype
information is known. Unfortunately, the number and exact locations of QTL
outside the testing interval are unknown. Hence, πl∣k and πr∣k are not
observable. Xu and Atchley [25] showed that when πl∣k and πr∣k are unknown, they can be estimated by some composite terms K(θlℒ,πℒ∣k) and K(θlℛ,πℛ∣k), where K(θlℒ,πℒ∣k) is a function
of the recombination fraction between the lth QTL and the
left marker ℒ as well as a
function of πℒ∣k, the IBD value for a pair of individuals at the left
marker ℒ. K(θlℛ,πℛ∣k) can be
similarly defined. Following Xu and Atchley [25], K(θlℒ,πℒ∣k) can be
expressed as a function of πℒ∣k multiplied by a
function of recombination frequency between the lth QTL and the
marker ℒ, f(θlℒ), that is, K(θlℒ,πℒ∣k)=πℒ∣kf(θlℒ). Similarly, K(θlℛ,πℛ∣k)=πℛ∣kf(θrℛ). When doing an interval scan, the covariance function
given in (2.34) between individuals i and j can be
reexpressed as(2.35)Cov(yik,yjk∣πℒ|k,π^imjm,π^im/jf,π^ifjf,πℛ|k)=∑l=1LK(θlℒ,πℒ∣k)σl2+π^imjmσm2+π^im/jf∣kσmf2+π^ifjfσf2+∑r=1RK(θlℛ,πℛ∣k)σr2+ϕijσg2+Iijσe2=πℒ∣k∑l=1Lf(θlℒ)σl2+π^imjmσm2+π^im/jfσmf2+π^ifjfσf2+πℛ∣k∑r=1Rf(θrℛ)σr2+ϕijσg2+Iijσe2=πℒ∣kσL2+π^imjmσm2+π^im/jf∣kσmf2+π^ififσf2+πℛ∣kσR2+ϕijσg2+Iijσe2.Instead of estimating individual
variance components σl2 and σr2, now we estimate the composite
terms
∑l=1Lf(θlℒ)σl2=σL2 and ∑r=1Rf(θrℛ)σr2=σR2. By conditioning the IBD sharing information for the
left and right markers ℒ and ℛ, the effects of those interference QTL are blocked. σL2 and σR2 absorb the
random effects of all QTL that are outside of the testing interval but are on
the same linkage group as the putative QTL. Estimation of the variance
components terms follows the same procedure as the single QTL analysis with
slight modification to consider multiple variance components.


## 2.1. Genetic Design

The dissection
of imprinting effects in line crosses depends on appropriate mating designs, where the allele parental origin can be traced and distinguished. Most commonly
used inbred line crosses are the backcross,F2, and recombinant
inbred line (RIL). Reciprocal backcross design has been proposed in iQTL
mapping [14, 16]. Considering parental origin of an allele, we use the
subscripts m and f to refer to an
allele inherited from the maternal and paternal parents, respectively. The
merit of a backcross design is that two reciprocal heterozygotes in offsprings, Amaf and amAf, can be distinguished and their mean effects can be
estimated and tested to assess imprinting [14, 16]. While all individuals in an
F2 segregation
population share the same parental information, theoretically it is impossible
to distinguish the phenotypic distribution of Amaf and amAf without extra
information. Considering sex-specific recombination rates, Cui et al. [15]
recently developed an imprinting model by incorporating this information into
an interval mapping framework. No study has been reported to use RILs for iQTL
mapping.The methods proposed in Cui [14] and Cui et al. [16]
are fixed-effects QTL models where the effects of an iQTL are considered as
fixed. While only four backcross families are considered, when extending to
multiple backcross families, the inference of iQTL variance calculation is less
efficient. The variance components method, initially proposed in human linkage
analysis [19], offers a powerful alternative in assessing genomic imprinting
[7]. In this paper, we will extend the variance components method to inbred
line populations by combining different backcross lines to map iQTL.A typical backcross design often starts with the cross
between one of the parental lines and theirF1 progeny to
create a segregation population. Then large number of offsprings are collected
for QTL mapping. When imprinting effect is considered, reciprocal backcrosses
are needed. A basic design framework is illustrated in Table 1 in Cui [14]. The
two reciprocal backcrosses are treated as the base mapping units. Multiple
backcross families are sampled based on these four crosses. For simplicity, we
sample equal number of families for each backcross category. For example, a
sample of 8 families would require two of each of the four backcrosses. Noted
that the variance components method assesses the degree of allele sharing among
siblings. When it is applied to inbred line crosses, each backcross population
is considered as one family and different families are considered as
independent. For fixed total sample size, one issue is to assess whether we
should sample large number of families each with small offspring size or small
number of families each with large offspring size. For example, to sample 400
individuals, shall we sample 4 backcross families each with 100 offsprings or
100 families each with 4 progenies or other sampling strategies? The choice of
optimal designs is intensively evaluated through simulations.Table 1
The IBD sharing coefficients for full-sib
pairs in a reciprocal backcross design considering allelic parental origin.
OffspringParent-specific
IBD sharingTotal
IBDBackcrossgenotypeπmmπffπm/fπQmQfQmqfQmQfQmqfQmQfQmqfQmQfQmqfQQ×QqQmQf0.50.50.5010.521Qmqf0.50.500.50.5011QmQfqmQfQmQfqmQfQmQfqmQfQmQfqmQfQq×QQQmQf0.500.50.510.521qmQf00.50.50.50.5011qmQfqmqfqmQfqmqfqmQfqmqfqmQfqmqfqq×QqqmQf0.50.50.5000.511qmqf0.50.500.50.5112QmqfqmqfQmqfqmqfQmqfqmqfQmqfqmqfQq×qqQmqf0.500.50.500.511qmqf00.50.50.50.5112


## 2.2. The Mixed-Effect Variance Components Model

Suppose there
is a putative QTL with two segregating allelesQ and q, located in an interval responsible for the variation
of a quantitative trait. The phenotype, yik, for individual i measured in
backcross family k(=1,…,K) can be written
as a linear function of QTL, polygene, and environmental
effects:(2.1)yik=μ+aik+Gik+eik,k=1,…,K;i=1,…,nk,where nk is the number
of offsprings in the kth backcross
family; μ denotes the
overall mean; aik is the random
additive effect of the major monogenic QTL assuming normal distribution with
mean zero; Gik is the
polygenic effect that reflects the effects of unlinked genes and is assumed to
be normally distributed with mean zero; eik~N(0,σe2) is the random
environmental error uncorrelated to other effects. The phenotypic variance-covariance
for the kth family can
be expressed as(2.2)Σk=Πkσa2+Φgσg2+Iσe2,where σa2 and σg2 are the
additive and polygene variances; Πk is a matrix
containing the proportion of marker alleles shared IBD for individuals in the kth backcross
family; Φg is a matrix of
the expected proportion of alleles shared IBD; I is the identity matrix. The calculation of the IBD sharing matrix with
inbred lines can be found in Xie et al. [18] which is based on the Malécot
coefficient of coancestry [20].Noted that a backcross offspring with genotypeQmqf may be obtained
by the QQ×Qq or the Qq×QQ cross. When
there is a significant maternal effect, the mean expression for genotype Qmqf may be
different depending on whether its maternal parents carrying QQ or Qq genotype. As
described in Cui [14], maternal effect is one source of potential confounding
factor for genomic imprinting. It should be appropriately modeled and adjusted
when testing imprinting. Here, we model the cytoplasmic maternal effects as
fixed effects, and the overall mean μ is replaced by μk which models
the maternal effect of the kth distinct
backcross family.To accommodate parent-of-origin effects, the QTL
additive effect (a) can be partitioned as two terms:(1)
a component
that reflects the influence of the QTL carried on the maternally derived
chromosome (am);(2)
a component that reflects the influence of the
QTL carried on the paternally derived chromosome (af).The model that accommodates the parent-specific
effects can be expressed as(2.3)yik=μk+aikm+aikf+Gik+eik,k=1,…,K;i=1,…,nk.
Note that the
proposed design contains three distinct maternal parent genotypes. Thus the k
maternal effects indexed by μk can be compressed to three distinct maternal
effects, instead of k terms. For data vector y in family k, the above model can be reexpressed
as(2.4)yk=Xkβ+akm+akf+Gk+ek,k=1,…,K,where Xk is an indicator
matrix corresponding to the kth backcross
family and β contains
parameters associated with the three maternal
effects; akm~N(0,Πm∣kσm2), akf~N(0,Πf∣kσf2), Gk~N(0,Φgσg2), ek~N(0,Iσe2), where Πm∣k and Πf∣k are matrices
containing the proportion of marker alleles shared IBD that are derived from
the mother and father, respectively; Φg is a matrix of
the expected proportion of alleles shared IBD; I is the identity matrix; σm2 and σf2 are the
variance of alleles inherited from the maternal and paternal parents, respectively.With noninbreeding mapping population, Hanson et al. [7] expressed the phenotypic variance-covariance for the kth family
as(2.5)Σk=Πm∣kσm2+Πf∣kσf2+Φgσg2+Iσe2.However, for an inbred mapping
population, this IBD-based variance partition method cannot be directly
applied. New method considering the inbreeding structure is needed.


## 2.3. Parent-Specific Allele Sharing and Covariances between Two Inbreeding Full-Sibs

Before we get
the phenotypic variance-covariance of a pair of individualsi and j, let us first consider the parent-specific allele
sharing status. Within each BC family, there are two alleles segregating at
each locus. Because of inbreeding, the IBD values between two backcross
individuals are different from those calculated from outbred full-sibs. Consider two sibs i and j in the kth backcross
family. Without considering allelic parental origin, Xie et al. [18] proposed
to calculate the IBD value at a QTL as(2.6)πij=2θij={2forQQ−QQ,1forQQ−QqorQq−Qqwith θij being the
Malécot coefficient of coancestry [20]. Thus, for an inbred population, πij is not the
actual IBD value between individuals i and j, rather interpreted as twice the coefficient of
coancestry [18, 21]. For individuals with itself, (2.7)πii=1+Fi={2forQQ−QQ,1forQq−Qq,where Fi is the
inbreeding coefficient for individual i at the QTL. The
elements in Φg matrix are just
the expected values of πij and πii which are ϕij=5/4 and ϕii=3/2 [18].When allelic parental origin is considered, the IBD
sharing matrix can also be calculated based on the coefficient of coancestry. By definition, the coefficient of coancestry is defined as the probability that
two randomly drawn alleles from individualsi and j are identical
by descent. Figure 1 displays possible alleles shared IBD for sibs drawn in
backcross families. Consider two backcross individuals i (with two
alleles Aim and Aif) and j (with two
alleles Ajm and Ajf). Define θij as the
coefficient of coancestry between individuals i and j. By definition, θij can be
calculated as(2.8)θij=14{Pr(Aim=Ajm)+Pr(Aim=Ajf)+Pr(Aif=Ajm)+Pr(Aif=Ajf)}=14(θimjm+θimjf+θifjm+θifjf),where
θi⋅j⋅
can be
interpreted as the allelic kinship coefficient, that is, the probability that a
randomly chosen allele from individual i is IBD to a
randomly chosen allele from individual j. Note that the two terms θimjf and θifjm are not
distinguishable. However, their sum is unique and therefore the two terms can
be combined as one single term, denoted as θim/jf(=θimjf+θifjm). After the manipulation, the coefficient of coancestry
for individuals i and j can be
expressed as θij=(1/4)(θimjm+θim/jf+θifjf) which is composed of three components.Figure 1
Possible alleles-shared IBD for
individualsi and j in inbreeding
backcross families. Lines indicate alleles-shared IBD.Following Xie et al. [18], the alleles shared IBD
between individuals i and j can be
expressed as(2.9)πij=2θij=12(θimjm+θim/jf+θifjf)=πimjm+πim/jf+πifjf,where πimjm=(1/2)θimjm and πifjf=(1/2)θifjf are the alleles
shared IBD derived from the mother and father, respectively; πim/jf=(1/2)θim/jf is the alleles
shared IBD due to alleles cross sharing, a special case for inbreeding sibs. Without inbreeding, πim/jf takes value of
zero.For completely inbreeding population, the inbreeding
coefficientFi is 1 if alleles
inherited from both parents are the same since these alleles can be traced back
to the same grandparent. For example, for an individual with genotype QmQf, Pr(Qm=Qf)=1 since both
alleles Qm and Qf are inherited
from the same grandparent. Therefore, for individuals with itself, πii=1+Fi would be the
same as πij(i≠j) when i and j carry the same
genotypes. The expected proportion of alleles shared IBD ϕij can also be
calculated.Thus, the proportion of alleles-shared IBD can be
partitioned as three components for inbreeding sibs, rather than two components
considering parent-of-origin effects proposed by Hanson et al. [7]. To further
illustrate the idea, we use one backcross family to demonstrate the derivation. A full list of possible IBD sharing values for the two reciprocal backcrosses is
given in Table 1. Considering a backcross family initiated with the Qq×QQ cross. Randomly
selecting two individuals i and j with genotype QmQf and QmQf, the Malécot coefficient of coancestry can be
calculated as(2.10)πij=2θij=12{Pr(Qim=Qjm)+Pr(Qim=Qjf)+Pr(Qif=Qjm)+Pr(Qif=Qjf)}=12[1+1+1+1]=2.Thus, πimjm=πifjf=0.5 and πim/jf=1. For sib pairs i (with genotype QmQf) and j (with genotype Qmqf), πimjm=0.5, πifjf=0 and πim/jf=0.5, and πij=1 which is the
same as given in (2.6) without considering parent-of-origin partition.Considering the allelic sharing status in a complete
inbreeding population, the relationship between the maternal and paternal
alleles is no longer independent if the two alleles are in identical form. There exists a covariance term (denoted asσmf2) due to
alleles cross sharing for two inbreeding full-sibs when calculating the phenotypic
variance. Corresponding to the partition of the IBD-sharing considering allelic
parental origin, the major QTL additive variance component can be partitioned
into three components, that is, σf2, σm2, and σmf2, in which σmf2 can be
interpreted as the covariance due to alleles cross sharing in inbreeding
families. Thus, the trait covariance between two individuals i and j can be
expressed as(2.11)Cov(yi,yj)=πimjmσm2+πifjfσf2+πim/jfσmf2+ϕijσg2+Iijσe2,where Iij is an indicator
variable taking value 1 if i=j and 0 if i≠j. The variance-covariance matrix for a phenotypic
vector in the kth backcross
family can then be expressed as(2.12)Σk=Πm∣kσm2+Πm/f∣kσmf2+Πf∣kσf2+Φgσg2+Iσe2,where the elements of Πm∣k, Πf∣k, and Πm/f∣k can be found in
Table 1.For noninbreeding sib pairs with random mating,πim/jf=0 and hence Cov(am,af)=0. Model (2.12) reduces to ∑k=Πm∣kσm2+Πf∣kσf2+Φgσg2+Iσe2, the same as the variance components partition model
considering parent-of-origin effects given by Hanson et al. [7].


## 2.4. Likelihood Function and Parameter Estimation

Assuming
multivariate normality, the density function of observing a particular vector
of datay for family k is given
by(2.13)f(yk;μk,Σk)=1(2π)nk/2|Σk|1/2exp[−12(yk−μk)TΣk−1(yk−μk)],where yk=(y1k,…,ynkk)T is an nk×1 vector of
phenotypes for the kth backcross
family, and nk is the kth backcross
family size. The overall log likelihood function for K independent
backcross families is give by(2.14)ℓ=∑k=1Klog[f(yk;μk,Σk)].Note that the maternal effectμk is the same for
families with the same maternal genotype. Thus, only three maternal effects
need to be estimated. Two commonly used methods can be applied to estimate
parameters in a mixed effects model, the ML method and the REML method. Both
methods have been applied in genetic linkage analysis in a variance components
model framework [19, 22]. In general, ML estimators tend to be downwardly
biased given that it does not account for the loss in degrees of freedom
resulted from estimation of the fixed effects [23]. The REML is based on a
linear transformation of the data such that the fixed effects are eliminated
from the model, hence it provides less-biased estimators. Even though standard
softwares such as SAS have standard procedures to estimate parameters for a
mixed effects model, the estimation for the proposed model cannot be directly
fitted into a standard software. The estimation procedures for the two methods
are detailed here.


### 2.4.1. The Ml Estimation

The phenotype
vector in thekth backcross
family follows a multivariate normal distribution, that is, yk~MVN(Xkβ,Σk). Parameters that need to be estimated are Ω=(β,σm2,σf2,σmf2,σg2,σe2) with β=(μ1,μ2,μ3)T.Defineσ2=σm2+σf2+σmf2+σg2+σe2, hm2=(σm2/σ2), hf2=(σf2/σ2), hmf2=(σmf2/σ2), hg2=(σg2/σ2), and he2=1−hm2−hf2−hmf2−hg2. σ2 is the total
phenotypic variance and hence hm2 and hf2 can be
considered as the heritability of maternal and paternal alleles, hm2+hf2+hmf2 is the total
genetic heritability due to the major QTL, hg2 is the polygene
heritability, and h2=hm2+hf2+hmf2+hg2 is the overall
heritability. The phenotypic variance-covariance between any two individuals i and j in the kth backcross
family can then be reexpressed as(2.15)Var(yikyjk)=σ2Hij∣k,  where(2.16)Hij∣k=(δiδijδijδj)with δi=πimimhm2+πififhf2+πim/ifhmf2+ϕiihg2+he2; δj is defined
similarly; δij=πimjmhm2+πifjfhf2+πim/jfhmf2+ϕijhg2.If there arenk sibs in each
backcross family, Hk={Hij∣k}nk×nk is simply an nk×nk matrix. Instead
of estimating Ω=(β,σm2,σf2,σmf2,σg2,σe2), we can estimate Ω=(β,σ2,hm2,hf2,hmf2,hg2) and solve
above equations to get the original variance estimates. Now, the log-likelihood
can be expressed as(2.17)ℓ(Ω)=∑k=1Klog[f(yk∣Ω)]∝−∑k=1K{nk2logσ2−12log|Hk|−12σ2(yk−Xkβ)′Hk−1(yk−Xkβ)}.Maximizing likelihood (2.14) is equivalent to maximize
(2.17). Here, we take an iterated estimation procedure to estimate the parameters
contained in Ω. For given values of hm2,hf2,hmf2,hg2, we can get the maximum likelihood estimates (MLEs)
of parameters (β,σ2) by setting the partial derivative of the
log-likelihood function (2.17) to zero, that is, (2.18)β^=∑k=1K(XkTHk−1Xk)−1(XkTHk−1yk),σ^2=1∑k=1Knk∑k=1K(yk−Xkβ^)THk−1(yk−Xkβ^).It can be seen thatβ^ and σ^2 are functions
of hm2, hf2, hmf2, and hg2. Plug the updated parameter values for β and σ2 into likelihood
equation (2.17), the log-likelihood function can be simplified as(2.19)ℓ(Ω)=∑k=1Klog[f(yk∣Ω)]∝−∑k=1Knk2logσ^2−12∑k=1Klog|Hk|.The simplex algorithm [24] can be applied to maximize
the function (2.19) with respect to parameters hm2, hf2, hmf2 and hg2 subject to the
the constraints that 0≤hm2,hf2,hmf2,hg2≤1 and 0≤h2≤1.To guarantee a positive definite covariance matrix
when searching for these heritability values over the constraint parameter
space, a reparameterization technique is adopted [25]. Taking δij=h2(πimjmγm2+πifjfγf2+πim/jfγmf2+ϕijγg2), where γm2=hm2/h2,γf2=hf2/h2,γmf2=hmf2/h2, γg2=hg2/h2, and h2=hm2+hf2+hmf2+hg2. We now have four new unknowns with the constraints:0≤h2≤1, γm2+γf2+γmf2+γg2=1, and γm2,γf2,γmf2,γg2≥0.The new constraints can be easily satisfied by a
reparameterization technique. Letu, vm, vf, vmf, and vg be any real
numbers. Estimating h2, γm2, γf2, γmf2, and γg2 can be done by
maximizing the likelihood function (2.19) via searching through the real domain
space with respect to u, vm, vf, vmf, and vg with the
reparameterization(2.20)h2=eu1+eu,γm2=evmevm+evf+evmf+evg,γf2=evfevm+evf+evmf+evg,γmf2=evmfevm+evf+evmf+evg,γg2=evgevm+evf+evmf+evg.MLEs ofh2, γm2, γf2, γmf2, and γg2 can be obtained
through the estimated values for u, vm, vf, vmf, and vg according to
the invariance property of MLEs. These estimated MLEs are used to update h2, hm2, hf2, hmf2, and hg2, and hence σ2 and β. The iteration steps continue until converge.


### 2.4.2. The Reml Estimation

The REML method
was first proposed by Patterson and Thompson [26]. This method has been broadly
applied to estimate variance components in a mixed-effect model framework. Taking Ω=(β,Θ), where Θ=(σm2,σf2,σmf2,σg2,σe2). The REML
method starts with maximizing the following likelihood
function:(2.21)ℓ*(Θ)=∑k=1Klog[f(yk∣Θ)]=−12∑k=1K{log|Σk|+log(|Xk′Σk−1Xk|)+yk′Pkyk},where Pk=Σk−1−Σk−1Xk(Xk′Σk−1Xk)−1Xk′Σk−1. We can combine all family data together as one N×1 vector denoted
as y, where N=∑k=1Knk. All the Xk and the
variance-covariance matrix ∑k corresponding
to each family can be combined. The log-likelihood function for the combined
data is expressed as(2.22)ℓ*(Θ)=log[f(y∣Θ)]=−12{log|Σ|+log(|X′Σ−1X|)+y'Py},where ∑ is a block
diagonal matrix with the kth diagonal
block ∑k corresponding
to the kth family and
off-diagonal blocks being zeros; P is also a block
diagonal matrix with block elements given by Pk. The dimension of ∑ is N×N. With this combination, we develop the following REML
estimation procedure.We apply the Fisher scoring algorithm to estimate the
unknowns, which has the form(2.23)Θ(t+1)=Θ(t)+I(Θ(t))−1∂ℓ*(Θ)∂Θ∣Θ(t),where I(Θ(t)) is the Fisher
information matrix evaluated at Θ(t) which can be
expressed as (2.24)I(σm2σf2σmf2σg2σe2)=12(tr(PΠmPΠm)tr(PΠmPΠf)tr(PΠmPΠm/f)tr(PΠmPΦg)tr(PΠmP)tr(PΠfPΠm)tr(PΠfPΠf)tr(PΠfPΠm/f)tr(PΠfPΦg)tr(PΠfP)tr(PΠm/fPΠm)tr(PΠm/fPΠf)tr(PΠm/fPΠm/f)tr(PΠm/fPΦg)tr(PΠm/fP)tr(PΦgPΠm)tr(PΦgPΠf)tr(PΦgPΠm/f)tr(PΦgPΦg)tr(PΦgP)tr(PΠmP)tr(PΠfP)tr(PΠm/fP)tr(PΦgP)tr(PP)).The first derivative of the log-likelihood functionℓ* with respective
to each variance components is given by(2.25)∂ℓ*∂σm2=−12(tr(PΠm)−yTPΠmPy),∂ℓ*∂σf2=−12(tr(PΠf)−yTPΠfPy),∂ℓ*∂σmf2=−12(tr(PΠm/f)−yTPΠm/fPy),∂ℓ*∂σg2=−12(tr(PΦg)−yTPΦgPy),∂ℓ*∂σe2=−12(tr(PIN)−yTPPy).The REML estimator of β is the
generalized least squares estimator, that is, (2.26)β^=(XTΣ^−1X)−1XTΣ^−1Y.Under the current mapping framework, the likelihood
function is very complex and no global maxima is guaranteed. Thus, initial
values are very important. In both ML and REML estimation procedures, the
estimated values under the null are set as the initial parameters for the
alternative model. The additional variance component(s) to be estimated under
the alternative model is(are) set to a small positive number. In this way, we
can guarantee that the alternative model always produces larger likelihood
values than the null model does, and hence positive likelihood ratio value.


## 2.4.1. The Ml Estimation

The phenotype
vector in thekth backcross
family follows a multivariate normal distribution, that is, yk~MVN(Xkβ,Σk). Parameters that need to be estimated are Ω=(β,σm2,σf2,σmf2,σg2,σe2) with β=(μ1,μ2,μ3)T.Defineσ2=σm2+σf2+σmf2+σg2+σe2, hm2=(σm2/σ2), hf2=(σf2/σ2), hmf2=(σmf2/σ2), hg2=(σg2/σ2), and he2=1−hm2−hf2−hmf2−hg2. σ2 is the total
phenotypic variance and hence hm2 and hf2 can be
considered as the heritability of maternal and paternal alleles, hm2+hf2+hmf2 is the total
genetic heritability due to the major QTL, hg2 is the polygene
heritability, and h2=hm2+hf2+hmf2+hg2 is the overall
heritability. The phenotypic variance-covariance between any two individuals i and j in the kth backcross
family can then be reexpressed as(2.15)Var(yikyjk)=σ2Hij∣k,  where(2.16)Hij∣k=(δiδijδijδj)with δi=πimimhm2+πififhf2+πim/ifhmf2+ϕiihg2+he2; δj is defined
similarly; δij=πimjmhm2+πifjfhf2+πim/jfhmf2+ϕijhg2.If there arenk sibs in each
backcross family, Hk={Hij∣k}nk×nk is simply an nk×nk matrix. Instead
of estimating Ω=(β,σm2,σf2,σmf2,σg2,σe2), we can estimate Ω=(β,σ2,hm2,hf2,hmf2,hg2) and solve
above equations to get the original variance estimates. Now, the log-likelihood
can be expressed as(2.17)ℓ(Ω)=∑k=1Klog[f(yk∣Ω)]∝−∑k=1K{nk2logσ2−12log|Hk|−12σ2(yk−Xkβ)′Hk−1(yk−Xkβ)}.Maximizing likelihood (2.14) is equivalent to maximize
(2.17). Here, we take an iterated estimation procedure to estimate the parameters
contained in Ω. For given values of hm2,hf2,hmf2,hg2, we can get the maximum likelihood estimates (MLEs)
of parameters (β,σ2) by setting the partial derivative of the
log-likelihood function (2.17) to zero, that is, (2.18)β^=∑k=1K(XkTHk−1Xk)−1(XkTHk−1yk),σ^2=1∑k=1Knk∑k=1K(yk−Xkβ^)THk−1(yk−Xkβ^).It can be seen thatβ^ and σ^2 are functions
of hm2, hf2, hmf2, and hg2. Plug the updated parameter values for β and σ2 into likelihood
equation (2.17), the log-likelihood function can be simplified as(2.19)ℓ(Ω)=∑k=1Klog[f(yk∣Ω)]∝−∑k=1Knk2logσ^2−12∑k=1Klog|Hk|.The simplex algorithm [24] can be applied to maximize
the function (2.19) with respect to parameters hm2, hf2, hmf2 and hg2 subject to the
the constraints that 0≤hm2,hf2,hmf2,hg2≤1 and 0≤h2≤1.To guarantee a positive definite covariance matrix
when searching for these heritability values over the constraint parameter
space, a reparameterization technique is adopted [25]. Taking δij=h2(πimjmγm2+πifjfγf2+πim/jfγmf2+ϕijγg2), where γm2=hm2/h2,γf2=hf2/h2,γmf2=hmf2/h2, γg2=hg2/h2, and h2=hm2+hf2+hmf2+hg2. We now have four new unknowns with the constraints:0≤h2≤1, γm2+γf2+γmf2+γg2=1, and γm2,γf2,γmf2,γg2≥0.The new constraints can be easily satisfied by a
reparameterization technique. Letu, vm, vf, vmf, and vg be any real
numbers. Estimating h2, γm2, γf2, γmf2, and γg2 can be done by
maximizing the likelihood function (2.19) via searching through the real domain
space with respect to u, vm, vf, vmf, and vg with the
reparameterization(2.20)h2=eu1+eu,γm2=evmevm+evf+evmf+evg,γf2=evfevm+evf+evmf+evg,γmf2=evmfevm+evf+evmf+evg,γg2=evgevm+evf+evmf+evg.MLEs ofh2, γm2, γf2, γmf2, and γg2 can be obtained
through the estimated values for u, vm, vf, vmf, and vg according to
the invariance property of MLEs. These estimated MLEs are used to update h2, hm2, hf2, hmf2, and hg2, and hence σ2 and β. The iteration steps continue until converge.


## 2.4.2. The Reml Estimation

The REML method
was first proposed by Patterson and Thompson [26]. This method has been broadly
applied to estimate variance components in a mixed-effect model framework. Taking Ω=(β,Θ), where Θ=(σm2,σf2,σmf2,σg2,σe2). The REML
method starts with maximizing the following likelihood
function:(2.21)ℓ*(Θ)=∑k=1Klog[f(yk∣Θ)]=−12∑k=1K{log|Σk|+log(|Xk′Σk−1Xk|)+yk′Pkyk},where Pk=Σk−1−Σk−1Xk(Xk′Σk−1Xk)−1Xk′Σk−1. We can combine all family data together as one N×1 vector denoted
as y, where N=∑k=1Knk. All the Xk and the
variance-covariance matrix ∑k corresponding
to each family can be combined. The log-likelihood function for the combined
data is expressed as(2.22)ℓ*(Θ)=log[f(y∣Θ)]=−12{log|Σ|+log(|X′Σ−1X|)+y'Py},where ∑ is a block
diagonal matrix with the kth diagonal
block ∑k corresponding
to the kth family and
off-diagonal blocks being zeros; P is also a block
diagonal matrix with block elements given by Pk. The dimension of ∑ is N×N. With this combination, we develop the following REML
estimation procedure.We apply the Fisher scoring algorithm to estimate the
unknowns, which has the form(2.23)Θ(t+1)=Θ(t)+I(Θ(t))−1∂ℓ*(Θ)∂Θ∣Θ(t),where I(Θ(t)) is the Fisher
information matrix evaluated at Θ(t) which can be
expressed as (2.24)I(σm2σf2σmf2σg2σe2)=12(tr(PΠmPΠm)tr(PΠmPΠf)tr(PΠmPΠm/f)tr(PΠmPΦg)tr(PΠmP)tr(PΠfPΠm)tr(PΠfPΠf)tr(PΠfPΠm/f)tr(PΠfPΦg)tr(PΠfP)tr(PΠm/fPΠm)tr(PΠm/fPΠf)tr(PΠm/fPΠm/f)tr(PΠm/fPΦg)tr(PΠm/fP)tr(PΦgPΠm)tr(PΦgPΠf)tr(PΦgPΠm/f)tr(PΦgPΦg)tr(PΦgP)tr(PΠmP)tr(PΠfP)tr(PΠm/fP)tr(PΦgP)tr(PP)).The first derivative of the log-likelihood functionℓ* with respective
to each variance components is given by(2.25)∂ℓ*∂σm2=−12(tr(PΠm)−yTPΠmPy),∂ℓ*∂σf2=−12(tr(PΠf)−yTPΠfPy),∂ℓ*∂σmf2=−12(tr(PΠm/f)−yTPΠm/fPy),∂ℓ*∂σg2=−12(tr(PΦg)−yTPΦgPy),∂ℓ*∂σe2=−12(tr(PIN)−yTPPy).The REML estimator of β is the
generalized least squares estimator, that is, (2.26)β^=(XTΣ^−1X)−1XTΣ^−1Y.Under the current mapping framework, the likelihood
function is very complex and no global maxima is guaranteed. Thus, initial
values are very important. In both ML and REML estimation procedures, the
estimated values under the null are set as the initial parameters for the
alternative model. The additional variance component(s) to be estimated under
the alternative model is(are) set to a small positive number. In this way, we
can guarantee that the alternative model always produces larger likelihood
values than the null model does, and hence positive likelihood ratio value.


## 2.5. QTL IBD Sharing and Genome-Wide Linkage Scan

The above IBD
computation procedure assumes that a putative QTL is located right on a marker. When a QTL is located within an interval, a more efficient approach would be to
do an interval scan and to test the imprinting property of QTLs at positions
across the entire linkage group. Under the proposed framework, essentially, we
need to estimate the proportion of putative QTL alleles shared IBD at every
genome position. Here, we propose a method to calculate QTL alleles-shared IBD
inside an interval conditional on the flanking markers. The so-called expected
conditional IBD values can be derived at each test position as a function of
recombination fraction between the two flanking markers, and the one between
one flanking marker and the QTL.We use one backcross initiated with the crossQQ×Qq as an example
to illustrate the idea. For a putative QTL with two alleles Q and q, four QTL genotype pairs QQ−QQ, QQ−Qq, Qq−QQ, and Qq−Qq can be formed. If the QTL genotype is observed, the corresponding QTL alleles-shared IBD can
be calculated (see Table 1). In general, the QTL genotype is unobservable, but
its conditional distribution can be calculated from the two flanking markers. For individuals i and j with flanking
marker genotypes gi and gj, let πv∣GiGj be the IBD
values calculated at the QTL position between individual i carrying QTL
genotype Gi(=1  or 2
corresponding to QQ or Qq, resp.) and individual j carrying
genotype Gj (similarly 1 or
2), where v=imjm,ifjf or im/jf. For example, πimjm∣GiGj is the
proportion of IBD sharing between individual i carrying QTL
genotype Gi and individual j carrying
genotype Gj for alleles
derived from the mother.LetφGi∣gi and φGj∣gj be the
conditional distribution of QTL genotype Gi and Gj for individuals i and j given on the
flanking markers gi and gj, respectively. This conditional probabilities can be
easily calculated and can be found at standard QTL mapping literature (see
[27]). The probability to observe πv∣GiGj is just φGi∣giφGj∣gj. Thus, the expected IBD values between
individuals i and j at the tested
QTL position conditioning on the flanking markers gi and gj can be
calculated as π^v = E(πv∣GiGj) = ∑Gi=12∑Gj=12πv∣GiGjφGi∣giφGj∣gj. For the above example, the IBD values derived from
the maternal and paternal parents can be calculated as π^imjm = E(πimjm∣GiGj) = 0.5 φ1∣giφ1∣gj + 0.5 φ1∣giφ2∣gj +0.5 φ2∣giφ1∣gj + 0.5 φ2∣giφ2∣gj and π^ifjf = E(πifjf∣GiGj) = 0.5 φ1∣giφ1∣gj + 0.5 φ2∣giφ2∣gj. Similarly, we can calculate the conditional
expectation of IBD sharing for other backcross families.SinceφGi∣gi and φGj∣gj are functions
of recombinations, the conditional QTL IBD values vary at different testing
positions. Once the estimated IBD matrix is calculated at every 1 or
2 cm on an interval bracketed by two markers
throughout the entire genome, a grid search can be done at all testing
positions. The amount of support for a QTL at a particular map position can be
displayed graphically through the use of likelihood ratio profiles, which plot
the likelihood ratio test statistic as a function of testing positions of
putative QTL (see details in Section 2.6). The peaks of the profile plot that
pass certain significant threshold correspond
to the positions of significant QTL.


## 2.6. Hypothesis Testing

With the
estimated parameters using either the ML or REML method, we are interested in
testing the existence of QTL across the genome and assess their imprinting
mechanism. The first hypothesis is to test the existence of major QTL, termed
overall QTL test, which can be formulated as(2.27)H0:σm2=σf2=σmf2=0,H1:atleastoneparameterisnotzero.Likelihood ratio (LR) test is applied which is computed between the full (there
is a QTL) and the reduced model (there is no QTL) corresponding toH1 and H0, respectively. Let Ω˜ and Ω^ be the
estimates of the unknown parameters under H0 and H1, respectively. The log-likelihood ratio can be
calculated as(2.28)LR1=−2[logL(Ω˜∣y)−logL(Ω^∣y)].When testing the hypothesis, the polygene and the
residual variances are nuisance parameters which are constrained to be nonnegative. The three tested genetic variance components under the null are lied on the
boundaries of their alternative parameter spaces. Following Self and Liang
[28], when the null is true, LR1 may
asymptotically follows a mixture of χ2 distribution on 0,…, 3 degrees of
freedom (df) with the mixture proportion for the χk2 components
being (3k)2−3. The theoretical distribution can be used to assess
significance in linkage scan. However, since there are many point tests across
the genome, the pointwise significance value may not guarantee an appropriate
genome-wide error rate. Another approach to assess significance is to use
nonparametric permutation tests in which the critical threshold value can be
empirically calculated on the basis of repeatedly shuffling the relationships
between marker genotypes and phenotypes [29]. In simulation studies, we also
simulate the null distribution and compare it with the theoretical
distribution.For those detected QTL, the next step is to assess
their imprinting property. An identified QTL can be imprinted, completely
imprinted, partially imprinted, or not imprinted at all. These can be tested
through the following sequential tests. The first imprinting test is to assess
whether a QTL shows imprinting effect, which can be done by formulating the
following hypotheses:(2.29)H0:σm2=σf2=σ2,H1:σm2≠σf2.Rejection of H0 provides
evidence of genomic imprinting and the QTL is called iQTL. Again, likelihood
ratio test can be applied in which the log-likelihood ratio test statistics
asymptotically follows an χ2 with one df
[7]. We denote the log-likelihood ratio test statistic as LRimp. If the null is rejected, one would be interested to
test if the detected iQTL is completely maternally or paternally imprinted. The
corresponding hypotheses can be formulated as(2.30)H0:σm2=0,H1:σm2≠0for testing completely maternal
imprinting and(2.31)H0:σf2=0,H1:σf2≠0for testing completely paternal
imprinting. The likelihood ratio test statistics for the above two tests
asymptotically follow a 50 : 50 mixture of χ02 and χ12
distributions
[28]. Rejection of complete imprinting indicates partial imprinting.


## 2.7. Multiple QTL Model

In reality, more than one QTL may contribute to the phenotypic variation located in one
chromosome region or across the whole genome. The polygenic effect in model (2.4)
absorbs the effects of multiple QTL located on other chromosomes. However, when
there are multiple QTL located on the same linkage group as the tested QTL, if
their effects are not properly adjusted, the estimation could be biased due to
interference caused by these QTL outside of the
testing interval [3, 30–32]. A multiple QTL model that can test the putative
QTL effect while adjusting the effects of interference QTL deserves more
attention.Zeng [32] previously showed that IBD variables share
the same property as the indicator variables in which the shared proportion of
alleles IBD for a QTL conditional on the IBD of one flanking marker is
independent of that of a QTL on the other side of that flanking marker. Thus, conditional on one flanking marker, the interference of QTL located on the other
side of the marker can be eliminated. By conditional on the IBD of the flanking
markers, the IBD sharing of a QTL is uncorrelated with that outside this
interval. Xu and Atchley [25] showed that one marker is enough to block the
interference caused by other QTL located on the same linkage group. The authors
derived the next-to-flanking markers structure to block additional QTL effects
from both sides of testing region in one chromosome. We derive a multiple QTL
model adopting a similar idea as Xu and Atchley [25]. Assume there are total S QTL located on
a linkage group. Considering parent-specific allelic effects, the multiple QTL
model can be expressed in general as(2.32)yik=μk+∑s=1Saikms+∑s=1Saikfs+Gik+eik,k=1,…,K;i=1,…,nk.In an interval-based linkage scan, only one putative
QTL is considered at each testing position conditioning on the effects of all
other QTL. Assuming there are totalL and R QTL located on
the left and right sides of the putative QTL on a
linkage group, model (2.32) can be modified as(2.33)yik=μk+∑l=1Laikl+(aikm+aikf)+∑r=1Raikr+Gik+eik,k=1,…,K;i=1,…,nk,where aikl and aikr are the lth and rth QTL random
effects on the left and right sides of the putative
QTL, respectively. When testing the putative QTL effect, we are only interested
in blocking the total effects of QTL outside of the tested interval. Therefore, in the modified model, the effects of QTL outside of the tested interval are
not partitioned. This, however, does not affect the inference of the tested
QTL.As shown by Zeng [32] and Jansen [31, 33], one marker
is enough to block the correlation between a locus on its left and a locus on
its right. Therefore, only two additional markers flanking the current interval
are needed to block interference caused by outside QTL [25]. Let ℳl and ℳr denote two
flanking markers for the tested interval, and ℒ and ℛ denote the two
markers next to ℳl and ℳl+1 with the marker
order ℒ - ℳl - ℳl+1 - ℛ. With the modified model given in (2.33), the
covariance of phenotypes between individuals i and j in the kth backcross
family can be expressed as(2.34)Cov(yik,yjk)=∑l=1LCov(aikl,ajkl)+Cov(aikm,ajkm)+Cov(aikf,ajkf)+Cov(aikm,ajkf)+∑r=1RCov(aikr,ajkr)+ϕijσg2+Iijσe2=∑l=1Lπl∣kσl2+πimjmσm2+πim/jfσmf2+πifjfσf2+∑r=1Rπr∣kσr2+ϕijσg2+Iijσe2,where πl∣k and πr∣k are the IBD
values for QTL located on the left and right sides
of the putative QTL in the kth backcross
family, and can be calculated following (2.6) and (2.7) if their genotype
information is known. Unfortunately, the number and exact locations of QTL
outside the testing interval are unknown. Hence, πl∣k and πr∣k are not
observable. Xu and Atchley [25] showed that when πl∣k and πr∣k are unknown, they can be estimated by some composite terms K(θlℒ,πℒ∣k) and K(θlℛ,πℛ∣k), where K(θlℒ,πℒ∣k) is a function
of the recombination fraction between the lth QTL and the
left marker ℒ as well as a
function of πℒ∣k, the IBD value for a pair of individuals at the left
marker ℒ. K(θlℛ,πℛ∣k) can be
similarly defined. Following Xu and Atchley [25], K(θlℒ,πℒ∣k) can be
expressed as a function of πℒ∣k multiplied by a
function of recombination frequency between the lth QTL and the
marker ℒ, f(θlℒ), that is, K(θlℒ,πℒ∣k)=πℒ∣kf(θlℒ). Similarly, K(θlℛ,πℛ∣k)=πℛ∣kf(θrℛ). When doing an interval scan, the covariance function
given in (2.34) between individuals i and j can be
reexpressed as(2.35)Cov(yik,yjk∣πℒ|k,π^imjm,π^im/jf,π^ifjf,πℛ|k)=∑l=1LK(θlℒ,πℒ∣k)σl2+π^imjmσm2+π^im/jf∣kσmf2+π^ifjfσf2+∑r=1RK(θlℛ,πℛ∣k)σr2+ϕijσg2+Iijσe2=πℒ∣k∑l=1Lf(θlℒ)σl2+π^imjmσm2+π^im/jfσmf2+π^ifjfσf2+πℛ∣k∑r=1Rf(θrℛ)σr2+ϕijσg2+Iijσe2=πℒ∣kσL2+π^imjmσm2+π^im/jf∣kσmf2+π^ififσf2+πℛ∣kσR2+ϕijσg2+Iijσe2.Instead of estimating individual
variance components σl2 and σr2, now we estimate the composite
terms
∑l=1Lf(θlℒ)σl2=σL2 and ∑r=1Rf(θrℛ)σr2=σR2. By conditioning the IBD sharing information for the
left and right markers ℒ and ℛ, the effects of those interference QTL are blocked. σL2 and σR2 absorb the
random effects of all QTL that are outside of the testing interval but are on
the same linkage group as the putative QTL. Estimation of the variance
components terms follows the same procedure as the single QTL analysis with
slight modification to consider multiple variance components.


## 3. Results

### 3.1. Simulation Design

To investigate
the performance of the proposed models and estimation methods, we conduct
intensive computer simulations. We start with the single-QTL simulation
followed by the multiple QTL analysis. Six evenly spaced markers (ℳ1−ℳ6) are simulated. The total length for the simulated
linkage group is 100 cm. We assume that all the backcross families share the
same linkage map constructed using Haldane map function. For simplicity, we
assume the sample size for all backcross families is the same (i.e., nk=n). The position
of the simulated QTL is assumed to be located 48 cm away from the first marker (ℳ1). The effect of the putative QTL is simulated by
assuming different imprinting mechanisms, that is, no imprinting, completely
imprinting, and partial imprinting. Once QTL genotypes are simulated, phenotypes can be simulated by randomly drawing multivariate normal
distribution with the covariance structure given in (2.12) with different
parameter combinations.To evaluate the effect of family and offspring size
combination on testing power and parameter estimation, we simulate data
assuming different sample size combinations. We fix the total sample size as
400 and vary the family and offspring size with different combinations, that
is,4×100, 8×50, 20×20, and 100×4. The first
number for each combination indicates the family size. For example, in the
combination 4×100, 4 families
each containing 100 offsprings are simulated. For each sib-pair, the IBD value
at a putative position at every 2 cm along the linkage group is calculated as
described in Section 2.7. For each simulation scenario, 100 simulation
replications are recorded, and the ML and the REML methods are used to estimate
the unknown parameters.


### 3.2. Simulation Results

#### 3.2.1. Single QTL Analysis

The single QTL
model assumes one QTL is located at the third interval in the simulated linkage
group, 48 cm away from the first marker. Results using both ML and REML
estimation methods are summarized in Table2. nF denotes the
number of families and nk denotes the
number of offsprings for each family. Without loss of generality, we assume
equal offspring size for all families in each simulation scenario. The
simulated parameter values are listed under each parameter. The root mean
square errors (RMSEs) are recorded for each parameter estimate to assess the
estimation precision. Overall, the fixed effects (three means) and most
variance components can be better estimated with large number of families. For
example, the RMSE of parameter μ1 is reduced from
1.869 (2.45) to 0.321 (0.305) when the number of families increases from 4 to
100 with the ML (REML) estimation method. The only exception is the two
variance components terms (σm2 and σf2) which are
better estimated with the 20×20 combination
design. Through the combination of different line crosses, the parameter
inference space is expanded, and as a result, better estimations are achieved
as expected. However, the QTL position is better estimated with the 8×50 and 20×20 designs than
the other two among the four simulation scenarios. The 100×4 design gives
the worst QTL position estimation with the largest RMSEs for both estimation
methods. Therefore, a balance of family and offspring size is needed. A
moderate family size with moderate offspring size would be necessary in order
to achieve reasonable parameter estimation for both QTL effects and position.Table 2
The power, MLEs, and REMLs of the QTL position, and effect parameters estimated based on 100
simulation replicates for a QTL with no imprinting effect under different
sampling designs. The square roots of the mean squared errors of the parameter
estimates are given in parentheses.
EstimationPositionμ1μ2μ3σm2σf2σmf2σg2σe2Type
InF×nkmethod48 cm10861.51.50.50.52Power1Power2error4×100ML48.789.858.085.800.910.730.780.172.120.910.990.42(17.87)(1.869)(1.246)(1.768)(1.296)(1.122)(1.242)(0.572)(0.229)REML47.29.838.136.231.511.421.521.062.020.920.960.17(11.95)(2.45)(1.39)(2.42)(1.65)(1.33)(2.88)(2.35)(0.57)8×50ML50.049.687.975.871.091.280.980.372.060.811.000.26(14.73)(1.373)(0.736)(1.157)(0.943)(0.990)(1.302)(0.701)(0.233)REML46.99.848.125.911.611.641.020.702.030.981.000.08(7.93)(1.16)(0.83)(1.10)(0.99)(1.26)(1.44)(1.10)(0.37)20×20ML48.969.997.956.051.341.250.600.322.010.980.980.12(10.70)(0.706)(0.495)(0.633)(0.768)(0.771)(0.810)(0.478)(0.182)REML49.7410.047.945.891.551.480.690.581.980.950.980.06(13.96)(0.597)(0.458)(0.675)(0.964)(0.786)(0.838)(0.608)(0.257)100×4ML51.3410.027.956.021.501.450.580.412.010.550.860.11(18.00)(0.321)(0.277)(0.345)(0.875)(0.878)(0.553)(0.444)(0.242)REML48.49.958.016.011.511.570.620.521.990.670.810.07(20.84)(0.305)(0.220)(0.319)(0.902)(0.892)(0.654)(0.485)(0.222)The locations of the  QTL are described by
the map distances (in cm) from the first marker of the linkage group (100 cm long). The true QTL is located at 48 cm; Powr1 is calculated using the empirical distribution
through simulation. Powr2 is calculated using the theoretical
distribution assuming mixture chi-square distribution. Type I error refers to
the imprinting type I error.Table2 also lists the results of power analysis under
different scenarios with two different estimation methods. Power1 denotes the
empirical power calculated from the simulated null distribution corresponding
to hypothesis test (2.27). We simulate the null distribution by simulating data
assuming no QTL effect (i.e., σm2 = σf2 = σmf2 = 0). The LR
test statistics is calculated for each simulation run, and the 95% cutoff is
reported as the threshold value. Power2 refers to the
theoretical power which is calculated assuming the mixture chi-square distribution, that is, (1/8)χ02+(3/8)χ12+(3/8)χ22+(1/8)χ32 [28]. Results
show that the threshold calculated from the theoretical distribution is smaller
than the one calculated from the simulation. Thus, the testing power based on
the theoretical cutoff is greater than the empirical power. The testing powers
under different sampling designs are very comparable except for the 100 × 4 design in
which the power is dramatically reduced compared to other designs. No
remarkable difference in power for both estimation methods is observed. Figure 2 shows the log-likelihood ratio test statistic calculated under the four
sampling designs across the simulated linkage group by using both ML and REML
estimation methods. The plotted LR curve is from averaged LR values out of 100
replications. It is clear that large offspring size always gives large test
statistics. As the family size increases from 4 to 100 and so decreased
offspring size, we observe a huge LR value decrease. Clearly, the 100 × 4 design is
less powerful than the others. The last column listed in Table 2 shows type I
error for testing genomic imprinting, that is, H0 : σm2=σf2. The simulated data assume no imprinting (σm2=σf2=1. 5). The imprinting test is only conducted at the
position where the overall QTL test shows significance. The imprinting test
statistic LRimp is compared
with a chi-square distribution with 1 df. Overall, the REML estimation method
results in smaller type I error rate than the ML method does. As the number of
families increase, type I error decreases. The 4 ×100 design yields
the largest type I error.The LR profile plot. The left and
right figures correspond to the LR profiles generated using the ML and REML
methods, respectively. The arrow indicates the true
QTL position.
(a)(b)In comparison of the ML and REML methods, the REML
method gives smaller estimation biases but larger RMSEs than the ML method
does. This reflects the large variability of the REML estimation. In terms of
computation speed, the ML method is faster than the REML method. For example, in a single-simulation run with the one-QTL model, the ML method takes about 9
minutes to scan the linkage group compared to 26 minutes with the REML method. The difference is more remarkable with the multiple QTL model (e.g., 10 minutes
for ML versus 43 minutes for REML). Even though the QTL position estimation is
better estimated by using the REML method when family size is small, as family
size increases, the REML method performs worse than the ML method (Table2). In
checking the LR profile plot in Figure 2 and the power analysis in Table 2, we
do not observe significant gain in power by using the REML method. The two
methods do no dominate each other and are very comparable in power analysis. With large sample size and limited computing resources, one might want to try
the ML method first. However, the REML method is suggested when testing
imprinting since it has small type I error.In a short summary of the results listed in Table2, the 8×50 and 20×20 designs give
better QTL position estimation and testing power. In terms of the type I error
for imprinting test, the 20×20 and 100×4 designs
provide reasonable type I error. Thus, a practical guidance is to choose the 20×20 design, and
one should always avoid designs with extremely large or extremely small family
size.To evaluate the proposed model under different
imprinting mechanisms, we simulated data assuming different
degrees of imprinting. Since the results in Table2
indicate that a 20×20 design
provides relatively reasonable parameter estimation, good power, and small type
I error rate for imprinting test, the evaluation of imprinting analysis is thus
focused on this design. The results for 100 simulation
replications are summarized in Table 3. Three
imprinting models are assumed complete maternal imprinting (σm2=0 and σf2=3), complete
paternal imprinting (σm2=3 and σf2=0), and partial
maternal imprinting (σm2=1 and σf2=2). Both ML
and REML estimators are reported. Overall, the two estimation methods produce
very comparable results with less-biased estimations by the REML method as we
expected. All the parameters can be properly estimated with reasonable
precision. Large imprinting power is observed when the variance difference
between the two parent-specific variance components is large. When the
difference between the two parent-specific variance components is reduced, the
power to detect imprinting is largely reduced. For example, when data are
simulated assuming complete paternal imprinting, the power is 0.91(0.86) by
using the ML(REML) estimation method. With partially imprinted data, the
imprinting power reduces to 0.24(0.09) by using the ML(REML) method, even
though it can be increased by increasing the offspring sample size (data not
shown).Table 3
The power, MLEs, and REMLs of the QTL position, and effect parameters estimated based on 100
simulation replicates for a QTL showing different imprinting effects under the20×20 sampling
design. The square roots of the mean squared errors of the parameters are given
in parentheses.
EstimationPositionμ1μ2μ3σm2σf2σmf2σg2σe2method48 cm10860.50.52Powr1Powr2ipower30ML48.9810.068.025.942.940.100.650.252.040.980.990.91(7.36)(0.716)(0.388)(0.722)(1.197)(0.217)(0.694)(0.426)(0.207)REML48.549.967.955.902.950.090.640.571.970.970.990.86(9.62)(0.705)(0.428)(0.747)(1.501)(0.180)(0.598)(0.609)(0.226)03ML46.6210.058.026.080.092.760.690.302.050.971.000.89(8.11)(0.680)(0.442)(0.584)(0.201)(1.052)(0.722)(0.430)(0.193)REML48.710.047.985.910.092.950.610.571.980.960.980.88(7.40)(0.560)(0.506)(0.615)(0.195)(1.279)(0.663)(0.544)(0.225)12ML47.810.098.006.110.842.040.660.322.020.971.000.24(7.83)(0.634)(0.468)(0.666)(0.615)(0.994)(0.811)(0.475)(0.203)REML49.6610.027.965.901.081.730.680.612.000.970.990.09(11.06)(0.606)(0.548)(0.668)(0.679)(0.958)(0.707)(0.616)(0.216)Por1 and Por2 correspond to the overall QTL effect test (2.27)
calculated using the empirical and theoretical cutoffs, respectively; ipower
refers to the imprinting test power corresponding to test (2.29). See Table 2 for
explanations of other parameters.In reality, whether a QTL is imprinted or not is an
unknown prior. When a QTL has Mendelian effect and is not imprinted, is there
any power loss by analyzing with the proposed imprinting model? Or when a QTL
is actually imprinted, is there any power loss by analyzing with regular
variance components approach? To answer these two questions, we simulated data
under different scenarios and analyzed with both Mendelian and imprinting
models. The first and second columns in Table4
refer to the simulation and analysis models, respectively. M refers to the
Mendelian model without variance components partition and I refers to the
imprinting model with allelic-specific partition of the variance components. For comparison purpose, heritabilities are recorded instead of original
variance components estimates. The polygene and residual variances are fixed as
0.5 (hg2=0.083) and 2, respectively for all the simulation
scenarios. We first simulated data with one additive genetic effect without
partitioning variance into allelic-specific components. This is equivalent to
simulate data assuming the Mendelian model. A single additive variance
component of 3.5 is assumed which corresponds to a heritability of ha2=0.583. The second scenario is to simulate data with three
allelic-specific variance components. Simulation models I1 and I2 correspond to a
complete maternal imprinting model (i.e., hm2=0 and hf2=0.5) and a partial
maternal imprinting model (i.e., hm2=0.083 and hf2=0.417), respectively. The variance component σmf2 is assumed to
be 0.5 (hmf2=0.083) for I1 and I2. In all the simulations, we use the 20×20 design to
make the comparison. Similar results are expected under the other sampling
designs. Since the true variance components values for the imprinting model
are unknown when data are simulated assuming
Mendelian effect and vice versa, only standard deviations for these parameter
estimates are recorded (listed as italic font in the parentheses).Table 4
The power, MLEs of the QTL position, and effect parameters estimated based on 100
simulation replicates for data simulated with Mendelian and imprinting models
based on the20×20 sampling
design.
SimulationAnalysisPositionμ1μ2μ3ha2hg2σe2modelmodel48 cm1086hm2hf2hmf20.0832Power1Power2MM48.0810.1237.8906.0230.5440.0931.9821.001.00(2.981)(0.887)(0.730)(1.046)(0.121)(0.137)(0.250)I48.2609.9228.0816.0740.2630.2720.2390.0491.9831.001.00(3.015)(1.010)(0.720)(0.949)(0.097)(0.110)(0.124)(0.080)(0.249)I1M48.94010.0797.9136.0230.3230.1112.0220.860.96(11.660)(0.588)(0.511)(0.712)(0.287)(0.134)(0.205)I47.8610.0777.9006.0180.0080.4730.0890.0912.0150.950.98(7.524)(0.579)(0.575)(0.665)(0.022)(0.151)(0.099)(0.053)(0.86)I2M48.64010.0707.9176.0280.3310.1032.0230.910.98(8.389)(0.613)(0.517)(0.706)(0.278)(0.129)(0.195)I49.20010.0597.9266.0080.0790.3960.0890.0892.0230.930.98(9.818)(0.617)(0.518)(0.673)(0.075)(0.156)(0.097)(0.050)(0.184)M and I refer to Mendelian and imprinting
models, respectively. Simulated parameters for model M:(ha2=0.583); I1:(hm2,hf2,hmf2)=(0, 0.5, 0.083); I2:(hm2,hf2,hmf2)=( 0.083, 0.417, 0.083). Power1 and Power2 correspond to the power calculated using the
empirical cutoff and the theoretical threshold, respectively. The numbers given
in the parenthesis with normal and italic fonts correspond to the RMSEs and
standard errors of the parameter estimates, respectively. See Table 2 for other
explanations.The simulation results are summarized in Table4
analyzed with the ML method. When the simulated model is Mendelian, QTL
position is better estimated with the Mendelian model than with the imprinting
model. No remarkable difference in power is observed for both models. The
estimated parent-specific variances due to maternal and paternal alleles are
almost identical and no imprinting is detected. When data are simulated
assuming imprinting (model I1 and I2), large power
is observed when analyzed with the imprinting model. For example, the power is
86% when analyze the I1 imprinting data
by the Mendelian model. The power is increased to 95% when data are analyzed by
the imprinting model. When imprinting data are analyzed with the Mendelian
model, the major QTL variance is underestimated, and the polygene variance is
slightly overestimated. No remarkable differences are observed for the
estimation of the three fixed mean effects and the residual variance under all
simulation cases. In any case, the imprinting model performs better or no worse
than the Mendelian model in terms of power. In checking type I error rate based
on the theoretical threshold, we find the imprinting model has slightly higher
type I error rate compared with the Mendelian model. In real data analysis, it
is more important to control the false negatives than the false positives. Thus, it is safe to apply the imprinting model for data with any inheritance
pattern in this regard.


#### 3.2.2. Multiple QTL Analysis

To see the
relative merit of multiple QTL analysis against single QTL analysis when
multiple QTL are located on the same linkage group, two QTL are simulated with
QTL 1 (denoted asQ1) located at
the second interval, 28cm away from the
first marker (ℳ1), and QTL 2 (denoted as Q2) located at
the fourth interval, 68 cm away from the first marker. Two simulation scenarios
are considered. The first scenario considers two nonimprinted QTL with equal
genetic effects. The second scenario assumes
Q1 is imprinted
and Q2 is not
imprinted. Simulated parameters for the two QTL are listed in Table 5. Data are
simulated assuming the 20×20 design. Parameters are estimated by the ML and REML approaches with 100 replicates.Table 5
The MLEs and
REMLs of the QTL position and effect parameters estimated based on 100
simulation replicates for data simulated with two QTL under the20×20 design. The
square roots of the mean squared errors are given in parentheses.
Q1Q2EstimationPositionσm2σf2σmf2Positionσm2σf2σmf2method28 cm0.2568 cm0.250.750.750.750.75ML31.060.9600.8190.38164.5600.8060.8910.414(12.825)(0.760)(0.693)(0.533)(12.126)(0.552)(0.679)(0.646)REML32.320.9490.9240.44064.600.9431.0020.461(13.951)(0.739)(0.686)(0.678)(12.169)(0.713)(0.711)(0.601)1.500.750.75ML29.521.5680.2100.60163.861.2100.5320.533(12.153)(0.872)(0.382)(0.836)(11.913)(0.932)(0.572)(0.847)REML31.221.6680.2820.43364.161.3550.6650.529(12.557)(0.881)(0.457)(0.598)(13.030)(1.128)(0.554)(0.650)Q1 and Q2 refer to two QTL located at 28 cm and 68 cm.Figure3 shows the LR profile plots for the single and
multiple QTL analyses. The single QTL model
indicates three major peaks. The highest peak for the single QTL analysis is
located at the wrong QTL interval where no QTL is assumed. The so-called
“ghost image” of QTL can be removed and the positions of the two QTL can
be precisely mapped on the chromosome by the multiple QTL model. Two clear
peaks indicating the correct QTL positions (arrow signs) are observed by the
multiple QTL analysis. However, we observe a remarkable reduction in LR values
by multiple QTL analysis compared to those by the single QTL analysis. Since
the threshold for multiple QTL analysis is unknown, we cannot make the
conclusion that multiple QTL analysis is less powerful than the single QTL
analysis. It is possible that we may gain accuracy in QTL position estimation
at the cost of power loss. Similar phenomenon and issues were also observed and
discussed in the literature [3, 25].The LR profile plot for single QTL
and multiple QTL analyses. The true QTL positions
are simulated at 28cm and 68 cm (see the arrow
sign). The dotted curve and the solid curve represent the LR profiles by single
QTL and multiple QTL analyses, respectively. The
left and right figures correspond to the LR profiles generated using the ML and
REML methods, respectively.
(a)(b)The results of the multiple QTL analysis are
summarized in Table5. The fixed mean effects, the polygene, and residual
variance components can be reasonably estimated with small RMSEs, similar
results shown in Table 2 for the 20×20 design and
hence are not reported here. Only the genetic factors for the two simulated
QTLs are reported. It can be seen that both ML and REML methods provide
reasonable parameter estimates and are very comparable. Under the first
simulation scenario in which both QTL are not imprinted, the genetic effects
are all slightly overestimated by both methods. This might be due to the
interference of the two QTL in the same linkage group. The multiple QTL model
may not completely block the effects of QTL outside of the tested interval. For
the second simulation scenario, an interesting pattern is observed. When one
QTL is imprinted (Q1), the maternal and paternal variance components for
the second one (Q2) tend to be estimated with bias in the direction as
the first imprinted QTL, that is, σm2 tends to be
overestimated and σf2 tends to be
underestimated. As we gain accuracy in QTL position estimation, we lose
precision for the parameter estimation. These effects are expected as described
in Zeng [3] and Xu and Atchley [25]. More investigations are needed in multiple
QTL analysis in order to maintain a good balance of QTL position and parameter
inference.


## 3.1. Simulation Design

To investigate
the performance of the proposed models and estimation methods, we conduct
intensive computer simulations. We start with the single-QTL simulation
followed by the multiple QTL analysis. Six evenly spaced markers (ℳ1−ℳ6) are simulated. The total length for the simulated
linkage group is 100 cm. We assume that all the backcross families share the
same linkage map constructed using Haldane map function. For simplicity, we
assume the sample size for all backcross families is the same (i.e., nk=n). The position
of the simulated QTL is assumed to be located 48 cm away from the first marker (ℳ1). The effect of the putative QTL is simulated by
assuming different imprinting mechanisms, that is, no imprinting, completely
imprinting, and partial imprinting. Once QTL genotypes are simulated, phenotypes can be simulated by randomly drawing multivariate normal
distribution with the covariance structure given in (2.12) with different
parameter combinations.To evaluate the effect of family and offspring size
combination on testing power and parameter estimation, we simulate data
assuming different sample size combinations. We fix the total sample size as
400 and vary the family and offspring size with different combinations, that
is,4×100, 8×50, 20×20, and 100×4. The first
number for each combination indicates the family size. For example, in the
combination 4×100, 4 families
each containing 100 offsprings are simulated. For each sib-pair, the IBD value
at a putative position at every 2 cm along the linkage group is calculated as
described in Section 2.7. For each simulation scenario, 100 simulation
replications are recorded, and the ML and the REML methods are used to estimate
the unknown parameters.


## 3.2. Simulation Results

### 3.2.1. Single QTL Analysis

The single QTL
model assumes one QTL is located at the third interval in the simulated linkage
group, 48 cm away from the first marker. Results using both ML and REML
estimation methods are summarized in Table2. nF denotes the
number of families and nk denotes the
number of offsprings for each family. Without loss of generality, we assume
equal offspring size for all families in each simulation scenario. The
simulated parameter values are listed under each parameter. The root mean
square errors (RMSEs) are recorded for each parameter estimate to assess the
estimation precision. Overall, the fixed effects (three means) and most
variance components can be better estimated with large number of families. For
example, the RMSE of parameter μ1 is reduced from
1.869 (2.45) to 0.321 (0.305) when the number of families increases from 4 to
100 with the ML (REML) estimation method. The only exception is the two
variance components terms (σm2 and σf2) which are
better estimated with the 20×20 combination
design. Through the combination of different line crosses, the parameter
inference space is expanded, and as a result, better estimations are achieved
as expected. However, the QTL position is better estimated with the 8×50 and 20×20 designs than
the other two among the four simulation scenarios. The 100×4 design gives
the worst QTL position estimation with the largest RMSEs for both estimation
methods. Therefore, a balance of family and offspring size is needed. A
moderate family size with moderate offspring size would be necessary in order
to achieve reasonable parameter estimation for both QTL effects and position.Table 2
The power, MLEs, and REMLs of the QTL position, and effect parameters estimated based on 100
simulation replicates for a QTL with no imprinting effect under different
sampling designs. The square roots of the mean squared errors of the parameter
estimates are given in parentheses.
EstimationPositionμ1μ2μ3σm2σf2σmf2σg2σe2Type
InF×nkmethod48 cm10861.51.50.50.52Power1Power2error4×100ML48.789.858.085.800.910.730.780.172.120.910.990.42(17.87)(1.869)(1.246)(1.768)(1.296)(1.122)(1.242)(0.572)(0.229)REML47.29.838.136.231.511.421.521.062.020.920.960.17(11.95)(2.45)(1.39)(2.42)(1.65)(1.33)(2.88)(2.35)(0.57)8×50ML50.049.687.975.871.091.280.980.372.060.811.000.26(14.73)(1.373)(0.736)(1.157)(0.943)(0.990)(1.302)(0.701)(0.233)REML46.99.848.125.911.611.641.020.702.030.981.000.08(7.93)(1.16)(0.83)(1.10)(0.99)(1.26)(1.44)(1.10)(0.37)20×20ML48.969.997.956.051.341.250.600.322.010.980.980.12(10.70)(0.706)(0.495)(0.633)(0.768)(0.771)(0.810)(0.478)(0.182)REML49.7410.047.945.891.551.480.690.581.980.950.980.06(13.96)(0.597)(0.458)(0.675)(0.964)(0.786)(0.838)(0.608)(0.257)100×4ML51.3410.027.956.021.501.450.580.412.010.550.860.11(18.00)(0.321)(0.277)(0.345)(0.875)(0.878)(0.553)(0.444)(0.242)REML48.49.958.016.011.511.570.620.521.990.670.810.07(20.84)(0.305)(0.220)(0.319)(0.902)(0.892)(0.654)(0.485)(0.222)The locations of the  QTL are described by
the map distances (in cm) from the first marker of the linkage group (100 cm long). The true QTL is located at 48 cm; Powr1 is calculated using the empirical distribution
through simulation. Powr2 is calculated using the theoretical
distribution assuming mixture chi-square distribution. Type I error refers to
the imprinting type I error.Table2 also lists the results of power analysis under
different scenarios with two different estimation methods. Power1 denotes the
empirical power calculated from the simulated null distribution corresponding
to hypothesis test (2.27). We simulate the null distribution by simulating data
assuming no QTL effect (i.e., σm2 = σf2 = σmf2 = 0). The LR
test statistics is calculated for each simulation run, and the 95% cutoff is
reported as the threshold value. Power2 refers to the
theoretical power which is calculated assuming the mixture chi-square distribution, that is, (1/8)χ02+(3/8)χ12+(3/8)χ22+(1/8)χ32 [28]. Results
show that the threshold calculated from the theoretical distribution is smaller
than the one calculated from the simulation. Thus, the testing power based on
the theoretical cutoff is greater than the empirical power. The testing powers
under different sampling designs are very comparable except for the 100 × 4 design in
which the power is dramatically reduced compared to other designs. No
remarkable difference in power for both estimation methods is observed. Figure 2 shows the log-likelihood ratio test statistic calculated under the four
sampling designs across the simulated linkage group by using both ML and REML
estimation methods. The plotted LR curve is from averaged LR values out of 100
replications. It is clear that large offspring size always gives large test
statistics. As the family size increases from 4 to 100 and so decreased
offspring size, we observe a huge LR value decrease. Clearly, the 100 × 4 design is
less powerful than the others. The last column listed in Table 2 shows type I
error for testing genomic imprinting, that is, H0 : σm2=σf2. The simulated data assume no imprinting (σm2=σf2=1. 5). The imprinting test is only conducted at the
position where the overall QTL test shows significance. The imprinting test
statistic LRimp is compared
with a chi-square distribution with 1 df. Overall, the REML estimation method
results in smaller type I error rate than the ML method does. As the number of
families increase, type I error decreases. The 4 ×100 design yields
the largest type I error.The LR profile plot. The left and
right figures correspond to the LR profiles generated using the ML and REML
methods, respectively. The arrow indicates the true
QTL position.
(a)(b)In comparison of the ML and REML methods, the REML
method gives smaller estimation biases but larger RMSEs than the ML method
does. This reflects the large variability of the REML estimation. In terms of
computation speed, the ML method is faster than the REML method. For example, in a single-simulation run with the one-QTL model, the ML method takes about 9
minutes to scan the linkage group compared to 26 minutes with the REML method. The difference is more remarkable with the multiple QTL model (e.g., 10 minutes
for ML versus 43 minutes for REML). Even though the QTL position estimation is
better estimated by using the REML method when family size is small, as family
size increases, the REML method performs worse than the ML method (Table2). In
checking the LR profile plot in Figure 2 and the power analysis in Table 2, we
do not observe significant gain in power by using the REML method. The two
methods do no dominate each other and are very comparable in power analysis. With large sample size and limited computing resources, one might want to try
the ML method first. However, the REML method is suggested when testing
imprinting since it has small type I error.In a short summary of the results listed in Table2, the 8×50 and 20×20 designs give
better QTL position estimation and testing power. In terms of the type I error
for imprinting test, the 20×20 and 100×4 designs
provide reasonable type I error. Thus, a practical guidance is to choose the 20×20 design, and
one should always avoid designs with extremely large or extremely small family
size.To evaluate the proposed model under different
imprinting mechanisms, we simulated data assuming different
degrees of imprinting. Since the results in Table2
indicate that a 20×20 design
provides relatively reasonable parameter estimation, good power, and small type
I error rate for imprinting test, the evaluation of imprinting analysis is thus
focused on this design. The results for 100 simulation
replications are summarized in Table 3. Three
imprinting models are assumed complete maternal imprinting (σm2=0 and σf2=3), complete
paternal imprinting (σm2=3 and σf2=0), and partial
maternal imprinting (σm2=1 and σf2=2). Both ML
and REML estimators are reported. Overall, the two estimation methods produce
very comparable results with less-biased estimations by the REML method as we
expected. All the parameters can be properly estimated with reasonable
precision. Large imprinting power is observed when the variance difference
between the two parent-specific variance components is large. When the
difference between the two parent-specific variance components is reduced, the
power to detect imprinting is largely reduced. For example, when data are
simulated assuming complete paternal imprinting, the power is 0.91(0.86) by
using the ML(REML) estimation method. With partially imprinted data, the
imprinting power reduces to 0.24(0.09) by using the ML(REML) method, even
though it can be increased by increasing the offspring sample size (data not
shown).Table 3
The power, MLEs, and REMLs of the QTL position, and effect parameters estimated based on 100
simulation replicates for a QTL showing different imprinting effects under the20×20 sampling
design. The square roots of the mean squared errors of the parameters are given
in parentheses.
EstimationPositionμ1μ2μ3σm2σf2σmf2σg2σe2method48 cm10860.50.52Powr1Powr2ipower30ML48.9810.068.025.942.940.100.650.252.040.980.990.91(7.36)(0.716)(0.388)(0.722)(1.197)(0.217)(0.694)(0.426)(0.207)REML48.549.967.955.902.950.090.640.571.970.970.990.86(9.62)(0.705)(0.428)(0.747)(1.501)(0.180)(0.598)(0.609)(0.226)03ML46.6210.058.026.080.092.760.690.302.050.971.000.89(8.11)(0.680)(0.442)(0.584)(0.201)(1.052)(0.722)(0.430)(0.193)REML48.710.047.985.910.092.950.610.571.980.960.980.88(7.40)(0.560)(0.506)(0.615)(0.195)(1.279)(0.663)(0.544)(0.225)12ML47.810.098.006.110.842.040.660.322.020.971.000.24(7.83)(0.634)(0.468)(0.666)(0.615)(0.994)(0.811)(0.475)(0.203)REML49.6610.027.965.901.081.730.680.612.000.970.990.09(11.06)(0.606)(0.548)(0.668)(0.679)(0.958)(0.707)(0.616)(0.216)Por1 and Por2 correspond to the overall QTL effect test (2.27)
calculated using the empirical and theoretical cutoffs, respectively; ipower
refers to the imprinting test power corresponding to test (2.29). See Table 2 for
explanations of other parameters.In reality, whether a QTL is imprinted or not is an
unknown prior. When a QTL has Mendelian effect and is not imprinted, is there
any power loss by analyzing with the proposed imprinting model? Or when a QTL
is actually imprinted, is there any power loss by analyzing with regular
variance components approach? To answer these two questions, we simulated data
under different scenarios and analyzed with both Mendelian and imprinting
models. The first and second columns in Table4
refer to the simulation and analysis models, respectively. M refers to the
Mendelian model without variance components partition and I refers to the
imprinting model with allelic-specific partition of the variance components. For comparison purpose, heritabilities are recorded instead of original
variance components estimates. The polygene and residual variances are fixed as
0.5 (hg2=0.083) and 2, respectively for all the simulation
scenarios. We first simulated data with one additive genetic effect without
partitioning variance into allelic-specific components. This is equivalent to
simulate data assuming the Mendelian model. A single additive variance
component of 3.5 is assumed which corresponds to a heritability of ha2=0.583. The second scenario is to simulate data with three
allelic-specific variance components. Simulation models I1 and I2 correspond to a
complete maternal imprinting model (i.e., hm2=0 and hf2=0.5) and a partial
maternal imprinting model (i.e., hm2=0.083 and hf2=0.417), respectively. The variance component σmf2 is assumed to
be 0.5 (hmf2=0.083) for I1 and I2. In all the simulations, we use the 20×20 design to
make the comparison. Similar results are expected under the other sampling
designs. Since the true variance components values for the imprinting model
are unknown when data are simulated assuming
Mendelian effect and vice versa, only standard deviations for these parameter
estimates are recorded (listed as italic font in the parentheses).Table 4
The power, MLEs of the QTL position, and effect parameters estimated based on 100
simulation replicates for data simulated with Mendelian and imprinting models
based on the20×20 sampling
design.
SimulationAnalysisPositionμ1μ2μ3ha2hg2σe2modelmodel48 cm1086hm2hf2hmf20.0832Power1Power2MM48.0810.1237.8906.0230.5440.0931.9821.001.00(2.981)(0.887)(0.730)(1.046)(0.121)(0.137)(0.250)I48.2609.9228.0816.0740.2630.2720.2390.0491.9831.001.00(3.015)(1.010)(0.720)(0.949)(0.097)(0.110)(0.124)(0.080)(0.249)I1M48.94010.0797.9136.0230.3230.1112.0220.860.96(11.660)(0.588)(0.511)(0.712)(0.287)(0.134)(0.205)I47.8610.0777.9006.0180.0080.4730.0890.0912.0150.950.98(7.524)(0.579)(0.575)(0.665)(0.022)(0.151)(0.099)(0.053)(0.86)I2M48.64010.0707.9176.0280.3310.1032.0230.910.98(8.389)(0.613)(0.517)(0.706)(0.278)(0.129)(0.195)I49.20010.0597.9266.0080.0790.3960.0890.0892.0230.930.98(9.818)(0.617)(0.518)(0.673)(0.075)(0.156)(0.097)(0.050)(0.184)M and I refer to Mendelian and imprinting
models, respectively. Simulated parameters for model M:(ha2=0.583); I1:(hm2,hf2,hmf2)=(0, 0.5, 0.083); I2:(hm2,hf2,hmf2)=( 0.083, 0.417, 0.083). Power1 and Power2 correspond to the power calculated using the
empirical cutoff and the theoretical threshold, respectively. The numbers given
in the parenthesis with normal and italic fonts correspond to the RMSEs and
standard errors of the parameter estimates, respectively. See Table 2 for other
explanations.The simulation results are summarized in Table4
analyzed with the ML method. When the simulated model is Mendelian, QTL
position is better estimated with the Mendelian model than with the imprinting
model. No remarkable difference in power is observed for both models. The
estimated parent-specific variances due to maternal and paternal alleles are
almost identical and no imprinting is detected. When data are simulated
assuming imprinting (model I1 and I2), large power
is observed when analyzed with the imprinting model. For example, the power is
86% when analyze the I1 imprinting data
by the Mendelian model. The power is increased to 95% when data are analyzed by
the imprinting model. When imprinting data are analyzed with the Mendelian
model, the major QTL variance is underestimated, and the polygene variance is
slightly overestimated. No remarkable differences are observed for the
estimation of the three fixed mean effects and the residual variance under all
simulation cases. In any case, the imprinting model performs better or no worse
than the Mendelian model in terms of power. In checking type I error rate based
on the theoretical threshold, we find the imprinting model has slightly higher
type I error rate compared with the Mendelian model. In real data analysis, it
is more important to control the false negatives than the false positives. Thus, it is safe to apply the imprinting model for data with any inheritance
pattern in this regard.


### 3.2.2. Multiple QTL Analysis

To see the
relative merit of multiple QTL analysis against single QTL analysis when
multiple QTL are located on the same linkage group, two QTL are simulated with
QTL 1 (denoted asQ1) located at
the second interval, 28cm away from the
first marker (ℳ1), and QTL 2 (denoted as Q2) located at
the fourth interval, 68 cm away from the first marker. Two simulation scenarios
are considered. The first scenario considers two nonimprinted QTL with equal
genetic effects. The second scenario assumes
Q1 is imprinted
and Q2 is not
imprinted. Simulated parameters for the two QTL are listed in Table 5. Data are
simulated assuming the 20×20 design. Parameters are estimated by the ML and REML approaches with 100 replicates.Table 5
The MLEs and
REMLs of the QTL position and effect parameters estimated based on 100
simulation replicates for data simulated with two QTL under the20×20 design. The
square roots of the mean squared errors are given in parentheses.
Q1Q2EstimationPositionσm2σf2σmf2Positionσm2σf2σmf2method28 cm0.2568 cm0.250.750.750.750.75ML31.060.9600.8190.38164.5600.8060.8910.414(12.825)(0.760)(0.693)(0.533)(12.126)(0.552)(0.679)(0.646)REML32.320.9490.9240.44064.600.9431.0020.461(13.951)(0.739)(0.686)(0.678)(12.169)(0.713)(0.711)(0.601)1.500.750.75ML29.521.5680.2100.60163.861.2100.5320.533(12.153)(0.872)(0.382)(0.836)(11.913)(0.932)(0.572)(0.847)REML31.221.6680.2820.43364.161.3550.6650.529(12.557)(0.881)(0.457)(0.598)(13.030)(1.128)(0.554)(0.650)Q1 and Q2 refer to two QTL located at 28 cm and 68 cm.Figure3 shows the LR profile plots for the single and
multiple QTL analyses. The single QTL model
indicates three major peaks. The highest peak for the single QTL analysis is
located at the wrong QTL interval where no QTL is assumed. The so-called
“ghost image” of QTL can be removed and the positions of the two QTL can
be precisely mapped on the chromosome by the multiple QTL model. Two clear
peaks indicating the correct QTL positions (arrow signs) are observed by the
multiple QTL analysis. However, we observe a remarkable reduction in LR values
by multiple QTL analysis compared to those by the single QTL analysis. Since
the threshold for multiple QTL analysis is unknown, we cannot make the
conclusion that multiple QTL analysis is less powerful than the single QTL
analysis. It is possible that we may gain accuracy in QTL position estimation
at the cost of power loss. Similar phenomenon and issues were also observed and
discussed in the literature [3, 25].The LR profile plot for single QTL
and multiple QTL analyses. The true QTL positions
are simulated at 28cm and 68 cm (see the arrow
sign). The dotted curve and the solid curve represent the LR profiles by single
QTL and multiple QTL analyses, respectively. The
left and right figures correspond to the LR profiles generated using the ML and
REML methods, respectively.
(a)(b)The results of the multiple QTL analysis are
summarized in Table5. The fixed mean effects, the polygene, and residual
variance components can be reasonably estimated with small RMSEs, similar
results shown in Table 2 for the 20×20 design and
hence are not reported here. Only the genetic factors for the two simulated
QTLs are reported. It can be seen that both ML and REML methods provide
reasonable parameter estimates and are very comparable. Under the first
simulation scenario in which both QTL are not imprinted, the genetic effects
are all slightly overestimated by both methods. This might be due to the
interference of the two QTL in the same linkage group. The multiple QTL model
may not completely block the effects of QTL outside of the tested interval. For
the second simulation scenario, an interesting pattern is observed. When one
QTL is imprinted (Q1), the maternal and paternal variance components for
the second one (Q2) tend to be estimated with bias in the direction as
the first imprinted QTL, that is, σm2 tends to be
overestimated and σf2 tends to be
underestimated. As we gain accuracy in QTL position estimation, we lose
precision for the parameter estimation. These effects are expected as described
in Zeng [3] and Xu and Atchley [25]. More investigations are needed in multiple
QTL analysis in order to maintain a good balance of QTL position and parameter
inference.


## 3.2.1. Single QTL Analysis

The single QTL
model assumes one QTL is located at the third interval in the simulated linkage
group, 48 cm away from the first marker. Results using both ML and REML
estimation methods are summarized in Table2. nF denotes the
number of families and nk denotes the
number of offsprings for each family. Without loss of generality, we assume
equal offspring size for all families in each simulation scenario. The
simulated parameter values are listed under each parameter. The root mean
square errors (RMSEs) are recorded for each parameter estimate to assess the
estimation precision. Overall, the fixed effects (three means) and most
variance components can be better estimated with large number of families. For
example, the RMSE of parameter μ1 is reduced from
1.869 (2.45) to 0.321 (0.305) when the number of families increases from 4 to
100 with the ML (REML) estimation method. The only exception is the two
variance components terms (σm2 and σf2) which are
better estimated with the 20×20 combination
design. Through the combination of different line crosses, the parameter
inference space is expanded, and as a result, better estimations are achieved
as expected. However, the QTL position is better estimated with the 8×50 and 20×20 designs than
the other two among the four simulation scenarios. The 100×4 design gives
the worst QTL position estimation with the largest RMSEs for both estimation
methods. Therefore, a balance of family and offspring size is needed. A
moderate family size with moderate offspring size would be necessary in order
to achieve reasonable parameter estimation for both QTL effects and position.Table 2
The power, MLEs, and REMLs of the QTL position, and effect parameters estimated based on 100
simulation replicates for a QTL with no imprinting effect under different
sampling designs. The square roots of the mean squared errors of the parameter
estimates are given in parentheses.
EstimationPositionμ1μ2μ3σm2σf2σmf2σg2σe2Type
InF×nkmethod48 cm10861.51.50.50.52Power1Power2error4×100ML48.789.858.085.800.910.730.780.172.120.910.990.42(17.87)(1.869)(1.246)(1.768)(1.296)(1.122)(1.242)(0.572)(0.229)REML47.29.838.136.231.511.421.521.062.020.920.960.17(11.95)(2.45)(1.39)(2.42)(1.65)(1.33)(2.88)(2.35)(0.57)8×50ML50.049.687.975.871.091.280.980.372.060.811.000.26(14.73)(1.373)(0.736)(1.157)(0.943)(0.990)(1.302)(0.701)(0.233)REML46.99.848.125.911.611.641.020.702.030.981.000.08(7.93)(1.16)(0.83)(1.10)(0.99)(1.26)(1.44)(1.10)(0.37)20×20ML48.969.997.956.051.341.250.600.322.010.980.980.12(10.70)(0.706)(0.495)(0.633)(0.768)(0.771)(0.810)(0.478)(0.182)REML49.7410.047.945.891.551.480.690.581.980.950.980.06(13.96)(0.597)(0.458)(0.675)(0.964)(0.786)(0.838)(0.608)(0.257)100×4ML51.3410.027.956.021.501.450.580.412.010.550.860.11(18.00)(0.321)(0.277)(0.345)(0.875)(0.878)(0.553)(0.444)(0.242)REML48.49.958.016.011.511.570.620.521.990.670.810.07(20.84)(0.305)(0.220)(0.319)(0.902)(0.892)(0.654)(0.485)(0.222)The locations of the  QTL are described by
the map distances (in cm) from the first marker of the linkage group (100 cm long). The true QTL is located at 48 cm; Powr1 is calculated using the empirical distribution
through simulation. Powr2 is calculated using the theoretical
distribution assuming mixture chi-square distribution. Type I error refers to
the imprinting type I error.Table2 also lists the results of power analysis under
different scenarios with two different estimation methods. Power1 denotes the
empirical power calculated from the simulated null distribution corresponding
to hypothesis test (2.27). We simulate the null distribution by simulating data
assuming no QTL effect (i.e., σm2 = σf2 = σmf2 = 0). The LR
test statistics is calculated for each simulation run, and the 95% cutoff is
reported as the threshold value. Power2 refers to the
theoretical power which is calculated assuming the mixture chi-square distribution, that is, (1/8)χ02+(3/8)χ12+(3/8)χ22+(1/8)χ32 [28]. Results
show that the threshold calculated from the theoretical distribution is smaller
than the one calculated from the simulation. Thus, the testing power based on
the theoretical cutoff is greater than the empirical power. The testing powers
under different sampling designs are very comparable except for the 100 × 4 design in
which the power is dramatically reduced compared to other designs. No
remarkable difference in power for both estimation methods is observed. Figure 2 shows the log-likelihood ratio test statistic calculated under the four
sampling designs across the simulated linkage group by using both ML and REML
estimation methods. The plotted LR curve is from averaged LR values out of 100
replications. It is clear that large offspring size always gives large test
statistics. As the family size increases from 4 to 100 and so decreased
offspring size, we observe a huge LR value decrease. Clearly, the 100 × 4 design is
less powerful than the others. The last column listed in Table 2 shows type I
error for testing genomic imprinting, that is, H0 : σm2=σf2. The simulated data assume no imprinting (σm2=σf2=1. 5). The imprinting test is only conducted at the
position where the overall QTL test shows significance. The imprinting test
statistic LRimp is compared
with a chi-square distribution with 1 df. Overall, the REML estimation method
results in smaller type I error rate than the ML method does. As the number of
families increase, type I error decreases. The 4 ×100 design yields
the largest type I error.The LR profile plot. The left and
right figures correspond to the LR profiles generated using the ML and REML
methods, respectively. The arrow indicates the true
QTL position.
(a)(b)In comparison of the ML and REML methods, the REML
method gives smaller estimation biases but larger RMSEs than the ML method
does. This reflects the large variability of the REML estimation. In terms of
computation speed, the ML method is faster than the REML method. For example, in a single-simulation run with the one-QTL model, the ML method takes about 9
minutes to scan the linkage group compared to 26 minutes with the REML method. The difference is more remarkable with the multiple QTL model (e.g., 10 minutes
for ML versus 43 minutes for REML). Even though the QTL position estimation is
better estimated by using the REML method when family size is small, as family
size increases, the REML method performs worse than the ML method (Table2). In
checking the LR profile plot in Figure 2 and the power analysis in Table 2, we
do not observe significant gain in power by using the REML method. The two
methods do no dominate each other and are very comparable in power analysis. With large sample size and limited computing resources, one might want to try
the ML method first. However, the REML method is suggested when testing
imprinting since it has small type I error.In a short summary of the results listed in Table2, the 8×50 and 20×20 designs give
better QTL position estimation and testing power. In terms of the type I error
for imprinting test, the 20×20 and 100×4 designs
provide reasonable type I error. Thus, a practical guidance is to choose the 20×20 design, and
one should always avoid designs with extremely large or extremely small family
size.To evaluate the proposed model under different
imprinting mechanisms, we simulated data assuming different
degrees of imprinting. Since the results in Table2
indicate that a 20×20 design
provides relatively reasonable parameter estimation, good power, and small type
I error rate for imprinting test, the evaluation of imprinting analysis is thus
focused on this design. The results for 100 simulation
replications are summarized in Table 3. Three
imprinting models are assumed complete maternal imprinting (σm2=0 and σf2=3), complete
paternal imprinting (σm2=3 and σf2=0), and partial
maternal imprinting (σm2=1 and σf2=2). Both ML
and REML estimators are reported. Overall, the two estimation methods produce
very comparable results with less-biased estimations by the REML method as we
expected. All the parameters can be properly estimated with reasonable
precision. Large imprinting power is observed when the variance difference
between the two parent-specific variance components is large. When the
difference between the two parent-specific variance components is reduced, the
power to detect imprinting is largely reduced. For example, when data are
simulated assuming complete paternal imprinting, the power is 0.91(0.86) by
using the ML(REML) estimation method. With partially imprinted data, the
imprinting power reduces to 0.24(0.09) by using the ML(REML) method, even
though it can be increased by increasing the offspring sample size (data not
shown).Table 3
The power, MLEs, and REMLs of the QTL position, and effect parameters estimated based on 100
simulation replicates for a QTL showing different imprinting effects under the20×20 sampling
design. The square roots of the mean squared errors of the parameters are given
in parentheses.
EstimationPositionμ1μ2μ3σm2σf2σmf2σg2σe2method48 cm10860.50.52Powr1Powr2ipower30ML48.9810.068.025.942.940.100.650.252.040.980.990.91(7.36)(0.716)(0.388)(0.722)(1.197)(0.217)(0.694)(0.426)(0.207)REML48.549.967.955.902.950.090.640.571.970.970.990.86(9.62)(0.705)(0.428)(0.747)(1.501)(0.180)(0.598)(0.609)(0.226)03ML46.6210.058.026.080.092.760.690.302.050.971.000.89(8.11)(0.680)(0.442)(0.584)(0.201)(1.052)(0.722)(0.430)(0.193)REML48.710.047.985.910.092.950.610.571.980.960.980.88(7.40)(0.560)(0.506)(0.615)(0.195)(1.279)(0.663)(0.544)(0.225)12ML47.810.098.006.110.842.040.660.322.020.971.000.24(7.83)(0.634)(0.468)(0.666)(0.615)(0.994)(0.811)(0.475)(0.203)REML49.6610.027.965.901.081.730.680.612.000.970.990.09(11.06)(0.606)(0.548)(0.668)(0.679)(0.958)(0.707)(0.616)(0.216)Por1 and Por2 correspond to the overall QTL effect test (2.27)
calculated using the empirical and theoretical cutoffs, respectively; ipower
refers to the imprinting test power corresponding to test (2.29). See Table 2 for
explanations of other parameters.In reality, whether a QTL is imprinted or not is an
unknown prior. When a QTL has Mendelian effect and is not imprinted, is there
any power loss by analyzing with the proposed imprinting model? Or when a QTL
is actually imprinted, is there any power loss by analyzing with regular
variance components approach? To answer these two questions, we simulated data
under different scenarios and analyzed with both Mendelian and imprinting
models. The first and second columns in Table4
refer to the simulation and analysis models, respectively. M refers to the
Mendelian model without variance components partition and I refers to the
imprinting model with allelic-specific partition of the variance components. For comparison purpose, heritabilities are recorded instead of original
variance components estimates. The polygene and residual variances are fixed as
0.5 (hg2=0.083) and 2, respectively for all the simulation
scenarios. We first simulated data with one additive genetic effect without
partitioning variance into allelic-specific components. This is equivalent to
simulate data assuming the Mendelian model. A single additive variance
component of 3.5 is assumed which corresponds to a heritability of ha2=0.583. The second scenario is to simulate data with three
allelic-specific variance components. Simulation models I1 and I2 correspond to a
complete maternal imprinting model (i.e., hm2=0 and hf2=0.5) and a partial
maternal imprinting model (i.e., hm2=0.083 and hf2=0.417), respectively. The variance component σmf2 is assumed to
be 0.5 (hmf2=0.083) for I1 and I2. In all the simulations, we use the 20×20 design to
make the comparison. Similar results are expected under the other sampling
designs. Since the true variance components values for the imprinting model
are unknown when data are simulated assuming
Mendelian effect and vice versa, only standard deviations for these parameter
estimates are recorded (listed as italic font in the parentheses).Table 4
The power, MLEs of the QTL position, and effect parameters estimated based on 100
simulation replicates for data simulated with Mendelian and imprinting models
based on the20×20 sampling
design.
SimulationAnalysisPositionμ1μ2μ3ha2hg2σe2modelmodel48 cm1086hm2hf2hmf20.0832Power1Power2MM48.0810.1237.8906.0230.5440.0931.9821.001.00(2.981)(0.887)(0.730)(1.046)(0.121)(0.137)(0.250)I48.2609.9228.0816.0740.2630.2720.2390.0491.9831.001.00(3.015)(1.010)(0.720)(0.949)(0.097)(0.110)(0.124)(0.080)(0.249)I1M48.94010.0797.9136.0230.3230.1112.0220.860.96(11.660)(0.588)(0.511)(0.712)(0.287)(0.134)(0.205)I47.8610.0777.9006.0180.0080.4730.0890.0912.0150.950.98(7.524)(0.579)(0.575)(0.665)(0.022)(0.151)(0.099)(0.053)(0.86)I2M48.64010.0707.9176.0280.3310.1032.0230.910.98(8.389)(0.613)(0.517)(0.706)(0.278)(0.129)(0.195)I49.20010.0597.9266.0080.0790.3960.0890.0892.0230.930.98(9.818)(0.617)(0.518)(0.673)(0.075)(0.156)(0.097)(0.050)(0.184)M and I refer to Mendelian and imprinting
models, respectively. Simulated parameters for model M:(ha2=0.583); I1:(hm2,hf2,hmf2)=(0, 0.5, 0.083); I2:(hm2,hf2,hmf2)=( 0.083, 0.417, 0.083). Power1 and Power2 correspond to the power calculated using the
empirical cutoff and the theoretical threshold, respectively. The numbers given
in the parenthesis with normal and italic fonts correspond to the RMSEs and
standard errors of the parameter estimates, respectively. See Table 2 for other
explanations.The simulation results are summarized in Table4
analyzed with the ML method. When the simulated model is Mendelian, QTL
position is better estimated with the Mendelian model than with the imprinting
model. No remarkable difference in power is observed for both models. The
estimated parent-specific variances due to maternal and paternal alleles are
almost identical and no imprinting is detected. When data are simulated
assuming imprinting (model I1 and I2), large power
is observed when analyzed with the imprinting model. For example, the power is
86% when analyze the I1 imprinting data
by the Mendelian model. The power is increased to 95% when data are analyzed by
the imprinting model. When imprinting data are analyzed with the Mendelian
model, the major QTL variance is underestimated, and the polygene variance is
slightly overestimated. No remarkable differences are observed for the
estimation of the three fixed mean effects and the residual variance under all
simulation cases. In any case, the imprinting model performs better or no worse
than the Mendelian model in terms of power. In checking type I error rate based
on the theoretical threshold, we find the imprinting model has slightly higher
type I error rate compared with the Mendelian model. In real data analysis, it
is more important to control the false negatives than the false positives. Thus, it is safe to apply the imprinting model for data with any inheritance
pattern in this regard.


## 3.2.2. Multiple QTL Analysis

To see the
relative merit of multiple QTL analysis against single QTL analysis when
multiple QTL are located on the same linkage group, two QTL are simulated with
QTL 1 (denoted asQ1) located at
the second interval, 28cm away from the
first marker (ℳ1), and QTL 2 (denoted as Q2) located at
the fourth interval, 68 cm away from the first marker. Two simulation scenarios
are considered. The first scenario considers two nonimprinted QTL with equal
genetic effects. The second scenario assumes
Q1 is imprinted
and Q2 is not
imprinted. Simulated parameters for the two QTL are listed in Table 5. Data are
simulated assuming the 20×20 design. Parameters are estimated by the ML and REML approaches with 100 replicates.Table 5
The MLEs and
REMLs of the QTL position and effect parameters estimated based on 100
simulation replicates for data simulated with two QTL under the20×20 design. The
square roots of the mean squared errors are given in parentheses.
Q1Q2EstimationPositionσm2σf2σmf2Positionσm2σf2σmf2method28 cm0.2568 cm0.250.750.750.750.75ML31.060.9600.8190.38164.5600.8060.8910.414(12.825)(0.760)(0.693)(0.533)(12.126)(0.552)(0.679)(0.646)REML32.320.9490.9240.44064.600.9431.0020.461(13.951)(0.739)(0.686)(0.678)(12.169)(0.713)(0.711)(0.601)1.500.750.75ML29.521.5680.2100.60163.861.2100.5320.533(12.153)(0.872)(0.382)(0.836)(11.913)(0.932)(0.572)(0.847)REML31.221.6680.2820.43364.161.3550.6650.529(12.557)(0.881)(0.457)(0.598)(13.030)(1.128)(0.554)(0.650)Q1 and Q2 refer to two QTL located at 28 cm and 68 cm.Figure3 shows the LR profile plots for the single and
multiple QTL analyses. The single QTL model
indicates three major peaks. The highest peak for the single QTL analysis is
located at the wrong QTL interval where no QTL is assumed. The so-called
“ghost image” of QTL can be removed and the positions of the two QTL can
be precisely mapped on the chromosome by the multiple QTL model. Two clear
peaks indicating the correct QTL positions (arrow signs) are observed by the
multiple QTL analysis. However, we observe a remarkable reduction in LR values
by multiple QTL analysis compared to those by the single QTL analysis. Since
the threshold for multiple QTL analysis is unknown, we cannot make the
conclusion that multiple QTL analysis is less powerful than the single QTL
analysis. It is possible that we may gain accuracy in QTL position estimation
at the cost of power loss. Similar phenomenon and issues were also observed and
discussed in the literature [3, 25].The LR profile plot for single QTL
and multiple QTL analyses. The true QTL positions
are simulated at 28cm and 68 cm (see the arrow
sign). The dotted curve and the solid curve represent the LR profiles by single
QTL and multiple QTL analyses, respectively. The
left and right figures correspond to the LR profiles generated using the ML and
REML methods, respectively.
(a)(b)The results of the multiple QTL analysis are
summarized in Table5. The fixed mean effects, the polygene, and residual
variance components can be reasonably estimated with small RMSEs, similar
results shown in Table 2 for the 20×20 design and
hence are not reported here. Only the genetic factors for the two simulated
QTLs are reported. It can be seen that both ML and REML methods provide
reasonable parameter estimates and are very comparable. Under the first
simulation scenario in which both QTL are not imprinted, the genetic effects
are all slightly overestimated by both methods. This might be due to the
interference of the two QTL in the same linkage group. The multiple QTL model
may not completely block the effects of QTL outside of the tested interval. For
the second simulation scenario, an interesting pattern is observed. When one
QTL is imprinted (Q1), the maternal and paternal variance components for
the second one (Q2) tend to be estimated with bias in the direction as
the first imprinted QTL, that is, σm2 tends to be
overestimated and σf2 tends to be
underestimated. As we gain accuracy in QTL position estimation, we lose
precision for the parameter estimation. These effects are expected as described
in Zeng [3] and Xu and Atchley [25]. More investigations are needed in multiple
QTL analysis in order to maintain a good balance of QTL position and parameter
inference.


## 4. Discussion

Statistical
methods assuming fixed effect models for iQTL mapping in controlled outbred and
inbred lines have been proposed (e.g., [11, 14–16]). Considering the limitation
of fixed-effect models, a random model that estimates the QTL variance by
extending single line cross to multiple line crosses should be more powerful in
QTL variance inference [18]. The IBD-based variance components method assuming
random genetic effect for iQTL mapping has been developed in human linkage
analysis [7]. However, no study has been proposed to map iQTL using variance
components method with inbred or partially inbred line cross. In this article, we have first time presented an IBD-based variance components framework to
search for the existence and distribution of iQTL throughout the entire genome
in multiple experimental line crosses. The idea of the method is demonstrated through
a backcross design. It can also be extended to multiple F2 line crosses
using the sex-specific recombination information as proposed by Cui et al. [15].The key point of the proposed iQTL variance components
analysis is to partition the additive genetic variance into parent-specific
components. We have proposed a new parent-specific allelic sharing method which
characterizes the relatedness of parent-specific alleles between pairs of
individuals in a backcross pedigree. The calculation of parent-specific allelic
sharing is based on the information of the coefficient of coancestry. More
complicated calculation of the coefficient of coancestry can be found at Harris
[21]. The quantification of the coefficient of the coancestry proposed by
Harris [21] can also be utilized to calculate the parent-specific IBD sharing
in an inbred human population, and thus for iQTL mapping in inbred human
populations.There have been extensive studies in the literature
about various methods in the estimation of variance components in a
mixed-effect model framework. The ML and REML are two commonly applied methods
in variance components estimation with less-biased estimation by the REML
method. Simulations show that the ML method yields high precision in parameter
estimation but with relatively large bias than the REML method. Power analysis
indicates that the ML method is a little more powerful than the REML method but
with large type I error when testing imprinting. In terms of computing speed, the ML method is faster than the REML method. Thus, no single method dominates
the other. In terms of overall QTL test, we suggest to use the ML method for
the genome-wide linkage scan and use the REML method for the imprinting test.The effect of sampling design is investigated by
extensive simulations. Results indicate that one can always achieve large power
with large offspring size when the total sample size is fixed. The LR value
differences under different sampling designs are shown in Figure2. However, the combination of small families each with large offsprings gives poor
parameter estimation and large type I error for imprinting test (Table 2). As
the number of families increase, we observe less-biased parameter estimates for
both fixed and random effects, but with poor QTL position estimation and small
power. This information implies that it is necessary to enlarge the number of
families to improve precision of parameter estimation. Meanwhile, a balance of
family and offspring size is needed to maintain good QTL detection power and
position estimation. Our simulations indicate that for a fixed total sample
size (n=400), both 8×50 and 20×20 designs
yield comparable results and both designs outperform the other two designs
(Table 2). Moreover, the 20×20 design
produces relatively small type I error in imprinting test. With the 20×20 design, results in Table 4 indicate that the imprinting model is better or as good as
the regular Mendelian analysis without considering imprinting. In real data
analysis, it should be safe to apply the proposed imprinting model for data
with any imprinting pattern.In this study, we have extended the single
marker-based analysis to an interval-based mapping for genome-wide scan and
testing of iQTL effects. Considering the interference of QTL located on the
same linkage group, we have extended the single QTL model to multiple QTL
analysis following the derivation of Xu and Atchley [25]. Simulation results
indicate the relative merits of the multiple QTL analysis with improved QTL
position inference, but with possible power loss (Figure 3). This, however, has
been a common issue in multiple QTL modeling (see [3, 25]). More investigations
are needed in deriving efficient and robust multiple QTL mapping models to
improve precision without suffering too much from power loss.The theoretical distribution for the likelihood ratio
test has been a challenging problem in QTL mapping. Dupuis and Siegmund [34]
first proposed theoretical properties for LR test statistics in a genome-wide
linkage scan for QTL in an interval mapping
framework with a fixed-effect model. Currently, most
linkage analysis using the variance components method assumes that the LR test
statistic follows a mixture of chi-square distribution [35]. The mixture
distribution is derived following Self and Liang [28]. With multiple testings
and multiple nuisance parameters in a genome-wide scan, the assumptions to get
the mixture chi-square distribution may not satisfy. Moreover, the multivariate
normal assumption for the phenotypic data required to get the mixture
distribution may not even valid. No theoretical work has been done to
investigate this in an IBD-based variance components linkage mapping. Our
simulations indicate that the theoretical threshold calculated from the mixture
chi-square distribution is smaller than the simulated cutoff. Thus, the power
calculated with the theoretical threshold is slightly inflated. A modified
mixture chi-quare distribution may be more appropriate. More theoretical
investigations are needed in this regard.



---

*Source: 689489-2009-03-26.xml*