R/simulate.R

In kbroman/qtl: Tools for Analyzing QTL Experiments

######################################################################
#
# simulate.R
#
# copyright (c) 2001-2019, Karl W Broman
# last modified Dec, 2019
# first written Apr, 2001
#
#     This program is free software; you can redistribute it and/or
#     modify it under the terms of the GNU General Public License,
#     version 3, as published by the Free Software Foundation.
#
#     This program is distributed in the hope that it will be useful,
#     but without any warranty; without even the implied warranty of
#     merchantability or fitness for a particular purpose.  See the GNU
#     General Public License, version 3, for more details.
#
#     A copy of the GNU General Public License, version 3, is available
#     at http://www.r-project.org/Licenses/GPL-3
#
# Part of the R/qtl package
# Contains: sim.map, sim.cross, sim.cross.bc, sim.cross.f2,
#           sim.cross.4way, sim.bcg
#
######################################################################

######################################################################
#
# sim.map: simulate a genetic map
#
######################################################################

sim.map <-
    function(len=rep(100,20), n.mar=10, anchor.tel=TRUE, include.x=TRUE,
             sex.sp=FALSE, eq.spacing=FALSE)
{
    if(length(len)!=length(n.mar) && length(len)!=1 && length(n.mar)!=1)
        stop("Lengths of vectors len and n.mar do not conform.")

    # make vectors the same length
    if(length(len) == 1) len <- rep(len,length(n.mar))
    else if(length(n.mar) == 1) n.mar <- rep(n.mar,length(len))

    n.chr <- length(n.mar)

    map <- vector("list",n.chr)
    names(map) <- as.character(1:n.chr)
    if(include.x) names(map)[n.chr] <- "X"

    for(i in 1:n.chr) {
        if(anchor.tel) {
            if(n.mar[i] < 2) n.mar[i] <- 2
            map[[i]] <- c(0,len[i])
            if(n.mar[i] > 2) {
                if(!eq.spacing)
                    map[[i]] <- sort(c(map[[i]],runif(n.mar[i]-2,0,len[i])))
                else # equal spacing
                    map[[i]] <- seq(0,len[i],length=n.mar[i])
            }
        }
        else {
            if(!eq.spacing) {
                map[[i]] <- sort(runif(n.mar[i],0,len[i]))
                map[[i]] <- map[[i]] - min(map[[i]])
            }
            else {  # equal spacing
                map[[i]] <- seq(0,len[i],length=n.mar[i]+1)
                map[[i]] <- map[[i]][-1] - map[[i]][2]/2
            }
        }
        names(map[[i]]) <- paste("D", names(map)[i], "M", 1:n.mar[i], sep="")
        class(map[[i]]) <- "A"
    }

    if(sex.sp) {
        for(i in 1:n.chr) {
            if(eq.spacing) tempmap <- map[[i]]
            else {
                if(anchor.tel) {
                    if(n.mar[i] < 2) n.mar[i] <- 2
                    tempmap <- c(0,len[i])
                    if(n.mar[i] > 2)
                        tempmap <- sort(c(tempmap,runif(n.mar[i]-2,0,len[i])))
                }
                else {
                    tempmap <- sort(runif(n.mar[i],0,len[i]))
                    tempmap <- tempmap - min(tempmap)
                }
            }
            map[[i]] <- rbind(map[[i]],tempmap)
            dimnames(map[[i]]) <- list(NULL,paste("D", names(map)[i], "M", 1:n.mar[i], sep=""))
            class(map[[i]]) <- "A"

            if(include.x && i==n.chr)  # if X chromosome, force no recombination in male
                map[[i]][2,] <- rep(0,ncol(map[[i]]))
        }
    }

    if(include.x) class(map[[n.chr]]) <- "X"

    class(map) <- "map"
    map
}

######################################################################
#
# sim.cross: Simulate an experimental cross
#
# Note: These functions are a bit of a mess.  I was in the "get it to
#       work without worrying about efficiency" mode while writing it.
#       Sorry!
#
######################################################################

