Assess estimated QTL effects by fitqtl

qtleff_riself.R

library(qtl)

# simulate data
set.seed(57173190)
data(map10)
map10 <- map10[1:2]
x <- sim.cross(map10, n.ind=128, type="riself")
g1 <- pull.geno(x)[,"D1M6"]
g2 <- pull.geno(x)[,"D2M6"]
x$pheno[,1] <- rnorm(nind(x)) + 0.5*(g1==2) + 1.5*(g2==2)
# additive effects should be about 0.25 and 0.75

# direct estimates of effects
diff(tapply(x$pheno[,1], g1, mean))/2 # 0.19
diff(tapply(x$pheno[,1], g2, mean))/2 # 0.66

# estimates fitting the two QTL jointly in a linear model
lm.out <- lm(x$pheno[,1] ~ g1 + g2)
lm.out$coef[-1]/2
# estimates are 0.299 and 0.705

# fitqtl by imputation
x <- sim.geno(x, n.draws=64)
qtl <- makeqtl(x, c(1,2), c(48.8, 51.8), what="draws")
out.fq <- fitqtl(x, qtl=qtl, get.ests=TRUE, dropone=FALSE)
summary(out.fq)
# estimates are 0.299 and 0.705

# fitqtl by Haley-Knott regression
x <- calc.genoprob(x)
qtl <- makeqtl(x, c(1,2), c(48.8, 51.8), what="prob")
out.fq.hk <- fitqtl(x, qtl=qtl, get.ests=TRUE, dropone=FALSE, method="hk")
summary(out.fq.hk)
# estimates are 0.299 and 0.705

qtleff_riself_selgeno.R

library(qtl)

# simulate data, selective genotyping
set.seed(91324919)
data(map10)
map10 <- map10[1:2]
n <- 128
frac <- 0.4
x <- sim.cross(map10, n.ind=n/frac, type="riself")
g1 <- pull.geno(x)[,"D1M6"]
g2 <- pull.geno(x)[,"D2M6"]
x$pheno[,1] <- rnorm(nind(x)) + 0.5*(g1==2) + 1.5*(g2==2)
# additive effects should be about 0.25 and 0.75

# omit genotypes for the middle 60%
o <- order(x$pheno[,1])
keep <- o[c(1:(n/2), nind(x)-((1:(n/2))-1))]
for(i in 1:nchr(x))
  x$geno[[i]]$data[-keep,] <- NA

# pull out genotypes again; 60% missing
g1 <- pull.geno(x)[,"D1M6"]
g2 <- pull.geno(x)[,"D2M6"]

# direct estimates of effects
diff(tapply(x$pheno[,1], g1, mean))/2 # 0.32
diff(tapply(x$pheno[,1], g2, mean))/2 # 1.57

# estimates fitting the two QTL jointly in a linear model
lm.out <- lm(x$pheno[,1] ~ g1 + g2)
lm.out$coef[-1]/2
# estimates are 0.084 and 1.55

# fitqtl by imputation
x <- sim.geno(x, n.draws=1024)
qtl <- makeqtl(x, c(1,2), c(48.8, 51.8), what="draws")
out.fq <- fitqtl(x, qtl=qtl, get.ests=TRUE, dropone=FALSE)
summary(out.fq)
# estimates are 0.083 and 0.686

# fitqtl by Haley-Knott regression (works poorly here)
x <- calc.genoprob(x)
qtl <- makeqtl(x, c(1,2), c(48.8, 51.8), what="prob")
out.fq.hk <- fitqtl(x, qtl=qtl, get.ests=TRUE, dropone=FALSE, method="hk")
summary(out.fq.hk)
# estimates are 0.092 and 1.54

# fitqtl by imputation, just using the genotyped individuals
x <- subset(x, ind=keep)
qtl <- makeqtl(x, c(1,2), c(48.8, 51.8), what="draws")
out.fq <- fitqtl(x, qtl=qtl, get.ests=TRUE, dropone=FALSE)
summary(out.fq)
# estimates are 0.084 and 1.55

# fitqtl by H-K regr, just using the genotyped individuals
qtl <- makeqtl(x, c(1,2), c(48.8, 51.8), what="prob")
out.fq.hk <- fitqtl(x, qtl=qtl, get.ests=TRUE, dropone=FALSE, method="hk")
summary(out.fq.hk)
# estimates are 0.084 and 1.55