######################################################################
#
# effectscan.R
#
# copyright (c) 2003-2020, Karl W. Broman
# [completely re-written in Sep, 2007, based partly on code from Hao Wu]
# last modified Feb, 2020
# first written Jan, 2003
#
# 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: effectscan
#
######################################################################
effectscan <-
function(cross, pheno.col=1, chr, get.se=FALSE, draw=TRUE,
gap=25, ylim, mtick=c("line","triangle"),
add.legend=TRUE, alternate.chrid=FALSE, ...)
{
type <- crosstype(cross)
mtick <- match.arg(mtick)
if(type == "4way")
stop("effect scan not working for 4-way cross yet.")
if(length(pheno.col) > 1) {
pheno.col <- pheno.col[1]
warning("effectscan 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.")
pheno <- cross$pheno[,pheno.col]
wh <- is.na(pheno)
if(any(wh)) {
pheno <- pheno[!wh]
cross <- subset.cross(cross, ind=(!wh))
}
if(!missing(chr)) cross <- subset.cross(cross, chr=chr)
chr_type <- sapply(cross$geno, chrtype)
n.ind <- length(pheno)
results <- NULL
for(i in 1:nchr(cross)) {
if(!("draws" %in% names(cross$geno[[i]])))
stop("You must first run sim.geno.")
draws <- cross$geno[[i]]$draws
# create map of positions
if("map" %in% names(attributes(cross$geno[[i]]$draws)))
map <- attr(cross$geno[[i]]$draws,"map")
else {
stp <- attr(cross$geno[[i]]$draws, "step")
oe <- attr(cross$geno[[i]]$draws, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[i]]$draws)))
stpw <- attr(cross$geno[[i]]$draws, "stepwidth")
else stpw <- "fixed"
map <- create.map(cross$geno[[i]]$map,stp,oe,stpw)
}
if(is.matrix(map)) {
marnam <- colnames(map)
map <- map[1,]
}
else marnam <- names(map)
if(type == "risib" || type=="riself" || type=="dh" || type=="haploid") {
mapping <- rbind(c(+1, -1),
c(+1, +1))
colnames(mapping) <- c("intercept","a")
dropcol <- 1
}
else if(type=="bc") {
if(chr_type[i] == "X") {
sexpgm <- getsex(cross)
draws <- reviseXdata(type, "full", sexpgm, draws=draws,
cross.attr=attributes(cross))
if(is.null(sexpgm$sex) || all(sexpgm$sex==0) || all(sexpgm$sex==1)) { # all one sex
mapping <- rbind(c(+1, -0.5),
c(+1, +0.5))
colnames(mapping) <- c("intercept", "a")
dropcol <- 1
}
else { # some of each sex
mapping <- rbind(c(+1, 0,-0.5, 0),
c(+1, 0,+0.5, 0),
c(+1,+1, 0, -0.5),
c(+1,+1, 0, +0.5))
colnames(mapping) <- c("intercept", "sex", "a.female", "a.male")
dropcol <- 1:2
}
} # end bc X chr
else {
mapping <- rbind(c(+1, -0.5),
c(+1, +0.5))
colnames(mapping) <- c("intercept", "a")
dropcol <- 1
} # end bc autosome
} # end bc
else { # intercross
if(chr_type[i] == "X") {
sexpgm <- getsex(cross)
draws <- reviseXdata(type, "full", sexpgm, draws=draws,
cross.attr=attributes(cross))
if(is.null(sexpgm$pgm) || all(sexpgm$pgm==0) || all(sexpgm$pgm==1)) { # all one direction
if(is.null(sexpgm$sex) || all(sexpgm$sex==0) || all(sexpgm$sex==1)) { # all one sex
mapping <- rbind(c(+1, -0.5),
c(+1, +0.5))
colnames(mapping) <- c("intercept", "a")
dropcol <- 1
}
else {
mapping <- rbind(c(+1, 0,-0.5, 0),
c(+1, 0,+0.5, 0),
c(+1,+1, 0, -0.5),
c(+1,+1, 0, +0.5))
colnames(mapping) <- c("intercept", "sex", "a.female", "a.male")
dropcol <- 1:2
}
}
else { # some of each direction
if(is.null(sexpgm$sex) || all(sexpgm$sex==0)) { # all female
mapping <- rbind(c(+1, 0,-0.5, 0),
c(+1, 0,+0.5, 0),
c(+1,+1, 0,-0.5),
c(+1,+1, 0,+0.5))
colnames(mapping) <- c("intercept","dir","a.forw","a.rev")
dropcol <- 1:2
}
else if(all(sexpgm$sex==1)) { # all male
mapping <- rbind(c(+1, -0.5),
c(+1, +0.5))
colnames(mapping) <- c("intercept", "a")
dropcol <- 1
}
else { # some of each sex
mapping <- rbind(c(+1, 0, 0, -0.5, 0, 0),
c(+1, 0, 0, +0.5, 0, 0),
c(+1,+1, 0, 0,-0.5, 0),
c(+1,+1, 0, 0,+0.5, 0),
c(+1, 0,+1, 0, 0,-0.5),
c(+1, 0,+1, 0, 0,+0.5))
colnames(mapping) <- c("intercept","dir","sex","a.femaleforw","a.femalerev","a.male")
dropcol <- 1:3
}
}
} # end f2 X chr
else {
mapping <- rbind(c(+1, -1, 0),
c(+1, 0, +1),
c(+1, +1, 0))
colnames(mapping) <- c("intercept","a","d")
dropcol <- 1
} # f2 autosome
} # end f2
n.gen <- ncol(mapping)
n.pos <- ncol(draws)
n.imp <- dim(draws)[3]
z <- .C("R_effectscan",
as.integer(n.ind),
as.integer(n.gen),
as.integer(n.imp),
as.integer(n.pos),
as.integer(draws-1),
as.double(pheno),
as.double(mapping),
beta=as.double(rep(0,n.pos*n.gen)),
se=as.double(rep(0,n.pos*n.gen)),
as.integer(get.se),
PACKAGE="qtl")
beta <- t(matrix(z$beta, ncol=n.pos))
colnames(beta) <- colnames(mapping)
if(get.se) {
se <- t(matrix(z$se, ncol=n.pos))
colnames(se) <- paste("se", colnames(mapping), sep=".")
beta <- cbind(beta, se[,-dropcol,drop=FALSE])
}
z <- beta[,-dropcol,drop=FALSE]
w <- marnam
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
z <- as.data.frame(z, stringsAsFactors=TRUE)
z <- cbind(chr=factor(rep(names(cross$geno)[i],length(map)),levels=names(cross$geno)),
pos=as.numeric(map), z)
rownames(z) <- w
if(i==1) results <- z
else {
w <- match(colnames(z), colnames(results))
if(any(is.na(w))) {
curnam <- colnames(results)
for(j in which(is.na(w)))
results <- cbind(results, rep(NA, nrow(results)))
colnames(results) <- c(curnam, colnames(z)[is.na(w)])
}
w <- match(colnames(results), colnames(z))
if(any(is.na(w))) {
curnam <- colnames(z)
for(j in which(is.na(w)))
z <- cbind(z, rep(NA, nrow(z)))
colnames(z) <- c(curnam, colnames(results)[is.na(w)])
}
results <- rbind(results, z)
}
} # end loop over chromosomes
class(results) <- c("effectscan", "scanone", "data.frame")
if(draw) { # make the figure
if(missing(ylim))
plot.effectscan(results, gap=gap, mtick=mtick, add.legend=add.legend,
alternate.chrid=alternate.chrid, ...)
else
plot.effectscan(results, gap=gap, mtick=mtick, add.legend=add.legend,
ylim=ylim, alternate.chrid=alternate.chrid, ...)
}
invisible(results)
}
# function to make the effectscan plot
plot.effectscan <-
function(x, gap=25, ylim, mtick=c("line","triangle"),
add.legend=TRUE, alternate.chrid=FALSE, ...)
{
col <- c("blue","red","darkorange","darkgreen","purple")
lightcol <- c("lightblue", "pink", "peachpuff1", "palegreen1", "thistle1")
results <- x
eff <- 3:ncol(results)
if(length(grep("^se", colnames(results)))>0) get.se <- TRUE
else get.se <- FALSE
if(get.se) {
se <- grep("^se", colnames(results)[eff])
eff <- eff[-se]
se <- se + 2
lo <- as.matrix(results[,eff]) - as.matrix(results[,se])
hi <- as.matrix(results[,eff]) + as.matrix(results[,se])
yl <- range(c(lo,hi), na.rm=TRUE)
}
else yl <- range(results[,eff], na.rm=TRUE)
if(!missing(ylim)) yl <- ylim
plot.scanone(results, lodcolumn=1, ylim=yl, gap=gap, mtick=mtick, alternate.chrid=alternate.chrid,
col=col[1], ...)
if(get.se) {
begend <- matrix(unlist(tapply(results[,2],results[,1],range)),ncol=2,byrow=TRUE)
rownames(begend) <- unique(results[,1])
chr <- unique(as.character(results[,1]))
begend <- begend[as.character(chr),,drop=FALSE]
len <- begend[,2]-begend[,1]
if(length(len)>1) start <- c(0,cumsum(len+gap))-c(begend[,1],0)
else start <- 0
x <- results[,2]
for(i in seq(along=chr))
x[results[,1]==chr[i]] <- results[results[,1]==chr[i],2]+start[i]
for(i in seq(along=chr)) {
wh <- results[,1]==chr[i]
for(j in 1:ncol(lo)) {
if(any(!is.na(lo[wh,j]))) {
xx <- c(x[wh], rev(x[wh]))
yy <- c(lo[wh,j], rev(hi[wh,j]))
polygon(xx, yy, col=lightcol[j], border=lightcol[j])
}
}
}
# go back and add lines at edges
for(i in seq(along=chr)) {
wh <- results[,1]==chr[i]
for(j in 1:ncol(lo)) {
if(any(!is.na(lo[wh,j]))) {
xx <- c(x[wh], rev(x[wh]))
yy <- c(lo[wh,j], rev(hi[wh,j]))
lines(xx, yy, col=lightcol[j])
}
}
}
plot.scanone(results, lodcolumn=1, add=TRUE, col=col[1])
}
if(length(eff) > 1) {
for(i in seq(along=eff)[-1])
plot.scanone(results, lodcolumn=eff[i]-2, gap=gap, add=TRUE, col=col[i])
}
if(add.legend)
legend("top", legend=names(results)[eff], col=col[1:length(eff)], lwd=2)
abline(h=0, lty=2)
}
# end of effectscan.R
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.