sim.cross <-
    function(map, model=NULL, n.ind=100,
             type=c("f2", "bc", "4way", "risib", "riself",
             "ri4sib", "ri4self", "ri8sib", "ri8self","bcsft"),
             error.prob=0, missing.prob=0, partial.missing.prob=0,
             keep.qtlgeno=TRUE, keep.errorind=TRUE, m=0, p=0,
             map.function=c("haldane","kosambi","c-f","morgan"),
             founderGeno, random.cross=TRUE, ...)
{
    type <- match.arg(type)
    map.function <- match.arg(map.function)

    # don't let error.prob be exactly zero (or >1)
    if(error.prob < 1e-50) error.prob <- 1e-50
    if(error.prob > 1) {
        error.prob <- 1-1e-50
        warning("error.prob shouldn't be > 1!")
    }

    # 2-way RIL by sibmating or selfing
    if(type=="risib" || type=="riself") {
        if(!is.null(model))
            warning('"model" argument currently ignored in simulating RILs')

        if(type=="risib") type <- "sibmating"
        else type <- "selfing"
        cross <- sim.ril(map, n.ind, type, "2", m=m, p=p,
                         error.prob=error.prob, missing.prob=missing.prob)
        cross$cross <- NULL

        return(cross)
    }
    # 4- or 8-way RIL by sibmating or selfing
    if(type=="ri4sib" || type=="ri4self" || type=="ri8sib" || type=="ri8self") {
        if(!is.null(model))
            warning('"model" argument currently ignored in simulating RILs')

        if(substr(type, 4, nchar(type))=="self") crosstype <- "selfing"
        else crosstype <- "sibmating"
        n.str <- substr(type, 3, 3)
        cross <- sim.ril(map, n.ind, crosstype, n.str, m=m, p=p,
                         random.cross=random.cross,
                         error.prob=0,
                         missing.prob=missing.prob)

        rcross <- convertMWril(cross, founderGeno, error.prob=error.prob)
        for(i in names(cross$geno))
            if(!("truegeno" %in% names(rcross$geno[[i]])))
                rcross$geno[[i]]$truegeno <- cross$geno[[i]]$data

        # remove "un" from cross type
        class(rcross) <- c(sub("un$", "", crosstype(cross)), "cross")

        fg <- t(founderGeno[[1]])
        if(length(founderGeno)>1)
            for(i in 2:length(founderGeno))
                fg <- cbind(fg, t(founderGeno[[i]]))
        colnames(fg) <- markernames(rcross)
        rcross$founderGeno <- fg
        return(rcross)
    }

    # sort the model matrix
    if(!is.null(model) && is.matrix(model))
        model <- model[order(model[,1],model[,2]),]

    if(type=="bc")
        cross <- sim.cross.bc(map,model,n.ind,error.prob,missing.prob,
                              keep.errorind,m,p,map.function)
    else if(type=="f2")
        cross <- sim.cross.f2(map,model,n.ind,error.prob,missing.prob,
                              partial.missing.prob,keep.errorind,
                              m,p,map.function)
    else if(type=="bcsft")
        cross <- sim.cross.bcsft(map,model,n.ind,error.prob,missing.prob,
                                 partial.missing.prob,keep.errorind,
                                 m,p,map.function, ...)
    else
        cross <- sim.cross.4way(map,model,n.ind,error.prob,missing.prob,
                                partial.missing.prob,keep.errorind,
                                m,p,map.function)

    # remove QTL genotypes from data and, if keep.qtlgeno=TRUE,
    #     place them in cross$qtlgeno
    qtlgeno <- NULL
    for(i in 1:nchr(cross)) {
        o <- grep("^QTL[0-9]+", colnames(cross$geno[[i]]$data))
        if(length(o) != 0) {
            qtlgeno <- cbind(qtlgeno, cross$geno[[i]]$data[,o,drop=FALSE])
            cross$geno[[i]]$data <- cross$geno[[i]]$data[,-o,drop=FALSE]
            if(is.matrix(cross$geno[[i]]$map))
                cross$geno[[i]]$map <- cross$geno[[i]]$map[,-o,drop=FALSE]
            else
                cross$geno[[i]]$map <- cross$geno[[i]]$map[-o]
        }
    }
    if(keep.qtlgeno) cross$qtlgeno <- qtlgeno

    # store genotype data as integers
    for(i in 1:nchr(cross))
        storage.mode(cross$geno[[i]]$data) <- "integer"

    if(is.null(names(cross$geno)))
        names(cross$geno) <- 1:length(cross$geno)

    cross
}


######################################################################
#
# sim.cross.bc
#
######################################################################

sim.cross.bc <-
    function(map,model,n.ind,error.prob,missing.prob,
             keep.errorind,m,p,map.function)
{
    if(map.function=="kosambi") mf <- mf.k
    else if(map.function=="c-f") mf <- mf.cf
    else if(map.function=="morgan") mf <- mf.m
    else mf <- mf.h

    if(any(sapply(map,is.matrix)))
        stop("Map must not be sex-specific.")

    chr.type <- sapply(map, chrtype)

    n.chr <- length(map)

    if(is.null(model)) n.qtl <- 0
    else {
        if(!((!is.matrix(model) && length(model) == 3) ||
             (is.matrix(model) && ncol(model) == 3)))
            stop("Model must be a matrix with 3 columns (chr, pos and effect).")
        if(!is.matrix(model)) model <- rbind(model)
        n.qtl <- nrow(model)
        if(any(model[,1] < 0 | model[,1] > n.chr))
            stop("Chromosome indicators in model matrix out of range.")
        model[,2] <- model[,2]+1e-14 # so QTL not on top of marker
    }

    # if any QTLs, place qtls on map
    if(n.qtl > 0) {
        for(i in 1:n.qtl) {
            temp <- map[[model[i,1]]]
            if(model[i,2] < min(temp)) {
                temp <- c(model[i,2],temp)
                names(temp)[1] <- paste("QTL",i,sep="")
            }
            else if(model[i,2] > max(temp)) {
                temp <- c(temp,model[i,2])
                names(temp)[length(temp)] <- paste("QTL",i,sep="")
            }
            else {
                j <- max((seq(along=temp))[temp < model[i,2]])
                temp <- c(temp[1:j],model[i,2],temp[(j+1):length(temp)])
                names(temp)[j+1] <- paste("QTL",i,sep="")
            }
            map[[model[i,1]]] <- temp
        }
    }

    geno <- vector("list", n.chr)
    names(geno) <- names(map)
    n.mar <- sapply(map,length)
    mar.names <- lapply(map,names)

    for(i in 1:n.chr) {
        # simulate genotype data
        thedata <- sim.bcg(n.ind, map[[i]], m, p, map.function)
        dimnames(thedata) <- list(NULL,mar.names[[i]])

        geno[[i]] <- list(data = thedata, map = map[[i]])
        class(geno[[i]]) <- chr.type[i]
        class(geno[[i]]$map) <- NULL

    } # end loop over chromosomes

    # simulate phenotypes
    pheno <- rnorm(n.ind,0,1)

    if(n.qtl > 0) {
        # find QTL positions in genotype data
        QTL.chr <- QTL.loc <- NULL
        for(i in 1:n.chr) {
            o <- grep("^QTL[0-9]+",mar.names[[i]])
            if(length(o)>0) {
                QTL.chr <- c(QTL.chr,rep(i,length(o)))
                QTL.loc <- c(QTL.loc,o)
            }
        }

        # incorporate QTL effects
        for(i in 1:n.qtl) {
            QTL.geno <- geno[[QTL.chr[i]]]$data[,QTL.loc[i]]
            pheno[QTL.geno==2] <- pheno[QTL.geno==2] + model[i,3]
        }

    } # end simulate phenotype

    n.mar <- sapply(geno, function(a) length(a$map))

    # add errors
    if(error.prob > 0) {
        for(i in 1:n.chr) {
            a <- sample(0:1,n.mar[i]*n.ind,replace=TRUE,
                        prob=c(1-error.prob,error.prob))
            geno[[i]]$data[a == 1] <- 3 - geno[[i]]$data[a == 1]
            if(keep.errorind) {
                errors <- matrix(0,n.ind,n.mar[i])
                errors[a==1] <- 1
                colnames(errors) <- colnames(geno[[i]]$data)
                geno[[i]]$errors <- errors
            }
        }
    }

    # add missing
    if(missing.prob > 0) {
        for(i in 1:n.chr) {
            o <- grep("^QTL[0-9]+",mar.names[[i]])
            if(length(o)>0)
                x <- geno[[i]]$data[,o]
            geno[[i]]$data[sample(c(TRUE,FALSE),n.mar[i]*n.ind,replace=TRUE,
                                  prob=c(missing.prob,1-missing.prob))] <- NA
            if(length(o)>0)
                geno[[i]]$data[,o] <- x
        }
    }

    pheno <- data.frame(phenotype=pheno, stringsAsFactors=TRUE)

    cross <- list(geno=geno,pheno=pheno)
    class(cross) <- c("bc","cross")

    cross
}

