Nothing
R/plot.R
In qtl: Tools for Analyzing QTL Experiments
Defines functions
plotPheno
plotPXG
plotInfo
plotGeno
plot.map
geno.image
plotMissing
Documented in
geno.image
plotGeno
plotInfo
plot.map
plotMissing
plotPheno
plotPXG
######################################################################
#
# plot.R
#
# copyright (c) 2000-2019, Karl W Broman
# [modifications of plot.cross from Brian Yandell]
# last modified Dec, 2019
# first written Mar, 2000
#
# 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: plotMissing, plotMap, plot.cross, plotGeno, plotInfo,
# plotPXG, plotPheno
#
######################################################################
plotMissing <-
function(x, chr, reorder=FALSE, main="Missing genotypes",
alternate.chrid=FALSE, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
if(!missing(chr)) cross <- subset(cross,chr=chr)
# get full genotype data into one matrix
Geno <- pull.geno(cross)
# reorder the individuals according to their phenotype
o <- 1:nrow(Geno)
if(reorder) {
# if reorder is a number, use the corresponding phenotype
if(is.numeric(reorder)) {
if(reorder < 1 || reorder > nphe(cross))
stop("reorder should be TRUE, FALSE, or an integer between 1 and", nphe(cross))
o <- order(cross$pheno[,reorder])
}
# otherwise, order according to the sum of the numeric phenotypes
else {
wh <- sapply(cross$pheno, is.numeric)
o <- order(apply(cross$pheno[,wh,drop=FALSE],1,sum))
}
}
# make matrix with 0 where genotype data is missing
# 1 where data is not missing
# 0.5 where data is partially missing
type <- crosstype(cross)
g <- t(Geno[o,])
g[is.na(g)] <- 0
if(type == "bc" || type=="risib" || type=="riself" || type=="bc")
g[g > 0] <- 1
else if(type=="f2") {
g[g > 0 & g < 4] <- 1
g[g > 3] <- 0.5
}
else if(type=="4way") {
g[g > 0 & g < 5] <- 1
g[g > 4] <- 0.5
}
else {
g[g > 0] <- 1
}
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
colors <- c("#000000", "gray80", "#FFFFFF")
# plot grid with black pixels where there is missing data
image(1:nrow(g),1:ncol(g),g,ylab="Individuals",xlab="Markers",col=colors,zlim=c(0,1))
# plot lines at the chromosome boundaries
n.mar <- nmar(cross)
n.chr <- nchr(cross)
a <- c(0.5,cumsum(n.mar)+0.5)
# the following makes the lines go slightly above the plotting region
b <- par("usr")
segments(a,b[3],a,b[4]+diff(b[3:4])*0.02)
# this line adds a line above the image
# (the image function seems to leave it out)
abline(h=0.5+c(0,ncol(g)),xpd=FALSE)
# add chromosome numbers
a <- par("usr")
wh <- cumsum(c(0.5,n.mar))
x <- 1:n.chr
for(i in 1:n.chr)
x[i] <- mean(wh[i+c(0,1)])
thechr <- names(cross$geno)
if(!alternate.chrid || length(thechr) < 2) {
for(i in seq(along=x))
axis(side=3, at=x[i], thechr[i], tick=FALSE, line=-0.5)
}
else {
odd <- seq(1, length(x), by=2)
even <- seq(2, length(x), by=2)
for(i in odd) {
axis(side=3, at=x[i], labels=thechr[i], line=-0.75, tick=FALSE)
}
for(i in even) {
axis(side=3, at=x[i], labels=thechr[i], line=+0, tick=FALSE)
}
}
title(main=main)
invisible()
}
geno.image <-
function(x, chr, reorder=FALSE, main="Genotype data",
alternate.chrid=FALSE, col=NULL, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
if(!missing(chr)) cross <- subset(cross,chr=chr)
type <- crosstype(cross)
# revise X chromosome data
if(type=="bc" || type=="f2") {
chr_type <- sapply(cross$geno, chrtype)
if(any(chr_type=="X")) {
for(i in which(chr_type=="X"))
cross$geno[[i]]$data <- reviseXdata(type, "simple", getsex(cross),
geno=cross$geno[[i]]$data,
cross.attr=attributes(cross))
}
}
# get full genotype data into one matrix
Geno <- pull.geno(cross)
# colors to use
maxgeno <- max(Geno, na.rm=TRUE)
if(type != "4way") {
thecolors <- c("white", "#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00")
thebreaks <- seq(-0.5, 5.5, by=1)
}
else {
if(maxgeno <= 5) {
thecolors <- c("white", "#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00")
thebreaks <- seq(-0.5, 5.5, by=1)
}
else {
thecolors <- c("white", "#8DD3C7", "#FFFFB3", "#BEBADA", "#FB8072",
"#80B1D3", "#FDB462", "#B3DE69", "#FCCDE5", "#D9D9D9", "#BC80BD")
thebreaks <- seq(-0.5, 10.5, by=1)
}
}
thecolors <- thecolors[1:(maxgeno+1)]
thebreaks <- thebreaks[1:(maxgeno+2)]
if(!is.null(col)) { # colors provided
if(length(col) < length(thecolors)) {
warning("col should have length ", maxgeno+1)
thecolors[seq_along(col)] <- col
}
else {
thecolors <- col[1:(maxgeno+1)]
}
}
# reorder the individuals according to their phenotype
o <- 1:nrow(Geno)
if(reorder) {
# if reorder is a number, use the corresponding phenotype
if(is.numeric(reorder)) {
if(reorder < 1 || reorder > nphe(cross))
stop("reorder should be TRUE, FALSE, or an integer between 1 and ", nphe(cross))
o <- order(cross$pheno[,reorder])
}
# otherwise, order according to the sum of the numeric phenotypes
else {
wh <- sapply(cross$pheno, is.numeric)
o <- order(apply(cross$pheno[,wh,drop=FALSE],1,sum))
}
}
g <- t(Geno[o,])
g[is.na(g)] <- 0
# make matrix with 0 where genotype data is missing
# 1 where data is not missing
# 0.5 where data is partially missing
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
# plot grid with black pixels where there is missing data
plot_image_sub <-
function(g, ylab="Individuals",xlab="Markers", col=thecolors, ...)
{
if(length(thebreaks) != length(col)+1)
stop("Must have one more break than color\n",
"length(breaks) = ", length(thebreaks),
"\nlength(col) = ", length(col))
image(1:nrow(g),1:ncol(g), g, col=col, xlab=xlab, ylab=ylab,
breaks=thebreaks, ...)
}
plot_image_sub(g, ...)
# plot lines at the chromosome boundaries
n.mar <- nmar(cross)
n.chr <- nchr(cross)
a <- c(0.5,cumsum(n.mar)+0.5)
# the following makes the lines go slightly above the plotting region
b <- par("usr")
segments(a,b[3],a,b[4]+diff(b[3:4])*0.02)
# this line adds a line above the image
# (the image function seems to leave it out)
abline(h=0.5+c(0,ncol(g)),xpd=FALSE)
# add chromosome numbers
a <- par("usr")
wh <- cumsum(c(0.5,n.mar))
x <- 1:n.chr
for(i in 1:n.chr)
x[i] <- mean(wh[i+c(0,1)])
thechr <- names(cross$geno)
if(!alternate.chrid || length(thechr) < 2) {
for(i in seq(along=x))
axis(side=3, at=x[i], thechr[i], tick=FALSE, line=-0.5)
}
else {
odd <- seq(1, length(x), by=2)
even <- seq(2, length(x), by=2)
for(i in odd)
axis(side=3, at=x[i], labels=thechr[i], line=-0.75, tick=FALSE)
for(i in even)
axis(side=3, at=x[i], labels=thechr[i], line=+0, tick=FALSE)
}
title(main=main)
invisible()
}
plotMap <- plot.map <-
function(x, map2, chr, horizontal=FALSE, shift=TRUE,
show.marker.names=FALSE, alternate.chrid=FALSE, ...)
{
dots <- list(...)
if("main" %in% names(dots)) {
themain <- dots$main
usemaindefault <- FALSE
}
else usemaindefault <- TRUE
if("xlim" %in% names(dots)) {
xlim <- dots$xlim
usexlimdefault <- FALSE
}
else usexlimdefault <- TRUE
if("ylim" %in% names(dots)) {
ylim <- dots$ylim
useylimdefault <- FALSE
}
else useylimdefault <- TRUE
if("xlab" %in% names(dots))
xlab <- dots$xlab
else {
if(horizontal)
xlab <- "Location (cM)"
else
xlab <- "Chromosome"
}
if("ylab" %in% names(dots))
ylab <- dots$ylab
else {
if(horizontal)
ylab <- "Chromosome"
else
ylab <- "Location (cM)"
}
map <- x
# figure out if the input is a cross (containing a map)
# or is the map itself
if(inherits(map, "cross"))
map <- pull.map(map)
if(!missing(map2) && inherits(map2, "cross"))
map2 <- pull.map(map2)
if(!inherits(map, "map") || (!missing(map2) && !inherits(map2, "map")))
warning("Input should have class \"cross\" or \"map\".")
if(!missing(map2) && is.matrix(map[[1]]) != is.matrix(map2[[1]]))
stop("Maps must be both sex-specific or neither sex-specific.")
if(!missing(chr)) {
map <- map[matchchr(chr, names(map))]
if(!missing(map2)) map2 <- map2[matchchr(chr, names(map2))]
}
sex.sp <- FALSE
if(is.matrix(map[[1]])) { # sex-specific map
one.map <- FALSE
sex.sp <- TRUE
if(!missing(map2)) {
if(is.logical(map2)) {
horizontal <- map2
map2 <- lapply(map,function(a) a[2,])
map <- lapply(map,function(a) a[1,])
}
else {
Map1 <- lapply(map,function(a) a[1,,drop=TRUE])
Map2 <- lapply(map,function(a) a[2,,drop=TRUE])
Map3 <- lapply(map2,function(a) a[1,,drop=TRUE])
Map4 <- lapply(map2,function(a) a[2,,drop=TRUE])
old.mfrow <- par("mfrow")
on.exit(par(mfrow=old.mfrow))
par(mfrow=c(2,1))
class(Map1) <- class(Map2) <- class(Map3) <- class(Map4) <- "map"
plotMap(Map1,Map3,horizontal=horizontal,shift=shift,
show.marker.names=show.marker.names,alternate.chrid=alternate.chrid)
plotMap(Map2,Map4,horizontal=horizontal,shift=shift,
show.marker.names=show.marker.names,alternate.chrid=alternate.chrid)
return(invisible(NULL))
}
}
else {
map2 <- lapply(map,function(a) a[2,])
map <- lapply(map,function(a) a[1,])
}
}
else { # single map
# determine whether a second map was given
if(!missing(map2))
one.map <- FALSE
else one.map <- TRUE
}
if(one.map) {
n.chr <- length(map)
if(!show.marker.names) { # locations of chromosomes
chrpos <- 1:n.chr
thelim <- range(chrpos)+c(-0.5, 0.5)
}
else {
chrpos <- seq(1, n.chr*2, by=2)
thelim <- range(chrpos)+c(-0.35, 2.35)
}
if(shift) map <- lapply(map, function(a) a-a[1])
maxlen <- max(unlist(lapply(map,max)))
if(horizontal) {
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE, las=1)
on.exit(par(xpd=old.xpd,las=old.las))
if(usexlimdefault) xlim <- c(0,maxlen)
if(useylimdefault) ylim <- rev(thelim)
plot(0,0,type="n",xlim=xlim, ylim=ylim,yaxs="i",
xlab=xlab, ylab=ylab, yaxt="n")
a <- par("usr")
for(i in 1:n.chr) {
segments(min(map[[i]]), chrpos[i], max(map[[i]]), chrpos[i])
segments(map[[i]], chrpos[i]-0.25, map[[i]], chrpos[i]+0.25)
if(show.marker.names)
text(map[[i]], chrpos[i]+0.35, names(map[[i]]), srt=90, adj=c(1,0.5))
}
# add chromosome labels
if(!alternate.chrid || length(chrpos) < 2) {
for(i in seq(along=chrpos))
axis(side=2, at=chrpos[i], labels=names(map)[i])
}
else {
odd <- seq(1, length(chrpos), by=2)
even <- seq(2, length(chrpos), by=2)
for(i in odd) {
axis(side=2, at=chrpos[i], labels="")
axis(side=2, at=chrpos[i], labels=names(map)[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=2, at=chrpos[i], labels="")
axis(side=2, at=chrpos[i], labels=names(map)[i], line=+0.4, tick=FALSE)
}
}
}
else {
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
if(usexlimdefault) xlim <- thelim
if(useylimdefault) ylim <- c(maxlen, 0)
plot(0,0,type="n",ylim=ylim,xlim=xlim,xaxs="i",
xlab=xlab, ylab=ylab, xaxt="n")
a <- par("usr")
for(i in 1:n.chr) {
segments(chrpos[i], min(map[[i]]), chrpos[i], max(map[[i]]))
segments(chrpos[i]-0.25, map[[i]], chrpos[i]+0.25, map[[i]])
if(show.marker.names)
text(chrpos[i]+0.35, map[[i]], names(map[[i]]), adj=c(0,0.5))
}
# add chromosome labels
if(!alternate.chrid || length(chrpos) < 2) {
for(i in seq(along=chrpos))
axis(side=1, at=chrpos[i], labels=names(map)[i])
}
else {
odd <- seq(1, length(chrpos), by=2)
even <- seq(2, length(chrpos), by=2)
for(i in odd) {
axis(side=1, at=chrpos[i], labels="")
axis(side=1, at=chrpos[i], labels=names(map)[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=1, at=chrpos[i], labels="")
axis(side=1, at=chrpos[i], labels=names(map)[i], line=+0.4, tick=FALSE)
}
}
}
if(usemaindefault)
title(main="Genetic map")
else if(themain != "")
title(main=themain)
}
else {
map1 <- map
# check that maps conform
if(is.matrix(map2[[1]]))
stop("Second map appears to be a sex-specific map.")
if(length(map1) != length(map2))
stop("Maps have different numbers of chromosomes.")
if(any(names(map1) != names(map2))) {
cat("Map1: ", names(map1), "\n")
cat("Map2: ", names(map2), "\n")
stop("Maps have different chromosome names.")
}
if(shift) {
map1 <- lapply(map1,function(a) a-a[1])
map2 <- lapply(map2,function(a) a-a[1])
}
n.mar1 <- sapply(map1, length)
n.mar2 <- sapply(map2, length)
markernames1 <- lapply(map1, names)
markernames2 <- lapply(map2, names)
if(any(n.mar1 != n.mar2)) {
if(show.marker.names) {
warning("Can't show marker names because of different numbers of markers.")
show.marker.names <- FALSE
}
}
else if(any(unlist(markernames1) != unlist(markernames2))) {
if(show.marker.names) {
warning("Can't show marker names because markers in different orders.")
show.marker.names <- FALSE
}
}
n.chr <- length(map1)
maxloc <- max(c(unlist(lapply(map1,max)),unlist(lapply(map2,max))))
if(!show.marker.names) { # locations of chromosomes
chrpos <- 1:n.chr
thelim <- range(chrpos)+c(-0.5, 0.5)
}
else {
chrpos <- seq(1, n.chr*2, by=2)
thelim <- range(chrpos)+c(-0.4, 2.4)
}
if(!horizontal) {
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
if(usexlimdefault) xlim <- thelim
if(useylimdefault) ylim <- c(maxloc, 0)
plot(0,0,type="n",ylim=ylim,xlim=xlim, xaxs="i",
xlab=xlab, ylab=ylab, xaxt="n")
a <- par("usr")
for(i in 1:n.chr) {
if(max(map2[[i]]) < max(map1[[i]]))
map2[[i]] <- map2[[i]] + (max(map1[[i]])-max(map2[[i]]))/2
else
map1[[i]] <- map1[[i]] + (max(map2[[i]])-max(map1[[i]]))/2
segments(chrpos[i]-0.3, min(map1[[i]]), chrpos[i]-0.3, max(map1[[i]]))
segments(chrpos[i]+0.3, min(map2[[i]]), chrpos[i]+0.3, max(map2[[i]]))
# lines between markers
wh <- match(markernames1[[i]], markernames2[[i]])
for(j in which(!is.na(wh)))
segments(chrpos[i]-0.3, map1[[i]][j], chrpos[i]+0.3, map2[[i]][wh[j]])
if(any(is.na(wh)))
segments(chrpos[i]-0.4, map1[[i]][is.na(wh)], chrpos[i]-0.2, map1[[i]][is.na(wh)])
wh <- match(markernames2[[i]], markernames1[[i]])
if(any(is.na(wh)))
segments(chrpos[i]+0.4, map2[[i]][is.na(wh)], chrpos[i]+0.2, map2[[i]][is.na(wh)])
if(show.marker.names)
text(chrpos[i]+0.35, map2[[i]], names(map2[[i]]), adj=c(0,0.5))
}
# add chromosome labels
if(!alternate.chrid || length(chrpos) < 2) {
for(i in seq(along=chrpos))
axis(side=1, at=chrpos[i], labels=names(map1)[i])
}
else {
odd <- seq(1, length(chrpos), by=2)
even <- seq(2, length(chrpos), by=2)
for(i in odd) {
axis(side=1, at=chrpos[i], labels="")
axis(side=1, at=chrpos[i], labels=names(map1)[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=1, at=chrpos[i], labels="")
axis(side=1, at=chrpos[i], labels=names(map1)[i], line=+0.4, tick=FALSE)
}
}
}
else {
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
if(usexlimdefault) xlim <- c(0,maxloc)
if(useylimdefault) ylim <- rev(thelim)
plot(0,0,type="n",xlim=xlim,ylim=ylim,
xlab=xlab, ylab=ylab, yaxt="n", yaxs="i")
a <- par("usr")
for(i in 1:n.chr) {
if(max(map2[[i]]) < max(map1[[i]]))
map2[[i]] <- map2[[i]] + (max(map1[[i]])-max(map2[[i]]))/2
else
map1[[i]] <- map1[[i]] + (max(map2[[i]])-max(map1[[i]]))/2
segments(min(map1[[i]]), chrpos[i]-0.3, max(map1[[i]]), chrpos[[i]]-0.3)
segments(min(map2[[i]]), chrpos[i]+0.3, max(map2[[i]]), chrpos[[i]]+0.3)
# lines between markers
wh <- match(markernames1[[i]], markernames2[[i]])
for(j in which(!is.na(wh)))
segments(map1[[i]][j], chrpos[i]-0.3, map2[[i]][wh[j]], chrpos[i]+0.3)
if(any(is.na(wh)))
segments(map1[[i]][is.na(wh)], chrpos[i]-0.4, map1[[i]][is.na(wh)], chrpos[i]-0.2)
wh <- match(markernames2[[i]], markernames1[[i]])
if(any(is.na(wh)))
segments(map2[[i]][is.na(wh)], chrpos[i]+0.4, map2[[i]][is.na(wh)], chrpos[i]+0.2)
if(show.marker.names)
text(map2[[i]], chrpos[i]+0.35, names(map2[[i]]), srt=90, adj=c(1,0.5))
}
# add chromosome labels
if(!alternate.chrid || length(chrpos) < 2) {
for(i in seq(along=chrpos))
axis(side=2, at=chrpos[i], labels=names(map1)[i])
}
else {
odd <- seq(1, length(chrpos), by=2)
even <- seq(2, length(chrpos), by=2)
for(i in odd) {
axis(side=2, at=chrpos[i], labels="")
axis(side=2, at=chrpos[i], labels=names(map1)[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=2, at=chrpos[i], labels="")
axis(side=2, at=chrpos[i], labels=names(map1)[i], line=+0.4, tick=FALSE)
}
}
}
if(usemaindefault) {
if(!sex.sp) title(main="Comparison of genetic maps")
else title(main="Genetic map")
}
else if(themain != "")
title(main=themain)
}
invisible()
}
plot.cross <-
function (x, auto.layout = TRUE, pheno.col,
alternate.chrid=TRUE, ...)
{
# look to see whether this should really be shipped to plotMap
if(inherits(auto.layout, "map") &&
(inherits(x, "map") || inherits(x, "cross"))) {
plotMap(x, auto.layout, alternate.chrid=alternate.chrid, ...)
return(invisible())
}
# make sure this is a cross
if(!inherits(x, "cross"))
stop("Input should have class \"cross\".")
old.yaxt <- par("yaxt")
old.mfrow <- par("mfrow")
on.exit(par(yaxt = old.yaxt, mfrow = old.mfrow))
n.phe <- nphe(x)
if(missing(pheno.col)) pheno.col <- 1:n.phe
if(is.character(pheno.col)) {
temp <- match(pheno.col, names(x$pheno))
if(any(is.na(temp)))
warning("Some phenotypes not found:",
paste(pheno.col[is.na(temp)], collapse=" "))
pheno.col <- temp[!is.na(temp)]
}
else pheno.col <- (1:nphe(x))[pheno.col]
n.plot = length(pheno.col) + 2
# automatically choose row/column structure for the plots
if(auto.layout) {
nr <- ceiling(sqrt(n.plot))
nc <- ceiling((n.plot)/nr)
par(mfrow = c(nr, nc))
}
plotMissing(x,alternate.chrid=alternate.chrid)
plotMap(x,alternate.chrid=alternate.chrid)
for(i in pheno.col) plotPheno(x, pheno.col=i)
invisible()
}
##################################################r####################
#
# plotGeno: Plot genotypes for a specified chromosome, with likely
# genotyping errors indicated.
#
######################################################################
plotGeno <-
function(x, chr, ind, include.xo=TRUE, horizontal=TRUE,
cutoff=4, min.sep=2, cex=1.2, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
if(missing(chr)) chr <- names(cross$geno)[1]
cross <- subset(cross,chr=chr)
if(nchr(cross) > 1)
cross <- subset(cross,chr=names(cross$geno)[1])
if(!missing(ind)) {
if(is.null(getid(cross))) cross$pheno$id <- 1:nind(cross)
if(!is.logical(ind)) ind <- unique(ind)
cross <- subset(cross, ind=ind)
}
id <- getid(cross)
if(is.null(id)) id <- 1:nind(cross)
use.id <- TRUE
type <- crosstype(cross)
old.las <- par("las")
on.exit(par(las=old.las))
par(las=1)
if(!("errorlod" %in% names(cross$geno[[1]]))) {
warning("First running calc.errorlod.")
cross <- calc.errorlod(cross,error.prob=0.01)
}
# indicators for apparent errors
errors <- matrix(0,ncol=ncol(cross$geno[[1]]$data),
nrow=nrow(cross$geno[[1]]$data))
dimnames(errors) <- dimnames(cross$geno[[1]]$data)
top <- top.errorlod(cross,names(cross$geno)[1],cutoff,FALSE)
if(length(top) > 0)
for(i in 1:nrow(top))
errors[match(top[i,2],id),as.character(top[i,3])] <- 1
# map, data, errors
map <- cross$geno[[1]]$map
if(is.matrix(map)) map <- map[1,] # if sex-specific map
L <- diff(range(map))
min.d <- L*min.sep/100
d <- diff(map)
d[d < min.d] <- min.d
map <- cumsum(c(0,d))
cross$geno[[1]]$map <- map
n.ind <- nrow(errors)
color <- c("white","gray60","black","green","orange","red")
# revise X chr data for backcross/intercross
data <- cross$geno[[1]]$data
chr_type <- chrtype(cross$geno[[1]])
if(chr_type=="X" && (type=="f2" || type=="bc"))
data <- reviseXdata(type, sexpgm=getsex(cross), geno=data, cross.attr=attributes(cross), force=TRUE)
if(include.xo) {
if(type != "4way") { # find crossover locations
xoloc <- locateXO(cross)
xoloc <- data.frame(ind=rep(1:length(xoloc),sapply(xoloc,length)),
loc=unlist(xoloc), stringsAsFactors=TRUE)
}
else { # 4-way cross
mcross <- dcross <- cross
class(mcross) <- class(dcross) <- c("bc", "cross")
mcross$geno[[1]]$data[!is.na(data) & data==1 | data==3 | data==5] <- 1
mcross$geno[[1]]$data[!is.na(data) & data==2 | data==4 | data==6] <- 2
mcross$geno[[1]]$data[!is.na(data) & data==7 | data==8 | data==9 | data==10] <- NA
dcross$geno[[1]]$data[!is.na(data) & data==1 | data==2 | data==7] <- 1
dcross$geno[[1]]$data[!is.na(data) & data==3 | data==4 | data==8] <- 2
dcross$geno[[1]]$data[!is.na(data) & data==5 | data==6 | data==9 | data==10] <- NA
mxoloc <- locateXO(mcross)
mxoloc <- data.frame(ind=rep(1:length(mxoloc),sapply(mxoloc,length)),
loc=unlist(mxoloc), stringsAsFactors=TRUE)
dxoloc <- locateXO(dcross)
dxoloc <- data.frame(ind=rep(1:length(dxoloc),sapply(dxoloc,length)),
loc=unlist(dxoloc), stringsAsFactors=TRUE)
}
}
# check for 'main' in the ...
args <- list(...)
if("main" %in% names(args))
themain <- args$main
else
themain <- paste("Chromosome",names(cross$geno)[1])
# check for 'xlim' and 'ylim'
if("xlim" %in% names(args)) thexlim <- args$xlim
else thexlim <- NULL
if("ylim" %in% names(args)) theylim <- args$ylim
else theylim <- NULL
if(type=="4way") {
jit <- 0.15
mdata <- data
ddata <- data
# mom's allele
mdata[!is.na(data) & (data==1 | data==3 | data==5)] <- 1
mdata[!is.na(data) & (data==2 | data==4 | data==6)] <- 2
mdata[!is.na(data) & (data==7 | data==8)] <- NA
# dad's allele
ddata[!is.na(data) & (data==1 | data==2 | data==7)] <- 1
ddata[!is.na(data) & (data==3 | data==4 | data==8)] <- 2
ddata[!is.na(data) & (data==5 | data==6)] <- NA
if(horizontal) {
if(is.null(thexlim)) thexlim <- c(0, max(map))
if(is.null(theylim)) theylim <- c(n.ind+1, 0)
plot(0,0,type="n",xlab="Location (cM)",ylab="Individual",
main=themain,
ylim=theylim,xlim=thexlim, yaxt="n", yaxs="i")
segments(0, 1:n.ind-jit, max(map), 1:n.ind-jit)
segments(0, 1:n.ind+jit, max(map), 1:n.ind+jit)
if(use.id) axis(side=2, at=1:n.ind, labels=id)
else axis(side=2, at=1:n.ind)
# A alleles
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=1] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind-jit,pch=21,col="black", bg=color[1],cex=cex)
# B alleles
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=2] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind-jit,pch=21,col="black", bg=color[3],cex=cex)
# 9/10 genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=9] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind-jit,pch=21,col="black", bg=color[4],cex=cex)
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=10] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind-jit,pch=21,col="black", bg=color[5],cex=cex)
# C alleles
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=1] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind+jit,pch=21,col="black", bg=color[1],cex=cex)
# D alleles
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=2] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind+jit,pch=21,col="black", bg=color[3],cex=cex)
# 9/10 genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=9] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind+jit,pch=21,col="black", bg=color[4],cex=cex)
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=10] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind+jit,pch=21,col="black", bg=color[5],cex=cex)
# plot map
u <- par("usr")
segments(map,u[3],map,u[3]-1/2)
segments(map,u[4],map,u[4]+1/2)
if(any(errors != 0)) {
ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
points(x,ind-jit,pch=0,col=color[6],cex=cex+0.4,lwd=2)
points(x,ind+jit,pch=0,col=color[6],cex=cex+0.4,lwd=2)
}
if(include.xo) {
points(mxoloc$loc,mxoloc$ind-jit,pch=4,col="blue",lwd=2)
points(dxoloc$loc,dxoloc$ind+jit,pch=4,col="blue",lwd=2)
}
}
else {
if(is.null(theylim)) theylim <- c(max(map), 0)
if(is.null(thexlim)) thexlim <- c(0, n.ind+1)
plot(0,0,type="n",ylab="Location (cM)",xlab="Individual",
main=themain,
xlim=thexlim,ylim=theylim, xaxt="n", xaxs="i")
segments(1:n.ind-jit, 0, 1:n.ind-jit, max(map))
segments(1:n.ind+jit, 0, 1:n.ind+jit, max(map))
if(use.id) axis(side=1, at=1:n.ind, labels=id)
else axis(side=1, at=1:n.ind)
# A alleles
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=1] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind-jit,y,pch=21,col="black",bg=color[1],cex=cex)
# B alleles
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=2] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind-jit,y,pch=21,col="black",bg=color[3],cex=cex)
# 9/10 genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=9] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind-jit,y,pch=21,col="black", bg=color[4],cex=cex)
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=10] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind-jit,y,pch=21,col="black", bg=color[5],cex=cex)
# C alleles
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=1] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind+jit,y,pch=21,col="black", bg=color[1],cex=cex)
# D alleles
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=2] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind+jit,y,pch=21,col="black", bg=color[3],cex=cex)
# 9/10 genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=9] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind+jit,y,pch=21,col="black", bg=color[4],cex=cex)
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=10] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind+jit,y,pch=21,col="black", bg=color[5],cex=cex)
# plot map
u <- par("usr")
segments(u[1],map,(u[1]+1)/2,map)
segments(u[2],map,(n.ind+u[2])/2,map)
if(any(errors != 0)) {
ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
points(ind-jit,y,pch=0,col=color[6],cex=cex+0.4,lwd=2)
points(ind+jit,y,pch=0,col=color[6],cex=cex+0.4,lwd=2)
}
if(include.xo) {
points(mxoloc$ind-jit,mxoloc$loc,pch=4,col="blue",lwd=2)
points(dxoloc$ind+jit,dxoloc$loc,pch=4,col="blue",lwd=2)
}
}
}
else {
if(horizontal) {
if(is.null(thexlim)) thexlim <- c(0, max(map))
if(is.null(theylim)) theylim <- c(n.ind+0.5,0.5)
plot(0,0,type="n",xlab="Location (cM)",ylab="Individual",
main=themain,
ylim=theylim,xlim=thexlim, yaxt="n")
segments(0, 1:n.ind, max(map), 1:n.ind)
if(use.id) axis(side=2, at=1:n.ind, labels=id)
else axis(side=2)
# AA genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(data)] <- NA
ind <- tind; ind[!is.na(data) & data!=1] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind,pch=21,col="black", bg=color[1],cex=cex)
# AB genotypes
ind <- tind; ind[!is.na(data) & data!=2] <- NA
if(type=="f2" || (type=="bc" && chr_type=="X"))
points(x,ind,pch=21,col="black", bg=color[2],cex=cex)
else points(x,ind,pch=21,col="black", bg=color[3],cex=cex)
if(type=="f2" || (type=="bc" && chr_type=="X")) {
# BB genotypes
ind <- tind; ind[!is.na(data) & data!=3] <- NA
points(x,ind,pch=21,col="black", bg=color[3],cex=cex)
}
if(type=="f2") {
# not BB (D in mapmaker/qtl) genotypes
ind <- tind; ind[!is.na(data) & data!=4] <- NA
points(x,ind,pch=21,col="black", bg=color[4],cex=cex)
# not AA (C in mapmaker/qtl) genotypes
ind <- tind; ind[!is.na(data) & data!=5] <- NA
points(x,ind,pch=21,col="black", bg=color[5],cex=cex)
}
# plot map
u <- par("usr")
segments(map,u[3],map,u[3]-1/2)
segments(map,u[4],map,u[4]+1/2)
if(any(errors != 0)) {
ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
points(x,ind,pch=0,col=color[6],cex=cex+0.4,lwd=2)
}
if(include.xo) points(xoloc$loc,xoloc$ind,pch=4,col="blue",lwd=2)
}
else {
if(is.null(theylim)) theylim <- c(max(map), 0)
if(is.null(thexlim)) thexlim <- c(0.5,n.ind+0.5)
plot(0,0,type="n",ylab="Location (cM)",xlab="Individual",
main=themain,
xlim=thexlim,ylim=theylim, xaxt="n")
segments(1:n.ind,0,1:n.ind,max(map))
if(use.id) axis(side=1, at=1:n.ind, labels=id)
else axis(side=1)
# AA genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(data)] <- NA
ind <- tind; ind[!is.na(data) & data!=1] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind,y,pch=21,col="black", bg="white",cex=cex)
# AB genotypes
ind <- tind; ind[!is.na(data) & data!=2] <- NA
if(type=="f2" || (type=="bc" && chr_type=="X"))
points(ind,y,pch=21,col="black", bg=color[2],cex=cex)
else points(ind,y,pch=21,col="black", bg=color[3],cex=cex)
if(type=="f2" || (type=="bc" && chr_type=="X")) {
# BB genotypes
ind <- tind; ind[!is.na(data) & data!=3] <- NA
points(ind,y,pch=21,col="black", bg=color[3],cex=cex)
}
if(type=="f2") {
# not BB genotypes
ind <- tind; ind[!is.na(data) & data!=4] <- NA
points(ind,y,pch=21,col="black", bg=color[4],cex=cex)
# not AA genotypes
ind <- tind; ind[!is.na(data) & data!=5] <- NA
points(ind,y,pch=21,col="black", bg=color[5],cex=cex)
}
# plot map
u <- par("usr")
segments(u[1],map,(u[1]+1)/2,map)
segments(u[2],map,(n.ind+u[2])/2,map)
if(any(errors != 0)) {
ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
points(ind,y,pch=0,col=color[6],cex=cex+0.4,lwd=2)
}
if(include.xo) points(xoloc$ind,xoloc$loc,pch=4,col="blue",lwd=2)
}
}
invisible()
}
######################################################################
#
# plotInfo: Plot the proportion of missing information in the
# genotype data.
#
######################################################################
plotInfo <-
function(x,chr,method=c("entropy","variance","both"), step=1,
off.end=0, error.prob=0.001,
map.function=c("haldane","kosambi","c-f","morgan"),
alternate.chrid=FALSE, fourwaycross=c("all", "AB", "CD"),
include.genofreq=FALSE, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
method <- match(match.arg(method),c("entropy","variance","both"))-1
map.function <- match.arg(map.function)
if(!missing(chr)) cross <- subset(cross,chr=chr)
n.chr <- nchr(cross)
results <- NULL
cross <- calc.genoprob(cross, step=step, error.prob=error.prob,
off.end=off.end, map.function=map.function)
gap <- 25
# 4-way cross: can consider just A/B or just C/D
fourwaycross <- match.arg(fourwaycross)
n.ind <- nind(cross)
if(include.genofreq) theprob <- NULL
for(i in 1:n.chr) {
n.gen <- dim(cross$geno[[i]]$prob)[3]
n.pos <- ncol(cross$geno[[i]]$prob)
# 4-way cross: can consider just A/B or just C/D
if(n.gen==4 && (fourwaycross!="all")) {
att <- attributes(cross$geno[[i]]$prob)
if(fourwaycross == "AB") {
cross$geno[[i]]$prob[,,1] <- cross$geno[[i]]$prob[,,1] + cross$geno[[i]]$prob[,,3]
cross$geno[[i]]$prob[,,2] <- cross$geno[[i]]$prob[,,2] + cross$geno[[i]]$prob[,,4]
}
else {
cross$geno[[i]]$prob[,,1] <- cross$geno[[i]]$prob[,,1] + cross$geno[[i]]$prob[,,2]
cross$geno[[i]]$prob[,,2] <- cross$geno[[i]]$prob[,,3] + cross$geno[[i]]$prob[,,4]
}
cross$geno[[i]]$prob <- cross$geno[[i]]$prob[,,1:2]
n.gen <- 2
for(j in seq(along=att))
if(names(att)[j] != "dim" && names(att)[j] != "dimnames")
attr(cross$geno[[i]]$prob, names(att)[j]) <- att[[j]]
}
# calculate information (between 0 and 1)
info <- .C("R_info",
as.integer(n.ind),
as.integer(n.pos),
as.integer(n.gen),
as.double(cross$geno[[i]]$prob),
info1=as.double(rep(0,n.pos)),
info2=as.double(rep(0,n.pos)),
as.integer(method),
PACKAGE="qtl")
if(method != 1) { # rescale entropy version
if(n.gen==3) maxent <- 1.5*log(2)
else maxent <- log(n.gen)
info$info1 <- -info$info1/maxent
}
if(method != 0) { # rescale variance version
maxvar <- c(0.25,0.5,1.25)[n.gen-1]
info$info2 <- info$info2/maxvar
}
# reconstruct map
if("map" %in% names(attributes(cross$geno[[i]]$prob)))
map <- attr(cross$geno[[i]]$prob,"map")
else {
stp <- attr(cross$geno[[i]]$prob, "step")
oe <- attr(cross$geno[[i]]$prob, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[i]]$prob)))
stpw <- attr(cross$geno[[i]]$prob, "stepwidth")
else stpw <- "fixed"
map <- create.map(cross$geno[[i]]$map,stp,oe,stpw)
}
if(is.matrix(map)) map <- map[1,]
z <- data.frame(chr=rep(names(cross$geno)[i],length(map)),
pos=as.numeric(map),
"Missing information"=info$info1,
"Missing information"=info$info2, stringsAsFactors=TRUE)
w <- names(map)
o <- grep("^loc-*[0-9]+",w)
if(length(o) > 0) # inter-marker locations cited as "c*.loc*"
w[o] <- paste("c",names(cross$geno)[i],".",w[o],sep="")
rownames(z) <- w
results <- rbind(results, z)
if(include.genofreq) {
p <- cross$geno[[i]]$prob
if(chrtype(cross$geno[[i]])=="X")
p <- reviseXdata(crosstype(cross), expandX="full", sexpgm=getsex(cross), prob=p,
cross.attr=attributes(cross))
p <- apply(p, 2:3, mean, na.rm=TRUE)
if(is.null(theprob))
theprob <- p
else {
m1 <- match(colnames(p), colnames(theprob))
m2 <- match(colnames(theprob), colnames(p))
if(!any(is.na(m1)) && !any(is.na(m2))) {
theprob <- rbind(theprob, p[,colnames(theprob)])
}
else { # some differences in column names
if(all(!is.na(m1))) { # need to add some new columns to p
noldcol <- ncol(p)
nnewcol <- sum(is.na(m2))
for(i in 1:nnewcol)
p <- cbind(p, 0)
colnames(p)[-(1:noldcol)] <- colnames(theprob)[is.na(m2)]
theprob <- rbind(theprob, p[,colnames(theprob)])
}
else if(all(!is.na(m2))) { # need to add some new columns to theprob
noldcol <- ncol(theprob)
nnewcol <- sum(is.na(m1))
for(i in 1:nnewcol)
theprob <- cbind(theprob, 0)
colnames(theprob)[-(1:noldcol)] <- colnames(p)[is.na(m1)]
theprob <- rbind(theprob, p[,colnames(theprob)])
}
else { # need to add some new columns to each
noldcol <- ncol(theprob)
nnewcol <- sum(is.na(m1))
oldcolnam <- colnames(theprob)
for(i in 1:nnewcol)
theprob <- cbind(theprob, 0)
colnames(theprob)[-(1:noldcol)] <- colnames(p)[is.na(m1)]
noldcol <- ncol(p)
nnewcol <- sum(is.na(m2))
for(i in 1:nnewcol)
p <- cbind(p, 0)
colnames(p)[-(1:noldcol)] <- oldcolnam[is.na(m2)]
theprob <- rbind(theprob, p[,colnames(theprob)])
}
}
}
} # end of if(include.genofreq)
}
colnames(results)[3:4] <- c("misinfo.entropy","misinfo.variance")
if(method==0) results <- results[,-4]
if(method==1) results <- results[,-3]
if(include.genofreq)
results <- cbind(results, theprob)
class(results) <- c("scanone","data.frame")
# check whether main was included as an argument
args <- list(...)
if("main" %in% names(args)) hasmain <- TRUE
else hasmain <- FALSE
# check whether gap was included as an argument
if(!("gap" %in% names(args))) {
if(method==0) {
if(hasmain)
plot.scanone(results,ylim=c(0,1),gap=gap,
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,ylim=c(0,1),gap=gap,
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
else if(method==1) {
if(hasmain)
plot.scanone(results,ylim=c(0,1),gap=gap,
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,ylim=c(0,1),gap=gap,
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
else if(method==2) {
if(hasmain)
plot.scanone(results,results,lodcolumn=1:2,ylim=c(0,1),gap=gap,
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,results,lodcolumn=1:2,ylim=c(0,1),gap=gap,
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
}
else { # gap was included in ...
if(method==0) {
if(hasmain)
plot.scanone(results,ylim=c(0,1),
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,ylim=c(0,1),
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
else if(method==1) {
if(hasmain)
plot.scanone(results,ylim=c(0,1),
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,ylim=c(0,1),
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
else if(method==2) {
if(hasmain)
plot.scanone(results,results,lodcolumn=1:2,ylim=c(0,1),
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,results,lodcolumn=1:2,ylim=c(0,1),
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
}
invisible(results)
}
# plot phenotypes against one or more markers
plotPXG <-
function(x, marker, pheno.col = 1, jitter = 1, infer = TRUE,
pch, ylab, main, col, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
type <- crosstype(cross)
if(LikePheVector(pheno.col, nind(cross), nphe(cross))) {
cross$pheno <- cbind(pheno.col, cross$pheno)
pheno.col <- 1
}
if(length(pheno.col) > 1) {
pheno.col <- pheno.col[1]
warning("plotPXG can take just one phenotype; only the first will be used")
}
if(is.character(pheno.col)) {
num <- find.pheno(cross, pheno.col)
if(is.na(num))
stop("Couldn't identify phenotype \"", pheno.col, "\"")
pheno.col <- num
}
if(pheno.col < 1 | pheno.col > nphe(cross))
stop("pheno.col values should be between 1 and the no. phenotypes")
if(!is.numeric(cross$pheno[,pheno.col]))
stop("phenotype \"", colnames(cross$pheno)[pheno.col], "\" is not numeric.")
if(missing(pch)) pch <- par("pch")
if(missing(ylab)) ylab <- colnames(cross$pheno)[pheno.col]
oldlas <- par("las")
on.exit(par(las = oldlas))
par(las = 1)
# find chromosomes containing the markers
o <- sapply(cross$geno, function(a, b) b %in% colnames(a$data),
marker)
if(length(marker)==1) o <- matrix(o,nrow=1)
if(!all(apply(o,1,any))) {
oo <- apply(o,1,any)
stop("Marker ", marker[!oo], " not found")
}
n.mark <- length(marker)
o <- apply(o, 1, which)
chr <- names(cross$geno)[o]
uchr <- unique(chr)
cross <- subset(cross, chr=uchr)
map <- pull.map(cross)
pos <- NULL
for(i in seq(length(chr))) pos[i] <- map[[chr[i]]][marker[i]]
chr_type <- sapply(cross$geno, chrtype)
names(chr_type) <- names(cross$geno)
chr_type <- chr_type[chr]
# if X chromosome and backcross or intercross, get sex/direction data
if(any(chr_type == "X") && (type == "bc" || type == "f2"))
sexpgm <- getsex(cross)
else sexpgm <- NULL
# number of possible genotypes
gen.names <- list()
for(i in seq(length(chr)))
gen.names[[i]] <- getgenonames(type, chr_type[i], "full", sexpgm, attributes(cross))
n.gen <- sapply(gen.names, length)
jitter <- jitter/10
if(any(n.gen == 2)) jitter <- jitter * 0.75
# function to determine whether genotype is fully known
tempf <-
function(x, type)
{
tmp <- is.na(x)
if(type=="f2") tmp[!is.na(x) & x>3] <- TRUE
if(type=="4way") tmp[!is.na(x) & x>4] <- TRUE
tmp
}
# if infer=TRUE, fill in genotype data by a single imputation
if(infer) {
which.missing <- tempf(cross$geno[[chr[1]]]$data[, marker[1]],type)
if(n.mark > 1)
for(i in 2:n.mark)
which.missing <- which.missing | tempf(cross$geno[[chr[i]]]$data[,marker[i]],type)
which.missing <- as.numeric(which.missing)
cross <- fill.geno(cross, method = "imp")
}
else which.missing <- rep(1,nind(cross))
# data to plot
x <- cross$geno[[chr[1]]]$data[, marker[1]]
if(n.mark > 1)
for(i in 2:n.mark)
x <- cbind(x, cross$geno[[chr[i]]]$data[, marker[i]])
else x <- as.matrix(x)
y <- cross$pheno[, pheno.col]
if(!infer) { # replace partially informative genotypes with NAs
if(type == "f2") x[x > 3] <- NA
if(type == "4way") x[x > 4] <- NA
if(sum(!is.na(x)) == 0)
stop("Can't use infer=FALSE as there are no fully informative genotypes")
}
# in case of X chromosome, recode some genotypes
if(any(chr_type == "X") && (type == "bc" || type == "f2")) {
ix = seq(n.mark)[chr_type == "X"]
for(i in ix)
x[, i] <- as.numeric(reviseXdata(type, "full", sexpgm,
geno = as.matrix(x[, i]),
cross.attr=attributes(cross)))
}
# save all of the data, returned invisibly
data <- as.data.frame(x, stringsAsFactors=TRUE)
names(data) <- marker
for(i in marker) data[[i]] <- ordered(data[[i]])
data$pheno <- y
data$inferred <- which.missing
# re-code the multi-marker genotypes
if(n.mark > 1) {
for(i in 2:n.mark)
x[, 1] <- n.gen[i] * (x[, 1] - 1) + x[, i]
}
x <- x[, 1]
observed <- sort(unique(x))
# amount of jitter
u <- runif(nind(cross), -jitter, jitter)
r <- (1 - 2 * jitter)/2
# genotype names
if(n.mark == 1)
gnames <- gen.names[[1]]
else {
gnames <- array(gen.names[[n.mark]], c(prod(n.gen), n.mark))
for(i in (n.mark - 1):1) {
tmpi <- rep(gen.names[[i]], rep(prod(n.gen[(i + 1):n.mark]),
n.gen[i]))
if(i > 1)
tmpi <- rep(tmpi, prod(n.gen[1:(i - 1)]))
gnames[, i] <- tmpi
}
gnames <- apply(gnames, 1, function(x) paste(x, collapse = "\n"))
}
# create plot
plot(x + u, y, xlab = "Genotype", ylab = ylab, type = "n",
main = "", xlim = c(1 - r + jitter, length(gnames) + r + jitter), xaxt = "n",
...)
# marker names at top
if(missing(main))
mtext(paste(marker, collapse = "\n"),
line=0.5, cex = ifelse(n.mark==1, 1.2, 0.8))
else
title(main=main)
abline(v = 1:prod(n.gen), col = "gray", lty = 3)
if(length(pch) == 1)
pch = rep(pch, length(x))
if(infer) {
points((x + u)[which.missing == 1], y[which.missing ==
1], col = "red", pch = pch[which.missing == 1])
points((x + u)[which.missing == 0], y[which.missing ==
0], pch = pch[which.missing == 0])
}
else points(x + u, y, pch = pch)
sux = sort(unique(x))
# add confidence intervals
me <- se <- array(NA, length(gnames))
me[sux] <- tapply(y, x, mean, na.rm = TRUE)
se[sux] <- tapply(y, x, function(a) sd(a, na.rm = TRUE)/sqrt(sum(!is.na(a))))
thecolors <- c("black", "blue", "red", "purple", "green", "orange")
if(missing(col)) {
col <- thecolors[1:n.gen[n.mark]]
if(n.gen[n.mark] == 3)
col <- c("blue", "purple", "red")
else if(n.gen[n.mark] == 2)
col <- c("blue", "red")
}
ng <- length(gnames)
segments(seq(ng) + jitter * 2, me, seq(ng) +
jitter * 4, me, lwd = 2, col = col)
segments(seq(ng) + jitter * 3, me - se, seq(ng) +
jitter * 3, me + se, lwd = 2, col = col)
segments(seq(ng) + jitter * 2.5, me - se, seq(ng) +
jitter * 3.5, me - se, lwd = 2, col = col)
segments(seq(ng) + jitter * 2.5, me + se, seq(ng) +
jitter * 3.5, me + se, lwd = 2, col = col)
# add genotypes below
u <- par("usr")
cxaxis <- par("cex.axis")
segments(seq(along=gnames), u[3], seq(along=gnames), u[3] - diff(u[3:4]) *
0.015, xpd = TRUE)
axis(side=1, at=seq(along=gnames), labels=gnames,
cex=ifelse(n.mark==1, cxaxis, cxaxis*0.8),
tick=FALSE, line = (length(marker)-1)/2)
invisible(data)
}
plotPheno <-
function(x, pheno.col=1, ...)
{
if(!inherits(x, "cross"))
stop("Input should have class \"cross\".")
if(LikePheVector(pheno.col, nind(x), nphe(x))) {
x$pheno <- cbind(pheno.col, x$pheno)
pheno.col <- 1
}
if(length(pheno.col) > 1) {
pheno.col <- pheno.col[1]
warning("Ignoring all but the first element in pheno.col.")
}
if(is.character(pheno.col)) {
num <- find.pheno(x, pheno.col)
if(is.na(num))
stop("Couldn't identify phenotype \"", pheno.col, "\"")
pheno.col <- num
}
if(pheno.col < 1 | pheno.col > nphe(x))
warning("pheno.col should be between 1 and ", nphe(x))
phe <- x$pheno[,pheno.col]
u <- length(unique(phe))
if(u==2 || (u < 10 && nind(x) > 50))
phe <- as.factor(phe)
plot_pheno_sub <-
function(phe, xlab=paste("phe", pheno.col),
main=colnames(x$pheno)[pheno.col], col="white",
breaks=ceiling(2*sqrt(nind(x))),
las=1,
...)
{
if(is.factor(phe)) {
barplot(table(phe), xlab=xlab, main=main, col=col, las=las, ...)
}
else {
phe <- as.numeric(phe)[1:nind(x)]
hist(phe, breaks = breaks,
xlab = xlab, main = main, las=las, ...)
}
}
plot_pheno_sub(phe, ...)
}
# end of plot.R
Try the qtl package in your browser
Any scripts or data that you put into this service are public.
qtl documentation built on Sept. 11, 2024, 5:43 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plotPheno plotPXG plotInfo plotGeno plot.map geno.image plotMissing
Documented in geno.image plotGeno plotInfo plot.map plotMissing plotPheno plotPXG
######################################################################
#
# plot.R
#
# copyright (c) 2000-2019, Karl W Broman
# [modifications of plot.cross from Brian Yandell]
# last modified Dec, 2019
# first written Mar, 2000
#
# 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: plotMissing, plotMap, plot.cross, plotGeno, plotInfo,
# plotPXG, plotPheno
#
######################################################################
plotMissing <-
function(x, chr, reorder=FALSE, main="Missing genotypes",
alternate.chrid=FALSE, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
if(!missing(chr)) cross <- subset(cross,chr=chr)
# get full genotype data into one matrix
Geno <- pull.geno(cross)
# reorder the individuals according to their phenotype
o <- 1:nrow(Geno)
if(reorder) {
# if reorder is a number, use the corresponding phenotype
if(is.numeric(reorder)) {
if(reorder < 1 || reorder > nphe(cross))
stop("reorder should be TRUE, FALSE, or an integer between 1 and", nphe(cross))
o <- order(cross$pheno[,reorder])
}
# otherwise, order according to the sum of the numeric phenotypes
else {
wh <- sapply(cross$pheno, is.numeric)
o <- order(apply(cross$pheno[,wh,drop=FALSE],1,sum))
}
}
# make matrix with 0 where genotype data is missing
# 1 where data is not missing
# 0.5 where data is partially missing
type <- crosstype(cross)
g <- t(Geno[o,])
g[is.na(g)] <- 0
if(type == "bc" || type=="risib" || type=="riself" || type=="bc")
g[g > 0] <- 1
else if(type=="f2") {
g[g > 0 & g < 4] <- 1
g[g > 3] <- 0.5
}
else if(type=="4way") {
g[g > 0 & g < 5] <- 1
g[g > 4] <- 0.5
}
else {
g[g > 0] <- 1
}
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
colors <- c("#000000", "gray80", "#FFFFFF")
# plot grid with black pixels where there is missing data
image(1:nrow(g),1:ncol(g),g,ylab="Individuals",xlab="Markers",col=colors,zlim=c(0,1))
# plot lines at the chromosome boundaries
n.mar <- nmar(cross)
n.chr <- nchr(cross)
a <- c(0.5,cumsum(n.mar)+0.5)
# the following makes the lines go slightly above the plotting region
b <- par("usr")
segments(a,b[3],a,b[4]+diff(b[3:4])*0.02)
# this line adds a line above the image
# (the image function seems to leave it out)
abline(h=0.5+c(0,ncol(g)),xpd=FALSE)
# add chromosome numbers
a <- par("usr")
wh <- cumsum(c(0.5,n.mar))
x <- 1:n.chr
for(i in 1:n.chr)
x[i] <- mean(wh[i+c(0,1)])
thechr <- names(cross$geno)
if(!alternate.chrid || length(thechr) < 2) {
for(i in seq(along=x))
axis(side=3, at=x[i], thechr[i], tick=FALSE, line=-0.5)
}
else {
odd <- seq(1, length(x), by=2)
even <- seq(2, length(x), by=2)
for(i in odd) {
axis(side=3, at=x[i], labels=thechr[i], line=-0.75, tick=FALSE)
}
for(i in even) {
axis(side=3, at=x[i], labels=thechr[i], line=+0, tick=FALSE)
}
}
title(main=main)
invisible()
}
geno.image <-
function(x, chr, reorder=FALSE, main="Genotype data",
alternate.chrid=FALSE, col=NULL, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
if(!missing(chr)) cross <- subset(cross,chr=chr)
type <- crosstype(cross)
# revise X chromosome data
if(type=="bc" || type=="f2") {
chr_type <- sapply(cross$geno, chrtype)
if(any(chr_type=="X")) {
for(i in which(chr_type=="X"))
cross$geno[[i]]$data <- reviseXdata(type, "simple", getsex(cross),
geno=cross$geno[[i]]$data,
cross.attr=attributes(cross))
}
}
# get full genotype data into one matrix
Geno <- pull.geno(cross)
# colors to use
maxgeno <- max(Geno, na.rm=TRUE)
if(type != "4way") {
thecolors <- c("white", "#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00")
thebreaks <- seq(-0.5, 5.5, by=1)
}
else {
if(maxgeno <= 5) {
thecolors <- c("white", "#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00")
thebreaks <- seq(-0.5, 5.5, by=1)
}
else {
thecolors <- c("white", "#8DD3C7", "#FFFFB3", "#BEBADA", "#FB8072",
"#80B1D3", "#FDB462", "#B3DE69", "#FCCDE5", "#D9D9D9", "#BC80BD")
thebreaks <- seq(-0.5, 10.5, by=1)
}
}
thecolors <- thecolors[1:(maxgeno+1)]
thebreaks <- thebreaks[1:(maxgeno+2)]
if(!is.null(col)) { # colors provided
if(length(col) < length(thecolors)) {
warning("col should have length ", maxgeno+1)
thecolors[seq_along(col)] <- col
}
else {
thecolors <- col[1:(maxgeno+1)]
}
}
# reorder the individuals according to their phenotype
o <- 1:nrow(Geno)
if(reorder) {
# if reorder is a number, use the corresponding phenotype
if(is.numeric(reorder)) {
if(reorder < 1 || reorder > nphe(cross))
stop("reorder should be TRUE, FALSE, or an integer between 1 and ", nphe(cross))
o <- order(cross$pheno[,reorder])
}
# otherwise, order according to the sum of the numeric phenotypes
else {
wh <- sapply(cross$pheno, is.numeric)
o <- order(apply(cross$pheno[,wh,drop=FALSE],1,sum))
}
}
g <- t(Geno[o,])
g[is.na(g)] <- 0
# make matrix with 0 where genotype data is missing
# 1 where data is not missing
# 0.5 where data is partially missing
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
# plot grid with black pixels where there is missing data
plot_image_sub <-
function(g, ylab="Individuals",xlab="Markers", col=thecolors, ...)
{
if(length(thebreaks) != length(col)+1)
stop("Must have one more break than color\n",
"length(breaks) = ", length(thebreaks),
"\nlength(col) = ", length(col))
image(1:nrow(g),1:ncol(g), g, col=col, xlab=xlab, ylab=ylab,
breaks=thebreaks, ...)
}
plot_image_sub(g, ...)
# plot lines at the chromosome boundaries
n.mar <- nmar(cross)
n.chr <- nchr(cross)
a <- c(0.5,cumsum(n.mar)+0.5)
# the following makes the lines go slightly above the plotting region
b <- par("usr")
segments(a,b[3],a,b[4]+diff(b[3:4])*0.02)
# this line adds a line above the image
# (the image function seems to leave it out)
abline(h=0.5+c(0,ncol(g)),xpd=FALSE)
# add chromosome numbers
a <- par("usr")
wh <- cumsum(c(0.5,n.mar))
x <- 1:n.chr
for(i in 1:n.chr)
x[i] <- mean(wh[i+c(0,1)])
thechr <- names(cross$geno)
if(!alternate.chrid || length(thechr) < 2) {
for(i in seq(along=x))
axis(side=3, at=x[i], thechr[i], tick=FALSE, line=-0.5)
}
else {
odd <- seq(1, length(x), by=2)
even <- seq(2, length(x), by=2)
for(i in odd)
axis(side=3, at=x[i], labels=thechr[i], line=-0.75, tick=FALSE)
for(i in even)
axis(side=3, at=x[i], labels=thechr[i], line=+0, tick=FALSE)
}
title(main=main)
invisible()
}
plotMap <- plot.map <-
function(x, map2, chr, horizontal=FALSE, shift=TRUE,
show.marker.names=FALSE, alternate.chrid=FALSE, ...)
{
dots <- list(...)
if("main" %in% names(dots)) {
themain <- dots$main
usemaindefault <- FALSE
}
else usemaindefault <- TRUE
if("xlim" %in% names(dots)) {
xlim <- dots$xlim
usexlimdefault <- FALSE
}
else usexlimdefault <- TRUE
if("ylim" %in% names(dots)) {
ylim <- dots$ylim
useylimdefault <- FALSE
}
else useylimdefault <- TRUE
if("xlab" %in% names(dots))
xlab <- dots$xlab
else {
if(horizontal)
xlab <- "Location (cM)"
else
xlab <- "Chromosome"
}
if("ylab" %in% names(dots))
ylab <- dots$ylab
else {
if(horizontal)
ylab <- "Chromosome"
else
ylab <- "Location (cM)"
}
map <- x
# figure out if the input is a cross (containing a map)
# or is the map itself
if(inherits(map, "cross"))
map <- pull.map(map)
if(!missing(map2) && inherits(map2, "cross"))
map2 <- pull.map(map2)
if(!inherits(map, "map") || (!missing(map2) && !inherits(map2, "map")))
warning("Input should have class \"cross\" or \"map\".")
if(!missing(map2) && is.matrix(map[[1]]) != is.matrix(map2[[1]]))
stop("Maps must be both sex-specific or neither sex-specific.")
if(!missing(chr)) {
map <- map[matchchr(chr, names(map))]
if(!missing(map2)) map2 <- map2[matchchr(chr, names(map2))]
}
sex.sp <- FALSE
if(is.matrix(map[[1]])) { # sex-specific map
one.map <- FALSE
sex.sp <- TRUE
if(!missing(map2)) {
if(is.logical(map2)) {
horizontal <- map2
map2 <- lapply(map,function(a) a[2,])
map <- lapply(map,function(a) a[1,])
}
else {
Map1 <- lapply(map,function(a) a[1,,drop=TRUE])
Map2 <- lapply(map,function(a) a[2,,drop=TRUE])
Map3 <- lapply(map2,function(a) a[1,,drop=TRUE])
Map4 <- lapply(map2,function(a) a[2,,drop=TRUE])
old.mfrow <- par("mfrow")
on.exit(par(mfrow=old.mfrow))
par(mfrow=c(2,1))
class(Map1) <- class(Map2) <- class(Map3) <- class(Map4) <- "map"
plotMap(Map1,Map3,horizontal=horizontal,shift=shift,
show.marker.names=show.marker.names,alternate.chrid=alternate.chrid)
plotMap(Map2,Map4,horizontal=horizontal,shift=shift,
show.marker.names=show.marker.names,alternate.chrid=alternate.chrid)
return(invisible(NULL))
}
}
else {
map2 <- lapply(map,function(a) a[2,])
map <- lapply(map,function(a) a[1,])
}
}
else { # single map
# determine whether a second map was given
if(!missing(map2))
one.map <- FALSE
else one.map <- TRUE
}
if(one.map) {
n.chr <- length(map)
if(!show.marker.names) { # locations of chromosomes
chrpos <- 1:n.chr
thelim <- range(chrpos)+c(-0.5, 0.5)
}
else {
chrpos <- seq(1, n.chr*2, by=2)
thelim <- range(chrpos)+c(-0.35, 2.35)
}
if(shift) map <- lapply(map, function(a) a-a[1])
maxlen <- max(unlist(lapply(map,max)))
if(horizontal) {
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE, las=1)
on.exit(par(xpd=old.xpd,las=old.las))
if(usexlimdefault) xlim <- c(0,maxlen)
if(useylimdefault) ylim <- rev(thelim)
plot(0,0,type="n",xlim=xlim, ylim=ylim,yaxs="i",
xlab=xlab, ylab=ylab, yaxt="n")
a <- par("usr")
for(i in 1:n.chr) {
segments(min(map[[i]]), chrpos[i], max(map[[i]]), chrpos[i])
segments(map[[i]], chrpos[i]-0.25, map[[i]], chrpos[i]+0.25)
if(show.marker.names)
text(map[[i]], chrpos[i]+0.35, names(map[[i]]), srt=90, adj=c(1,0.5))
}
# add chromosome labels
if(!alternate.chrid || length(chrpos) < 2) {
for(i in seq(along=chrpos))
axis(side=2, at=chrpos[i], labels=names(map)[i])
}
else {
odd <- seq(1, length(chrpos), by=2)
even <- seq(2, length(chrpos), by=2)
for(i in odd) {
axis(side=2, at=chrpos[i], labels="")
axis(side=2, at=chrpos[i], labels=names(map)[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=2, at=chrpos[i], labels="")
axis(side=2, at=chrpos[i], labels=names(map)[i], line=+0.4, tick=FALSE)
}
}
}
else {
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
if(usexlimdefault) xlim <- thelim
if(useylimdefault) ylim <- c(maxlen, 0)
plot(0,0,type="n",ylim=ylim,xlim=xlim,xaxs="i",
xlab=xlab, ylab=ylab, xaxt="n")
a <- par("usr")
for(i in 1:n.chr) {
segments(chrpos[i], min(map[[i]]), chrpos[i], max(map[[i]]))
segments(chrpos[i]-0.25, map[[i]], chrpos[i]+0.25, map[[i]])
if(show.marker.names)
text(chrpos[i]+0.35, map[[i]], names(map[[i]]), adj=c(0,0.5))
}
# add chromosome labels
if(!alternate.chrid || length(chrpos) < 2) {
for(i in seq(along=chrpos))
axis(side=1, at=chrpos[i], labels=names(map)[i])
}
else {
odd <- seq(1, length(chrpos), by=2)
even <- seq(2, length(chrpos), by=2)
for(i in odd) {
axis(side=1, at=chrpos[i], labels="")
axis(side=1, at=chrpos[i], labels=names(map)[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=1, at=chrpos[i], labels="")
axis(side=1, at=chrpos[i], labels=names(map)[i], line=+0.4, tick=FALSE)
}
}
}
if(usemaindefault)
title(main="Genetic map")
else if(themain != "")
title(main=themain)
}
else {
map1 <- map
# check that maps conform
if(is.matrix(map2[[1]]))
stop("Second map appears to be a sex-specific map.")
if(length(map1) != length(map2))
stop("Maps have different numbers of chromosomes.")
if(any(names(map1) != names(map2))) {
cat("Map1: ", names(map1), "\n")
cat("Map2: ", names(map2), "\n")
stop("Maps have different chromosome names.")
}
if(shift) {
map1 <- lapply(map1,function(a) a-a[1])
map2 <- lapply(map2,function(a) a-a[1])
}
n.mar1 <- sapply(map1, length)
n.mar2 <- sapply(map2, length)
markernames1 <- lapply(map1, names)
markernames2 <- lapply(map2, names)
if(any(n.mar1 != n.mar2)) {
if(show.marker.names) {
warning("Can't show marker names because of different numbers of markers.")
show.marker.names <- FALSE
}
}
else if(any(unlist(markernames1) != unlist(markernames2))) {
if(show.marker.names) {
warning("Can't show marker names because markers in different orders.")
show.marker.names <- FALSE
}
}
n.chr <- length(map1)
maxloc <- max(c(unlist(lapply(map1,max)),unlist(lapply(map2,max))))
if(!show.marker.names) { # locations of chromosomes
chrpos <- 1:n.chr
thelim <- range(chrpos)+c(-0.5, 0.5)
}
else {
chrpos <- seq(1, n.chr*2, by=2)
thelim <- range(chrpos)+c(-0.4, 2.4)
}
if(!horizontal) {
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
if(usexlimdefault) xlim <- thelim
if(useylimdefault) ylim <- c(maxloc, 0)
plot(0,0,type="n",ylim=ylim,xlim=xlim, xaxs="i",
xlab=xlab, ylab=ylab, xaxt="n")
a <- par("usr")
for(i in 1:n.chr) {
if(max(map2[[i]]) < max(map1[[i]]))
map2[[i]] <- map2[[i]] + (max(map1[[i]])-max(map2[[i]]))/2
else
map1[[i]] <- map1[[i]] + (max(map2[[i]])-max(map1[[i]]))/2
segments(chrpos[i]-0.3, min(map1[[i]]), chrpos[i]-0.3, max(map1[[i]]))
segments(chrpos[i]+0.3, min(map2[[i]]), chrpos[i]+0.3, max(map2[[i]]))
# lines between markers
wh <- match(markernames1[[i]], markernames2[[i]])
for(j in which(!is.na(wh)))
segments(chrpos[i]-0.3, map1[[i]][j], chrpos[i]+0.3, map2[[i]][wh[j]])
if(any(is.na(wh)))
segments(chrpos[i]-0.4, map1[[i]][is.na(wh)], chrpos[i]-0.2, map1[[i]][is.na(wh)])
wh <- match(markernames2[[i]], markernames1[[i]])
if(any(is.na(wh)))
segments(chrpos[i]+0.4, map2[[i]][is.na(wh)], chrpos[i]+0.2, map2[[i]][is.na(wh)])
if(show.marker.names)
text(chrpos[i]+0.35, map2[[i]], names(map2[[i]]), adj=c(0,0.5))
}
# add chromosome labels
if(!alternate.chrid || length(chrpos) < 2) {
for(i in seq(along=chrpos))
axis(side=1, at=chrpos[i], labels=names(map1)[i])
}
else {
odd <- seq(1, length(chrpos), by=2)
even <- seq(2, length(chrpos), by=2)
for(i in odd) {
axis(side=1, at=chrpos[i], labels="")
axis(side=1, at=chrpos[i], labels=names(map1)[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=1, at=chrpos[i], labels="")
axis(side=1, at=chrpos[i], labels=names(map1)[i], line=+0.4, tick=FALSE)
}
}
}
else {
old.xpd <- par("xpd")
old.las <- par("las")
par(xpd=TRUE,las=1)
on.exit(par(xpd=old.xpd,las=old.las))
if(usexlimdefault) xlim <- c(0,maxloc)
if(useylimdefault) ylim <- rev(thelim)
plot(0,0,type="n",xlim=xlim,ylim=ylim,
xlab=xlab, ylab=ylab, yaxt="n", yaxs="i")
a <- par("usr")
for(i in 1:n.chr) {
if(max(map2[[i]]) < max(map1[[i]]))
map2[[i]] <- map2[[i]] + (max(map1[[i]])-max(map2[[i]]))/2
else
map1[[i]] <- map1[[i]] + (max(map2[[i]])-max(map1[[i]]))/2
segments(min(map1[[i]]), chrpos[i]-0.3, max(map1[[i]]), chrpos[[i]]-0.3)
segments(min(map2[[i]]), chrpos[i]+0.3, max(map2[[i]]), chrpos[[i]]+0.3)
# lines between markers
wh <- match(markernames1[[i]], markernames2[[i]])
for(j in which(!is.na(wh)))
segments(map1[[i]][j], chrpos[i]-0.3, map2[[i]][wh[j]], chrpos[i]+0.3)
if(any(is.na(wh)))
segments(map1[[i]][is.na(wh)], chrpos[i]-0.4, map1[[i]][is.na(wh)], chrpos[i]-0.2)
wh <- match(markernames2[[i]], markernames1[[i]])
if(any(is.na(wh)))
segments(map2[[i]][is.na(wh)], chrpos[i]+0.4, map2[[i]][is.na(wh)], chrpos[i]+0.2)
if(show.marker.names)
text(map2[[i]], chrpos[i]+0.35, names(map2[[i]]), srt=90, adj=c(1,0.5))
}
# add chromosome labels
if(!alternate.chrid || length(chrpos) < 2) {
for(i in seq(along=chrpos))
axis(side=2, at=chrpos[i], labels=names(map1)[i])
}
else {
odd <- seq(1, length(chrpos), by=2)
even <- seq(2, length(chrpos), by=2)
for(i in odd) {
axis(side=2, at=chrpos[i], labels="")
axis(side=2, at=chrpos[i], labels=names(map1)[i], line=-0.4, tick=FALSE)
}
for(i in even) {
axis(side=2, at=chrpos[i], labels="")
axis(side=2, at=chrpos[i], labels=names(map1)[i], line=+0.4, tick=FALSE)
}
}
}
if(usemaindefault) {
if(!sex.sp) title(main="Comparison of genetic maps")
else title(main="Genetic map")
}
else if(themain != "")
title(main=themain)
}
invisible()
}
plot.cross <-
function (x, auto.layout = TRUE, pheno.col,
alternate.chrid=TRUE, ...)
{
# look to see whether this should really be shipped to plotMap
if(inherits(auto.layout, "map") &&
(inherits(x, "map") || inherits(x, "cross"))) {
plotMap(x, auto.layout, alternate.chrid=alternate.chrid, ...)
return(invisible())
}
# make sure this is a cross
if(!inherits(x, "cross"))
stop("Input should have class \"cross\".")
old.yaxt <- par("yaxt")
old.mfrow <- par("mfrow")
on.exit(par(yaxt = old.yaxt, mfrow = old.mfrow))
n.phe <- nphe(x)
if(missing(pheno.col)) pheno.col <- 1:n.phe
if(is.character(pheno.col)) {
temp <- match(pheno.col, names(x$pheno))
if(any(is.na(temp)))
warning("Some phenotypes not found:",
paste(pheno.col[is.na(temp)], collapse=" "))
pheno.col <- temp[!is.na(temp)]
}
else pheno.col <- (1:nphe(x))[pheno.col]
n.plot = length(pheno.col) + 2
# automatically choose row/column structure for the plots
if(auto.layout) {
nr <- ceiling(sqrt(n.plot))
nc <- ceiling((n.plot)/nr)
par(mfrow = c(nr, nc))
}
plotMissing(x,alternate.chrid=alternate.chrid)
plotMap(x,alternate.chrid=alternate.chrid)
for(i in pheno.col) plotPheno(x, pheno.col=i)
invisible()
}
##################################################r####################
#
# plotGeno: Plot genotypes for a specified chromosome, with likely
# genotyping errors indicated.
#
######################################################################
plotGeno <-
function(x, chr, ind, include.xo=TRUE, horizontal=TRUE,
cutoff=4, min.sep=2, cex=1.2, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
if(missing(chr)) chr <- names(cross$geno)[1]
cross <- subset(cross,chr=chr)
if(nchr(cross) > 1)
cross <- subset(cross,chr=names(cross$geno)[1])
if(!missing(ind)) {
if(is.null(getid(cross))) cross$pheno$id <- 1:nind(cross)
if(!is.logical(ind)) ind <- unique(ind)
cross <- subset(cross, ind=ind)
}
id <- getid(cross)
if(is.null(id)) id <- 1:nind(cross)
use.id <- TRUE
type <- crosstype(cross)
old.las <- par("las")
on.exit(par(las=old.las))
par(las=1)
if(!("errorlod" %in% names(cross$geno[[1]]))) {
warning("First running calc.errorlod.")
cross <- calc.errorlod(cross,error.prob=0.01)
}
# indicators for apparent errors
errors <- matrix(0,ncol=ncol(cross$geno[[1]]$data),
nrow=nrow(cross$geno[[1]]$data))
dimnames(errors) <- dimnames(cross$geno[[1]]$data)
top <- top.errorlod(cross,names(cross$geno)[1],cutoff,FALSE)
if(length(top) > 0)
for(i in 1:nrow(top))
errors[match(top[i,2],id),as.character(top[i,3])] <- 1
# map, data, errors
map <- cross$geno[[1]]$map
if(is.matrix(map)) map <- map[1,] # if sex-specific map
L <- diff(range(map))
min.d <- L*min.sep/100
d <- diff(map)
d[d < min.d] <- min.d
map <- cumsum(c(0,d))
cross$geno[[1]]$map <- map
n.ind <- nrow(errors)
color <- c("white","gray60","black","green","orange","red")
# revise X chr data for backcross/intercross
data <- cross$geno[[1]]$data
chr_type <- chrtype(cross$geno[[1]])
if(chr_type=="X" && (type=="f2" || type=="bc"))
data <- reviseXdata(type, sexpgm=getsex(cross), geno=data, cross.attr=attributes(cross), force=TRUE)
if(include.xo) {
if(type != "4way") { # find crossover locations
xoloc <- locateXO(cross)
xoloc <- data.frame(ind=rep(1:length(xoloc),sapply(xoloc,length)),
loc=unlist(xoloc), stringsAsFactors=TRUE)
}
else { # 4-way cross
mcross <- dcross <- cross
class(mcross) <- class(dcross) <- c("bc", "cross")
mcross$geno[[1]]$data[!is.na(data) & data==1 | data==3 | data==5] <- 1
mcross$geno[[1]]$data[!is.na(data) & data==2 | data==4 | data==6] <- 2
mcross$geno[[1]]$data[!is.na(data) & data==7 | data==8 | data==9 | data==10] <- NA
dcross$geno[[1]]$data[!is.na(data) & data==1 | data==2 | data==7] <- 1
dcross$geno[[1]]$data[!is.na(data) & data==3 | data==4 | data==8] <- 2
dcross$geno[[1]]$data[!is.na(data) & data==5 | data==6 | data==9 | data==10] <- NA
mxoloc <- locateXO(mcross)
mxoloc <- data.frame(ind=rep(1:length(mxoloc),sapply(mxoloc,length)),
loc=unlist(mxoloc), stringsAsFactors=TRUE)
dxoloc <- locateXO(dcross)
dxoloc <- data.frame(ind=rep(1:length(dxoloc),sapply(dxoloc,length)),
loc=unlist(dxoloc), stringsAsFactors=TRUE)
}
}
# check for 'main' in the ...
args <- list(...)
if("main" %in% names(args))
themain <- args$main
else
themain <- paste("Chromosome",names(cross$geno)[1])
# check for 'xlim' and 'ylim'
if("xlim" %in% names(args)) thexlim <- args$xlim
else thexlim <- NULL
if("ylim" %in% names(args)) theylim <- args$ylim
else theylim <- NULL
if(type=="4way") {
jit <- 0.15
mdata <- data
ddata <- data
# mom's allele
mdata[!is.na(data) & (data==1 | data==3 | data==5)] <- 1
mdata[!is.na(data) & (data==2 | data==4 | data==6)] <- 2
mdata[!is.na(data) & (data==7 | data==8)] <- NA
# dad's allele
ddata[!is.na(data) & (data==1 | data==2 | data==7)] <- 1
ddata[!is.na(data) & (data==3 | data==4 | data==8)] <- 2
ddata[!is.na(data) & (data==5 | data==6)] <- NA
if(horizontal) {
if(is.null(thexlim)) thexlim <- c(0, max(map))
if(is.null(theylim)) theylim <- c(n.ind+1, 0)
plot(0,0,type="n",xlab="Location (cM)",ylab="Individual",
main=themain,
ylim=theylim,xlim=thexlim, yaxt="n", yaxs="i")
segments(0, 1:n.ind-jit, max(map), 1:n.ind-jit)
segments(0, 1:n.ind+jit, max(map), 1:n.ind+jit)
if(use.id) axis(side=2, at=1:n.ind, labels=id)
else axis(side=2, at=1:n.ind)
# A alleles
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=1] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind-jit,pch=21,col="black", bg=color[1],cex=cex)
# B alleles
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=2] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind-jit,pch=21,col="black", bg=color[3],cex=cex)
# 9/10 genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=9] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind-jit,pch=21,col="black", bg=color[4],cex=cex)
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=10] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind-jit,pch=21,col="black", bg=color[5],cex=cex)
# C alleles
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=1] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind+jit,pch=21,col="black", bg=color[1],cex=cex)
# D alleles
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=2] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind+jit,pch=21,col="black", bg=color[3],cex=cex)
# 9/10 genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=9] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind+jit,pch=21,col="black", bg=color[4],cex=cex)
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=10] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind+jit,pch=21,col="black", bg=color[5],cex=cex)
# plot map
u <- par("usr")
segments(map,u[3],map,u[3]-1/2)
segments(map,u[4],map,u[4]+1/2)
if(any(errors != 0)) {
ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
points(x,ind-jit,pch=0,col=color[6],cex=cex+0.4,lwd=2)
points(x,ind+jit,pch=0,col=color[6],cex=cex+0.4,lwd=2)
}
if(include.xo) {
points(mxoloc$loc,mxoloc$ind-jit,pch=4,col="blue",lwd=2)
points(dxoloc$loc,dxoloc$ind+jit,pch=4,col="blue",lwd=2)
}
}
else {
if(is.null(theylim)) theylim <- c(max(map), 0)
if(is.null(thexlim)) thexlim <- c(0, n.ind+1)
plot(0,0,type="n",ylab="Location (cM)",xlab="Individual",
main=themain,
xlim=thexlim,ylim=theylim, xaxt="n", xaxs="i")
segments(1:n.ind-jit, 0, 1:n.ind-jit, max(map))
segments(1:n.ind+jit, 0, 1:n.ind+jit, max(map))
if(use.id) axis(side=1, at=1:n.ind, labels=id)
else axis(side=1, at=1:n.ind)
# A alleles
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=1] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind-jit,y,pch=21,col="black",bg=color[1],cex=cex)
# B alleles
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=2] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind-jit,y,pch=21,col="black",bg=color[3],cex=cex)
# 9/10 genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=9] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind-jit,y,pch=21,col="black", bg=color[4],cex=cex)
tind <- rep(1:n.ind,length(map));tind[is.na(mdata)] <- NA
ind <- tind; ind[!is.na(mdata) & mdata!=10] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind-jit,y,pch=21,col="black", bg=color[5],cex=cex)
# C alleles
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=1] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind+jit,y,pch=21,col="black", bg=color[1],cex=cex)
# D alleles
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=2] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind+jit,y,pch=21,col="black", bg=color[3],cex=cex)
# 9/10 genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=9] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind+jit,y,pch=21,col="black", bg=color[4],cex=cex)
tind <- rep(1:n.ind,length(map));tind[is.na(ddata)] <- NA
ind <- tind; ind[!is.na(ddata) & ddata!=10] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind+jit,y,pch=21,col="black", bg=color[5],cex=cex)
# plot map
u <- par("usr")
segments(u[1],map,(u[1]+1)/2,map)
segments(u[2],map,(n.ind+u[2])/2,map)
if(any(errors != 0)) {
ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
points(ind-jit,y,pch=0,col=color[6],cex=cex+0.4,lwd=2)
points(ind+jit,y,pch=0,col=color[6],cex=cex+0.4,lwd=2)
}
if(include.xo) {
points(mxoloc$ind-jit,mxoloc$loc,pch=4,col="blue",lwd=2)
points(dxoloc$ind+jit,dxoloc$loc,pch=4,col="blue",lwd=2)
}
}
}
else {
if(horizontal) {
if(is.null(thexlim)) thexlim <- c(0, max(map))
if(is.null(theylim)) theylim <- c(n.ind+0.5,0.5)
plot(0,0,type="n",xlab="Location (cM)",ylab="Individual",
main=themain,
ylim=theylim,xlim=thexlim, yaxt="n")
segments(0, 1:n.ind, max(map), 1:n.ind)
if(use.id) axis(side=2, at=1:n.ind, labels=id)
else axis(side=2)
# AA genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(data)] <- NA
ind <- tind; ind[!is.na(data) & data!=1] <- NA
x <- rep(map,rep(n.ind,length(map)))
points(x,ind,pch=21,col="black", bg=color[1],cex=cex)
# AB genotypes
ind <- tind; ind[!is.na(data) & data!=2] <- NA
if(type=="f2" || (type=="bc" && chr_type=="X"))
points(x,ind,pch=21,col="black", bg=color[2],cex=cex)
else points(x,ind,pch=21,col="black", bg=color[3],cex=cex)
if(type=="f2" || (type=="bc" && chr_type=="X")) {
# BB genotypes
ind <- tind; ind[!is.na(data) & data!=3] <- NA
points(x,ind,pch=21,col="black", bg=color[3],cex=cex)
}
if(type=="f2") {
# not BB (D in mapmaker/qtl) genotypes
ind <- tind; ind[!is.na(data) & data!=4] <- NA
points(x,ind,pch=21,col="black", bg=color[4],cex=cex)
# not AA (C in mapmaker/qtl) genotypes
ind <- tind; ind[!is.na(data) & data!=5] <- NA
points(x,ind,pch=21,col="black", bg=color[5],cex=cex)
}
# plot map
u <- par("usr")
segments(map,u[3],map,u[3]-1/2)
segments(map,u[4],map,u[4]+1/2)
if(any(errors != 0)) {
ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
points(x,ind,pch=0,col=color[6],cex=cex+0.4,lwd=2)
}
if(include.xo) points(xoloc$loc,xoloc$ind,pch=4,col="blue",lwd=2)
}
else {
if(is.null(theylim)) theylim <- c(max(map), 0)
if(is.null(thexlim)) thexlim <- c(0.5,n.ind+0.5)
plot(0,0,type="n",ylab="Location (cM)",xlab="Individual",
main=themain,
xlim=thexlim,ylim=theylim, xaxt="n")
segments(1:n.ind,0,1:n.ind,max(map))
if(use.id) axis(side=1, at=1:n.ind, labels=id)
else axis(side=1)
# AA genotypes
tind <- rep(1:n.ind,length(map));tind[is.na(data)] <- NA
ind <- tind; ind[!is.na(data) & data!=1] <- NA
y <- rep(map,rep(n.ind,length(map)))
points(ind,y,pch=21,col="black", bg="white",cex=cex)
# AB genotypes
ind <- tind; ind[!is.na(data) & data!=2] <- NA
if(type=="f2" || (type=="bc" && chr_type=="X"))
points(ind,y,pch=21,col="black", bg=color[2],cex=cex)
else points(ind,y,pch=21,col="black", bg=color[3],cex=cex)
if(type=="f2" || (type=="bc" && chr_type=="X")) {
# BB genotypes
ind <- tind; ind[!is.na(data) & data!=3] <- NA
points(ind,y,pch=21,col="black", bg=color[3],cex=cex)
}
if(type=="f2") {
# not BB genotypes
ind <- tind; ind[!is.na(data) & data!=4] <- NA
points(ind,y,pch=21,col="black", bg=color[4],cex=cex)
# not AA genotypes
ind <- tind; ind[!is.na(data) & data!=5] <- NA
points(ind,y,pch=21,col="black", bg=color[5],cex=cex)
}
# plot map
u <- par("usr")
segments(u[1],map,(u[1]+1)/2,map)
segments(u[2],map,(n.ind+u[2])/2,map)
if(any(errors != 0)) {
ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
points(ind,y,pch=0,col=color[6],cex=cex+0.4,lwd=2)
}
if(include.xo) points(xoloc$ind,xoloc$loc,pch=4,col="blue",lwd=2)
}
}
invisible()
}
######################################################################
#
# plotInfo: Plot the proportion of missing information in the
# genotype data.
#
######################################################################
plotInfo <-
function(x,chr,method=c("entropy","variance","both"), step=1,
off.end=0, error.prob=0.001,
map.function=c("haldane","kosambi","c-f","morgan"),
alternate.chrid=FALSE, fourwaycross=c("all", "AB", "CD"),
include.genofreq=FALSE, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
method <- match(match.arg(method),c("entropy","variance","both"))-1
map.function <- match.arg(map.function)
if(!missing(chr)) cross <- subset(cross,chr=chr)
n.chr <- nchr(cross)
results <- NULL
cross <- calc.genoprob(cross, step=step, error.prob=error.prob,
off.end=off.end, map.function=map.function)
gap <- 25
# 4-way cross: can consider just A/B or just C/D
fourwaycross <- match.arg(fourwaycross)
n.ind <- nind(cross)
if(include.genofreq) theprob <- NULL
for(i in 1:n.chr) {
n.gen <- dim(cross$geno[[i]]$prob)[3]
n.pos <- ncol(cross$geno[[i]]$prob)
# 4-way cross: can consider just A/B or just C/D
if(n.gen==4 && (fourwaycross!="all")) {
att <- attributes(cross$geno[[i]]$prob)
if(fourwaycross == "AB") {
cross$geno[[i]]$prob[,,1] <- cross$geno[[i]]$prob[,,1] + cross$geno[[i]]$prob[,,3]
cross$geno[[i]]$prob[,,2] <- cross$geno[[i]]$prob[,,2] + cross$geno[[i]]$prob[,,4]
}
else {
cross$geno[[i]]$prob[,,1] <- cross$geno[[i]]$prob[,,1] + cross$geno[[i]]$prob[,,2]
cross$geno[[i]]$prob[,,2] <- cross$geno[[i]]$prob[,,3] + cross$geno[[i]]$prob[,,4]
}
cross$geno[[i]]$prob <- cross$geno[[i]]$prob[,,1:2]
n.gen <- 2
for(j in seq(along=att))
if(names(att)[j] != "dim" && names(att)[j] != "dimnames")
attr(cross$geno[[i]]$prob, names(att)[j]) <- att[[j]]
}
# calculate information (between 0 and 1)
info <- .C("R_info",
as.integer(n.ind),
as.integer(n.pos),
as.integer(n.gen),
as.double(cross$geno[[i]]$prob),
info1=as.double(rep(0,n.pos)),
info2=as.double(rep(0,n.pos)),
as.integer(method),
PACKAGE="qtl")
if(method != 1) { # rescale entropy version
if(n.gen==3) maxent <- 1.5*log(2)
else maxent <- log(n.gen)
info$info1 <- -info$info1/maxent
}
if(method != 0) { # rescale variance version
maxvar <- c(0.25,0.5,1.25)[n.gen-1]
info$info2 <- info$info2/maxvar
}
# reconstruct map
if("map" %in% names(attributes(cross$geno[[i]]$prob)))
map <- attr(cross$geno[[i]]$prob,"map")
else {
stp <- attr(cross$geno[[i]]$prob, "step")
oe <- attr(cross$geno[[i]]$prob, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[i]]$prob)))
stpw <- attr(cross$geno[[i]]$prob, "stepwidth")
else stpw <- "fixed"
map <- create.map(cross$geno[[i]]$map,stp,oe,stpw)
}
if(is.matrix(map)) map <- map[1,]
z <- data.frame(chr=rep(names(cross$geno)[i],length(map)),
pos=as.numeric(map),
"Missing information"=info$info1,
"Missing information"=info$info2, stringsAsFactors=TRUE)
w <- names(map)
o <- grep("^loc-*[0-9]+",w)
if(length(o) > 0) # inter-marker locations cited as "c*.loc*"
w[o] <- paste("c",names(cross$geno)[i],".",w[o],sep="")
rownames(z) <- w
results <- rbind(results, z)
if(include.genofreq) {
p <- cross$geno[[i]]$prob
if(chrtype(cross$geno[[i]])=="X")
p <- reviseXdata(crosstype(cross), expandX="full", sexpgm=getsex(cross), prob=p,
cross.attr=attributes(cross))
p <- apply(p, 2:3, mean, na.rm=TRUE)
if(is.null(theprob))
theprob <- p
else {
m1 <- match(colnames(p), colnames(theprob))
m2 <- match(colnames(theprob), colnames(p))
if(!any(is.na(m1)) && !any(is.na(m2))) {
theprob <- rbind(theprob, p[,colnames(theprob)])
}
else { # some differences in column names
if(all(!is.na(m1))) { # need to add some new columns to p
noldcol <- ncol(p)
nnewcol <- sum(is.na(m2))
for(i in 1:nnewcol)
p <- cbind(p, 0)
colnames(p)[-(1:noldcol)] <- colnames(theprob)[is.na(m2)]
theprob <- rbind(theprob, p[,colnames(theprob)])
}
else if(all(!is.na(m2))) { # need to add some new columns to theprob
noldcol <- ncol(theprob)
nnewcol <- sum(is.na(m1))
for(i in 1:nnewcol)
theprob <- cbind(theprob, 0)
colnames(theprob)[-(1:noldcol)] <- colnames(p)[is.na(m1)]
theprob <- rbind(theprob, p[,colnames(theprob)])
}
else { # need to add some new columns to each
noldcol <- ncol(theprob)
nnewcol <- sum(is.na(m1))
oldcolnam <- colnames(theprob)
for(i in 1:nnewcol)
theprob <- cbind(theprob, 0)
colnames(theprob)[-(1:noldcol)] <- colnames(p)[is.na(m1)]
noldcol <- ncol(p)
nnewcol <- sum(is.na(m2))
for(i in 1:nnewcol)
p <- cbind(p, 0)
colnames(p)[-(1:noldcol)] <- oldcolnam[is.na(m2)]
theprob <- rbind(theprob, p[,colnames(theprob)])
}
}
}
} # end of if(include.genofreq)
}
colnames(results)[3:4] <- c("misinfo.entropy","misinfo.variance")
if(method==0) results <- results[,-4]
if(method==1) results <- results[,-3]
if(include.genofreq)
results <- cbind(results, theprob)
class(results) <- c("scanone","data.frame")
# check whether main was included as an argument
args <- list(...)
if("main" %in% names(args)) hasmain <- TRUE
else hasmain <- FALSE
# check whether gap was included as an argument
if(!("gap" %in% names(args))) {
if(method==0) {
if(hasmain)
plot.scanone(results,ylim=c(0,1),gap=gap,
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,ylim=c(0,1),gap=gap,
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
else if(method==1) {
if(hasmain)
plot.scanone(results,ylim=c(0,1),gap=gap,
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,ylim=c(0,1),gap=gap,
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
else if(method==2) {
if(hasmain)
plot.scanone(results,results,lodcolumn=1:2,ylim=c(0,1),gap=gap,
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,results,lodcolumn=1:2,ylim=c(0,1),gap=gap,
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
}
else { # gap was included in ...
if(method==0) {
if(hasmain)
plot.scanone(results,ylim=c(0,1),
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,ylim=c(0,1),
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
else if(method==1) {
if(hasmain)
plot.scanone(results,ylim=c(0,1),
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,ylim=c(0,1),
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
else if(method==2) {
if(hasmain)
plot.scanone(results,results,lodcolumn=1:2,ylim=c(0,1),
alternate.chrid=alternate.chrid,...)
else
plot.scanone(results,results,lodcolumn=1:2,ylim=c(0,1),
main="Missing information",
alternate.chrid=alternate.chrid,...)
}
}
invisible(results)
}
# plot phenotypes against one or more markers
plotPXG <-
function(x, marker, pheno.col = 1, jitter = 1, infer = TRUE,
pch, ylab, main, col, ...)
{
cross <- x
if(!inherits(cross, "cross"))
stop("Input should have class \"cross\".")
type <- crosstype(cross)
if(LikePheVector(pheno.col, nind(cross), nphe(cross))) {
cross$pheno <- cbind(pheno.col, cross$pheno)
pheno.col <- 1
}
if(length(pheno.col) > 1) {
pheno.col <- pheno.col[1]
warning("plotPXG can take just one phenotype; only the first will be used")
}
if(is.character(pheno.col)) {
num <- find.pheno(cross, pheno.col)
if(is.na(num))
stop("Couldn't identify phenotype \"", pheno.col, "\"")
pheno.col <- num
}
if(pheno.col < 1 | pheno.col > nphe(cross))
stop("pheno.col values should be between 1 and the no. phenotypes")
if(!is.numeric(cross$pheno[,pheno.col]))
stop("phenotype \"", colnames(cross$pheno)[pheno.col], "\" is not numeric.")
if(missing(pch)) pch <- par("pch")
if(missing(ylab)) ylab <- colnames(cross$pheno)[pheno.col]
oldlas <- par("las")
on.exit(par(las = oldlas))
par(las = 1)
# find chromosomes containing the markers
o <- sapply(cross$geno, function(a, b) b %in% colnames(a$data),
marker)
if(length(marker)==1) o <- matrix(o,nrow=1)
if(!all(apply(o,1,any))) {
oo <- apply(o,1,any)
stop("Marker ", marker[!oo], " not found")
}
n.mark <- length(marker)
o <- apply(o, 1, which)
chr <- names(cross$geno)[o]
uchr <- unique(chr)
cross <- subset(cross, chr=uchr)
map <- pull.map(cross)
pos <- NULL
for(i in seq(length(chr))) pos[i] <- map[[chr[i]]][marker[i]]
chr_type <- sapply(cross$geno, chrtype)
names(chr_type) <- names(cross$geno)
chr_type <- chr_type[chr]
# if X chromosome and backcross or intercross, get sex/direction data
if(any(chr_type == "X") && (type == "bc" || type == "f2"))
sexpgm <- getsex(cross)
else sexpgm <- NULL
# number of possible genotypes
gen.names <- list()
for(i in seq(length(chr)))
gen.names[[i]] <- getgenonames(type, chr_type[i], "full", sexpgm, attributes(cross))
n.gen <- sapply(gen.names, length)
jitter <- jitter/10
if(any(n.gen == 2)) jitter <- jitter * 0.75
# function to determine whether genotype is fully known
tempf <-
function(x, type)
{
tmp <- is.na(x)
if(type=="f2") tmp[!is.na(x) & x>3] <- TRUE
if(type=="4way") tmp[!is.na(x) & x>4] <- TRUE
tmp
}
# if infer=TRUE, fill in genotype data by a single imputation
if(infer) {
which.missing <- tempf(cross$geno[[chr[1]]]$data[, marker[1]],type)
if(n.mark > 1)
for(i in 2:n.mark)
which.missing <- which.missing | tempf(cross$geno[[chr[i]]]$data[,marker[i]],type)
which.missing <- as.numeric(which.missing)
cross <- fill.geno(cross, method = "imp")
}
else which.missing <- rep(1,nind(cross))
# data to plot
x <- cross$geno[[chr[1]]]$data[, marker[1]]
if(n.mark > 1)
for(i in 2:n.mark)
x <- cbind(x, cross$geno[[chr[i]]]$data[, marker[i]])
else x <- as.matrix(x)
y <- cross$pheno[, pheno.col]
if(!infer) { # replace partially informative genotypes with NAs
if(type == "f2") x[x > 3] <- NA
if(type == "4way") x[x > 4] <- NA
if(sum(!is.na(x)) == 0)
stop("Can't use infer=FALSE as there are no fully informative genotypes")
}
# in case of X chromosome, recode some genotypes
if(any(chr_type == "X") && (type == "bc" || type == "f2")) {
ix = seq(n.mark)[chr_type == "X"]
for(i in ix)
x[, i] <- as.numeric(reviseXdata(type, "full", sexpgm,
geno = as.matrix(x[, i]),
cross.attr=attributes(cross)))
}
# save all of the data, returned invisibly
data <- as.data.frame(x, stringsAsFactors=TRUE)
names(data) <- marker
for(i in marker) data[[i]] <- ordered(data[[i]])
data$pheno <- y
data$inferred <- which.missing
# re-code the multi-marker genotypes
if(n.mark > 1) {
for(i in 2:n.mark)
x[, 1] <- n.gen[i] * (x[, 1] - 1) + x[, i]
}
x <- x[, 1]
observed <- sort(unique(x))
# amount of jitter
u <- runif(nind(cross), -jitter, jitter)
r <- (1 - 2 * jitter)/2
# genotype names
if(n.mark == 1)
gnames <- gen.names[[1]]
else {
gnames <- array(gen.names[[n.mark]], c(prod(n.gen), n.mark))
for(i in (n.mark - 1):1) {
tmpi <- rep(gen.names[[i]], rep(prod(n.gen[(i + 1):n.mark]),
n.gen[i]))
if(i > 1)
tmpi <- rep(tmpi, prod(n.gen[1:(i - 1)]))
gnames[, i] <- tmpi
}
gnames <- apply(gnames, 1, function(x) paste(x, collapse = "\n"))
}
# create plot
plot(x + u, y, xlab = "Genotype", ylab = ylab, type = "n",
main = "", xlim = c(1 - r + jitter, length(gnames) + r + jitter), xaxt = "n",
...)
# marker names at top
if(missing(main))
mtext(paste(marker, collapse = "\n"),
line=0.5, cex = ifelse(n.mark==1, 1.2, 0.8))
else
title(main=main)
abline(v = 1:prod(n.gen), col = "gray", lty = 3)
if(length(pch) == 1)
pch = rep(pch, length(x))
if(infer) {
points((x + u)[which.missing == 1], y[which.missing ==
1], col = "red", pch = pch[which.missing == 1])
points((x + u)[which.missing == 0], y[which.missing ==
0], pch = pch[which.missing == 0])
}
else points(x + u, y, pch = pch)
sux = sort(unique(x))
# add confidence intervals
me <- se <- array(NA, length(gnames))
me[sux] <- tapply(y, x, mean, na.rm = TRUE)
se[sux] <- tapply(y, x, function(a) sd(a, na.rm = TRUE)/sqrt(sum(!is.na(a))))
thecolors <- c("black", "blue", "red", "purple", "green", "orange")
if(missing(col)) {
col <- thecolors[1:n.gen[n.mark]]
if(n.gen[n.mark] == 3)
col <- c("blue", "purple", "red")
else if(n.gen[n.mark] == 2)
col <- c("blue", "red")
}
ng <- length(gnames)
segments(seq(ng) + jitter * 2, me, seq(ng) +
jitter * 4, me, lwd = 2, col = col)
segments(seq(ng) + jitter * 3, me - se, seq(ng) +
jitter * 3, me + se, lwd = 2, col = col)
segments(seq(ng) + jitter * 2.5, me - se, seq(ng) +
jitter * 3.5, me - se, lwd = 2, col = col)
segments(seq(ng) + jitter * 2.5, me + se, seq(ng) +
jitter * 3.5, me + se, lwd = 2, col = col)
# add genotypes below
u <- par("usr")
cxaxis <- par("cex.axis")
segments(seq(along=gnames), u[3], seq(along=gnames), u[3] - diff(u[3:4]) *
0.015, xpd = TRUE)
axis(side=1, at=seq(along=gnames), labels=gnames,
cex=ifelse(n.mark==1, cxaxis, cxaxis*0.8),
tick=FALSE, line = (length(marker)-1)/2)
invisible(data)
}
plotPheno <-
function(x, pheno.col=1, ...)
{
if(!inherits(x, "cross"))
stop("Input should have class \"cross\".")
if(LikePheVector(pheno.col, nind(x), nphe(x))) {
x$pheno <- cbind(pheno.col, x$pheno)
pheno.col <- 1
}
if(length(pheno.col) > 1) {
pheno.col <- pheno.col[1]
warning("Ignoring all but the first element in pheno.col.")
}
if(is.character(pheno.col)) {
num <- find.pheno(x, pheno.col)
if(is.na(num))
stop("Couldn't identify phenotype \"", pheno.col, "\"")
pheno.col <- num
}
if(pheno.col < 1 | pheno.col > nphe(x))
warning("pheno.col should be between 1 and ", nphe(x))
phe <- x$pheno[,pheno.col]
u <- length(unique(phe))
if(u==2 || (u < 10 && nind(x) > 50))
phe <- as.factor(phe)
plot_pheno_sub <-
function(phe, xlab=paste("phe", pheno.col),
main=colnames(x$pheno)[pheno.col], col="white",
breaks=ceiling(2*sqrt(nind(x))),
las=1,
...)
{
if(is.factor(phe)) {
barplot(table(phe), xlab=xlab, main=main, col=col, las=las, ...)
}
else {
phe <- as.numeric(phe)[1:nind(x)]
hist(phe, breaks = breaks,
xlab = xlab, main = main, las=las, ...)
}
}
plot_pheno_sub(phe, ...)
}
# end of plot.R
Try the qtl package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.