######################################################################
#
# plot.scantwo.R
#
# copyright (c) 2001-2019, Karl W Broman, Hao Wu and Brian Yandell
# 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
#
# Hao Wu (The Jackson Lab) wrote the initial code
#
# Part of the R/qtl package
# Contains: plot.scantwo
#
######################################################################
plot.scantwo <-
function(x, chr, incl.markers = FALSE, zlim, lodcolumn=1,
lower = c("full", "add", "cond-int", "cond-add", "int"),
upper = c("int", "cond-add", "cond-int", "add", "full"),
nodiag = TRUE,
contours = FALSE, main, zscale = TRUE, point.at.max=FALSE,
col.scheme = c("viridis", "redblue","cm","gray","heat","terrain","topo"),
gamma=0.6, allow.neg=FALSE, alternate.chrid=FALSE, ...)
{
if(!inherits(x, "scantwo"))
stop("Input should have class \"scantwo\".")
col.scheme <- match.arg(col.scheme)
if(length(dim(x$lod)) > 2) { # results from multiple phenotypes
if(length(lodcolumn) > 1) {
warning("Argument lodcolumn should be of length 1.")
lodcolumn <- lodcolumn[1]
}
if(lodcolumn < 0 || lodcolumn > dim(x$lod)[3])
stop("Argument lodcolumn misspecified.")
x$lod <- x$lod[,,lodcolumn]
}
if(!missing(chr)) x <- subset(x, chr=chr)
if(nrow(x$lod)==0) {
warning("Empty scantwo object.")
return(invisible(NULL))
}
chr <- as.character(unique(x$map[,1]))
addpair <- attr(x, "addpair")
if(!is.null(addpair) && addpair) {
lower <- "full"
upper <- "add"
if(missing(zlim)) {
if(allow.neg)
zlim <- rep(max(abs(x$lod), na.rm=TRUE), 2)
else
zlim <- rep(max(x$lod, na.rm=TRUE), 2)
}
addpair <- TRUE
}
else addpair <- FALSE
if(length(lower)==1 && lower == "fv1") lower <- "cond-int"
if(length(lower)==1 && lower == "av1") lower <- "cond-add"
if(length(upper)==1 && upper == "fv1") upper <- "cond-int"
if(length(upper)==1 && upper == "av1") upper <- "cond-add"
lower <- match.arg(lower)
upper <- match.arg(upper)
if(!any(class(x) == "scantwo"))
stop("Input variable is not an object of class scantwo!")
lod <- x$lod
map <- x$map
# backward compatibility for previous version of R/qtl
if(!("scanoneX" %in% names(x))) {
warning("It would be best to re-run scantwo() with the R/qtl version 0.98 or later.\n")
scanoneX <- NULL
}
else scanoneX <- x$scanoneX
# if incl.markers is FALSE, drop positions
# for which third column of map is 0
if(!incl.markers && any(map[, 3] == 0)) {
o <- (map[, 3] == 1)
lod <- lod[o, o]
map <- map[o, ]
if(!is.null(scanoneX)) scanoneX <- scanoneX[o]
}
if(all(diag(lod) < 1e-14) && (lower == "cond-int" || lower=="cond-add") )
stop("Need to run scantwo with run.scanone=TRUE.")
oldlod <- lod
lo <- lower.tri(lod)
up <- upper.tri(lod)
# grab the interaction LOD scores
if(upper=="int")
lod[up] <- t(oldlod)[up] - oldlod[up]
if(lower=="int")
lod[lo] <- oldlod[lo] - t(oldlod)[lo]
if(lower=="add")
lod[lo] <- t(oldlod)[lo]
if(upper=="full")
lod[up] <- t(oldlod)[up]
# get conditional LOD scores
if(lower=="cond-int" || lower=="cond-add") {
if(lower=="cond-add")
lod[lo] <- t(oldlod)[lo]
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])
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]
}
n.chr <- length(chr)
for(i in 1:n.chr) {
pi <- which(thechr==uchr[i])
for(j in i:n.chr) {
pj <- which(thechr==uchr[j])
temp <- lod[pj,pi] - max(maxo[c(i,j)])
temp[!is.na(temp) & temp<0] <- 0
if(i==j) lod[pj,pi][lower.tri(temp)] <- temp[lower.tri(temp)]
else lod[pj,pi] <- temp
}
}
}
if(upper=="cond-int" || upper=="cond-add") {
if(upper=="cond-int")
lod[up] <- t(oldlod)[up]
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])
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]
}
n.chr <- length(chr)
for(i in 1:n.chr) {
pi <- which(thechr==uchr[i])
for(j in i:n.chr) {
pj <- which(thechr==uchr[j])
temp <- lod[pi,pj] - max(maxo[c(i,j)])
temp[!is.na(temp) & temp<0] <- 0
if(i==j) lod[pi,pj][upper.tri(temp)] <- temp[upper.tri(temp)]
else lod[pi,pj] <- temp
}
}
}
if(nodiag) diag(lod) <- 0
# deal with bad LOD score values
if(any(is.na(lod))) {
u <- is.na(lod)
n <- sum(u)
warning(n, " LOD scores NA, set to 0")
lod[u] <- 0
}
if(!allow.neg && any(!is.na(lod) & lod < -1e-6)) {
u <- !is.na(lod) & lod < 0
n <- sum(u)
warning(n, " LOD scores <0, set to 0")
lod[u] <- 0
}
if(any(!is.na(lod) & lod == Inf)) {
u <- !is.na(lod) & lod == Inf
n <- sum(u)
warning(n, " LOD scores =Inf, set to maximum observed value")
lod[u] <- max(lod[!is.na(lod) & lod < Inf])
}
if(missing(zlim)) { # no given zlim
# calculate the zlim for interactive and full LODs
if(allow.neg) {
zlim.int <- max(abs(lod[row(lod) < col(lod)]))
zlim.jnt <- max(abs(lod[row(lod) >= col(lod)]))
}
else {
zlim.int <- max(lod[row(lod) < col(lod)])
zlim.jnt <- max(lod[row(lod) >= col(lod)])
}
}
else {
zlim.jnt <- zlim[1]
if(length(zlim) < 2)
zlim.int <- zlim[1]
else
zlim.int <- zlim[2]
}
# rescale the data in upper triangle based on zlims.jnt
lod[row(lod) < col(lod)] <- lod[row(lod) < col(lod)] * zlim.jnt/zlim.int
# make sure LOD values are below (0,zlim.jnt) or update zlim.jnt
# if(max(lod) > zlim.jnt) {
# warning("LOD values out of range; updating zlim.")
# temp <- max(lod)
# zlim.int <- zlim.int * temp/zlim.jnt
# zlim.jnt <- temp
# }
# save old par parameters, to restore them on exit
if(zscale) {
old.mar <- par("mar")
old.mfrow <- par("mfrow")
old.las <- par("las")
on.exit(par(las = old.las, mar = old.mar, mfrow = old.mfrow))
}
else {
old.las <- par("las")
on.exit(par(las = old.las))
}
par(las = 1)
dots <- list(...)
if(zscale) {
if("layout" %in% names(dots))
layout(dots[["layout"]][[1]],dots[["layout"]][[2]])
else
layout(cbind(1, 2), c(6, 1))
if("mar1" %in% names(dots))
par(mar=dots[["mar1"]])
else
par(mar = c(5, 4, 4, 2) + 0.1)
}
if( gamma < 0 && col.scheme == "redblue")
stop( "gamma must be non-negative" )
cols <- switch(col.scheme,
gray = if( gamma <= 0) rev(gray(seq(0,1,len=256)))
else rev(gray(log(seq(1,exp(gamma),len=256))/gamma)),
heat = heat.colors(256),
terrain = terrain.colors(256),
topo = topo.colors(256),
cm = cm.colors(256),
redblue = rev(rainbow(256, start = 0, end = 2/3)),
viridis = viridis_qtl(256)
)
if(col.scheme=="redblue") {
# convert colors using gamma=0.6 (which will no longer be available in R)
rgbval <- (col2rgb(cols)/255)^0.6
cols <- rgb(rgbval[1,], rgbval[2,], rgbval[3,])
}
if(allow.neg) {
lo <- -zlim.jnt
lo.int <- -zlim.int
}
else lo.int <- lo <- 0
if("xlab" %in% names(dots)) {
xlab <- dots$xlab
if("ylab" %in% names(dots))
ylab <- dots$ylab
else ylab <- xlab
}
else {
if("ylab" %in% names(dots))
xlab <- ylab <- dots$ylab
else {
if(length(chr) > 1)
xlab <- ylab <- "Chromosome"
else
xlab <- ylab <- "Location (cM)"
}
}
if(length(chr) > 1)
image(1:ncol(lod), 1:nrow(lod), lod, ylab = ylab,
xlab = xlab, zlim = c(lo, zlim.jnt), col = cols,
xaxt = "n", yaxt = "n")
else
image(map[,2], map[,2], lod, ylab = ylab,
xlab = xlab, zlim = c(lo, zlim.jnt), col = cols)
# plot point at maximum, if requested
if(point.at.max) {
temp <- lod
temp[upper.tri(temp)] <- -50
temp[diag(temp)] <- -50
wh <- which(temp == max(temp), arr.ind=TRUE)
if(length(chr) > 1)
points(wh,rev(wh),pch=4,lwd=2)
else {
points(map[wh[,1],2],map[wh[,2],2],pch=4,lwd=2,col="blue")
points(map[wh[,2],2],map[wh[,1],2],pch=4,lwd=2,col="blue")
}
}
# add contours if requested
if(any(contours > 0)) {
if(is.logical(contours))
contours = 1.5
tmp = lod
tmp[row(lod) < col(lod)] <- NA
if(length(chr) > 1)
thepos <- 1:ncol(lod)
else thepos <- map[,2]
contour(thepos, thepos, tmp, add = TRUE,drawlabels=FALSE,
levels = max(tmp,na.rm=TRUE) - contours, col = "blue",
lwd = 2)
tmp = lod
tmp[row(lod) > col(lod)] <- NA
contour(thepos, thepos, tmp, add = TRUE,drawlabels=FALSE,
levels = max(tmp,na.rm=TRUE) - contours * zlim.jnt/zlim.int,
col = "blue", lwd = 2)
}
if(length(chr) > 1) {
# calculate how many markers in each chromesome
n.mar <- NULL
for(i in 1:length(chr)) n.mar[i] <- sum(map[, 1] == chr[i])
# plot lines at the chromosome boundaries
if(length(chr) > 1)
wh <- c(0.5, cumsum(n.mar) + 0.5)
abline(v = wh, xpd = FALSE)
abline(h = wh, xpd = FALSE)
# add chromesome numbers
a <- par("usr")
placement <- (wh[-1] + wh[-length(wh)])/2
if(!alternate.chrid || length(chr)<2) {
for(i in 1:length(chr)) {
axis(side=1, at=placement[i], labels=chr[i])
axis(side=2, at=placement[i], labels=chr[i])
}
}
else {
odd <- seq(1, length(chr), by=2)
even <- seq(2, length(chr), by=2)
for(i in odd) {
axis(side=1, at=placement[i], labels="")
axis(side=2, at=placement[i], labels="")
axis(side=1, at=placement[i], labels=chr[i], line=-0.4, tick=FALSE)
axis(side=2, at=placement[i], labels=chr[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=1, at=placement[i], labels="")
axis(side=2, at=placement[i], labels="")
axis(side=1, at=placement[i], labels=chr[i], line=+0.4, tick=FALSE)
axis(side=2, at=placement[i], labels=chr[i], line=+0.4, tick=FALSE)
}
}
}
else {
u <- par("usr")
abline(v=u[1:2], h=u[3:4])
}
# add title
if(!missing(main))
title(main = main)
if(zscale) {
# plot the colormap
dots <- list(...)
if("mar2" %in% names(dots))
par(mar=dots[["mar2"]])
else
par(mar = c(5, 2, 4, 2) + 0.1)
colorstep <- (zlim.jnt-lo)/255
image(x = 1:1, y = seq(lo, zlim.jnt, colorstep), z = matrix(1:256, 1, 256),
zlim = c(1, 256), ylab = "", xlab = "",
xaxt = "n", yaxt = "n", col = cols)
# make sure there's a box around it
u <- par("usr")
abline(v = u[1:2], xpd = FALSE)
abline(h = u[3:4], xpd = FALSE)
if(any(contours) > 0) {
for(i in seq(length(contours))) {
segments(mean(u[1:2]),
max(lod[row(lod) > col(lod)]) - contours[i],
u[2], max(lod[row(lod) > col(lod)]) - contours[i],
xpd = FALSE, col = "blue", lwd = 2)
segments(u[1], max(lod[row(lod) < col(lod)]) - contours[i] *
zlim.jnt/zlim.int, mean(u[1:2]),
max(lod[row(lod) < col(lod)]) - contours[i] * zlim.jnt / zlim.int,
xpd = FALSE, col = "blue", lwd = 2)
}
}
# figure out how big the axis labels should be
fin <- par("fin")[1] # figure width in inches
pin <- par("pin")[1] # plot width in inches
mai <- par("mai")[2] # margin width in inches
# note: pin + 2*mai = fin
xlen.mar <- mai/pin * diff(u[1:2])
# axis for full LODs
yloc <- pretty(c(lo, zlim.jnt), 4)
yloc <- yloc[yloc >= u[3] & yloc <= u[4]]
# segments(u[2], yloc, u[2] + xlen.mar/4, yloc, xpd = TRUE)
# text(u[2] + xlen.mar/3, yloc, as.character(yloc), xpd = TRUE, adj = 0)
axis(side=4, at=yloc, labels=yloc)
# axis for int've LODs
yloc <- pretty(c(lo.int, zlim.int), 4)
yloc.rev <- yloc * zlim.jnt/zlim.int
yloc <- yloc[yloc.rev >= u[3] & yloc.rev <= u[4]]
yloc.rev <- yloc.rev[yloc.rev >= u[3] & yloc.rev <= u[4]]
# segments(u[1], yloc.rev, u[1] - xlen.mar/4, yloc.rev, xpd = TRUE)
# text(u[1] - xlen.mar/3, yloc.rev, as.character(yloc), xpd = TRUE, adj = 1)
if(!addpair)
axis(side=2, at=yloc.rev, labels=yloc)
}
invisible()
}
# end of plot.scantwo.R
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.