######################################################################
#
# sim.cross.f2
#
######################################################################

sim.cross.f2 <-
    function(map,model,n.ind,error.prob,missing.prob,partial.missing.prob,
             keep.errorind,m,p,map.function)
{
    if(map.function=="kosambi") mf <- mf.k
    else if(map.function=="c-f") mf <- mf.cf
    else if(map.function=="morgan") mf <- mf.m
    else mf <- mf.h

    if(any(sapply(map,is.matrix)))
        stop("Map must not be sex-specific.")

    # chromosome types
    chr.type <- sapply(map, chrtype)

    n.chr <- length(map)
    if(is.null(model)) n.qtl <- 0
    else {
        if(!((!is.matrix(model) && length(model) == 4) ||
             (is.matrix(model) && ncol(model) == 4))) {
            stop("Model must be a matrix with 4 columns (chr, pos and effects).")
        }
        if(!is.matrix(model)) model <- rbind(model)
        n.qtl <- nrow(model)
        if(any(model[,1] < 0 | model[,1] > n.chr))
            stop("Chromosome indicators in model matrix out of range.")
        model[,2] <- model[,2]+1e-14 # so QTL not on top of marker
    }

    # if any QTLs, place qtls on map
    if(n.qtl > 0) {
        for(i in 1:n.qtl) {
            temp <- map[[model[i,1]]]
            if(model[i,2] < min(temp)) {
                temp <- c(model[i,2],temp)
                names(temp)[1] <- paste("QTL",i,sep="")
            }
            else if(model[i,2] > max(temp)) {
                temp <- c(temp,model[i,2])
                names(temp)[length(temp)] <- paste("QTL",i,sep="")
            }
            else {
                j <- max((seq(along=temp))[temp < model[i,2]])
                temp <- c(temp[1:j],model[i,2],temp[(j+1):length(temp)])
                names(temp)[j+1] <- paste("QTL",i,sep="")
            }
            map[[model[i,1]]] <- temp
        }
    }

    geno <- vector("list", n.chr)
    names(geno) <- names(map)
    n.mar <- sapply(map,length)
    mar.names <- lapply(map,names)

    for(i in 1:n.chr) {

        # simulate genotype data
        thedata <- sim.bcg(n.ind, map[[i]], m, p, map.function)
        dimnames(thedata) <- list(NULL,mar.names[[i]])

        if(chr.type[i] != "X")
            thedata <- thedata + sim.bcg(n.ind, map[[i]], m, p, map.function) - 1

        geno[[i]] <- list(data = thedata, map = map[[i]])
        class(geno[[i]]) <- chr.type[i]
        class(geno[[i]]$map) <- NULL

    } # end loop over chromosomes

    # simulate phenotypes
    pheno <- rnorm(n.ind,0,1)

    if(n.qtl > 0) {
        # find QTL positions in genotype data
        QTL.chr <- QTL.loc <- NULL
        for(i in 1:n.chr) {
            o <- grep("^QTL[0-9]+",mar.names[[i]])
            if(length(o)>0) {
                QTL.chr <- c(QTL.chr,rep(i,length(o)))
                QTL.loc <- c(QTL.loc,o)
            }
        }

        # incorporate QTL effects
        for(i in 1:n.qtl) {
            QTL.geno <- geno[[QTL.chr[i]]]$data[,QTL.loc[i]]
            pheno[QTL.geno==1] <- pheno[QTL.geno==1] - model[i,3]
            pheno[QTL.geno==2] <- pheno[QTL.geno==2] + model[i,4]
            pheno[QTL.geno==3] <- pheno[QTL.geno==3] + model[i,3]
        }

    } # end simulate phenotype

    n.mar <- sapply(geno, function(a) length(a$map))

    # add errors
    if(error.prob > 0) {
        for(i in 1:n.chr) {
            if(chr.type[i]=="X") {
                a <- sample(0:1,n.mar[i]*n.ind,replace=TRUE,
                            prob=c(1-error.prob,error.prob))
                geno[[i]]$data[a == 1] <- 3 - geno[[i]]$data[a == 1]
            }
            else {
                a <- sample(0:2,n.mar[i]*n.ind,replace=TRUE,
                            prob=c(1-error.prob,error.prob/2,error.prob/2))
                if(any(a>0 & geno[[i]]$data==1))
                    geno[[i]]$data[a>0 & geno[[i]]$data==1] <-
                        (geno[[i]]$data+a)[a>0 & geno[[i]]$data==1]
                if(any(a>0 & geno[[i]]$data==2)) {
                    geno[[i]]$data[a>0 & geno[[i]]$data==2] <-
                        (geno[[i]]$data+a)[a>0 & geno[[i]]$data==2]
                    geno[[i]]$data[geno[[i]]$data>3] <- 1
                }
                if(any(a>0 & geno[[i]]$data==3))
                    geno[[i]]$data[a>0 & geno[[i]]$data==3] <-
                        (geno[[i]]$data-a)[a>0 & geno[[i]]$data==3]
            }

            if(keep.errorind) {
                errors <- matrix(0,n.ind,n.mar[i])
                errors[a>0] <- 1
                colnames(errors) <- colnames(geno[[i]]$data)
                geno[[i]]$errors <- errors
            }

        } # end loop over chromosomes
    } # end simulate genotyping errors

    # add partial missing
    if(partial.missing.prob > 0) {
        for(i in 1:n.chr) {
            if(chr.type[i] != "X") {
                o <- sample(c(TRUE,FALSE),n.mar[i],replace=TRUE,
                            prob=c(partial.missing.prob,1-partial.missing.prob))
                if(any(o)) {
                    o2 <- grep("^QTL[0-9]+",mar.names[[i]])
                    if(length(o2)>0)
                        x <- geno[[i]]$data[,o2]
                    m <- (1:n.mar[i])[o]
                    for(j in m) {
                        if(runif(1) < 0.5)
                            geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==2,j] <- 4
                        else
                            geno[[i]]$data[geno[[i]]$data[,j]==3 | geno[[i]]$data[,j]==2,j] <- 5
                    }
                    if(length(o2)>0)
                        geno[[i]]$data[,o2] <- x
                }
            }

        } # end loop over chromosomes
    } # end simulate partially missing data

    # add missing
    if(missing.prob > 0) {
        for(i in 1:n.chr) {
            o <- grep("^QTL[0-9]+",mar.names[[i]])
            if(length(o)>0)
                x <- geno[[i]]$data[,o]
            geno[[i]]$data[sample(c(TRUE,FALSE),n.mar[i]*n.ind,replace=TRUE,
                                  prob=c(missing.prob,1-missing.prob))] <- NA
            if(length(o)>0)
                geno[[i]]$data[,o] <- x
        }
    }

    pheno <- data.frame(phenotype=pheno, stringsAsFactors=TRUE)

    cross <- list(geno=geno,pheno=pheno)
    class(cross) <- c("f2","cross")

    cross
}

