R/iplotScantwo.R

In kbroman/qtlcharts: Interactive Graphics for QTL Experiments

## iplotScantwo
## Karl W Broman

#' Interactive plot of 2d genome scan
#'
#' Creates an interactive plot of the results of
#' [qtl::scantwo()], for a two-dimensional, two-QTL genome
#' scan.
#'
#' @param scantwoOutput Output of [qtl::scantwo()]
#' @param cross (Optional) Object of class `"cross"`, see
#'   [qtl::read.cross()].
#' @param lodcolumn Numeric value indicating LOD score column to plot.
#' @param pheno.col (Optional) Phenotype column in cross object.
#' @param chr (Optional) Optional vector indicating the chromosomes
#'   for which LOD scores should be calculated. This should be a vector
#'   of character strings referring to chromosomes by name; numeric
#'   values are converted to strings. Refer to chromosomes with a
#'   preceding - to have all chromosomes but those considered. A logical
#'   (TRUE/FALSE) vector may also be used.
#' @param chartOpts A list of options for configuring the chart.  Each
#'   element must be named using the corresponding option.
#' @param digits Round data to this number of significant digits
#'     before passing to the chart function. (Use NULL to not round.)
#'
#' @return An object of class `htmlwidget` that will
#' intelligently print itself into HTML in a variety of contexts
#' including the R console, within R Markdown documents, and within
#' Shiny output bindings.
#'
#' @details The estimated QTL effects, and the genotypes in the
#' phenotype x genotype plot, in the right-hand panels, are derived
#' following a single imputation to fill in missing data, and so are a
#' bit crude.
#'
#' Note that the usual `height` and `width` options in
#' `chartOpts` are ignored here. Instead, you may provide
#' `pixelPerCell` (number of pixels per cell in the heat map),
#' `chrGap` (gaps between chr in heat map), `wright` (width
#' in pixels of right panels), and `hbot` (height in pixels of
#' each of the lower panels)
#'
#' @keywords hplot
#' @seealso [iplotScanone()]
#'
#' @examples
#' library(qtl)
#' data(fake.f2)
#' \dontshow{fake.f2 <- fake.f2[c(1, 13, "X"),]}
#' fake.f2 <- calc.genoprob(fake.f2, step=5)
#' out <- scantwo(fake.f2, method="hk", verbose=FALSE)
#' \donttest{
#' iplotScantwo(out, fake.f2)}
#'
#' @export
iplotScantwo <-
function(scantwoOutput, cross=NULL, lodcolumn=1, pheno.col=1, chr=NULL,
         chartOpts=NULL, digits=5)
{
    if(!inherits(scantwoOutput, "scantwo"))
        stop('"scantwoOutput" should have class "scantwo".')

    if(!is.null(chr)) {
        scantwoOutput <- subset(scantwoOutput, chr=chr)
        if(!is.null(cross)) cross <- subset(cross, chr=chr)
    }

    if(length(lodcolumn) > 1) {
        lodcolumn <- lodcolumn[1]
        warning("lodcolumn should have length 1; using first value")
    }
    if(length(dim(scantwoOutput)) < 3) {
        if(lodcolumn != 1) {
            warning("scantwoOutput contains just one set of lod scores; using lodcolumn=1")
            lodcolumn <- 1
        }
    }
    else {
        d <- dim(scantwoOutput)[3]
        if(lodcolumn < 1 || lodcolumn > 3)
            stop("lodcolumn should be between 1 and ", d)
        scantwoOutput$lod <- scantwoOutput$lod[,,lodcolumn]
    }

    if(length(pheno.col) > 1) {
        pheno.col <- pheno.col[1]
        warning("pheno.col should have length 1; using first value")
    }
    if(!is.null(cross))
        pheno <- qtl::pull.pheno(cross, pheno.col)
    else cross <- pheno <- NULL

    scantwo_list <- data4iplotScantwo(scantwoOutput)

    phenogeno_list <- cross4iplotScantwo(scantwoOutput, cross, pheno)

    defaultAspect <- 1 # width/height
    browsersize <- getPlotSize(defaultAspect)

    x <- list(scantwo_data=scantwo_list,
              phenogeno_data=phenogeno_list,
              chartOpts=chartOpts)
    if(!is.null(digits))
        attr(x, "TOJSON_ARGS") <- list(digits=digits)

    htmlwidgets::createWidget("iplotScantwo", x,
                              width=chartOpts$width,
                              height=chartOpts$height,
                              sizingPolicy=htmlwidgets::sizingPolicy(
                                  browser.defaultWidth=browsersize$width,
                                  browser.defaultHeight=browsersize$height,
                                  knitr.defaultWidth=1000,
                                  knitr.defaultHeight=1000/defaultAspect
                              ),
                              package="qtlcharts")
}

