Created
May 20, 2016 02:08
-
-
Save kbroman/b1805bc027bb558e692f1f7191355f7c to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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