######################################################################
#
# sim.cross.4way
#
######################################################################

sim.cross.4way <-
    function(map,model,n.ind,error.prob,missing.prob,partial.missing.prob,
             keep.errorind,m,p,map.function)
{
    if(map.function=="kosambi") mf <- mf.k
    else if(map.function=="c-f") mf <- mf.cf
    else if(map.function=="morgan") mf <- mf.m
    else mf <- mf.h

    if(!all(sapply(map,is.matrix)))
        stop("Map must be sex-specific.")

    n.chr <- length(map)
    if(is.null(model)) n.qtl <- 0
    else {
        if(!((!is.matrix(model) && length(model) == 5) ||
             (is.matrix(model) && ncol(model) == 5))) {
            stop("Model must be a matrix with 5 columns (chr, pos and effects).")
        }
        if(!is.matrix(model)) model <- rbind(model)
        n.qtl <- nrow(model)
        if(any(model[,1] < 0 | model[,1] > n.chr))
            stop("Chromosome indicators in model matrix out of range.")
        model[,2] <- model[,2]+1e-14 # so QTL not on top of marker
    }

    chr.type <- sapply(map, chrtype)

    # if any QTLs, place qtls on map
    if(n.qtl > 0) {
        for(i in 1:n.qtl) {
            temp <- map[[model[i,1]]]
            temp1 <- temp[1,]
            temp2 <- temp[2,]
            qtlloc <- model[i,2]

            if(qtlloc < min(temp1)) {
                temp1 <- c(qtlloc,temp1)
                temp2 <- min(temp2) - (min(temp1)-qtlloc)/diff(range(temp1))*diff(range(temp2))
                temp1 <- temp1-min(temp1)
                temp2 <- temp2-min(temp2)
                n <- c(paste("QTL",i,sep=""),colnames(temp))
            }
            else if(qtlloc > max(temp1)) {
                temp1 <- c(temp1,qtlloc)
                temp2 <- (qtlloc-max(temp1))/diff(range(temp1))*diff(range(temp2))+max(temp2)
                n <- c(colnames(temp),paste("QTL",i,sep=""))
            }
            else {
                temp1 <- c(temp1,qtlloc)
                o <- order(temp1)
                wh <- (seq(along=temp1))[order(temp1)==length(temp1)]
                temp2 <- c(temp2[1:(wh-1)],NA,temp2[-(1:(wh-1))])
                temp2[wh] <- temp2[wh-1] + (temp1[wh]-temp1[wh-1])/(temp1[wh+1]-temp1[wh-1]) *
                    (temp2[wh+1]-temp2[wh-1])
                temp1 <- sort(temp1)
                n <- c(colnames(temp),paste("QTL",i,sep=""))[o]
            }
            map[[model[i,1]]] <- rbind(temp1,temp2)
            dimnames(map[[model[i,1]]]) <- list(NULL, n)
        }
    }

    geno <- vector("list", n.chr)
    names(geno) <- names(map)
    n.mar <- sapply(map,ncol)
    mar.names <- lapply(map,function(a) colnames(a))

    for(i in 1:n.chr) {

        # simulate sex
        sex <- NULL
        if(chr.type[i]=="X")
            sex <- rep(0,n.ind)

        # simulate genotype data
        thedata <- sim.bcg(n.ind, map[[i]], m, p, map.function)
        dimnames(thedata) <- list(NULL,mar.names[[i]])

        if(chr.type[i] != "X")
            thedata <- thedata + 2*sim.bcg(n.ind, map[[i]][2:1,], m, p, map.function) - 2

        dimnames(thedata) <- list(NULL,mar.names[[i]])

        geno[[i]] <- list(data = thedata, map = map[[i]])

        class(geno[[i]]) <- chr.type[i]
        class(geno[[i]]$map) <- NULL
    } # end loop over chromosomes

    # simulate phenotypes
    pheno <- rnorm(n.ind,0,1)

    if(n.qtl > 0) {
        # find QTL positions
        QTL.chr <- QTL.loc <- NULL
        for(i in 1:n.chr) {
            o <- grep("^QTL[0-9]+",mar.names[[i]])
            if(length(o)>0) {
                QTL.chr <- c(QTL.chr,rep(i,length(o)))
                QTL.loc <- c(QTL.loc,o)
            }
        }

        # incorporate QTL effects
        for(i in 1:n.qtl) {
            QTL.geno <- geno[[QTL.chr[i]]]$data[,QTL.loc[i]]
            pheno[QTL.geno==1] <- pheno[QTL.geno==1] + model[i,3]
            pheno[QTL.geno==2] <- pheno[QTL.geno==2] + model[i,4]
            pheno[QTL.geno==3] <- pheno[QTL.geno==3] + model[i,5]
        }

    } # end simulate phenotype

    n.mar <- sapply(geno, function(a) ncol(a$map))

    # add errors
    if(error.prob > 0) {
        for(i in 1:n.chr) {
            if(chr.type[i] != "X") { # 4-way cross; autosomal
                a <- sample(0:3,n.mar[i]*n.ind,replace=TRUE,
                            prob=c(1-error.prob,rep(error.prob/3,3)))
                if(any(a>0 & geno[[i]]$data==1))
                    geno[[i]]$data[a>0 & geno[[i]]$data==1] <-
                        geno[[i]]$data[a>0 & geno[[i]]$data==1] + a[a>0 & geno[[i]]$data==1]
                if(any(a>0 & geno[[i]]$data==2))
                    geno[[i]]$data[a>0 & geno[[i]]$data==2] <-
                        geno[[i]]$data[a>0 & geno[[i]]$data==2] + c(-1,1,2)[a[a>0 & geno[[i]]$data==2]]
                if(any(a>0 & geno[[i]]$data==3))
                    geno[[i]]$data[a>0 & geno[[i]]$data==3] <-
                        geno[[i]]$data[a>0 & geno[[i]]$data==3] + c(-2,-1,1)[a[a>0 & geno[[i]]$data==3]]
                if(any(a>0 & geno[[i]]$data==4))
                    geno[[i]]$data[a>0 & geno[[i]]$data==4] <-
                        geno[[i]]$data[a>0 & geno[[i]]$data==4] - a[a>0 & geno[[i]]$data==4]
            }
            else {
                a <- sample(0:1,n.mar[i]*n.ind,replace=TRUE,
                            prob=c(1-error.prob,error.prob))
                if(any(a>0 & geno[[i]]$data==1))
                    geno[[i]]$data[a>0 & geno[[i]]$data==1] <-
                        geno[[i]]$data[a>0 & geno[[i]]$data==1] + 1
                if(any(a>0 & geno[[i]]$data==2))
                    geno[[i]]$data[a>0 & geno[[i]]$data==2] <-
                        geno[[i]]$data[a>0 & geno[[i]]$data==2] - 1
                if(any(a>0 & geno[[i]]$data==3))
                    geno[[i]]$data[a>0 & geno[[i]]$data==3] <-
                        geno[[i]]$data[a>0 & geno[[i]]$data==3] + 1
                if(any(a>0 & geno[[i]]$data==4))
                    geno[[i]]$data[a>0 & geno[[i]]$data==4] <-
                        geno[[i]]$data[a>0 & geno[[i]]$data==4] - 1
            }

            if(keep.errorind) {
                errors <- matrix(0,n.ind,n.mar[i])
                errors[a>0] <- 1
                colnames(errors) <- colnames(geno[[i]]$data)
                geno[[i]]$errors <- errors
            }

        } # end loop over chromosomes
    } # end simulate genotyping errors

    # add partial missing
    if(partial.missing.prob > 0) {
        for(i in 1:n.chr) {
            if(chr.type[i] != "X") {
                o <- sample(c(TRUE,FALSE),n.mar[i],replace=TRUE,
                            prob=c(partial.missing.prob,1-partial.missing.prob))

                if(any(o)) {
                    o2 <- grep("^QTL[0-9]+",mar.names[[i]])
                    if(length(o2)>0)
                        x <- geno[[i]]$data[,o2]
                    m <- (1:n.mar[i])[o]
                    for(j in m) {
                        a <- sample(1:4,1)
                        if(a==1) { # AB:AA marker
                            geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==3,j] <- 5
                            geno[[i]]$data[geno[[i]]$data[,j]==2 | geno[[i]]$data[,j]==4,j] <- 6
                        }
                        else if(a==2) { # AA:AB marker
                            geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==2,j] <- 7
                            geno[[i]]$data[geno[[i]]$data[,j]==3 | geno[[i]]$data[,j]==4,j] <- 8
                        }
                        else if(a==3)  # AB:AB marker
                            geno[[i]]$data[geno[[i]]$data[,j]==2 | geno[[i]]$data[,j]==3,j] <- 10
                        else  # AB:BA marker
                            geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==4,j] <- 9
                    }
                    if(length(o2) > 0)
                        geno[[i]]$data[,o2] <- x
                }
            }

        } # end loop over chromosomes
    } # end simulate partially missing data

    # add missing
    if(missing.prob > 0) {
        for(i in 1:n.chr) {
            o <- grep("^QTL[0-9]+",mar.names[[i]])
            if(length(o)>0)
                x <- geno[[i]]$data[,o]
            geno[[i]]$data[sample(c(TRUE,FALSE),n.mar[i]*n.ind,replace=TRUE,
                                  prob=c(missing.prob,1-missing.prob))] <- NA
            if(length(o)>0)
                geno[[i]]$data[,o] <- x
        }
    }

    if(!is.null(sex)) {
        pheno <- cbind(pheno,sex)
        dimnames(pheno) <- list(NULL, c("phenotype", "sex"))
    }
    else {
        pheno <- cbind(pheno)
        dimnames(pheno) <- list(NULL, "phenotype")
    }

    pheno <- as.data.frame(pheno, stringsAsFactors=TRUE)
    cross <- list(geno=geno,pheno=pheno)
    class(cross) <- c("4way","cross")

    cross
}


