sim_dof1.R

# simulating DOF1 data
# cross between DO individuals and some 9th strain

library(simcross)
library(qtl2geno)

# load some DO data, for the founder genotypes and map
file <- paste0("https://raw.githubusercontent.com/rqtl/",
               "qtl2data/master/DO_Recla/recla.zip")
do_recla <- read_cross2(file)

# simulate DOF1 pedigree (n=144)
ped <- sim_dof1_pedigree(20, nkids_per=1)

# length of chr 1
L <- max(do_recla$gmap[[1]])

# simulate genotypes along chr 1, omit all but DOF1 individuals
xodata <- sim_from_pedigree(ped, L=L)[ped$dof1]
xodata <- collapse_do_alleles(xodata) # collapse alleles 9-16 to 1-8

# simulate founder genotypes for the backcross strain
# (here just doing all SNPs with 50% MAF)
founder_geno <- do_recla$founder_geno[[1]]
founder_geno <- rbind(founder_geno,
                      I=sample(c(1,3), ncol(founder_geno), replace=TRUE))

# drop markers with missing founder genotypes
todrop <- colSums(is.na(founder_geno) | founder_geno==0)>0

# convert DOF1 data to SNP genotypes
geno <- convert2geno(xodata, do_recla$gmap[[1]][!todrop], founder_geno[,!todrop,drop=FALSE])

# number of individuals
n <- sum(ped$dof1)

# create a cross2 object
x <- list(geno=list("1"=geno),
          founder_geno=list("1"=founder_geno[,!todrop,drop=FALSE]),
          gmap=list("1"=do_recla$gmap[[1]][!todrop]),
          crosstype="dof1",
          is_female=(ped$sex[ped$dof1]==0),
          is_x_chr=c("1"=FALSE),
          cross_info=cbind(as.integer(ped$gen[ped$dof1])))
rownames(x$geno[[1]]) <- names(x$is_female) <- rownames(x$cross_info) <- 1:n
class(x) <- c("cross2", "list")

# simulate a QTL
# grab genotype at position 35
qtl <- get_geno(xodata, 35)[,1] # first column is 1-8; second column is == 9

# very large effect; alleles 1-4 vs 5-8
pheno <- rnorm(n, (qtl<4)*2, 1)
names(pheno) <- 1:n

x$pheno <- cbind("pheno"=pheno)

# QTL analysis
# calc genotype probabilities
pr <- calc_genoprob(x, step=1)

# plain Haley-Knott
library(qtl2scan)
out <- scan1(pr, x$pheno)

# plot LOD curve
library(qtl2plot)
plot(out)

# estimate effects
coef <- scan1coef(pr, x$pheno)
plot_coefCC(coef) # use CC colors

# plot effects + LOD curve
plot_coefCC(coef, scan1_output=out)

# estimate effects with blup
blup <- scan1blup(pr, x$pheno)
plot_coefCC(blup, scan1_output=out)