R/scan1perm_gen.R

In qtl2: Quantitative Trait Locus Mapping in Experimental Crosses

# scan1perm with general function but no covariates, kinship, weights
scan1perm_gen_simple <-
    function(genoprobs, pheno, n_perm=1, perm_strata=NULL,
             cores=1, scan_func, ind2keep, ..., max_batch=NULL)
{
    dotargs <- list(...)
    tol <- grab_dots(dotargs, "tol", 1e-12)
    quiet <- grab_dots(dotargs, "quiet", TRUE)

    # set up parallel analysis
    cores <- setup_cluster(cores)
    if(!quiet && n_cores(cores)>1) {
        message(" - Using ", n_cores(cores), " cores")
        quiet <- TRUE # make the rest quiet
    }

    # max batch size
    if(is.null(max_batch))
        max_batch <- min(1000, ceiling(n_perm*length(genoprobs)*ncol(pheno)/n_cores(cores)))

    # generate permutations
    perms <- gen_strat_perm(n_perm, ind2keep, perm_strata)

    # batch permutations
    phe_batches <- batch_vec( rep(seq_len(ncol(pheno)), n_perm), max_batch)
    perm_batches <- batch_vec( rep(seq_len(n_perm), each=ncol(pheno)), max_batch)

    # drop cols in genotype probs that are all 0 (just looking at the X chromosome)
    genoprob_Xcol2drop <- genoprobs_col2drop(genoprobs)
    is_x_chr <- attr(genoprobs, "is_x_chr")
    if(is.null(is_x_chr)) is_x_chr <- rep(FALSE, length(genoprobs))

    # batches for analysis, to allow parallel analysis
    run_batches <- data.frame(chr=rep(seq_len(length(genoprobs)), length(phe_batches)),
                              phe_batch=rep(seq_along(phe_batches), each=length(genoprobs)))

    run_indexes <- seq_len(length(genoprobs)*length(phe_batches))

    # subset the phenotypes
    pheno <- pheno[ind2keep,,drop=FALSE]
    nind <- nrow(pheno)

    # the function that does the work
    by_group_func <- function(i) {
        # deal with batch information, including individuals to drop due to missing phenotypes
        chr <- run_batches$chr[i]
        chrnam <- names(genoprobs)[chr]
        phebatch <- phe_batches[[run_batches$phe_batch[i]]]
        permbatch <- perm_batches[[run_batches$phe_batch[i]]]

        pr <- genoprobs[[chr]][ind2keep,,,drop=FALSE]
        pr <- setNames(list(pr), chrnam)
        attr(pr, "is_x_chr") <- attr(genoprobs, "is_x_chr")[chr]
        attr(pr, "crosstype") <- attr(genoprobs, "crosstype")
        attr(pr, "alleles") <- attr(genoprobs, "alleles")
        attr(pr, "alleleprobs") <- attr(genoprobs, "alleleprobs")
        class(pr) <- class(genoprobs)

        ph <- pheno[,phebatch,drop=FALSE]
        for(col in seq_len(ncol(ph))) # permute columns
            ph[,col] <- ph[perms[,permbatch[col]] , col]

        lod <- scan_func(genoprobs=pr, pheno=ph,
                         addcovar=NULL, Xcovar=NULL,
                         intcovar=NULL, kinship=NULL,
                         cores=1, ...)
        lod <- fix_output_if_snps(lod)

        # return column maxima
        apply(lod, 2, max, na.rm=TRUE)
    }

    # calculations in parallel
    list_result <- cluster_lapply(cores, run_indexes, by_group_func)

    # check for problems (if clusters run out of memory, they'll return NULL)
    result_is_null <- vapply(list_result, is.null, TRUE)
    if(any(result_is_null))
        stop("cluster problem: returned ", sum(result_is_null), " NULLs.")

    # reorganize results
    result <- array(dim=c(length(genoprobs), n_perm, ncol(pheno)))

    for(i in run_indexes) {
        chr <- run_batches$chr[i]
        phebatch <- phe_batches[[run_batches$phe_batch[i]]]
        permbatch <- perm_batches[[run_batches$phe_batch[i]]]

        for(j in seq_along(phebatch))
            result[chr,permbatch[j], phebatch[j]] <- list_result[[i]][j]

    }

    result <- apply(result, c(2,3), max)
    colnames(result) <- colnames(pheno)

    class(result) <- c("scan1perm", "matrix")
    result

}