######################################################################
# sim.bcg
#
# simulate backcross genotype data for a single chromosome;
# output is a matrix of 1's and 0's
######################################################################

sim.bcg <-
    function(n.ind, map, m, p,
             map.function=c("haldane","kosambi","c-f","morgan"))
{
    map.function <- match.arg(map.function)

    if(map.function=="kosambi") mf <- mf.k
    else if(map.function=="c-f") mf <- mf.cf
    else if(map.function=="morgan") mf <- mf.m
    else mf <- mf.h

    if(m < 0 || p < 0 || p > 1)
        stop("Must have m >= 0 and 0 <= p <= 1")

    if(is.matrix(map)) map <- map[1,]
    map <- map-map[1]
    n.mar <- length(map)

    if(m==0 || p==1) { # no interference
        g <- .C("R_sim_bc_ni",
                as.integer(n.mar),
                as.integer(n.ind),
                as.double(mf(diff(map))),
                g=as.integer(rep(0, n.mar*n.ind)),
                PACKAGE="qtl")$g
    }
    else {
        g <- .C("R_sim_bc",
                as.integer(n.mar),
                as.integer(n.ind),
                as.double(map),
                as.integer(m),
                as.double(p),
                g=as.integer(rep(0, n.mar*n.ind)),
                PACKAGE="qtl")$g
    }
    matrix(g, ncol=n.mar)
}



# end of simulate.R