Skip to content

Instantly share code, notes, and snippets.

@kbroman
Created May 20, 2016 02:08
Show Gist options
  • Save kbroman/b1805bc027bb558e692f1f7191355f7c to your computer and use it in GitHub Desktop.
Save kbroman/b1805bc027bb558e692f1f7191355f7c to your computer and use it in GitHub Desktop.
# simulate RIL
#
# QTL on chr 6 at 19.9 cM
# a = -6.85, hsq = 0.1768
# load R/qtl
library(qtl)
# QTL effect and residual SD
a <- -6.85
hsq <- 0.1768
sigma <- sqrt( a^2 * (1-hsq) / hsq)
# an example map
data(map10)
# add QTL position to map
map10[[6]] <- sort(c(map10[[6]], "qtl"=19.9))
# simulate N=100 RIL by selfing
N <- 100
simRIL <- sim.cross(map10, type="riself", n.ind=N)
# grab the genotypes at the qtl
qtlg <- pull.geno(simRIL)[,"qtl"]
# drop the QTL from the data
simRIL <- drop.markers(simRIL, "qtl")
# add 20% missing data
for(chr in chrnames(simRIL)) {
g <- simRIL$geno[[chr]]$data
ng <- prod(dim(g))
nmis <- rbinom(1, ng, 0.2)
g[sample(1:ng, nmis)] <- NA
simRIL$geno[[chr]]$data <- g
}
# simulate the phenotype
qtlg <- (qtlg-1.5)*2 # convert from 1/2, to -1/+1
simRIL$pheno[,1] <- qtlg * a + rnorm(nind(simRIL), 0, sigma)
# genome scan
simRIL <- calc.genoprob(simRIL, step=1)
out <- scanone(simRIL, method="hk")
# fitqtl to get estimated effect
mx <- max(out, chr=6) # biggest peak on chr 6
qtl <- makeqtl(simRIL, chr=mx[[1]], pos=mx[[2]], what="prob")
out.fq <- fitqtl(simRIL, qtl=qtl, method="hk")
summary(out.fq)
# estimated hsq
hsq_est <- out.fq[[1]][1,5]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment