#############################################################################################################
## Legacy routines to be purged.
#############################################################################################################
hotspot.scan <- function(cross, highobj, lod.thr = NULL, quant.level = NULL, window = NULL,
verbose = FALSE)
{
hotsize(highobj, lod.thr, window, quant.level)
}
#############################################################################################################
hotspot.plot <- function(hotobj, quant.thr = NULL, maps = NULL, main = "")
{
plot(hotobj, ..., by.chr = TRUE, quant.thr = quant.thr, maps = maps, title = main)
}
#############################################################################################################
##############################################################################
OGetCisCandReg <- function(cand.reg, scan, lod.thr, drop = 1.5) {
all.trait.nms <- cand.reg[, 1]
cand.reg <- cand.reg[cand.reg[, 2] == cand.reg[, 4],]
n <- nrow(cand.reg)
trait.nms <- names(scan)[-c(1, 2)]
is.cis <- rep(FALSE, n)
peak.pos.lower <- rep(NA, n)
peak.pos.upper <- rep(NA, n)
for (i in 1 : n) {
phys.pos <- cand.reg[i, 3]
trait.index <- which(trait.nms == cand.reg[i, 1]) + 2
sscan <- scan[scan[, 1] == cand.reg[i, 4], c(1, 2, trait.index)]
lod.interval <- lodint(sscan, chr = cand.reg[i, 4], drop = drop)
peak.pos.lower[i] <- lod.interval[1, 2]
peak.pos.upper[i] <- lod.interval[nrow(lod.interval), 2]
lod.phys.pos <- sscan[which.min(abs(sscan[, 2] - cand.reg[i, 3])), 3]
if (phys.pos >= peak.pos.lower[i] & phys.pos <= peak.pos.upper[i] &
lod.phys.pos >= lod.thr) {
is.cis[i] <- TRUE
}
}
out <- data.frame(cand.reg, peak.pos.lower, peak.pos.upper)
index <- match(out[is.cis, 1], all.trait.nms)
list(cis.reg = out[is.cis,], cis.index = index)
}
##############################################################################
OGetCandReg <- function(scan, annot, traits, lod.thr, drop = 1.5,
nms = names(scan)[-c(1,2)])
{
## want to use highlod instead of scan below.
## need to decode highlod.
## currently this only gets max over genome; want max by chr, yes?
traits <- unique(traits)
n <- length(traits)
out <- data.frame(matrix(NA, n, 6))
names(out) <- c("gene", "phys.chr", "phys.pos", "peak.chr", "peak.pos",
"peak.lod")
for (i in 1:n) {
ii <- match(traits[i], annot[, 1])
out[i, 1:3] <- annot[ii, c(1, 3, 5)]
trait.chr <- annot[ii, 3]
trait.pos <- annot[ii, 5]
if (!is.na(trait.pos)) {
peak.index <- which.max(scan[, traits[i]])
peak.lod <- scan[peak.index, traits[i]]
peak.chr <- scan[peak.index, 1]
peak.pos <- scan[peak.index, 2]
if (peak.lod >= lod.thr) {
out[i, 4] <- peak.chr
out[i, 5] <- peak.pos
out[i, 6] <- peak.lod
}
}
cat(" ", i, "\n")
}
subset(out, !is.na(out[, 4]))
}
##############################################################################
OGetCoMappingTraits <- function(traits, scan, lod.thr, drop = 1.5) {
n <- nrow(traits)
out <- vector(mode = "list", length = n)
names(out) <- traits[, 1]
nms1 <- names(scan)[-c(1, 2)]
for (i in 1:n) {
peak.pos <- traits[i, 5]
sub.scan <- scan[scan[, 1] == traits[i, 4], ]
peak.lods <- apply(sub.scan[, -c(1, 2)], 2, max)
aux <- which(peak.lods >= lod.thr)
aux <- aux[-match(traits[i, 1], names(aux))]
nms2 <- names(aux)
is.target <- rep(FALSE, length(nms2))
for (j in 1:length(nms2)) {
trait.index <- match(nms2[j], nms1)
sscan <- sub.scan[, c(1, 2, trait.index + 2)]
lod.interval <- lodint(sscan, chr = traits[i, 4], drop = drop)
lb <- lod.interval[1, 2]
ub <- lod.interval[nrow(lod.interval), 2] ## was 3; error if there are tied LOD scores.
if(peak.pos >= lb & peak.pos <= ub)
is.target[j] <- TRUE
}
out[[i]] <- nms2[is.target]
cat(" ", i, "\n")
}
out
}
##############################################################################
## Run permuations
lod.quantile.permutations.2 <- function(cross, pheno.col, N, n.perm, lod.thr,
batch.effect, window, seed = 123456789, verbose = FALSE)
{
## Set up pseudo-random number generation seeds.
set.seed(seed[[1]])
all.seeds <- sample(c(98765:987654), n.perm, replace=FALSE)
## Set up matrices to record values.
nphe.cross <- nphe(cross)
N <- N[N <= nphe.cross]
quants <- 1 - (N - 1)/nphe.cross
l.N <- length(N)
max.lod.quant <- matrix(NA, n.perm, l.N)
dimnames(max.lod.quant)[[2]] <- paste(paste(round(quants,4)*100, "%", sep=""),
N, sep="_")
l.lod.thr <- length(lod.thr)
max.N <- matrix(NA, n.perm, l.lod.thr)
max.lod.quant <- matrix(NA, n.perm, l.N)
max.N.window <- matrix(NA, n.perm, l.lod.thr)
if(!is.null(batch.effect))
cross <- subset(cross, ind = !is.na(cross$pheno[[batch.effect]]))
covars <- sexbatch.covar(cross, batch.effect)
n.ind <- nind(cross)
for(i in 1:n.perm){
perm.cross <- cross
set.seed(all.seeds[i])
tmp <- sample(c(1:n.ind), n.ind, replace=FALSE)
perm.cross$pheno <- cross$pheno[tmp,]
per.scan <- scanone(perm.cross, pheno.col=pheno.col, method="hk",
addcovar=covars$addcovar, intcovar=covars$intcovar)
## Elias' quantiles.
quant <- apply(per.scan[,-(1:2)], 1, quantile, quants)
max.lod.quant[i,] <- apply(quant,1,max)
## Jansen's count.
max.N[i,] <- apply(count.thr(per.scan, lod.thr, droptwo = TRUE), 1, sum)
## Smoothed count.
maxlod <- apply(per.scan[,-(1:2)], 2, tapply, per.scan[,1], max)
chrs <- dimnames(maxlod)[[1]]
chr <- factor(per.scan[,1], levels = chrs)
pos <- as.numeric(per.scan[,2])
maxlod.pos <- maxlod.thr.pos <- vector("list", length(chrs))
names(maxlod.pos) <- names(maxlod.thr.pos) <- chrs
for(k in seq(length(chrs))) {
scan.out.bychr <- per.scan[per.scan[,1] == chrs[k], ]
maxlod.pos[[k]] <- apply(scan.out.bychr[,-(1:2)], 2, function(a,b)
mean(b[!is.nan(a) & a==max(a, na.rm=TRUE)]), scan.out.bychr[,2])
}
for(j in seq(along = lod.thr)){
for(k in chrs)
maxlod.thr.pos[[k]] <- maxlod.pos[[k]][maxlod[k,] >= lod.thr[j]]
neqtl.pos <- smoothall(maxlod.thr.pos, thechr = chr, thepos = pos, window = window)
max.N.window[i,j] <- max(neqtl.pos[,3])
}
print(i)
}
list(max.lod.quant=max.lod.quant,
max.N=max.N,
max.N.window=max.N.window)
}
##############################################################################
OWW.perm <- function(scanmat, lod.thrs, n.perm, verbose = FALSE)
{
## Permute all but first column.
nlod <- length(lod.thrs)
max.WW <- matrix(NA, n.perm, nlod)
dimnames(max.WW) <- list(NULL, as.character(lod.thrs))
for(i in 1:n.perm) {
if(verbose)
cat("\n", i)
## Separately permute columns of the scan object separately by trait.
mycat("sample...", verbose, last = "")
scans <- apply(scanmat, 2, sample)
## Original code had reshuffle of same values. Better to start over each time.
##scanmat <- apply(scanmat, 2, sample)
## Get the maximum spurious hotspot size (N-method) across genome
## for different QTL mapping significance levels.
mycat("sum...", verbose, last = "")
max.WW[i,] <- apply(count.thr(scans, lod.thrs, FALSE), 1, max)
}
max.WW
}
##############################################################################
ONL.N.perm <- function(cross, Nmax, n.perm, lod.thrs, drop = 1.5,
verbose = FALSE, init.seed = 0)
{
set.seed(init.seed)
n.phe <- nphe(cross)
n.ind <- nind(cross)
Nmax <- min(Nmax, n.phe)
Ns <- seq(Nmax)
quants <- 1 - (Ns - 1) / n.phe
n.lod <- length(lod.thrs)
max.N <- matrix(NA, n.perm, n.lod)
dimnames(max.N) <- list(NULL, as.character(lod.thrs))
max.lod.quant <- matrix(NA, n.perm, Nmax)
dimnames(max.lod.quant) <- list(NULL, as.character(Ns))
for(i in 1:n.perm){
if(verbose)
cat("\n", i)
## permute rows of the phenotype data matrix
perm.cross <- cross
tmp <- sample(c(1:n.ind), n.ind, replace=FALSE)
perm.cross$pheno <- cross$pheno[tmp,]
## perform mapping analysis in the permuted data
mycat("scanone...", verbose, last = "")
## NB: scanone groups phenos in batches based on missing data patterns.
scanmat <- scanone(perm.cross, pheno.col = c(1:n.phe), method = "hk")
chr <- scanmat[, 1]
scanmat <- as.matrix(scanmat[, -(1:2), drop = FALSE])
## apply LOD drop interval to per.scan
mycat("LOD drop...", verbose, last = "")
scanmat <- set.to.zero.beyond.drop.int(chr, scanmat, min(lod.thrs), drop)
## get lod quantiles at each genomic location
## rows indexes the quants associated with the Ns
## columns indexes the genomic positions
mycat("quantile...", verbose, last = "")
quant <- apply(scanmat, 1, quantile, quants)
## get the maximum lod-quantile across the genome
## rows indexes the permutations
## columns indexes the Ns
mycat("quant...", verbose, last = "")
max.lod.quant[i,] <- apply(quant, 1, max)
## Get the maximum spurious hotspot size (N-method) across genome
## for different QTL mapping significance levels.
mycat("max...", verbose, last = "")
max.N[i,] <- apply(count.thr(scanmat, lod.thrs, FALSE), 1, max)
}
list(max.lod.quant = max.lod.quant, max.N = max.N)
}
######################################################################
## This function takes a scanone object (scan) and for each trait and each
## chromosome it: (1) finds out the max LOD peak; (2) if the max LOD value
## is smaller than the map threshold (thr) it set to zero all LOD values;
## (3) if the max LOD value is higher than thr, it computes the LOD drop
## interval around the peak and set to zero all LOD values outside the LOD
## LOD interval.
##
set.to.zero.beyond.drop.int <- function(chr, scanmat, thr, drop = 1.5)
{
mylodint <- function(x, drop = 1.5)
{
drops <- x < (max(x) - drop)
if(any(drops))
x[drops] <- 0
x
}
mychr <- levels(chr)
for(i in mychr) {
pointer <- which(chr == i)
if(length(pointer)) {
if(max(scanmat[pointer,]) < thr)
## Zero out if no peak above threshold.
scanmat[pointer,] <- 0
else
## Zero out below drop from peak.
scanmat[pointer,] <-
apply(scanmat[pointer,, drop = FALSE], 2, mylodint, drop)
}
}
scanmat
}
######################################################################
qtlhot.scan <- function(...) quant.scanone(...)
quant.plot <- function(...)
{
x <- quant.lod(...)
thr.level <- x$thr.level
lod.thr <- x$lod.thr
n.probs <- length(probs)
tmp.plot <- function(x.vals, quant, x.crit, probs, level, is.quantile = FALSE, main = "",
add.level = FALSE)
{
n.probs <- length(probs)
thr.level <- min(which(probs >= level))
quant.thr <- rev(quant[thr.level,])[thr.level]
xlabs <- "single trait LOD threshold"
if(is.quantile)
xlabs <- paste(xlabs, "quantile")
plot(range(x.vals), c(0, max(quant)), type = "n", xlab = "", ylab = "")
mtext(xlabs, 1, 2)
mtext("hotspot size", 2, 2)
abline(v = x.crit, col = "darkgray", lty = 2)
abline(h = quant.thr, col = "darkgray", lty = 2)
mtext(ceiling(quant.thr), 2, at = quant.thr, las = 2, cex = 0.5)
for(i in seq(along = probs)) {
lines(rev(sort(x.vals)), quant[i,],
lwd = 1 + 2 * (round(probs[i] - level, 2) == 0))
}
text(x.crit, rev(quant[n.probs,])[thr.level] + 5, 1 - max(probs), adj = 0)
text(x.crit, rev(quant[1,])[thr.level] - 5, 1 - min(probs), adj = 1)
text(x.vals[n.probs], quant.thr + 5, 1 - level, adj = 1)
if(add.level)
main <- paste(main, "\n hotspot size significance level =", 1 - max(probs), "to", 1 - min(probs))
mtext(main, 3, 0.5)
quant.thr
}
tmpar <- par(mfrow = c(2,2), mar = c(3.1,3.1,2.1,0.1))
## Jansen method, smoothing.
x$quant.thr <- tmp.plot(lod.thrs, x$quant.N.window, lod.thr, probs, level, FALSE,
"Jansen method 5cM window")
tmp.plot(probs, x$quant.N.window, level, probs, level, TRUE,
"Jansen method 5cM window")
tmp.plot(lod.thrs, x$quant.N, lod.thr, probs, level, FALSE,
"Jansen method per locus")
tmp.plot(probs, x$quant.N, level, probs, level, TRUE,
"Jansen method per locus")
par(tmpar)
## Chaibub Neto method.
plot(c(1,x$max.quant), c(min(x$quant, na.rm = TRUE), max(x$quant, na.rm = TRUE)), type = "n",
xlab = "significant hotspot size with given threshold",
ylab = "hotspot LOD score threshold",
log = "xy")
abline(h = lod.thr, col = "darkgray", lty = 2)
abline(v = x$quant.thr, col = "darkgray", lty = 2)
mtext(ceiling(x$quant.thr), 1, at = x$quant.thr, las = 2)
for(i in seq(along = probs)) {
tmp <- (x$quant[i,seq(x$max.quant)] >= lod.thr)
if(any(tmp))
lines(seq(x$max.quant)[tmp], x$quant[i,seq(x$max.quant)][tmp],
lwd = 1 + 2 * (round(probs[i] - 0.95, 2) == 0))
if(any(!tmp))
lines(seq(x$max.quant)[!tmp], x$quant[i,seq(x$max.quant)][!tmp],
lwd = 1 + 2 * (round(probs[i] - 0.95, 2) == 0),
col = "darkgray")
}
n.thr2 <- length(lod.thrs) / 2
text(n.thr2 + 1, x$quant[n.probs, n.thr2], 1 - max(probs), adj = 0)
text(n.thr2 - 1, x$quant[1,n.thr2], 1 - min(probs), adj = 1)
title(paste("hotspot LOD threshold by hotspot size\nsignificance level =",
1 - max(probs), "to", 1 - min(probs)))
invisible(x)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.