R/calc.pairprob.R

In qtl: Tools for Analyzing QTL Experiments

######################################################################
#
# calc.pairprob.R
#
# copyright (c) 2001-2019, Karl W Broman
# last modified Dec, 2019
# first written Nov, 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: calc.pairprob
#
######################################################################

######################################################################
#
# calc.pairprob: calculate joint genotype probabilities for all pairs
#                of putative QTLs, conditional on the observed marker
#                data
#
# This is an *internal* function, not to be called by the user.
#
# The input argument cross is assumed to have just one chromosome.
#
######################################################################

calc.pairprob <-
    function(cross, step=0, off.end=0, error.prob=0.0001,
             map.function=c("haldane","kosambi","c-f","morgan"),
             map, assumeCondIndep=FALSE)
{
    # which type of cross is this?
    type <- crosstype(cross)

    if(assumeCondIndep) { # assume conditional independence of QTL given markers
        if(!("prob" %in% names(cross$geno[[1]]))) {
            cross <- calc.genoprob(subset(cross, chr=1), step=step, off.end=off.end,
                                   error.prob=error.prob, map.function=map.function)
        }
        prob <- cross$geno[[1]]$prob
        n.ind <- dim(prob)[1]
        n.pos <- dim(prob)[2]
        n.gen <- dim(prob)[3]

        if(n.pos < 2) return(NULL)

        z <- .C("R_calc_pairprob_condindep",
                as.integer(n.ind),
                as.integer(n.pos),
                as.integer(n.gen),
                as.double(prob),
                pairprob=as.double(rep(0,n.ind*choose(n.pos, 2)*n.gen*n.gen)),
                PACKAGE="qtl")

        pairprob <- array(z$pairprob, dim=c(n.ind,n.pos*(n.pos-1)/2,n.gen,n.gen))

        return(pairprob)
    }

    if(step==0 && off.end > 0) step <- off.end*2

    # map function
    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

    # 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!")
    }
    n.ind <- nind(cross)
    n.chr <- nchr(cross)

    # type of chromosome?
    chrtype <- chrtype(cross$geno[[1]])
    if(chrtype=="X") xchr <- TRUE
    else xchr <- FALSE

    if(type == "f2") {
        one.map <- TRUE
        if(!xchr) { # autosome
            cfunc <- "calc_pairprob_f2"
            n.gen <- 3
            gen.names <- getgenonames("f2", "A", cross.attr=attributes(cross))
        }
        else {                             # X chromsome
            cfunc <- "calc_pairprob_bc"
            n.gen <- 2
            gen.names <- c("g1","g2")
        }
    }
    else if(type == "bc") {
        cfunc <- "calc_pairprob_bc"
        n.gen <- 2
        if(!xchr) # autosome
            gen.names <- getgenonames("bc", "A", cross.attr=attributes(cross))
        else gen.names <- c("g1","g2")
        one.map <- TRUE
    }
    else if(type == "riself" || type=="risib" || type=="dh" || type=="haploid") {
        cfunc <- "calc_pairprob_bc"
        n.gen <- 2
        gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
        one.map <- TRUE
    }
    else if(type == "4way") {
        cfunc <- "calc_pairprob_4way"
        n.gen <- 4
        one.map <- FALSE
        gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
    }
    else if(type == "ri4self" || type=="ri4sib" || type=="ri8self" || type=="ri8sib" || type=="magic16") {
        cfunc <- paste("calc_pairprob_", type, sep="")
        if(type=="magic16") n.gen <- 16
        else n.gen <- as.numeric(substr(type, 3, 3))
        one.map <- TRUE
        gen.names <- LETTERS[1:n.gen]
        if(xchr)
            warning("calc.pairprob not working properly for the X chromosome for 4- or 8-way RIL.")
    }
    else if(type == "bcsft") {
        one.map <- TRUE
        cfunc <- "calc_pairprob_bcsft"
        cross.scheme <- attr(cross, "scheme") ## c(s,t) for BC(s)F(t)
        if(!xchr) { # autosome
            gen.names <- getgenonames("bcsft", "A", cross.attr=attributes(cross))
            n.gen <- 2 + (cross.scheme[2] > 0)
        }
        else { ## X chromsome
            cross.scheme[1] <- cross.scheme[1] + cross.scheme[2] - (cross.scheme[1] == 0)
            cross.scheme[2] <- 0
            gen.names <- c("g1","g2")
            n.gen <- 2
        }
    }
    else
        stop("calc.pairprob not available for cross type ", type, ".")

    # genotype data
    gen <- cross$geno[[1]]$data
    gen[is.na(gen)] <- 0

    # get recombination fractions
    if(one.map) {
        #    map <- create.map(cross$geno[[1]]$map,step,off.end)
        rf <- mf(diff(map))
        if(type=="risib" || type=="riself")
            rf <- adjust.rf.ri(rf,substr(type,3,nchar(type)),chrtype(cross$geno[[1]]))
        rf[rf < 1e-14] <- 1e-14

        # new genotype matrix with pseudomarkers filled in
        newgen <- matrix(ncol=length(map),nrow=nrow(gen))
        colnames(newgen) <- names(map)
        newgen[,colnames(gen)] <- gen
        newgen[is.na(newgen)] <- 0
        n.pos <- ncol(newgen)
        marnames <- names(map)
    }
    else {
        #    map <- create.map(cross$geno[[1]]$map,step,off.end)
        rf <- mf(diff(map[1,]))
        rf[rf < 1e-14] <- 1e-14
        rf2 <- mf(diff(map[2,]))
        rf2[rf2 < 1e-14] <- 1e-14

        # new genotype matrix with pseudomarkers filled in
        newgen <- matrix(ncol=ncol(map),nrow=nrow(gen))
        colnames(newgen) <- colnames(map)
        newgen[,colnames(gen)] <- gen
        newgen[is.na(newgen)] <- 0
        n.pos <- ncol(newgen)
        marnames <- colnames(map)
    }

    if(n.pos < 2) return(NULL)

    # below: at least two positions
    # call the C function
    if(one.map) {
        ## Hide cross scheme in genoprob to pass to routine. BY
        temp <- as.double(rep(0,n.gen*n.ind*n.pos))
        if(type == "bcsft")
            temp[1:2] <- cross.scheme

        z <- .C(cfunc,
                as.integer(n.ind),         # number of individuals
                as.integer(n.pos),         # number of markers
                as.integer(newgen),        # genotype data
                as.double(rf),             # recombination fractions
                as.double(error.prob),     #
                as.double(temp),
                pairprob=as.double(rep(0,n.ind*n.pos*(n.pos-1)/2*n.gen^2)),
                PACKAGE="qtl")
    }
    else {
        z <- .C(cfunc,
                as.integer(n.ind),         # number of individuals
                as.integer(n.pos),         # number of markers
                as.integer(newgen),        # genotype data
                as.double(rf),             # recombination fractions
                as.double(rf2),            # recombination fractions
                as.double(error.prob),     #
                as.double(rep(0,n.gen*n.ind*n.pos)),
                pairprob=as.double(rep(0,n.ind*n.pos*(n.pos-1)/2*n.gen^2)),
                PACKAGE="qtl")
    }

    pairprob <- array(z$pairprob, dim=c(n.ind,n.pos*(n.pos-1)/2,n.gen,n.gen))

    # 4- and 8-way RIL: reorganize the results
    if(type=="ri4self" || type=="ri4sib" || type=="ri8self" || type=="ri8sib" || type=="bgmagic16")
        pairprob <- reorgRIpairprob(cross, pairprob)

    pairprob
}

# end of calc.pairprob.R