###################
# Normal Model
###################
data(hyper)
# Genotype probabilities for EM and H-K
hyper <- calc.genoprob(hyper, step=2.5)
out.em <- scanone(hyper, method="em")
out.hk <- scanone(hyper, method="hk")
# Summarize results: peaks above 3
summary(out.em, 3)
summary(out.hk, 3)
# Plot the results
plot(out.hk, out.em)
plot(out.hk, out.em, chr=c(1,4), lty=1, col=c("blue","black"))
# Imputation; first need to run sim.geno
# Do just chromosomes 1 and 4, to save time
hyper.c1n4 <- sim.geno(subset(hyper, chr=c(1,4)),
step=2.5, n.draws=64)
<testonly>hyper.c1n4 <- sim.geno(subset(hyper, chr=c(1,4)),
step=2.5, n.draws=16)</testonly>
out.imp <- scanone(hyper.c1n4, method="imp")
summary(out.imp, 3)
# Plot all three results
plot(out.imp, out.hk, out.em, chr=c(1,4), lty=1,
col=c("red","blue","black"))
# Permutation tests
permo <- scanone(hyper, method="hk", n.perm=1000)<testonly>permo <- scanone(hyper, method="hk", n.perm=3)</testonly>quantile(permo, 0.95)
###################
# Non-parametric
###################
out.np <- scanone(hyper, model="np")
summary(out.np, 3)
# Plot with previous results
plot(out.np, chr=c(1,4), lty=1, col="green")
plot(out.imp, out.hk, out.em, chr=c(1,4), lty=1,
col=c("red","blue","black"), add=TRUE)
###################
# Two-part Model
###################
data(listeria)
listeria <- calc.genoprob(listeria,step=2.5)
out.2p <- scanone(listeria, model="2part", upper=TRUE)
summary(out.2p, 5)
# Plot all three LOD scores together
plot(out.2p[,-3], out.2p[,-(3:4)], out.2p, lty=1, chr=c(1,5,13),
col=c("red","blue","black"))
# Permutation test
permo <- scanone(listeria, model="2part", upper=TRUE,
n.perm=1000)<testonly>permo <- scanone(listeria, model="2part", upper=TRUE,
n.perm=3)</testonly>apply(permo, 2, quantile, 0.95)
###################
# Binary model
###################
listeria <- subset(listeria, ind=!is.na(listeria$pheno[,1]))
listeria$pheno[,2] <- rep(0,nind(listeria))
listeria$pheno[listeria$pheno[,1]==264,2] <- 1
out.bin <- scanone(listeria, pheno.col=2, model="binary")
summary(out.bin, 3)
# Plot LOD for binary model with LOD(p) from 2-part model
plot(out.bin, out.2p[,-3], lty=1, col=c("black", "red"),
chr=c(1,5,13))
# Permutation test
permo <- scanone(listeria, pheno.col=2, model="binary",
n.perm=1000)<testonly>permo <- scanone(listeria, pheno.col=2, model="binary",
n.perm=3)</testonly>quantile(permo, 0.95)
###################
# Covariates
###################
data(fake.bc)
plot(fake.bc)
fake.bc <- calc.genoprob(fake.bc, step=2.5)
# genome scans without covariates
out.nocovar <- scanone(fake.bc)
# genome scans with covariates
ac <- fake.bc$pheno[,c("sex","age")]
ic <- fake.bc$pheno[,"sex"]
out.covar <- scanone(fake.bc, pheno.col=1, addcov=ac, intcov=ic)
summary(out.nocovar,3)
summary(out.covar,3)
plot(out.covar,out.nocovar,chr=c(2,5,10))
Run the code above in your browser using DataLab