scan1perm_gen <-
    function(genoprobs, pheno, kinship=NULL, addcovar=NULL, Xcovar=NULL,
             intcovar=NULL, weights=NULL,
             n_perm=1, perm_strata=NULL,
             cores=1, scan_func, ind2keep, ..., max_batch=NULL)
{
    dotargs <- list(...)
    tol <- grab_dots(dotargs, "tol", 1e-12)
    quiet <- grab_dots(dotargs, "quiet", TRUE)

    # set up parallel analysis
    cores <- setup_cluster(cores)
    if(!quiet && n_cores(cores)>1) {
        message(" - Using ", n_cores(cores), " cores")
        quiet <- TRUE # make the rest quiet
    }

    # max batch size
    if(is.null(max_batch))
        max_batch <- min(1000, ceiling(n_perm*length(genoprobs)*ncol(pheno)/n_cores(cores)))

    # generate permutations
    perms <- gen_strat_perm(n_perm, ind2keep, perm_strata)

    # batch permutations
    phe_batches <- batch_cols(pheno[ind2keep,,drop=FALSE], max_batch)

    # batches for analysis, to allow parallel analysis
    run_batches <- data.frame(chr=rep(seq_len(length(genoprobs)), length(phe_batches)*n_perm),
                              phe_batch=rep(seq_along(phe_batches), each=length(genoprobs)*n_perm),
                              perm_batch=rep(rep(seq_len(n_perm), each=length(genoprobs), length(phe_batches))))

    run_indexes <- seq_len(length(genoprobs)*length(phe_batches)*n_perm)

    # the function that does the work
    by_group_func <- function(i) {
        # deal with batch information, including individuals to drop due to missing phenotypes
        chr <- run_batches$chr[i]
        chrnam <- names(genoprobs)[chr]
        phebatch <- phe_batches[[run_batches$phe_batch[i]]]
        permbatch <- run_batches$perm_batch[i]
        phecol <- phebatch$cols
        omit <- phebatch$omit
        these2keep <- ind2keep # individuals 2 keep for this batch
        if(length(omit) > 0) these2keep <- ind2keep[-omit]
        if(length(these2keep)<=2) return(NULL) # not enough individuals

        # apply permutation to the probs
        #     (I initially thought we'd need to apply these just to the rownames, but
        #     actually we've already aligned probs and pheno and so forth
        #     via ind2keep/these2keep, so now we need to just permute the rows)
        #
        # (also need some contortions here to keep the sizes the same)
        pr <- genoprobs[[chr]][ind2keep,,,drop=FALSE]
        pr <- pr[perms[,permbatch],,,drop=FALSE]
        rownames(pr) <- ind2keep

        # make the probabilities back into a genoprobs object
        pr <- setNames(list(pr), chrnam)
        attr(pr, "is_x_chr") <- attr(genoprobs, "is_x_chr")[chr]
        attr(pr, "crosstype") <- attr(genoprobs, "crosstype")
        attr(pr, "alleles") <- attr(genoprobs, "alleles")
        attr(pr, "alleleprobs") <- attr(genoprobs, "alleleprobs")
        class(pr) <- class(genoprobs)

        # subset kinship and pheno
        k <- kinship; if(is_kinship_list(k)) k <- k[[chr]]
        ph <- pheno[, phecol, drop=FALSE]

        # scan1 function taking clean data (with no missing values)
        lod <- scan_func(genoprobs=pr, pheno=ph, kinship=k, addcovar=addcovar,
                         intcovar=intcovar, Xcovar=Xcovar, weights=weights, ...)
        lod <- fix_output_if_snps(lod)

        # return column maxima
        apply(lod, 2, max, na.rm=TRUE)
    }

    # calculations in parallel
    list_result <- cluster_lapply(cores, run_indexes, by_group_func)

    # check for problems (if clusters run out of memory, they'll return NULL)
    result_is_null <- vapply(list_result, is.null, TRUE)
    if(any(result_is_null))
        stop("cluster problem: returned ", sum(result_is_null), " NULLs.")

    # reorganize results
    result <- array(dim=c(length(genoprobs), n_perm, ncol(pheno)))

    for(i in run_indexes) {
        chr <- run_batches$chr[i]
        phebatch <- phe_batches[[run_batches$phe_batch[i]]]
        permbatch <- run_batches$perm_batch[i]

        result[chr, permbatch, phebatch$cols] <- list_result[[i]]
    }

    result <- apply(result, c(2,3), max)
    colnames(result) <- colnames(pheno)

    class(result) <- c("scan1perm", "matrix")
    result
}



# generate (potentially) stratified permutations
#   (actual work is done in c++, see src/random.cpp)
gen_strat_perm <-
    function(n_perm, ind2keep, perm_strata=NULL)
{
    if(!is.null(perm_strata)) {
        perm_strata <- perm_strata[ind2keep]
        u <- unique(perm_strata)

        if(length(u) > 1) {

            strat_numeric <- match(perm_strata, u)-1

            return(permute_nvector_stratified(n_perm,
                                              seq_along(ind2keep),
                                              strat_numeric,
                                              length(u)))
        }
    }

    permute_nvector(n_perm, seq_along(ind2keep))
}

# detect if scan1snps out (list with "lod" and "snpinfo")
fix_output_if_snps <-
    function(scan_output)
{
    if(is.list(scan_output) && length(scan_output)==2
       && all(names(scan_output) == c("lod", "snpinfo")))
        return(scan_output$lod)

    scan_output
}