Nothing
######################################################################
#
# summary.scantwo.R
#
# copyright (c) 2001-2020, Karl W Broman, Hao Wu, and Brian Yandell
# last modified Feb, 2020
# first written Nov, 2001
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License,
# version 3, as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but without any warranty; without even the implied warranty of
# merchantability or fitness for a particular purpose. See the GNU
# General Public License, version 3, for more details.
#
# A copy of the GNU General Public License, version 3, is available
# at http://www.r-project.org/Licenses/GPL-3
#
# Part of the R/qtl package
# Contains: summary.scantwo, print.summary.scantwo,
# max.scantwo, clean.scantwo, print.scantwo, subset.scantwo
# summary.scantwoperm, print.summary.scantwoperm
# condense.scantwo, summary.scantwocondensed
# max.scantwocondensed, print.summary.addpair
# rbind.scantwoperm, c.scantwoperm, subset.scantwoperm
# [.scantwoperm
#
######################################################################
#
# summarize the result from scantwo
#
######################################################################
summary.scantwo <-
function(object, thresholds,
what=c("best", "full", "add", "int"),
perms, alphas, lodcolumn=1, pvalues=FALSE,
allpairs=TRUE, ...)
{
if(!inherits(object, "scantwo") &&
!inherits(object, "scantwocondensed"))
stop("Input should have class \"scantwo\".")
addpair <- attr(object, "addpair")
if(!is.null(addpair) && addpair) { # results from addpair() that need special treatment
if("lod.minus1" %in% names(attributes(object))) { # asymmetric formula
attr(object, "addpair") <- NULL
x <- summary.scantwo(object, allpairs=allpairs)
class(x) <- "data.frame"
x <- x[,c(1,2,8,9,10,3,4,5),drop=FALSE]
mlod.minus1 <- tapply(attr(object, "lod.minus1"), object$map$chr, max, na.rm=TRUE)
mlod.minus2 <- tapply(attr(object, "lod.minus2"), object$map$chr, max, na.rm=TRUE)
w <- which(x[,1] == x[,2])
for(i in seq(along=w))
if(x[w[i],5] < x[w[i],8]) x[w[i],3:5] <- x[w[i],c(7,6,8)]
if(any(x[,1] != x[,2])) {
w <- which(x[,1] != x[,2])
y <- x[w,,drop=FALSE]
for(i in 1:nrow(y))
y[,1:5] <- y[,c(2,1,7,6,8)]
neworder <- cbind(1:nrow(x), rep(NA, nrow(x)))
neworder[w,2] <- nrow(x) + 1:nrow(y)
neworder <- as.numeric(t(neworder))
x <- rbind(x, y)[neworder[!is.na(neworder)],,drop=FALSE]
}
x <- cbind(x[,1:5], lod.2v1b=rep(NA,nrow(x)), lod.2v1a=rep(NA,nrow(x)))
names(x)[5] <- "lod.2v0"
x[,6] <- x[,5] - mlod.minus1[as.character(x[,2])]
x[,7] <- x[,5] - mlod.minus2[as.character(x[,1])]
if(!missing(thresholds)) {
if(length(thresholds) > 2)
warning("Only the first two values in thresholds are used.")
if(length(thresholds) == 1) thresholds <- c(thresholds, 0)
x <- x[!is.na(x[,5]) & x[,5] >= thresholds[1] &
((!is.na(x[,6]) & x[,6] >= thresholds[2]) |
(!is.na(x[,7]) & x[,7] >= thresholds[2])),,drop=FALSE]
}
class(x) <- c("summary.addpair", "data.frame")
return(x)
}
else { # symmetric formula
attr(object, "addpair") <- NULL
if(missing(thresholds))
x <- summary.scantwo(object, allpairs=allpairs)
else
x <- summary.scantwo(object, thresholds=thresholds, allpairs=allpairs)
x <- x[,1:6]
colnames(x)[5:6] <- c("lod.2v0", "lod.2v1")
if(!missing(thresholds)) {
if(length(thresholds) > 2)
warning("Only the first two values in thresholds are used.")
if(length(thresholds) == 1) thresholds <- c(thresholds, 0)
x <- x[!is.na(x[,5]) & x[,5] >= thresholds[1] &
!is.na(x[,6]) & x[,6] >= thresholds[2],,drop=FALSE]
}
class(x) <- c("summary.addpair", "data.frame")
return(x)
}
}
what <- match.arg(what)
if(!missing(thresholds)) {
if(length(thresholds) != 5)
stop("If thresholds are given, there must be 5 of them.")
else if(!is.numeric(thresholds))
stop("thresholds should be a numeric vector")
}
if(!missing(alphas)) {
if(length(alphas) != 5) {
if(length(alphas)==1)
alphas <- rep(alphas, 5)
else
stop("If alphas are given, there must be 5 of them.")
}
else if(!is.numeric(alphas))
stop("alphas should be a numeric vector")
}
if(!missing(perms) && !inherits(perms, "scantwoperm"))
stop("perms must be in scantwoperm format.")
# subset object and permutations, if necessary
if(inherits(object, "scantwo")) {
d <- dim(object$lod)
if(length(d)==3) {
if(!missing(perms)) {
if("AA" %in% names(perms)) { # contains X-specific results
ncp <- sapply(perms$AA, ncol)
}
else {
ncp <- sapply(perms, ncol)
}
if(all(ncp==1)) onepermcol <- TRUE
else onepermcol <- FALSE
if(any(ncp != d[3])) {
if(onepermcol) {
if(lodcolumn > 1)
warning("Just one column of permutation results; assuming they apply to all LOD score columns.")
}
else
stop("perms have different numbers of columns as object input.\n")
}
}
if(lodcolumn < 1 || lodcolumn > d[3])
stop("lodcolumn must be between 1 and ", d[3])
object$lod <- object$lod[,,lodcolumn]
if(!missing(perms) && !onepermcol) {
if("AA" %in% names(perms)) {
for(i in seq(along=perms))
perms[[i]] <- lapply(perms[[i]], function(a, b) a[,b,drop=FALSE], lodcolumn)
}
else {
perms <- lapply(perms, function(a, b) a[,b,drop=FALSE], lodcolumn)
}
}
}
}
else { # condensed version
if(is.matrix(object$pos1.jnt)) {
d <- ncol(object$pos1.jnt)
if(!missing(perms)) {
ncp <- sapply(perms, ncol)
if(all(ncp==1)) onepermcol <- TRUE
else onepermcol <- FALSE
if(any(ncp != d[3])) {
if(onepermcol)
warning("Just one column of permutation results; reusing for all LOD score columns.")
else
stop("perms have different numbers of columns as object input.\n")
}
}
if(lodcolumn < 1 || lodcolumn > d)
stop("lodcolumn must be between 1 and ", d)
for(i in 3:length(object))
object[[i]] <- object[[i]][,lodcolumn]
if(!missing(perms) && !onepermcol)
perms <- lapply(perms, function(a, b) a[,b,drop=FALSE], lodcolumn)
}
}
# check input
if(missing(perms) && !missing(alphas))
stop("If alphas are to be used, permutation results must be provided.")
if(!missing(thresholds) && !missing(alphas))
stop("Only one of threshold and alpha should be specified.")
if(pvalues && what != "best") {
pvalues <- FALSE
warning("pvalues shown only with what=\"best\".")
}
if(pvalues && missing(perms)) {
pvalues <- FALSE
warning("p-values may be calculated only if perms are provided.")
}
if(inherits(object, "scantwo"))
out <- subrousummaryscantwo(object, for.perm=FALSE)
else
out <- as.data.frame(object, stringsAsFactors=TRUE)
if(!allpairs) # only look at self-self cases
out <- out[out$chr1==out$chr2,]
if(what=="best") {
p1.f <- out$pos1.jnt
p2.f <- out$pos2.jnt
p1.a <- out$pos1.add
p2.a <- out$pos2.add
lf <- out$jnt.lod.full
li <- out$jnt.lod.full - out$add.lod.add
lfv1 <- out$jnt.lod.full - out$lod.1qtl
la <- out$add.lod.add
lav1 <- out$add.lod.add - out$lod.1qtl
}
else if(what=="full") {
p1.f <- p1.a <- out$pos1.jnt
p2.f <- p2.a <- out$pos2.jnt
lf <- out$jnt.lod.full
li <- out$jnt.lod.full - out$jnt.lod.add
lfv1 <- out$jnt.lod.full - out$lod.1qtl
la <- out$jnt.lod.add
lav1 <- out$jnt.lod.add - out$lod.1qtl
}
else if(what=="add") {
p1.f <- p1.a <- out$pos1.add
p2.f <- p2.a <- out$pos2.add
lf <- out$add.lod.full
li <- out$add.lod.full - out$add.lod.add
lfv1 <- out$add.lod.full - out$lod.1qtl
la <- out$add.lod.add
lav1 <- out$add.lod.add - out$lod.1qtl
}
else { # what == "int"
p1.f <- p1.a <- out$pos1.int
p2.f <- p2.a <- out$pos2.int
lf <- out$int.lod.full
li <- out$int.lod.full - out$int.lod.add
lfv1 <- out$int.lod.full - out$lod.1qtl
la <- out$int.lod.add
lav1 <- out$int.lod.add - out$lod.1qtl
}
out <- data.frame(chr1=out$chr1, chr2=out$chr2,
pos1f=p1.f,
pos2f=p2.f,
lod.full=lf, lod.fv1=lfv1, lod.int=li,
pos1a=p1.a,
pos2a=p2.a,
lod.add=la, lod.av1=lav1, stringsAsFactors=TRUE)
if(what != "best") {
out <- out[,-(8:9)]
names(out)[3:4] <- c("pos1", "pos2")
}
if(!missing(perms) && "AA" %in% names(perms)) {
xchr_specific <- TRUE
chr_type <- attr(perms, "chrtype")
chrpair_type <- paste0(chr_type[out$chr1], chr_type[out$chr2])
if(all(chrpair_type=="AA")) {
perms <- perms$AA
xchr_specific <- FALSE
}
else if(all(chrpair_type=="XX")) {
perms <- perms$XX
xchr_specific <- FALSE
}
else if(all(chrpair_type=="AX")) {
perms <- perms$AX
xchr_specific <- FALSE
}
}
else xchr_specific <- FALSE
if(!missing(alphas)) { # get thresholds
if(!xchr_specific) {
thresholds <- rep(0,5)
for(i in 1:5)
thresholds[i] <- quantile(perms[[i]], 1-alphas[i])
thresholds[alphas==1] <- 0
thresholds[alphas==0] <- Inf
}
else { # x-chr-specific
thresholds <- matrix(nrow=3, ncol=5)
for(j in 1:3) {
for(i in 1:5)
thresholds[j,i] <- quantile(perms[[j]][[i]], 1-alphas[i])
}
thresholds[alphas==1] <- 0
thresholds[alphas==0] <- Inf
}
}
if(!missing(thresholds) || !missing(alphas)) { # apply thresholds
if(!xchr_specific) {
out <- out[(out$lod.full >= thresholds[1] & (out$lod.fv1 >= thresholds[2] |
out$lod.int >= thresholds[3])) |
(out$lod.add >= thresholds[4] & out$lod.av1 >= thresholds[5]),,drop=FALSE]
}
else {
tokeep <- NULL
for(i in 1:nrow(out)) {
this_chrpair_type <- match(chrpair_type[i], c("AA", "AX", "XX"))
if((out$lod.full[i] >= thresholds[this_chrpair_type, 1] &
(out$lod.fv1[i] >= thresholds[this_chrpair_type, 2] |
out$lod.int[i] >= thresholds[this_chrpair_type, 3])) |
(out$lod.add[i] >= thresholds[this_chrpair_type, 4] &
out$lod.av1[i] >= thresholds[this_chrpair_type, 5]))
tokeep <- c(tokeep, i)
}
out <- out[tokeep, , drop=FALSE]
}
}
if(pvalues && nrow(out) > 0) {
result <- as.data.frame(matrix(ncol=11+5, nrow=nrow(out)), stringsAsFactors=TRUE)
wh <- c(1,2,3,4,5,7,9,11,12,13,15)
wh2 <- (1:16)[-wh]
result[,wh] <- out
names(result)[wh] <- names(out)
colnames(result)[wh2] <- rep("pval",5)
if(!xchr_specific) {
for(i in 1:5) {
for(j in 1:nrow(out))
result[j,wh2[i]] <- mean(perms[[i]] >= result[j,wh2[i]-1], na.rm=TRUE)
}
}
else { # X-chr-specific
L <- attr(perms, "L")
sumL <- sum(L)
for(j in 1:nrow(out)) {
this_chrpair_type <- match(chrpair_type[j], c("AA", "AX", "XX"))
pow <- sumL/L[this_chrpair_type]
for(i in 1:5) {
nominal_p <- mean(perms[[this_chrpair_type]][[i]] >= result[j,wh2[i]-1], na.rm=TRUE)
result[j,wh2[i]] <- 1 - (1-nominal_p)^pow # adjusted P-value
}
}
}
out <- result
}
class(out) <- c("summary.scantwo", "data.frame")
out
}
# subroutine for summary.scantwo; pulls out the key info
# on each pair of chromosomes
subrousummaryscantwo <-
function(object, for.perm=FALSE)
{
lod <- object$lod
lod[is.na(lod) | lod == Inf | lod == -Inf] <- 0
map <- object$map
pos <- map[,2]
chr <- as.factor(map[,1])
tchr <- as.numeric(chr)
n.chr <- max(tchr)
xchr <- tapply(map[,4], map[,1], function(a) a[1])
xchr <- xchr[!is.na(xchr)]
n.phe <- 1
if(length(dim(lod)) == 3) n.phe <- dim(lod)[3]
if(!("scanoneX" %in% names(object)) ||
is.null(object$scanoneX) || length(object$scanoneX)==0) {
if(n.phe==1) scanoneX <- diag(lod)
else {
if(nrow(lod)==1) scanoneX <- lod[1,1,1]
else scanoneX <- diag(lod[,,1])
for(i in 2:n.phe) {
if(nrow(lod)==1) scanoneX <- cbind(scanoneX, lod[1,1,i])
else scanoneX <- cbind(scanoneX, diag(lod[,,i]))
}
}
}
else scanoneX <- object$scanoneX
if((is.matrix(scanoneX) && nrow(scanoneX) != nrow(lod)) ||
(!is.matrix(scanoneX) && length(scanoneX) != nrow(lod)))
stop("scanoneX component has length ", length(scanoneX), " but should have length ", nrow(lod))
n.chrpair <- n.chr*(n.chr+1)/2
fill <- matrix(0, nrow=n.chrpair, ncol=n.phe)
out <- .C("R_summary_scantwo",
as.integer(nrow(map)),
as.integer(n.phe),
as.double(lod),
as.integer(n.chr),
as.integer(tchr),
as.double(pos),
as.integer(xchr),
as.double(scanoneX),
as.integer(n.chrpair),
chr1=as.integer(rep(0,n.chrpair)),
chr2=as.integer(rep(0,n.chrpair)),
as.integer(rep(0,n.chr*n.chr)),
pos1.jnt=as.double(fill),
pos2.jnt=as.double(fill),
pos1.add=as.double(fill),
pos2.add=as.double(fill),
pos1.int=as.double(fill),
pos2.int=as.double(fill),
jnt.lod.full=as.double(fill),
jnt.lod.add=as.double(fill),
add.lod.full=as.double(fill),
add.lod.add=as.double(fill),
int.lod.full=as.double(fill),
int.lod.add=as.double(fill),
lod.1qtl=as.double(fill),
PACKAGE="qtl")
chr1 <- factor(levels(chr)[out$chr1+1], levels=levels(chr))
chr2 <- factor(levels(chr)[out$chr2+1], levels=levels(chr))
if(n.phe == 1) {
out <- data.frame(chr1=chr1, chr2=chr2,
pos1.jnt=out$pos1.jnt,
pos2.jnt=out$pos2.jnt,
jnt.lod.full=out$jnt.lod.full,
jnt.lod.add=out$jnt.lod.add,
pos1.add=out$pos1.add,
pos2.add=out$pos2.add,
add.lod.full=out$add.lod.full,
add.lod.add=out$add.lod.add,
pos1.int=out$pos1.int,
pos2.int=out$pos2.int,
int.lod.full=out$int.lod.full,
int.lod.add=out$int.lod.add,
lod.1qtl=out$lod.1qtl,
stringsAsFactors=TRUE)
}
else
out <- list(chr1=chr1, chr2=chr2,
pos1.jnt=matrix(out$pos1.jnt, ncol=n.phe),
pos2.jnt=matrix(out$pos2.jnt, ncol=n.phe),
jnt.lod.full=matrix(out$jnt.lod.full, ncol=n.phe),
jnt.lod.add=matrix(out$jnt.lod.add, ncol=n.phe),
pos1.add=matrix(out$pos1.add, ncol=n.phe),
pos2.add=matrix(out$pos2.add, ncol=n.phe),
add.lod.full=matrix(out$add.lod.full, ncol=n.phe),
add.lod.add=matrix(out$add.lod.add, ncol=n.phe),
pos1.int=matrix(out$pos1.int, ncol=n.phe),
pos2.int=matrix(out$pos2.int, ncol=n.phe),
int.lod.full=matrix(out$int.lod.full, ncol=n.phe),
int.lod.add=matrix(out$int.lod.add, ncol=n.phe),
lod.1qtl=matrix(out$lod.1qtl, ncol=n.phe))
if(for.perm) {
if(n.phe==1) {
out <- c("full"=max(out$jnt.lod.full,na.rm=TRUE),
"fv1"=max(out$jnt.lod.full - out$lod.1qtl, na.rm=TRUE),
"int"=max(out$jnt.lod.full - out$add.lod.add, na.rm=TRUE),
"add"=max(out$add.lod.add, na.rm=TRUE),
"av1"=max(out$add.lod.add - out$lod.1qtl, na.rm=TRUE),
"one"=max(out$lod.1qtl, na.rm=TRUE))
}
else {
out <- list("full"=apply(out$jnt.lod.full, 2, max, na.rm=TRUE),
"fv1"=apply(out$jnt.lod.full - out$lod.1qtl, 2, max, na.rm=TRUE),
"int"=apply(out$jnt.lod.full - out$add.lod.add, 2, max, na.rm=TRUE),
"add"=apply(out$add.lod.add, 2, max, na.rm=TRUE),
"av1"=apply(out$add.lod.add - out$lod.1qtl, 2, max, na.rm=TRUE),
"one"=apply(out$lod.1qtl, 2, max, na.rm=TRUE))
}
}
out
}
print.summary.scantwo <-
function(x, ...)
{
if(nrow(x)==0) {
cat(" There were no pairs of loci meeting the criteria.\n")
return(invisible(NULL))
}
z <- as.character(unlist(x[,1]))
if(max(nchar(z)) == 1)
rownames(x) <- apply(x[,1:2], 1, function(a)
paste("c", a, collapse=":", sep=""))
else
rownames(x) <- apply(x[,1:2], 1, function(a)
paste(sprintf("c%-2s", a), collapse=":"))
x <- x[,-(1:2)]
cn <- colnames(x)
if(any(cn=="pos1a"))
cn[cn=="pos1a"] <- " pos1a"
wh <- grep("^pval", cn)
if(length(wh) > 0)
cn[wh] <- "pval"
colnames(x) <- cn
print.data.frame(x, digits=3)
}
print.summary.addpair <-
function(x, ...)
{
if(nrow(x)==0) {
cat(" There were no pairs of loci meeting the criteria.\n")
return(invisible(NULL))
}
z <- as.character(unlist(x[,1]))
if(max(nchar(z)) == 1)
rownames(x) <- apply(x[,1:2], 1, function(a)
paste("c", a, collapse=":", sep=""))
else
rownames(x) <- apply(x[,1:2], 1, function(a)
paste(sprintf("c%-2s", a), collapse=":"))
x <- x[,-(1:2), drop=FALSE]
print.data.frame(x, digits=3)
}
print.scantwo <-
function(x, ...)
{
d <- dim(x$lod)
dn <- dimnames(x$lod)[[3]]
if(is.null(dn)) dn <- attr(x, "phenotypes")
if(nrow(x$lod) == 0)
cat("Empty scantwo object.\n")
else {
if(length(d)==2)
print(summary(x))
else {
for(i in 1:d[3]) {
if(is.null(dn)) cat("Phenotype", i, "\n")
else cat(dn[i], "\n")
print(summary(x, lod=i))
}
}
}
}
######################################################################
#
# max.scantwo: Give pair of chromosome with maximum 2-locus LOD score
#
######################################################################
max.scantwo <-
function(object, lodcolumn=1,
what=c("best", "full", "add", "int"),
na.rm=TRUE, ...)
{
if(!inherits(object, "scantwo") &&
!inherits(object, "scantwocondensed"))
stop("Input must have class \"scantwo\".")
addpair <- attr(object, "addpair")
if(!is.null(addpair) && addpair) { # special treatment for output for addpair
temp <- summary(object)
mx <- max(temp[,5],na.rm=TRUE)
return(temp[!is.na(temp[,5]) & temp[,5]==mx,])
}
what <- match.arg(what)
if(inherits(object, "scantwo")) {
d <- dim(object$lod)
if(length(d)==3) {
if(lodcolumn < 1 || lodcolumn > d[3])
stop("lodcolumn must be between 1 and ", d[3])
object$lod <- object$lod[,,lodcolumn]
}
out <- subrousummaryscantwo(object, for.perm=FALSE)
}
else { # condensed version
if(is.matrix(object$pos1.jnt)) {
d <- ncol(object$pos1.jnt)
if(lodcolumn < 1 || lodcolumn > d)
stop("lodcolumn must be between 1 and ", d)
for(i in 3:length(object))
object[[i]] <- object[[i]][,lodcolumn]
}
out <- as.data.frame(object, stringsAsFactors=TRUE)
}
if(what=="best") {
wh <- which(!is.na(out$jnt.lod.full) & out$jnt.lod.full==max(out$jnt.lod.full, na.rm=TRUE))
p1.f <- out$pos1.jnt[wh]
p2.f <- out$pos2.jnt[wh]
p1.a <- out$pos1.add[wh]
p2.a <- out$pos2.add[wh]
lf <- out$jnt.lod.full[wh]
li <- out$jnt.lod.full[wh] - out$add.lod.add[wh]
lfv1 <- out$jnt.lod.full[wh] - out$lod.1qtl[wh]
la <- out$add.lod.add[wh]
lav1 <- out$add.lod.add[wh] - out$lod.1qtl[wh]
}
else if(what=="full") {
wh <- which(!is.na(out$jnt.lod.full) & out$jnt.lod.full==max(out$jnt.lod.full, na.rm=TRUE))
p1.f <- p1.a <- out$pos1.jnt[wh]
p2.f <- p2.a <- out$pos2.jnt[wh]
lf <- out$jnt.lod.full[wh]
li <- out$jnt.lod.full[wh] - out$jnt.lod.add[wh]
lfv1 <- out$jnt.lod.full[wh] - out$lod.1qtl[wh]
la <- out$jnt.lod.add[wh]
lav1 <- out$jnt.lod.add[wh] - out$lod.1qtl[wh]
}
else if(what=="add") {
wh <- which(!is.na(out$add.lod.add) & out$add.lod.add==max(out$add.lod.add, na.rm=TRUE))
p1.f <- p1.a <- out$pos1.add[wh]
p2.f <- p2.a <- out$pos2.add[wh]
lf <- out$add.lod.full[wh]
li <- out$add.lod.full[wh] - out$add.lod.add[wh]
lfv1 <- out$add.lod.full[wh] - out$lod.1qtl[wh]
la <- out$add.lod.add[wh]
lav1 <- out$add.lod.add[wh] - out$lod.1qtl[wh]
}
else { # what == "int"
lod.int <- out$int.lod.full - out$int.lod.add
wh <- which(!is.na(lod.int) & lod.int == max(lod.int, na.rm=TRUE))
p1.f <- p1.a <- out$pos1.int[wh]
p2.f <- p2.a <- out$pos2.int[wh]
lf <- out$int.lod.full[wh]
li <- out$int.lod.full[wh] - out$int.lod.add[wh]
lfv1 <- out$int.lod.full[wh] - out$lod.1qtl[wh]
la <- out$int.lod.add[wh]
lav1 <- out$int.lod.add[wh] - out$lod.1qtl[wh]
}
out <- data.frame(chr1=out$chr1[wh], chr2=out$chr2[wh],
pos1f=p1.f,
pos2f=p2.f,
lod.full=lf, lod.fv1=lfv1, lod.int=li,
pos1a=p1.a,
pos2a=p2.a,
lod.add=la, lod.av1=lav1, stringsAsFactors=TRUE)
if(what != "best") {
out <- out[,-(8:9)]
names(out)[3:4] <- c("pos1", "pos2")
}
class(out) <- c("summary.scantwo", "data.frame")
rownames(out) <- what
out
}
######################################################################
# clean.scantwo
#
# sets LOD scores that are missing or < 0 to 0
# If full LOD < add've LOD, set full = add've
# sets LOD scores, for pairs of positions that are not separated
# by n.mar markers and distance cM, to 0
######################################################################
clean.scantwo <-
function(object, n.mar=1, distance=0, ...)
{
if(!inherits(object, "scantwo"))
stop("Input should have class \"scantwo\".")
addpair <- attr(object, "addpair")
if(is.null(addpair)) addpair <- FALSE
lod <- object$lod
map <- object$map
if(!("fullmap" %in% names(attributes(object))))
stop("clean.scantwo only works on scantwo objects created with R/qtl ver >= 1.04-38.\n")
fmap <- attr(object, "fullmap")
lod[is.na(lod) | lod < 0] <- 0
subrou <- function(x,y,z)
{
out <- x
for(i in seq(along=x))
out[i] <- sum((z > x[i] & z < y[i]) | (z < x[i] & z > y[i])) +
(x[i] == y[i])
out
}
last <- 0
for(i in seq(along=fmap)) {
m <- map[map[,1]==names(fmap)[i],2]
idx <- 1:length(m)+last
last <- last + length(m)
toclean <- (outer(m, m, subrou, fmap[[i]]) < n.mar) | (abs(outer(m, m, "-")) <distance)
diag(toclean) <- FALSE
if(length(dim(lod))==3) # case of multiple phenotypes
lod[idx,idx,][toclean] <- 0
else
lod[idx,idx][toclean] <- 0
}
# if full LOD < add've LOD, set full = add've
if(!addpair) {
if(length(dim(lod)) == 2) {
lo <- lower.tri(lod)
wh <- (lod[lo] < t(lod)[lo])
if(any(wh))
lod[lo][wh] <- t(lod)[lo][wh]
}
else {
lo <- lower.tri(lod[,,1])
for(i in 1:dim(lod)[3]) {
wh <- (lod[,,i][lo] < t(lod[,,i])[lo])
if(any(wh))
lod[,,i][lo][wh] <- t(lod[,,i])[lo][wh]
}
}
}
object$lod <- lod
object
}
######################################################################
#
# subset.scantwo
#
######################################################################
subset.scantwo <-
function(x, chr, lodcolumn, ...)
{
if(!inherits(x, "scantwo"))
stop("Input should have class \"scantwo\".")
if((missing(chr) || length(chr)==0) && missing(lodcolumn)) return(x)
if(!missing(lodcolumn)) {
if(any(lodcolumn>0) && any(lodcolumn<0))
stop("lodcolumn values can't be both >0 and <0.")
if(any(lodcolumn<0) || is.logical(lodcolumn))
lodcolumn <- (1:(dim(x$lod)[3]))[lodcolumn]
if(length(lodcolumn)==0)
stop("You must retain at least one LOD column.")
if(any(lodcolumn < 1 || lodcolumn > dim(x$lod)[3]))
stop("lodcolumn values must be >=1 and <=",dim(x$lod)[3])
x$lod <- x$lod[,,lodcolumn]
if("scanoneX" %in% names(x))
x$scanoneX <- x$scanoneX[,lodcolumn]
}
if(!missing(chr)) {
chr <- matchchr(chr, unique(x$map[,1]))
wh <- x$map[,1] %in% chr
x$map <- x$map[wh,,drop=FALSE]
x$map[,1] <- droplevels(x$map[,1])
if(length(dim(x$lod))==2)
x$lod <- x$lod[wh,wh,drop=FALSE]
else
x$lod <- x$lod[wh,wh,,drop=FALSE]
if(!is.null(x$scanoneX))
x$scanoneX <- x$scanoneX[wh]
if("fullmap" %in% names(attributes(x))) {
fmap <- attr(x, "fullmap")
fmap <- fmap[chr]
attr(x, "fullmap") <- fmap
}
}
x
}
######################################################################
# summary.scantwoperm
#
# Give genome-wide LOD thresholds on the basis of results of
# scantwo permutation test (from scantwo with n.perm > 0)
######################################################################
summary.scantwoperm <-
function(object, alpha=c(0.05, 0.10), ...)
{
if(!inherits(object, "scantwoperm"))
stop("Input should have class \"scantwoperm\".")
if("AA" %in% names(object)) { # X-chr-specific version
# get region-specific (nominal) signif levels
L <- attr(object, "L")
LL <- attr(object, "LL")
if(is.null(LL)) stop("LL attribute not found in input object")
if(length(alpha)==1) {
tmp <- L/sum(L); tmp <- c(tmp[1], 1, tmp[2])
one_minus_alpha_onechr <- cbind((1-alpha)^tmp)
one_minus_alpha <- cbind((1-alpha)^(LL/sum(LL)))
}
else {
tmp <- L/sum(L); tmp <- c(tmp[1], 1, tmp[2])
one_minus_alpha_onechr <- vapply(alpha, function(a,b) (1-a)^b, rep(0, length(tmp)), tmp)
one_minus_alpha <- vapply(alpha, function(a,b) (1-a)^b, rep(0, length(LL)), LL/sum(LL))
}
# get quantiles
out <- vector("list", length(object))
names(out) <- names(object)
for(i in seq(along=out)) {
f <- function(a, qu) {
b <- apply(a, 2, quantile, qu)
if(!is.matrix(b)) {
nam <- names(b)
b <- matrix(b, nrow=1)
colnames(b) <- nam
}
rownames(b) <- paste0(alpha*100, "%")
b
}
out[[i]] <- lapply(unclass(object[[i]])[1:5], f, one_minus_alpha[i,,drop=FALSE])
qu_one <- f(object[[i]][[6]], one_minus_alpha_onechr[i,,drop=FALSE])
out[[i]] <- c(out[[i]], "one"=list(qu_one))
}
attr(out, "n.perm") <- vapply(object, function(a) nrow(a[[1]]), 0)
class(out) <- c("summary.scantwoperm", "list")
return(out)
}
out <- lapply(object, apply, 2, quantile, 1-alpha)
for(i in 1:length(out)) {
if(!is.matrix(out[[i]]))
out[[i]] <- matrix(out[[i]], nrow=length(alpha))
rownames(out[[i]]) <- paste0(alpha*100, "%")
}
attr(out, "n.perm") <- nrow(object[[1]])
class(out) <- c("summary.scantwoperm", "list")
out
}
print.summary.scantwoperm <-
function(x, ...)
{
if("AA" %in% names(x)) { # X-sp perms
nperm <- attr(x, "n.perm")
rn <- rownames(x[[1]][[1]])
nc <- rapply(x, ncol)
if(length(unique(nc)) != 1)
stop("The components shouldn't have varying numbers of columns.\n")
nc <- nc[1]
phe <- colnames(x[[1]][[1]])
convert2matrix <-
function(a, which_col=1) {
a <- lapply(a, "[", , which_col, drop=FALSE)
b <- matrix(unlist(a), ncol=6)
rownames(b) <- rn
colnames(b) <- names(a)
b
}
for(i in seq(along=phe)) {
cat(phe[i], ":\n", sep="")
y <- lapply(x, convert2matrix, i)
lab <- c(AA="A:A", AX="A:X", XX="X:X")
for(j in c("AA", "AX", "XX")) {
cat(" ", lab[j], " (", nperm[j], " permutations)\n", sep="")
print(y[[j]], digits=3)
cat("\n")
}
if(i != length(phe)) cat("\n")
}
invisible(return(x))
}
nam <- names(x)
rn <- rownames(x[[1]])
n.perm <- attr(x, "n.perm")
nc <- sapply(x, ncol)
if(length(unique(nc)) != 1)
stop("The components shouldn't have varying numbers of columns.\n")
nc <- nc[1]
if(nc==1) {
phe <- colnames(x[[1]])
if(is.null(phe)) phe <- ""
x <- matrix(unlist(x), nrow=length(rn))
dimnames(x) <- list(rn, nam)
if(is.null(phe))
cat("(", n.perm, " permutations)\n", sep="")
else
cat(phe, " (", n.perm, " permutations)\n", sep="")
print(x, digits=3)
}
else {
phe <- colnames(x[[1]])
if(is.null(phe)) phe <- paste("pheno", 1:nc)
for(i in 1:nc) {
y <- matrix(ncol=length(x), nrow=nrow(x[[1]]))
for(j in 1:length(x))
y[,j] <- x[[j]][,i]
dimnames(y) <- list(rn, nam)
cat(phe[i], " (", n.perm, " permutations)\n", sep="")
print(y, digits=3)
if(i != nc) cat("\n")
}
}
}
######################################################################
# combine scantwo results ... paste the phenotype columns together
######################################################################
cbind.scantwo <- c.scantwo <-
function(...)
{
dots <- list(...)
if(length(dots)==1 && is.list(dots[[1]])) dots <- dots[[1]]
if(length(dots)==1) return(dots[[1]])
for(i in seq(along=dots)) {
if(!inherits(dots[[i]], "scantwo"))
stop("Input should have class \"scantwo\".")
}
# check dimensions of LODs
nr <- sapply(dots, function(a) nrow(a$lod))
nc <- sapply(dots, function(a) ncol(a$lod))
nd3 <- sapply(dots, function(a) dim(a$lod)[3])
if(any(nr[-1]!=nr[1]) || any(nc[-1] != nc[1]))
stop("Input objects are not the same dimensions.")
# check maps
map1 <- dots$map[[1]]
for(i in 2:length(dots)) {
if(!all(dots$map[[i]] != map1))
stop("Maps are not all the same.")
}
output <- dots[[1]]
nd3[is.na(nd3)] <- 1
output$lod <- array(dim=c(nr[1], nc[1], sum(nd3)))
start <- cumsum(c(0,nd3))[-length(nd3)-1]
end <- start+nd3
for(i in seq(along=dots))
output$lod[,,(start[i]+1):end[i]] <- dots[[i]]$lod
attr(output, "phenotypes") <- unlist(lapply(dots, attr, "phenotypes"))
dimnames(output$lod) <- list(NULL, NULL, attr(output, "phenotypes"))
# check scanoneX
if(is.null(dots[[1]]$scanoneX)) {
for(i in 2:length(dots)) {
if(!is.null(dots[[i]]$scanoneX))
stop("Some but not all input objects have null scanoneX.")
}
}
else {
nrx <- sapply(dots, function(a) nrow(a$scanoneX))
if(any(nrx[-1]!=nrx[1]))
stop("Mismatch in scanoneX dimensions.")
for(i in 2:length(dots))
output$scanoneX <- cbind(output$scanoneX, dots[[i]]$scanoneX)
colnames(output$scanoneX) <- attr(output, "phenotypes")
}
methods <- sapply(dots, attr, "method")
if(length(unique(methods)) != 1)
attr(output, "method") <- rep(sapply(dots, attr, "method"), nd3)
fm <- lapply(dots, attr, "fullmap")
for(i in 2:length(fm)) {
if(length(fm[[1]]) != length(fm[[i]]) || !all(names(fm[[1]]) == names(fm[[i]])))
stop("Mismatch in \"fullmap\" attributes (1).")
for(j in 1:length(fm[[1]])) {
if(length(fm[[1]][[j]]) != length(fm[[i]][[j]]) ||
!all(names(fm[[1]][[j]]) == names(fm[[i]][[j]])) ||
max(abs(fm[[1]][[j]] - fm[[i]][[j]])) > 0.001)
stop("Mismatch in \"fullmap\" attributes (2).")
}
}
output
}
######################################################################
# combine scantwoperm results ... paste the rows together
######################################################################
rbind.scantwoperm <- c.scantwoperm <-
function(...)
{
dots <- list(...)
if(length(dots)==1 && is.list(dots[[1]])) dots <- dots[[1]]
if(length(dots)==1) return(dots[[1]])
xchrsp <- vapply(dots, function(a) "AA" %in% names(a), TRUE)
if(any(xchrsp)) {
if(!all(xchrsp)) stop("Some but not all inputs are X-chr specific")
for(i in 2:length(dots)) {
for(j in seq(along=dots[[1]])) {
for(k in seq(along=dots[[1]][[j]])) {
if(ncol(dots[[1]][[j]][[k]]) != ncol(dots[[i]][[j]][[k]]))
stop("Mismatch in no. columns")
if(any(colnames(dots[[1]][[j]][[k]]) != colnames(dots[[1]][[j]][[k]])))
warning("Mismatch in column names")
dots[[1]][[j]][[k]] <- rbind(dots[[1]][[j]][[k]], dots[[i]][[j]][[k]])
}
}
}
return(dots[[1]])
}
for(i in seq(along=dots)) {
if(!inherits(dots[[i]], "scantwoperm"))
stop("Input should have class \"scantwoperm\".")
}
nc <- sapply(dots, function(a) ncol(a[[1]]))
if(length(unique(nc)) != 1)
stop("Number of LOD columns in the input objects must be constant.\n")
flag <- 0
for(j in 1:length(dots[[1]])) {
for(i in 2:length(dots)) {
dots[[1]][[j]] <- rbind(dots[[1]][[j]], dots[[i]][[j]])
if(any(colnames(dots[[i]][[j]]) != colnames(dots[[1]][[j]]))) flag <- 1
}
}
if(flag) warning("Mismatch in column names; input may not be consistent.\n")
dots[[1]]
}
# paste columns together
cbind.scantwoperm <-
function(...)
{
dots <- list(...)
if(length(dots)==1 && is.list(dots[[1]])) dots <- dots[[1]]
if(length(dots)==1) return(dots)
xchrsp <- vapply(dots, function(a) "AA" %in% names(a), TRUE)
if(any(xchrsp)) {
if(!all(xchrsp)) stop("Some but not all inputs are X-chr specific")
for(i in 2:length(dots)) {
for(j in seq(along=dots[[1]])) {
for(k in seq(along=dots[[1]][[j]])) {
if(nrow(dots[[1]][[j]][[k]]) != nrow(dots[[i]][[j]][[k]]))
stop("Mismatch in no. permutations")
dots[[1]][[j]][[k]] <- cbind(dots[[1]][[j]][[k]], dots[[i]][[j]][[k]])
}
}
}
return(dots[[1]])
}
for(i in 2:length(dots)) {
for(j in seq(along=dots[[1]])) {
if(nrow(dots[[1]][[j]]) != nrow(dots[[i]][[j]]))
stop("Mismatch in no. permutations")
dots[[1]][[j]] <- cbind(dots[[1]][[j]], dots[[i]][[j]])
}
}
dots[[1]]
}
######################################################################
# condensed scantwo output
condense <-
function(object)
UseMethod("condense")
condense.scantwo <-
function(object)
{
out <- subrousummaryscantwo(object, for.perm=FALSE)
class(out) <- c("scantwocondensed", "list")
out
}
summary.scantwocondensed <- summary.scantwo
max.scantwocondensed <- max.scantwo
##############################
# subset.scantwoperm: pull out a set of lodcolumns
##############################
subset.scantwoperm <-
function(x, repl, lodcolumn, ...)
{
if("AA" %in% names(x)) { # x chr specific
for(j in seq(along=x)) {
if(missing(lodcolumn)) lodcolumn <- 1:ncol(x[[j]][[1]])
else if(!check_colindex(lodcolumn, x[[j]][[1]]))
stop("lodcolumn misspecified.")
repl <- 1:nrow(x[[j]][[1]])
x[[j]] <- lapply(x[[j]], function(a,b,d) unclass(a)[b,d,drop=FALSE], repl, lodcolumn)
}
return(x)
}
att <- attributes(x)
if(any(!sapply(x, is.matrix)))
x <- lapply(x, as.matrix)
if(missing(lodcolumn)) lodcolumn <- 1:ncol(x[[1]])
else if(!check_colindex(lodcolumn, x[[1]]))
stop("lodcolumn misspecified.")
if(missing(repl)) repl <- 1:nrow(x[[1]])
else if(!check_rowindex(repl, x[[1]]))
stop("repl misspecified.")
cl <- class(x)
x <- lapply(x, function(a,b,d) unclass(a)[b,d,drop=FALSE], repl, lodcolumn)
class(x) <- cl
for(i in seq(along=att)) {
if(names(att)[i] == "dim" || length(grep("names", names(att)[i]))>0) next
attr(x, names(att)[i]) <- att[[i]]
}
x
}
# subset.scantwoperm using [,]
`[.scantwoperm` <-
function(x, repl, lodcolumn)
subset.scantwoperm(x, repl, lodcolumn)
# end of summary.scantwo.R
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.