###################
# Normal Model
###################
data(hyper)
# Genotype probabilities for EM and H-K
## Not run: 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, thr=3)
summary(out.hk, thr=3)
# An alternate method of summarizing:
# patch them together and then summarize
out <- c(out.em, out.hk)
summary(out, thr=3, format="allpeaks")
# 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
## Not run: hyper.c1n4 <- sim.geno(subset(hyper, chr=c(1,4)),
# step=2.5, n.draws=8)
# ## End(Not run)
out.imp <- scanone(hyper.c1n4, method="imp")
summary(out.imp, thr=3)
# Plot all three results
plot(out.imp, out.hk, out.em, chr=c(1,4), lty=1,
col=c("red","blue","black"))
# extended Haley-Knott
out.ehk <- scanone(hyper, method="ehk")
plot(out.hk, out.em, out.ehk, chr=c(1,4))
# Permutation tests
## Not run: permo <- scanone(hyper, method="hk", n.perm=1000)
# Threshold from the permutation test
summary(permo, alpha=c(0.05, 0.10))
# Results above the 0.05 threshold
summary(out.hk, perms=permo, alpha=0.05)
####################
# scan with square-root of phenotype
# (Note that pheno.col can be a vector of phenotype values)
####################
out.sqrt <- scanone(hyper, pheno.col=sqrt(pull.pheno(hyper, 1)))
plot(out.em - out.sqrt, ylim=c(-0.1,0.1),
ylab="Difference in LOD")
abline(h=0, lty=2, col="gray")
####################
# Stratified permutations
####################
extremes <- (nmissing(hyper)/totmar(hyper) < 0.5)
## Not run: operm.strat <- scanone(hyper, method="hk", n.perm=1000,
# perm.strata=extremes)
# ## End(Not run)
summary(operm.strat)
####################
# X-specific permutations
####################
data(fake.f2)
## Not run: fake.f2 <- calc.genoprob(fake.f2, step=2.5)
# genome scan
out <- scanone(fake.f2, method="hk")
# X-chr-specific permutations
## Not run: operm <- scanone(fake.f2, method="hk", n.perm=1000, perm.Xsp=TRUE)
# thresholds
summary(operm)
# scanone summary with p-values
summary(out, perms=operm, alpha=0.05, pvalues=TRUE)
###################
# Non-parametric
###################
out.np <- scanone(hyper, model="np")
summary(out.np, thr=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)
## Not run: listeria <- calc.genoprob(listeria,step=2.5)
out.2p <- scanone(listeria, model="2part", upper=TRUE)
summary(out.2p, thr=c(5,3,3), format="allpeaks")
# Plot all three LOD scores together
plot(out.2p, out.2p, out.2p, lodcolumn=c(2,3,1), lty=1, chr=c(1,5,13),
col=c("red","blue","black"))
# Permutation test
## Not run: permo <- scanone(listeria, model="2part", upper=TRUE,
# n.perm=1000)
# ## End(Not run)
# Thresholds
summary(permo)
###################
# Binary model
###################
binphe <- as.numeric(pull.pheno(listeria,1)==264)
out.bin <- scanone(listeria, pheno.col=binphe, model="binary")
summary(out.bin, thr=3)
# Plot LOD for binary model with LOD(p) from 2-part model
plot(out.bin, out.2p, lodcolumn=c(1,2), lty=1, col=c("black", "red"),
chr=c(1,5,13))
# Permutation test
## Not run: permo <- scanone(listeria, pheno.col=binphe, model="binary",
# n.perm=1000)
# ## End(Not run)
# Thresholds
summary(permo)
###################
# Covariates
###################
data(fake.bc)
## Not run: fake.bc <- calc.genoprob(fake.bc, step=2.5)
# genome scans without covariates
out.nocovar <- scanone(fake.bc)
# genome scans with covariates
ac <- pull.pheno(fake.bc, c("sex","age"))
ic <- pull.pheno(fake.bc, "sex")
out.covar <- scanone(fake.bc, pheno.col=1,
addcovar=ac, intcovar=ic)
summary(out.nocovar, thr=3)
summary(out.covar, thr=3)
plot(out.covar, out.nocovar, chr=c(2,5,10))
Run the code above in your browser using DataLab