#' @rdname qtlcharts-shiny
#' @export
iplotScantwo_output <- function(outputId, width="100%", height="1000") {
    htmlwidgets::shinyWidgetOutput(outputId, "iplotScantwo", width, height, package="qtlcharts")
}

#' @rdname qtlcharts-shiny
#' @export
iplotScantwo_render <- function(expr, env=parent.frame(), quoted=FALSE) {
    if(!quoted) { expr <- substitute(expr) } # force quoted
    htmlwidgets::shinyRenderWidget(expr, iplotScantwo_output, env, quoted=TRUE)
}


# convert scantwo output to JSON format
data4iplotScantwo <-
    function(scantwoOutput)
{
    lod <- scantwoOutput$lod
    map <- scantwoOutput$map
    lod <- lod[map$eq.spacing==1, map$eq.spacing==1,drop=FALSE]
    map <- map[map$eq.spacing==1,,drop=FALSE]

    lodv1 <- get_lodv1(lod, map, scantwoOutput$scanoneX)

    chr <- as.character(map$chr)
    chrnam <- unique(chr)
    n.mar <- tapply(chr, chr, length)[chrnam]
    names(n.mar) <- NULL

    labels <- revisePmarNames(rownames(map))

    dimnames(lod) <- NULL
    dimnames(lodv1) <- NULL
    list(lod=lod, lodv1=lodv1,
         nmar=n.mar, chrname=chrnam,
         marker=labels, chr=chr, pos=map$pos)
}


# convert pseudomarker names, e.g. "c5.loc25" -> "5@25"
revisePmarNames <-
    function(labels)
{
    wh.pmar <- grep("^c.+\\.loc-*[0-9]+", labels)
    pmar <- labels[wh.pmar]
    pmar.spl <- strsplit(pmar, "\\.loc")
    pmar_chr <- vapply(pmar.spl, "[", "", 1)
    pmar_chr <- substr(pmar_chr, 2, nchar(pmar_chr))
    pmar_pos <- vapply(pmar.spl, "[", "", 2)
    labels[wh.pmar] <- paste0(pmar_chr, "@", pmar_pos)
    labels
 }


# get fv1/av1 LOD scores from the full/add LOD scores
get_lodv1 <-
    function(lod, map, scanoneX=NULL)
{
    thechr <- map$chr
    uchr <- unique(thechr)
    thechr <- factor(as.character(thechr), levels=as.character(uchr))
    uchr <- factor(as.character(uchr), levels=levels(thechr))
    xchr <- tapply(map$xchr, thechr, function(a) a[1])

    # maximum 1-d LOD score on each chromosome
    maxo <- tapply(diag(lod), thechr, max, na.rm=TRUE)
    if(any(xchr) && !is.null(scanoneX)) {
        maxox <- tapply(scanoneX, thechr, max, na.rm=TRUE)
        maxo[xchr] <- maxox[xchr]
    }

    # subtract max(i,j) from each chr pair (i,j)
    n.chr <- length(uchr)
    for(i in 1:n.chr) {
        pi <- which(thechr==uchr[i])
        for(j in 1:n.chr) {
            pj <- which(thechr==uchr[j])
            lod[pj,pi] <- lod[pj,pi] - max(maxo[c(i,j)])
        }
    }
    # if negative, replace with 0
    lod[is.na(lod) | lod < 0] <- 0

    lod
}


