#setOldClass(c("haplotype", "genotype"), test=TRUE)
#content should be raw - TODO: put restriction
setClass("SnpMatrix", contains="matrix")
setClass("XSnpMatrix", representation("SnpMatrix", diploid="logical"),
contains="SnpMatrix")
# constructor
setMethod("initialize",
"SnpMatrix",
function(.Object, ...) {
if(is.null(.Object)) {
stop("object is null to constructor");
}
.Object <- callNextMethod()
# Please don't remove the tests - they are to do with promises and not redundant.
# The accessor methods transparently passes to .Data, the assignment methods are different...
if(!is.null(dim(.Object)) && (dim(.Object)[1] > 0) && (dim(.Object)[2] > 0)) {
if (mode(.Object) != "raw") {
cat("coercing object of mode ", mode(.Object), " to SnpMatrix\n")
mode(.Object@.Data) <- "raw"
}
if (is.null(dimnames(.Object))) {
cat("object has no names - using numeric order for row/column names\n")
dimnames(.Object@.Data) <- c(list(as.character(1:(dim(.Object)[1]))), list(as.character(1:(dim(.Object)[2]))))
}
}
.Object
})
setMethod("initialize",
"XSnpMatrix",
function(.Object, ...) {
.Object <- callNextMethod()
if (length(.Object@diploid)==0){
warning("diploid slot not supplied; it has been estimated from heterozygosity")
guess <- .guessPloidy(.Object)
.Object@diploid <- guess$diploid
het <- guess$Heterozygosity
}
else {
lend <- length(.Object@diploid)
if (lend==1)
.Object@diploid <- rep(.Object@diploid, nrow(.Object))
else if (lend!=nrow(.Object))
stop("diploid argument incorrect length")
dpna <- is.na(.Object@diploid)
if (any(dpna)) {
warning("diploid argument contains NAs; these have been replaced by guesses based on heterozygosity")
guess <- .guessPloidy(.Object, .Object@diploid)
.Object@diploid <- guess$diploid
het <- guess$Heterozygosity
}
else
het <- row.summary(.Object)$Heterozygosity
}
hetm <- !.Object@diploid & (!is.na(het)&(het>0))
if (any(hetm)){
warning(sum(hetm, na.rm=TRUE),
" heterozygous calls for haploid genotypes; these were set to NA")
.Object@.Data <- .forceHom(.Object@.Data, .Object@diploid)
}
.Object
})
setMethod("[", signature(x="SnpMatrix",i="ANY",j="ANY",drop="ANY"),
function(x, i, j, drop=FALSE) {
if (drop!=FALSE)
stop("dimensions cannot be dropped from a SnpMatrix object")
if (missing(i)) {
if (missing(j))
return(x)
else
x <- x@.Data[,j,drop=FALSE]
}
else {
if (missing(j))
x <- x@.Data[i,,drop=FALSE]
else
x <- x@.Data[i,j,drop=FALSE]
}
cl <- "SnpMatrix"
attr(cl, "package") <- "snpStats"
class(x) <- cl
# setting an S4 class doesn't not automatically
# set the object's internal tag; do it manually
x <- asS4(x)
x
}
)
setMethod("[", signature(x="XSnpMatrix",i="ANY",j="ANY",drop="ANY"),
function(x, i, j, drop=FALSE) {
if (drop!=FALSE)
stop("dimensions cannot be dropped from an XSnpMatrix object")
diploid <- x@diploid
if (missing(i)) {
if (missing(j))
return(x)
else
x <- x@.Data[,j, drop=FALSE]
}
else {
if (is.character(i))
i <- match(i, rownames(x))
diploid <- diploid[i]
if (missing(j))
x <- x@.Data[i,,drop=FALSE]
else
x <- x@.Data[i,j,drop=FALSE]
}
new("XSnpMatrix", x, diploid=diploid)
})
# The sub-assignment methods
setMethod("[<-", signature(x="XSnpMatrix",i="ANY",j="ANY",value="XSnpMatrix"),
function(x, i, j, value) {
# The raw sub assignment should *just work*:
if (!missing(i) & !missing(j)) {
x@.Data[i,j] <- value
} else if (missing(i) & missing(j)) {
x@.Data[,] <- value
} else if (missing(i)) {
x@.Data[,j] <- value
} else if (missing(j)) {
x@.Data[i,] <- value
}
# All we want to do is to shuffle the rows
# if a row index is passed
if(!missing(i)) {
slot(x, "diploid")[i] <- slot(value, "diploid")
}
x
})
setAs("SnpMatrix", "numeric",
function(from) {
.Call("asnum", from, PACKAGE="snpStats")
})
setAs("SnpMatrix", "character",
function(from) {
df <- dim(from)
dnames <- dimnames(from)
from <- 1+as.integer(from)
to <- ifelse(from<5, c("NA", "A/A", "A/B", "B/B")[from], "Uncertain")
dim(to) <- df
dimnames(to) <- dnames
to
})
setAs("SnpMatrix", "XSnpMatrix",
function(from) {
new("XSnpMatrix", from)
})
setAs("XSnpMatrix", "character",
function(from) {
df <- dim(from)
ifr <- 1 + as.integer(from)
offset <- 4*rep(from@diploid, ncol(from))
to <- ifelse(ifr<5, c("NA", "A", "Error", "B",
"NA", "A/A", "A/B", "B/B")[ifr + offset],
"Uncertain")
dim(to) <- df
dimnames(to) <- dimnames(from)
to
})
setAs("matrix", "SnpMatrix",
function(from) {
if (is.numeric(from)) {
valid <- (from==0 | from==1 | from==2)
if(!all(is.na(from) | valid))
warning("values other than 0, 1 or 2 set to NA")
to <- from+1
to[is.na(from)|!valid] <- 0
mode(to) <- "raw"
}
else if (is.raw(from)) {
not.valid <- (from>3)
if (any(not.valid))
warning("values other than 01, 02 or 03 set to NA (00)")
to <- from
to[not.valid] <- as.raw(0)
}
else
stop("Can only convert `numeric' or `raw' matrices")
new("SnpMatrix", to)
})
# Summary of a SnpMatrix
setGeneric("summary")
setMethod("summary", "SnpMatrix",
function(object) {
list(rows = summary(row.summary(object)),
cols = summary(col.summary(object, uncertain=TRUE)))
})
setMethod("summary", "XSnpMatrix",
function(object) {
tab <- table(object@diploid)
names(tab) <- ifelse(as.logical(names(tab)), "diploid", "haploid")
list(sex = tab,
rows = summary(row.summary(object)),
cols = summary(col.summary(object, uncertain=TRUE)))
})
# Imputation
setClass("ImputationRules", contains="list")
setMethod("summary", "ImputationRules",
function(object) {
info <- .Call("r2_impute", object, PACKAGE="snpStats")
i2 <- info[,2]
levs <- sort(unique(i2))
labs <- paste(levs, "tags")
size <- factor(i2, levels= levs, labels=labs)
r2 <- info[,1]
r2[r2>1] <- 1
byr2 <- cut(r2, c((0:9)/10, 0.95, 0.99, 1.0), right=FALSE,
include.lowest=TRUE)
table(byr2, size, dnn=c("R-squared", "SNPs used"), useNA="ifany")
})
setMethod("[",
signature(x="ImputationRules", i="ANY", j="missing", drop="missing"),
function(x, i){
if (is.character(i))
i <- match(i, names(x), nomatch=0)
ss <- x@.Data[i]
names(ss) <- names(x)[i]
res <- new("ImputationRules", ss)
attr(res, "Max.predictors") <-
max(sapply(res, function(rule) length(rule$snps)))
res
})
# Show and plot methods
setMethod("show", "ImputationRules",
function(object) {
to <- names(object)
for (i in 1:length(object)) {
if (is.null(object[[i]]$snps))
cat(to[i], "~ No imputation available\n")
else {
if (is.null(object[[i]]$hap.probs))
stop("Old imputation rule - not supported by this version")
else
cat(to[i], " ~ ", paste(object[[i]]$snps, collapse="+"))
cat(" (MAF = ", object[[i]]$maf,
", R-squared = ", object[[i]]$r.squared,
")\n", sep="")
}
}
})
setMethod("plot", signature(x="ImputationRules", y="missing"),
function(x, y, ...) {
mat <- summary(x)
n <- nrow(mat)
m <- ncol(mat)
if (is.na(rownames(mat)[n])) {
mat <- mat[-n,-m]
n <- n-1
m <- m-1
}
val <- barplot(t(mat[n:1,]), legend.text=TRUE,
beside=F, col=heat.colors(m),
xlab=expression(r^2), ylab="Number of imputed SNPs",
names.arg=rep("", n), cex.lab=1.3, ...)
mtext(rownames(mat)[n:1], at=val, side=1, las=2, line=0.5, cex=0.8)
})
setMethod("show", "XSnpMatrix",
function(object) {
nr <- nrow(object)
nc <- ncol(object)
dip <- object@diploid
cat("An XSnpMatrix with ", nr, " rows (", sum(!dip), " haploid and ",
sum(dip), " diploid), and ", nc, " columns\n", sep="")
if (nr>1)
cat("Row names: ", rownames(object)[1],"...", rownames(object)[nr],"\n")
else
cat("Row name: ", rownames(object)[1], "\n")
if (nc>1)
cat("Col names: ", colnames(object)[1],"...", colnames(object)[nc],"\n")
else
cat("Col name: ", colnames(object)[1],"\n")
})
setMethod("show", "SnpMatrix",
function(object) {
nr <- nrow(object)
nc <- ncol(object)
cat("A SnpMatrix with ", nr, "rows and ", nc,
"columns\n")
if (nr>1)
cat("Row names: ", rownames(object)[1],"...", rownames(object)[nr],"\n")
else
cat("Row name: ", rownames(object)[1],"\n")
if (nc>1)
cat("Col names: ", colnames(object)[1],"...", colnames(object)[nc],"\n")
else
cat("Col name: ", colnames(object)[1],"\n")
})
setMethod("is.na", "SnpMatrix", function(x){ x==0})
snp.rbind <- function(...){
.External("snp_rbind", ..., PACKAGE="snpStats")
}
snp.cbind <- function(...){
.External("snp_cbind", ..., PACKAGE="snpStats")
}
setMethod("rbind", signature(...="SnpMatrix"), snp.rbind)
setMethod("cbind", signature(...="SnpMatrix"), snp.cbind)
# Tests
setClass("SingleSnpTests",
representation(snp.names="character", chisq="matrix",
N="integer", N.r2="numeric"))
setClass("SingleSnpTestsScore",
representation("SingleSnpTests", U="matrix", V="matrix"),
contains="SingleSnpTests")
setClass("GlmTests",
representation(snp.names="ANY", var.names="character",
chisq="numeric", df="integer",
N="integer"))
setClass("GlmTestsScore", representation("GlmTests", score="list"),
contains="GlmTests")
setMethod("[",
signature(x="SingleSnpTests", i="ANY",
j="missing", drop="missing"),
function(x, i, j, drop) {
if (is.character(i)) {
i <- match(i, x@snp.names)
if (any(is.na(i))) {
warning(sum(is.na(i)),
" SNPs couldn't be found in SingleSnpTests object\n")
i <- i[!is.na(i)]
}
} else if (is.logical(i)) {
if (length(i)!=length(x@snp.names))
stop("logical selection array has incorrect length")
i <- (1:length(i))[i]
} else if (is.numeric(i)) {
if (min(i)<0 || max(i)>length(x@snp.names))
stop("selection index out of range")
} else {
stop("illegal selection array")
}
if (length(x@N.r2)>0)
N.r2 <- x@N.r2[i]
else
N.r2 <- numeric(0)
new("SingleSnpTests",
snp.names=x@snp.names[i], chisq=x@chisq[i,, drop=FALSE],
N=x@N[i], N.r2=N.r2)
})
setMethod("[",
signature(x="SingleSnpTestsScore", i="ANY",
j="missing", drop="missing"),
function(x, i, j, drop) {
if (is.character(i)) {
i <- match(i, x@snp.names)
if (any(is.na(i))) {
warning(sum(is.na(i)),
" SNPs couldn't be found in SingleSnpTests object\n")
i <- i[!is.na(i)]
}
} else if (is.logical(i)) {
if (length(i)!=length(x@snp.names))
stop("logical selection array has incorrect length")
i <- (1:length(i))[i]
} else if (is.numeric(i)) {
if (min(i)<0 || max(i)>length(x@snp.names))
stop("selection index out of range")
} else {
stop("illegal selection array")
}
if (length(x@N.r2)>0)
N.r2 <- x@N.r2[i]
else
N.r2 <- numeric(0)
new("SingleSnpTestsScore",
snp.names=x@snp.names[i], chisq=x@chisq[i,,drop=FALSE],
N=x@N[i], N.r2=N.r2,
U=x@U[i,,drop=FALSE], V=x@V[i,,drop=FALSE])
})
setMethod("[",
signature(x="GlmTests", i="ANY",
j="missing", drop="missing"),
function(x, i, j, drop) {
if (is.character(i))
i <- match(i, names(x), nomatch=0)
new("GlmTests",
snp.names = x@snp.names[i],
var.names = x@var.names,
chisq = x@chisq[i],
df = x@df[i],
N = x@N[i])
})
setMethod("[",
signature(x="GlmTestsScore", i="ANY",
j="missing", drop="missing"),
function(x, i, j, drop) {
if (is.character(i))
i <- match(i, names(x), nomatch=0)
new("GlmTestsScore",
snp.names = x@snp.names[i],
var.names = x@var.names,
chisq = x@chisq[i],
df = x@df[i],
N = x@N[i],
score = x@score[i])
})
setAs("SingleSnpTests", "data.frame",
function(from) {
if (length(from@N.r2)>0)
data.frame(N=from@N, N.r2=from@N.r2,
Chi.squared=from@chisq,
P.1df=pchisq(from@chisq[,1], df=1, lower.tail=FALSE),
P.2df=pchisq(from@chisq[,2], df=2, lower.tail=FALSE),
row.names=from@snp.names)
else
data.frame(N=from@N,
Chi.squared=from@chisq,
P.1df=pchisq(from@chisq[,1], df=1, lower.tail=FALSE),
P.2df=pchisq(from@chisq[,2], df=2, lower.tail=FALSE),
row.names=from@snp.names)
})
setMethod("summary", "SingleSnpTests",
function(object) {
summary(as(object, 'data.frame'))
})
setMethod("show", "SingleSnpTests",
function(object) {
print(as(object, 'data.frame'))
})
setAs("GlmTests", "data.frame",
function(from) {
chi2 <- from@chisq
df <- from@df
p <- pchisq(chi2, df=df, lower.tail=FALSE)
N <- from@N
data.frame(row.names=names(from),
Chi.squared=chi2, Df=df, p.value=p)
})
setMethod("summary", "GlmTests",
function(object) {
summary(as(object, "data.frame"))
})
setMethod("show", "GlmTests",
function(object) {
print(as(object, "data.frame"))
})
# There are no standard generics for these, but the call to standardGeneric
# avoids a warning message
setGeneric("p.value", function(x, df) standardGeneric("p.value"),
useAsDefault=FALSE)
setGeneric("chi.squared", function(x, df) standardGeneric("chi.squared"),
useAsDefault=FALSE)
setGeneric("deg.freedom", function(x) standardGeneric("deg.freedom"),
useAsDefault=FALSE)
setGeneric("effect.sign", function(x, simplify) standardGeneric("effect.sign"),
useAsDefault=FALSE)
setGeneric("sample.size", function(x) standardGeneric("sample.size"),
useAsDefault=FALSE)
setGeneric("effective.sample.size",
function(x) standardGeneric("effective.sample.size"),
useAsDefault=FALSE)
setGeneric("pool2", function(x, y, score) standardGeneric("pool2"),
useAsDefault=FALSE)
setGeneric("names")
setMethod("p.value", signature(x="SingleSnpTests", df="numeric"),
function(x, df) {
if (df>2)
stop("df must be 1 or 2")
p <- pchisq(x@chisq[,df], df=df, lower.tail=FALSE)
names(p) <- x@snp.names
p
})
setMethod("p.value", signature(x="GlmTests", df="missing"),
function(x, df) {
p <- pchisq(q=x@chisq, x@df, lower.tail=FALSE)
p
})
setMethod("chi.squared", signature(x="SingleSnpTests", df="numeric"),
function(x, df) {
if (df>2)
stop("df must be 1 or 2")
chi2 <- x@chisq[,df]
names(chi2) <- x@snp.names
chi2
})
setMethod("chi.squared", signature(x="GlmTests", df="missing"),
function(x, df) {
chi2 <- x@chisq
chi2
})
setMethod("deg.freedom", signature(x="GlmTests"),
function(x) {
df <- x@df
df
})
setMethod("effect.sign", signature(x="GlmTests", simplify="logical"),
function(x, simplify=TRUE) {
res <- sapply(x@score, function(x) sign(x$U))
res
})
setMethod("effect.sign",
signature(x="SingleSnpTestsScore", simplify="missing"),
function(x, simplify) {
res <- sign(x@U[,1])
names(res) <- x@snp.names
res
})
setMethod("sample.size", signature(x="GlmTests"),
function(x) {
res <- x@N
res
})
setMethod("sample.size", signature(x="SingleSnpTests"),
function(x) {
res <- x@N
names(res) <- x@snp.names
res
})
setMethod("effective.sample.size", signature(x="SingleSnpTests"),
function(x) {
if (length(x@N.r2)==0)
res <- x@N
else
res <- x@N.r2
names(res) <- x@snp.names
res
})
setMethod("names", signature(x="SingleSnpTests"), function(x) x@snp.names)
setMethod("names", "GlmTests",
function(x) {
objn <- x@snp.names
if (is.list(objn)) {
nobjn <- attr(objn, "names")
if (!is.null(nobjn))
nobjn
else
vapply(objn, function(x){
paste(x[1],x[length(x)], sep="...")}, "character")
}
else
as.character(objn)
})
setMethod("pool2",
signature(x="SingleSnpTestsScore",y="SingleSnpTestsScore",
score="logical"),
function(x, y, score) {
all.snps <- union(x@snp.names, y@snp.names)
can.pool <- intersect(x@snp.names, y@snp.names)
x.only <- setdiff(x@snp.names, can.pool)
y.only <- setdiff(y@snp.names, can.pool)
if (length(can.pool)>0) {
xsel <- match(can.pool, x@snp.names)
ysel <- match(can.pool, y@snp.names)
N <- x@N[xsel] + y@N[ysel]
U <- x@U[xsel,,drop=FALSE] + y@U[ysel,,drop=FALSE]
V <- x@V[xsel,,drop=FALSE] + y@V[ysel,,drop=FALSE]
if (length(x@N.r2)>0) {
if (length(y@N.r2)>0)
Nr2 <- x@N.r2[xsel] + y@N.r2[ysel]
else
Nr2 <- x@N.r2[xsel] + y@N[ysel]
} else {
if (length(y@N.r2)>0)
Nr2 <- x@N[xsel]+ y@N.r2[ysel]
else
Nr2 <- numeric(0)
}
} else {
N <- NULL
U <- NULL
V <- NULL
Nr2 <- numeric(0)
}
if (length(x.only)>0) {
xsel <- match(x.only, x@snp.names)
U <- rbind(U, x@U[xsel,,drop=FALSE])
V <- rbind(V, x@V[xsel,,drop=FALSE])
if (length(Nr2>0)) {
if (length(x@N.r2)>0)
Nr2 <- c(Nr2, x@N.r2[xsel])
else
Nr2 <- c(Nr2, x@N[xsel])
} else {
if (length(x@N.r2)>0)
Nr2 <- c(N, x@N.r2[xsel])
}
N <- c(N, x@N[xsel])
}
if (length(y.only)>0) {
ysel <- match(y.only, y@snp.names)
U <- rbind(U, y@U[ysel,,drop=FALSE])
V <- rbind(V, y@V[ysel,,drop=FALSE])
if (length(Nr2>0)) {
if (length(y@N.r2)>0)
Nr2 <- c(Nr2, y@N.r2[ysel])
else
Nr2 <- c(Nr2, y@N[ysel])
} else {
if (length(y@N.r2)>0)
Nr2 <- c(N, y@N.r2[ysel])
}
N <- c(N, y@N[ysel])
}
chisq <- .Call("chisq_single", list(U=U, V=V, N=N),
PACKAGE="snpStats")
if (score)
res <- new("SingleSnpTestsScore",
snp.names=c(can.pool, x.only, y.only),
chisq=chisq, N=N, U=U, V=V, N.r2=Nr2)
else
res <- new("SingleSnpTests",
snp.names=c(can.pool, x.only, y.only),
chisq=chisq, N=N, N.r2=Nr2)
res
})
# need to edit this
setMethod("pool2",
signature(x="GlmTestsScore",y="GlmTestsScore",
score="logical"),
function(x, y, score) {
if (!is.null(x@var.names) && !is.null(y@var.names) &&
(any(x@var.names!=y@var.names)))
warning("tests to be pooled have non-matching var.names slots")
nm.x <- names(x)
nm.y <- names(y)
if (any(duplicated(nm.x))||any(duplicated(nm.y)))
stop("At least one object has duplicated test names")
to.pool <- intersect(nm.x, nm.y)
if (length(to.pool)>0) {
res <- .Call("pool2_glm",
x[to.pool],
y[to.pool], score,
PACKAGE="snpStats")
}
else {
res <- NULL
}
x.only <- setdiff(nm.x, to.pool)
y.only <- setdiff(nm.y, to.pool)
if (length(x.only)==0 && length(y.only)==0)
return(res);
ix <- match(x.only, nm.x, nomatch=0)
iy <- match(y.only, nm.y, nomatch=0)
if (is.null(res)) {
res.chisq <- NULL
res.df <- NULL
res.N <- NULL
res.score <- NULL
}
else {
res.chisq <- res@chisq
res.df <- res@df
res.N <- res@N
if (score)
res.score <- res@score
}
if (score)
res <- new("GlmTestsScore",
snp.names=c(to.pool, x.only, y.only),
var.names=x@var.names,
chisq=c(res.chisq, x@chisq[ix], y@chisq[iy]),
df=c(res.df, x@df[ix], y@df[iy]),
N=c(res.N, x@N[ix], y@N[iy]),
score=append(res.score,
append(x@score[ix], y@score[iy])))
else
res <- new("GlmTests",
snp.names=c(to.pool, x.only, y.only),
var.names=x@var.names,
chisq=c(res.chisq, x@chisq[ix], y@chisq[iy]),
df=c(res.df, x@df[ix], y@df[iy]),
N=c(res.N, x@N[ix], y@N[iy]))
res
})
# Switch allele methods
setGeneric("switch.alleles",
function(x, snps) standardGeneric("switch.alleles"),
useAsDefault=FALSE)
setMethod("switch.alleles", signature(x="SnpMatrix", snps="ANY"),
function(x, snps) {
if (is.character(snps))
snps <- match(snps, colnames(x))
else if (is.logical(snps)) {
if (length(snps)!=ncol(x))
stop("logical snp selection vector is wrong length")
snps <- (1:ncol(x))[snps]
}
else if (!is.integer(snps))
stop("snp selection must be character, logical or integer")
if (any(is.na(snps) | snps>ncol(x) | snps<1))
stop("illegal snp selection")
.Call("smat_switch", x, snps, PACKAGE="snpStats")
})
setMethod("switch.alleles", signature(x="SingleSnpTestsScore", snps="ANY"),
function(x, snps) {
if (is.character(snps)) {
snps <- match(snps, x@snp.names)
if (any(is.na(snps)))
warning(sum(is.na(snps)),
" SNP names were not found in tests object")
snps <- snps[!is.na(snps)]
}
ntest <- length(x@snp.names)
if (is.logical(snps)) {
if (length(snps)!=ntest)
stop("incompatible arguments")
if (sum(snps)==0)
return(x)
} else if (is.numeric(snps)) {
if (length(snps)==0)
return(x)
if (max(snps)>ntest || min(snps)<1)
stop("incompatible arguments")
} else {
stop("SNPs to be switched must be indicated by name, position, or by a logical vector")
}
res <- x
res@U[snps,1] <- -x@U[snps,1]
if (ncol(x@U)==3) {
res@U[snps,2] <- -x@U[snps,2]
res@U[snps,3] <- x@U[snps,3] - x@U[snps,2]
res@V[snps,4] <- x@V[snps,4] - 2*x@V[snps,3] +
x@V[snps,2]
res@V[snps,3] <- x@V[snps,3] - x@V[snps,2]
} else {
res@U[snps,2] <- x@U[snps,2] - x@U[snps,1]
res@V[snps,3] <- x@V[snps,3] - 2*x@V[snps,2] +
x@V[snps,1]
res@V[snps,2] <- x@V[snps,2] - x@V[snps,1]
}
res
})
## This next will be slow but needed rarely
setMethod("switch.alleles", signature(x="GlmTestsScore",
snps="character"),
function(x, snps) {
if (is.list(x@snp.names))
switch <- lapply(x@snp.names, function(x, y) x%in%y, snps)
else
switch <- x@snp.names %in% snps
N <- length(x@score)
new.score <- vector("list", N)
for (i in 1:N) {
Ui <- x$U
Vi <- x$V
Si <- switch[i]
lu <- length(Ui)
ls <- length(Si)
if (lu!=ls) {
if (lu%%ls)
stop("error in GlmTest object")
Si <- rep(Si, rep(lu/ls, ls))
}
new.score[[i]] <- list(U=ifelse(Si, -Ui, Ui), V=Vi)
}
new("GlmTestsScore", snp.names=x@snp.names, var.names=x@var.names,
chisq=x@chisq, df=x@df, N=x@N,
score=new.score)
})
pool <- function(..., score=FALSE) {
argl <- list(...)
na <- length(argl)
while (na > 0 && is.null(argl[[na]])) {
argl <- argl[-na]
na <- na - 1
}
if (na<2)
stop("need at least two test sets to pool")
if (na>2) {
p2 <- pool2(..1, ..2, score=TRUE)
r <- do.call(pool, c(p2, argl[3:na], score=score))
}
else
r <- pool2(..1, ..2, score=score)
r
}
## Parameter estimates
setClass("GlmEstimates", contains="list")
setMethod("[", signature("GlmEstimates", i="ANY", j="missing",
drop="missing"),
function(x, i) {
if (is.character(i))
i <- match(i, names(x))
res <- x@.Data[i, drop=FALSE]
if (length(res)==0)
return(NULL)
names(res) <- names(x)[i]
new("GlmEstimates", res)
})
setMethod("show", "GlmEstimates",
function(object) {
len0 <- max(5, nchar(names(object)))
len1 <- max(10,
sapply(object, function(x)
max(if(is.null(x)) 0 else nchar(x$Y.var))))
len2 <- max(9,
sapply(object, function(x)
max(if(is.null(x)) 0 else nchar(names(x$beta)))))
space <- 2
nwidth <- 10
di <- 5
divide <- rep("-", len0+len1+len2+3*nwidth+5*space)
cat("\n", format("Model", justify="right", width=len0),
format("Y-variable", justify="right", width=len1+space),
format("Parameter", justify="right", width=len2+space),
format("Estimate", justify="right", width=nwidth+space),
format("S.E.", justify="right", width=nwidth+space),
format("z-value", justify="right", width=nwidth+space),"\n",
divide, "\n", sep="")
labs <- names(object)
for (j in 1:length(object)) {
labj <- labs[j]
x <- object[[j]]
if (is.null(x)) {
cat(format(labj, justify="right", width=len1),
"- No estimates available\n")
}
else {
nyj <- x$Y.var
nb <- length(x$beta)
namep <- names(x$beta)
ii <- 0
for (i in 1:nb){
ii <- ii+i
bi <- x$beta[i]
si <- sqrt(x$Var.beta[ii])
zi <- bi/si
cat(format(labj, justify="right", width=len0), sep="")
cat(format(nyj, justify="right", width=len1+space), sep="")
labj <- ""
nyj <- ""
cat(format(namep[i], justify="right", width=len2+space),
format(bi, justify="right", width=nwidth+space,
digits=di),
format(si, justify="right", width=nwidth+space,
digits=di),
formatC(zi, format="f", width=nwidth+space, digits=3),
"\n", sep="")
}
}
cat(divide, "\n", sep="")
}
}
)
## Wald test
setAs("GlmEstimates", "GlmTests",
function(from) {
.Call("wald", from, package="snpStats")
}
)
## To do
## names method for snp and X.snp
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.