R/bcsft.R

In qtl: Tools for Analyzing QTL Experiments

convert2bcsft <- function(cross, BC.gen = 0, F.gen = 0, estimate.map = TRUE,
                          error.prob=0.0001, map.function=c("haldane","kosambi","c-f","morgan"),
                          verbose=FALSE)
{
    cross.class <- crosstype(cross)
    if(cross.class=="bcsft") cross.class <- "f2"

    if((cross.class %in% c("bc","f2"))) {
        class(cross) <- c("bcsft", "cross")
        ## If BC.gen = 0 and F.gen = 0, then set to BC1F0 (bc) or BC0F2 (f2).
        if(cross.class == "bc" & F.gen > 0) {
            stop("input cross has only 2 genotypes--cannot have F.gen > 0")
            if(BC.gen == 0)
                BC.gen <- 1
        }
        if(cross.class == "f2") {
            if(F.gen == 0) {
                if(BC.gen == 0)
                    F.gen <- 2
                else
                    stop("input cross has 3 genotypes--cannot have F.gen = 0")
            }
        }
        if(BC.gen < 0 | F.gen < 0)
            stop("BC.gen and F.gen cannot be negative")

        attr(cross, "scheme") <- c(BC.gen, F.gen)
        cross
    }
    else stop("cross object has to be of class bc or f2 to be converted to bcsft")

    # re-estimate map?
    if(estimate.map) {
        cat(" --Estimating genetic map\n")
        newmap <- est.map(cross, error.prob=error.prob, map.function=map.function, verbose=verbose)
        cross <- replace.map(cross, newmap)
    }

    cross
}

read.cross.bcsft <- function(..., BC.gen = 0, F.gen = 0, cross = NULL, force.bcsft = FALSE,
                             estimate.map=TRUE)
{
    ## Must specify s = BC.gen and t = F.gen.
    ## Later: Could import in clever way from qtlcart? See qtlcart_io.R and their software.

    ## Make sure we only estimate map once!
    if(is.null(cross)) # Estimate map at end of this routine (called read.cross.bcsft directly).
        cross <- read.cross(..., estimate.map = FALSE)
    else # Estimate map in parent read.cross() call (read.cross.bcsft is pass-through from read.cross).
        estimate.map <- FALSE

    force.bcsft <- force.bcsft | (BC.gen > 0 | F.gen > 0)
    if((crosstype(cross) %in% c("bc","f2","bcsft")) && force.bcsft) {

        # deal with ... args
        dots <- list(...)

        if("verbose" %in% names(dots)) verbose <- dots$verbose
        else verbose <- TRUE

        if("error.prob" %in% names(dots)) error.prob <- dots$error.prob
        else error.prob <- 0.0001

        if("map.function" %in% names(dots)) map.function <- dots$map.function
        else map.function <- "haldane"

        cross <- convert2bcsft(cross, BC.gen, F.gen, estimate.map = estimate.map,
                               error.prob=error.prob, map.function=map.function,
                               verbose=verbose)
    }

    cross
}

######################################################################

sim.cross.bcsft <- function(map,model,n.ind,error.prob,missing.prob,
                            partial.missing.prob,keep.errorind,
                            m,p,map.function,
                            cross.scheme)
{
    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.")

    ## cross.scheme = c(s,t) for bcsft.
    if(missing(cross.scheme))
        stop("must specify cross.scheme for bcsft")
    if(length(cross.scheme) != 2)
        stop("cross.scheme for bcsft must have 2 values")
    cross.scheme <- round(cross.scheme)
    if(min(cross.scheme) < 0)
        stop("cross.scheme for bcsft must have 2 non-negative integers")
    n.eff <- 3 + (cross.scheme[2] > 0)

    ## 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) == n.eff) ||
             (is.matrix(model) && ncol(model) == n.eff))) {
            stop(paste("Model must be a matrix with ", n.eff, " columns (chr, pos and effect",
                       ifelse(n.eff == 4, "s", ""), ").", sep = ""))
        }
        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)
    BC.gen <- cross.scheme[1]
    F.gen <- cross.scheme[2] - (BC.gen == 0)

    for(i in 1:n.chr) {

        # simulate genotype data
        bcallele1 <- sim.bcg(n.ind, map[[i]], m, p, map.function) - 1
        ## BCs: multiply independent instances of meiosis together.
        if(BC.gen > 0) {
            if(BC.gen > 1) for(j in seq(2, BC.gen))
                bcallele1 <- bcallele1 * (sim.bcg(n.ind, map[[i]], m, p, map.function) - 1)
        }

        if(F.gen == 0) ## BCs only.
            thedata <- bcallele1 + 1
        else {
            if(chr.type[i] == "X") {
                if(F.gen > 1) for(j in seq(F.gen))
                    bcallele1 <- bcallele1 * (sim.bcg(n.ind, map[[i]], m, p, map.function) - 1)

                thedata <- bcallele1 + 1
            }
            else { ## chr.type[i] != "X"
                if(BC.gen > 0) { ## Two unique alleles from BC(s).
                    bcallele2 <- bcallele1 * (sim.bcg(n.ind, map[[i]], m, p, map.function) - 1)
                    bcallele1 <- bcallele1 * (sim.bcg(n.ind, map[[i]], m, p, map.function) - 1)
                }
                else ## Starting from F(1) with two unique alleles.
                    bcallele2 <- sim.bcg(n.ind, map[[i]], m, p, map.function) - 1

                if(F.gen > 1) for(j in seq(2, F.gen)) {
                    ## need two instances.
                    allelemask1 <- sim.bcg(n.ind, map[[i]], m, p, map.function) - 1
                    allelemask2 <- sim.bcg(n.ind, map[[i]], m, p, map.function) - 1
                    bcallele1 <- bcallele1 * allelemask1 + bcallele2 * (1 - allelemask1)
                    bcallele2 <- bcallele2 * allelemask2 + bcallele1 * (1 - allelemask2)
                }
                thedata <- bcallele1 + bcallele2 + 1
            }
        }

        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,n.eff]
            if(n.eff == 4) {
                pheno[QTL.geno==1] <- pheno[QTL.geno==1] - model[i,3]
                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)

    cross <- list(geno=geno,pheno=pheno)
    class(cross) <- c("bcsft","cross")
    attr(cross, "scheme") <- cross.scheme

    cross
}

## End bcsft.R