# convert genotype/phenotype information to JSON format
cross4iplotScantwo <-
    function(scantwoOutput, cross=NULL, pheno=NULL)
{
    # if no cross or phenotype, just return null
    if(is.null(cross) || is.null(pheno))
        return(NULL)

    # pull out locations of LOD calculations, on grid
    map <- scantwoOutput$map
    map <- map[map$eq.spacing==1,,drop=FALSE]

    # local function to get pseudomarker names
    getPmarNames <-
        function(cross)
        {
            nam <- lapply(cross$geno, function(a) colnames(a$draws))
            chrnam <- names(cross$geno)
            for(i in seq(along=nam)) {
                g <- grep("^loc", nam[[i]])
                if(length(g) > 0)
                    nam[[i]][g] <- paste0("c", chrnam[i], ".", nam[[i]][g])
            }
            nam
        }

    # attempt to get imputations at pseudomarker locations in scantwo object
    needImp <- TRUE
    if("draws" %in% names(cross$geno[[1]])) { # contains imputations; attempt to use these
        nam <- getPmarNames(cross)
        if(all(rownames(scantwoOutput) %in% unlist(nam))) { # same pseudomarkers as in scantwo
            needImp <- FALSE
        }
    }
    if(needImp) {
        if("prob" %in% names(cross$geno[[1]])) { # contains genotype probabilities
            # do imputation with positions in calc.genoprob
            cross <- qtl::sim.geno(cross,
                                   error.prob=attr(cross$geno[[1]]$prob, "error.prob"),
                                   step=attr(cross$geno[[1]]$prob, "step"),
                                   off.end=attr(cross$geno[[1]]$prob, "off.end"),
                                   map.function=attr(cross$geno[[1]]$prob, "map.function"),
                                   stepwidth=attr(cross$geno[[1]]$prob, "stepwidth"),
                                   n.draws=1)

            nam <- getPmarNames(cross)
            if(all(rownames(scantwoOutput) %in% unlist(nam))) {
                needImp <- FALSE
            }
        }

        if(needImp) # give up
            stop('cross object needs imputed genotypes at pseudomarkers\n',
                 'Run sim.geno with step/stepwidth as used for scantwo')
    }

    # X chr imputations: 1/2 -> AA/AB/BB/AY/BY
    cross_type <- crosstype(cross)
    chr_type <- sapply(cross$geno, chrtype)
    sexpgm <- qtl::getsex(cross)
    cross.attr <- attributes(cross)
    if(cross_type %in% c("f2", "bc", "bcsft") && any(chr_type=="X")) {
        for(i in which(chr_type=="X")) {
            cross$geno[[i]]$draws <- qtl::reviseXdata(cross_type, "standard", sexpgm,
                                                      draws=cross$geno[[i]]$draws,
                                                      cross.attr=cross.attr)
        }
    }

    # pull out imputed genotypes as a list
    geno <- vector("list", nrow(map))
    names(geno) <- rownames(map)
    for(i in seq(along=cross$geno)) {
        m <- which(nam[[i]] %in% rownames(map))
        for(j in m)
            geno[[nam[[i]][j]]] <- cross$geno[[i]]$draws[,j,1]
    }
    names(geno) <- revisePmarNames(rownames(map))

    # names of the possible genotypes
    genonames <- vector("list", length(cross$geno))
    names(genonames) <- names(cross$geno)
    for(i in seq(along=genonames))
        genonames[[i]] <- qtl::getgenonames(cross_type, class(cross$geno[[i]]),
                                            "standard", sexpgm, cross.attr)

    X_geno_by_sex <- NULL
    chr_type <- vapply(cross$geno, chrtype, "")
    if(any(chr_type=="X")) {
        f_sexpgm <- m_sexpgm <- sexpgm
        f_sexpgm$sex[f_sexpgm$sex==1] <- 0
        m_sexpgm$sex[m_sexpgm$sex==0] <- 1
        X_geno_by_sex <- list(
            qtl::getgenonames(cross_type, class(cross$geno[[i]]), "standard", f_sexpgm, cross.attr),
            qtl::getgenonames(cross_type, class(cross$geno[[i]]), "standard", m_sexpgm, cross.attr)
        )
    }


    # chr for each marker
    chr <- as.character(map$chr)
    names(chr) <- names(geno) # the revised pseudomarker names

    # individual IDs
    indID <- qtl::getid(cross)
    if(is.null(indID)) indID <- 1:qtl::nind(cross)
    indID <- as.character(indID)

    list(geno=geno, chr=as.list(chr), genonames=genonames, X_geno_by_sex=X_geno_by_sex,
         pheno=pheno, indID=indID, chrtype=as.list(chr_type))
}