R/util.R

In kbroman/qtlcharts: Interactive Graphics for QTL Experiments

## generic utilities
## Karl W Broman

# grab argument from a list
#
# for example:
#   grabarg(list(method="argmax", map.function="c-f"), "method", "imp")
grabarg <-
function(arguments, argname, default)
    ifelse(argname %in% names(arguments), arguments[[argname]], default)


# return selected phenotype columns as a character vector
getPhename <-
function(cross, pheno.col)
{
    if(is.character(pheno.col)) return(pheno.col)
    names(cross$pheno)[pheno.col]
}

# turn a selection of matrix columns into a numeric vector
selectMatrixColumns <-
function(matrix, cols)
{
    stopifnot(is.matrix(matrix))

    origcols <- cols

    if(is.character(cols)) {
        cols <- match(cols, colnames(matrix))
        if(any(is.na(cols)))
            stop("Unmatched columns: ", paste(origcols[is.na(cols)], collapse=" "))
    }

    (1:ncol(matrix))[cols]
}

# extract phenotypes
extractPheno <-
function(cross, pheno.col)
{
    if(is.character(pheno.col)) {
        pheindex <- qtl::find.pheno(cross, pheno.col)
        if(any(is.na(pheindex)))
            stop("Some phenotypes not found: ",
                 paste(pheno.col[is.na(pheindex)], collapse=" "))
        pheno.col <- pheindex
    }
    if(is.matrix(pheno.col) && nrow(pheno.col) == qtl::nind(cross))
        return(pheno.col) # treat as phenotype matrix
    if(is.numeric(pheno.col) && length(pheno.col) == qtl::nind(cross))
        return(cbind("phenotype"=pheno.col))

    if(is.numeric(pheno.col)) { # look for problem indices
        if(any(pheno.col == 0))
            stop("Cannot have pheno.col == 0")
        if(any(pheno.col < 0) && !all(pheno.col < 0))
            stop("Cannot give a mixture of positive and negative indices")
    }

    # handle negative indices and logical values
    pheno.col <- (1:qtl::nphe(cross))[pheno.col]

    phe <- cross$pheno[,pheno.col,drop=FALSE]
    isnum <- vapply(phe, is.numeric, TRUE)
    if(!all(isnum))
        stop("Some phenotypes not numeric: ",
             paste(qtl::phenames(cross)[pheno.col[!isnum]], collapse=" "))

    as.matrix(phe)
}

# signed LOD scores
#
# If columns==1 and the first column is not "a", we don't change signs
# Otherwise, take the average across these columns to determine the sign
calcSignedLOD <-
function(scanoneOutput, effects, columns=1)
{
    stopifnot(length(effects) == nrow(scanoneOutput))
    stopifnot(all(vapply(effects, nrow, 1) == ncol(scanoneOutput)-2))

    signs <- t(vapply(effects, function(a) {
        if(length(columns)==1 & columns==1 && colnames(a)[1]!="a") return(rep(1, nrow(a)))
        (rowMeans(a[,columns,drop=FALSE], na.rm=TRUE)>=0)*2-1
    }, effects[[1]][,1]))

    scanoneOutput[,-(1:2)] <- signs * scanoneOutput[,-(1:2)]

    scanoneOutput
}

# test if a vector is a set of equally-spaced values
is_equally_spaced <-
function(vec, tol=1e-5)
{
    if(length(vec) < 2) {
        warning("vector has length < 2")
        return(FALSE)
    }

    if(!is.numeric(vec)) {
        warning("vector is not numeric")
        return(FALSE)
    }

    if(any(is.na(vec))) {
        warning("vector contains missing values")
        return(FALSE)
    }

    d <- diff(vec)
    if(!(all(d >= 0) || all(d <= 0))) {
        warning("vector is not monotonic")
        return(FALSE)
    }

    return(stats::sd(d)/abs(stats::median(d)) < tol) # if TRUE, looks equally spaced
}

# check that sets of individual IDs are all the same
# and then return them
# argument n = expected length
get_indID <-
    function(n, ...)
{
    ids <- list(...)

    # get rid of NULLs
    nulls <- vapply(ids, is.null, TRUE)
    if(sum(!nulls)==0)
        return(as.character(1:n))
    ids <- ids[!nulls]

    # there's just one
    if(length(ids)==1)
        ids <- ids[[1]]
    else {
        flag_incompat <- FALSE
        for(i in 2:length(ids)) {
            if(length(ids[[2]]) != length(ids[[1]]) ||
               any(ids[[2]] != ids[[1]]))
                flag_incompat <- TRUE
        }
        if(flag_incompat)
            warning("Incompatibilities among possible individual IDs")
        ids <- ids[[1]]
    }

    if(length(ids) != n)
        warning("Number of ids (", length(ids), ") != expected length (", n, ")")

    ids
}

# ensure that a "group" vector is really the numbers 1, 2, ..., k
# if missing values, put in the last group
group2numeric <-
function(group, preserveNA=FALSE)
{
    if(is.null(group)) return(NULL)

    if(is.factor(group))
        group <- as.numeric(group)
    else if(!is.numeric(group)) # if numeric, leave as it is
        group <-  match(group, sort(unique(group)))

    # NAs -> last group
    if(!preserveNA && any(is.na(group)))
        group[is.na(group)] <- max(group, na.rm=TRUE)+1

    group
}

# determine cross type
crosstype <-
    function(cross)
    {
        type <- class(cross)
        type <- type[type != "cross" & type != "list"]
        if(length(type) > 1) {
            warning("cross has multiple classes")
        }
        type[1]
    }

# determine chromosome type
chrtype <-
    function(object)
    {
        if(inherits(object, "X")) return("X")
        "A"
    }