R/util.R
In byandell/qtl: Tools for analyzing QTL experiments
#####################################################################
#
# util.R
#
# copyright (c) 2001-2012, Karl W Broman
# [find.pheno, find.flanking, and a modification to create.map
# from Brian Yandell]
# last modified Nov, 2012
# first written Feb, 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: markernames, c.cross, create.map,
# clean, clean.cross, drop.nullmarkers, nullmarkers
# drop.markers, pull.markers, drop.dupmarkers
# geno.table, genotab.em
# mf.k, mf.h, imf.k, imf.h, mf.cf, imf.cf, mf.m, imf.m,
# mf.stahl, imf.stahl
# switch.order, makeSSmap,
# subset.cross, fill.geno, checkcovar, find.marker,
# find.pseudomarker,
# adjust.rf.ri, lodint, bayesint,
# comparecrosses, movemarker, summary.map (aka summaryMap),
# print.summary.map, find.pheno,
# convert, convert.scanone, convert.scantwo
# find.flanking, strip.partials, comparegeno
# qtlversion, locateXO, jittermap, getid,
# find.markerpos, geno.crosstab, LikePheVector,
# matchchr, convert2sa, charround, testchr,
# scantwoperm2scanoneperm, subset.map, [.map, [.cross,
# findDupMarkers, convert2riself, convert2risib,
# switchAlleles, nqrank, cleanGeno, typingGap,
# calcPermPval, phenames
#
######################################################################
######################################################################
#
# markernames
#
# pull out the marker names for selected chromosomes as one big vector
#
######################################################################
markernames <-
function(cross, chr)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
if(!missing(chr)) cross <- subset(cross, chr=chr)
temp <- unlist(lapply(cross$geno, function(a) colnames(a$data)))
names(temp) <- NULL
temp
}
######################################################################
#
# chrnames
#
# pull out the chrnames for a cross
#
######################################################################
chrnames <-
function(cross)
{
names(cross$geno)
}
######################################################################
#
# create.map
#
# create a new map with inserted inter-marker locations
#
# Note: map is a vector or a matrix with 2 rows
#
# stepwidth = "fixed" is the original R/qtl version;
# stepwidth="variable" is for Brian Yandell and the qtlbim package
# stepwidth="max" creates the minimal number of inserted pseudomarkers
# to have the maximum stepwidth = step
######################################################################
create.map <-
function(map, step, off.end, stepwidth = c("fixed", "variable", "max"))
{
stepwidth <- match.arg(stepwidth)
if(step<0 || off.end<0) stop("step and off.end must be > 0.")
if(!is.matrix(map)) { # sex-ave map
if(stepwidth == "variable") {
if(off.end > 0) {
tmp <- names(map)
## Append and prepend by off.end value (exact here, no roundoff).
map <- c(map[1] - off.end, map, map[length(map)] + off.end)
names(map) <- c("loc000", tmp, "loc999")
}
if(step == 0)
return(unclass(map))
## Determine differences and expansion vector.
dif <- diff(map)
expand <- pmax(1, floor(dif / step))
## Create pseudomarker map.
a <- min(map) + cumsum(c(0, rep(dif / expand, expand)))
## Names are marker names or locNNN.
namesa <- paste("loc", seq(length(a)), sep = "")
namesa[cumsum(c(1, expand))] <- names(map)
names(a) <- namesa
return(unclass(a))
}
if(stepwidth == "max") {
if(off.end > 0) {
toadd <- c(map[1] - off.end, map[length(map)]+off.end)
if(step==0) {
names(toadd) <- paste("loc", 1:2, sep="")
map <- sort(c(map, toadd))
return(unclass(map))
}
nmap <- c(map[1] - off.end, map, map[length(map)]+off.end)
}
else {
nmap <- map
toadd <- NULL
}
if(step==0 || (length(map)==1 && off.end==0)) return(unclass(map))
d <- diff(nmap)
nadd <- ceiling(d/step)-1
if(sum(nadd) > 0) {
for(j in 1:(length(nmap)-1)) {
if(nadd[j]>0)
toadd <- c(toadd, seq(nmap[j], nmap[j+1], len=nadd[j]+2)[-c(1,nadd[j]+2)])
}
}
if(length(toadd) > 0) {
names(toadd) <- paste("loc", 1:length(toadd), sep="")
map <- sort(c(map, toadd))
}
return(unclass(map))
}
if(length(map) == 1) { # just one marker!
if(off.end==0) {
if(step == 0) step <- 1
nam <- names(map)
map <- c(map,map+step)
names(map) <- c(nam,paste("loc",step,sep=""))
}
else {
if(step==0) m <- c(-off.end,off.end)
else m <- seq(-off.end,off.end,by=step)
m <- m[m!=0]
names(m) <- paste("loc",m,sep="")
map <- sort(c(m+map,map))
}
return(map)
}
minloc <- min(map)
map <- map-minloc
if(step==0 && off.end==0) return(map+minloc)
else if(step==0 && off.end > 0) {
a <- c(floor(min(map)-off.end),ceiling(max(map)+off.end))
names(a) <- paste("loc", a, sep="")
return(sort(c(a,map))+minloc)
}
else if(step>0 && off.end == 0) {
a <- seq(floor(min(map)),max(map),
by = step)
if(any(is.na(match(a, map)))) {
a <- a[is.na(match(a,map))]
names(a) <- paste("loc",a,sep="")
return(sort(c(a,map))+minloc)
}
else return(map+minloc)
}
else {
a <- seq(floor(min(map)-off.end),ceiling(max(map)+off.end+step),
by = step)
a <- a[is.na(match(a,map))]
# no more than one point above max(map)+off.end
z <- (seq(along=a))[a >= max(map)+off.end]
if(length(z) > 1) a <- a[-z[-1]]
names(a) <- paste("loc",a,sep="")
return(sort(c(a,map))+minloc)
}
} # end sex-ave map
else { # sex-specific map
if(stepwidth == "variable") {
if(off.end > 0) {
tmp <- colnames(map)
map <- cbind(map[, 1] - off.end, map, map[, ncol(map)] + off.end)
dimnames(map) <- list(NULL, c("loc000", tmp, "loc999"))
}
if(step == 0)
return(unclass(map))
## Determine differences and expansion vector.
dif <- diff(map[1, ])
expand <- pmax(1, floor(dif / step))
## Create pseudomarker map.
a <- min(map[1, ]) + cumsum(c(0, rep(dif / expand, expand)))
b <- min(map[2, ]) + cumsum(c(0, rep(diff(map[2, ]) / expand, expand)))
namesa <- paste("loc", seq(length(a)), sep = "")
namesa[cumsum(c(1, expand))] <- dimnames(map)[[2]]
map <- rbind(a,b)
dimnames(map) <- list(NULL, namesa)
return(unclass(map))
}
if(stepwidth == "max") {
if(step==0 && off.end==0) return(unclass(map))
if(step==0 && off.end>0) {
if(ncol(map)==1) { # only one marker; assume equal recomb in sexes
L1 <- L2 <- 1
}
else {
L1 <- diff(range(map[1,]))
L2 <- diff(range(map[2,]))
}
nam <- colnames(map)
nmap1 <- c(map[1,1]-off.end, map[1,], map[1,ncol(map)]+off.end)
nmap2 <- c(map[2,1]-off.end*L2/L1, map[2,], map[2,ncol(map)]+off.end*L2/L1)
map <- rbind(nmap1, nmap2)
colnames(map) <- c("loc1", nam, "loc2")
return(unclass(map))
}
if(ncol(map)==1) L1 <- L2 <- 1
else {
L1 <- diff(range(map[1,]))
L2 <- diff(range(map[2,]))
}
nam <- colnames(map)
if(off.end > 0) {
toadd1 <- c(map[1,1] - off.end, map[1,ncol(map)]+off.end)
toadd2 <- c(map[2,1] + off.end*L2/L1, map[2,ncol(map)]+off.end*L2/L1)
neword <- order(c(map[1,], toadd1))
nmap1 <- c(map[1,], toadd1)[neword]
nmap2 <- c(map[2,], toadd2)[neword]
}
else {
nmap1 <- map[1,]
nmap2 <- map[2,]
toadd1 <- toadd2 <- NULL
}
d <- diff(nmap1)
nadd <- ceiling(d/step)-1
if(sum(nadd) > 0) {
for(j in 1:(length(nmap1)-1)) {
if(nadd[j]>0) {
toadd1 <- c(toadd1, seq(nmap1[j], nmap1[j+1], len=nadd[j]+2)[-c(1,nadd[j]+2)])
toadd2 <- c(toadd2, seq(nmap2[j], nmap2[j+1], len=nadd[j]+2)[-c(1,nadd[j]+2)])
}
}
}
newnam <- paste("loc", 1:length(toadd1), sep="")
toadd1 <- sort(toadd1)
toadd2 <- sort(toadd2)
neword <- order(c(map[1,], toadd1))
nmap1 <- c(map[1,], toadd1)[neword]
nmap2 <- c(map[2,], toadd2)[neword]
map <- rbind(nmap1, nmap2)
colnames(map) <- c(nam, newnam)[neword]
return(unclass(map))
}
minloc <- c(min(map[1,]),min(map[2,]))
map <- unclass(map-minloc)
markernames <- colnames(map)
if(step==0 && off.end==0) return(map+minloc)
else if(step==0 && off.end > 0) {
map <- map+minloc
if(ncol(map)==1) { # only one marker; assume equal recomb in sexes
L1 <- L2 <- 1
}
else {
L1 <- diff(range(map[1,]))
L2 <- diff(range(map[2,]))
}
nam <- colnames(map)
nmap1 <- c(map[1,1]-off.end, map[1,], map[1,ncol(map)]+off.end)
nmap2 <- c(map[2,1]-off.end*L2/L1, map[2,], map[2,ncol(map)]+off.end*L2/L1)
map <- rbind(nmap1, nmap2)
colnames(map) <- c("loc1", nam, "loc2")
return(map)
}
else if(step>0 && off.end == 0) {
if(ncol(map)==1) return(map+minloc)
a <- seq(floor(min(map[1,])),max(map[1,]),
by = step)
a <- a[is.na(match(a,map[1,]))]
if(length(a)==0) return(map+minloc)
b <- sapply(a,function(x,y,z) {
ZZ <- min((seq(along=y))[y > x])
(x-y[ZZ-1])/(y[ZZ]-y[ZZ-1])*(z[ZZ]-z[ZZ-1])+z[ZZ-1] }, map[1,],map[2,])
m1 <- c(a,map[1,])
m2 <- c(b,map[2,])
names(m1) <- names(m2) <- c(paste("loc",a,sep=""),markernames)
return(rbind(sort(m1),sort(m2))+minloc)
}
else {
a <- seq(floor(min(map[1,])-off.end),ceiling(max(map[1,])+off.end+step),
by = step)
a <- a[is.na(match(a,map[1,]))]
# no more than one point above max(map)+off.end
z <- (seq(along=a))[a >= max(map[1,])+off.end]
if(length(z) > 1) a <- a[-z[-1]]
b <- sapply(a,function(x,y,z,ml) {
if(x < min(y)) {
return(min(z) - (min(y)-x)/diff(range(y))*diff(range(z)) - ml)
}
else if(x > max(y)) {
return(max(z) + (x - max(y))/diff(range(y))*diff(range(z)) - ml)
}
else {
ZZ <- min((seq(along=y))[y > x])
(x-y[ZZ-1])/(y[ZZ]-y[ZZ-1])*(z[ZZ]-z[ZZ-1])+z[ZZ-1]
}
}, map[1,],map[2,], minloc[2])
m1 <- c(a,map[1,])
m2 <- c(b,map[2,])
names(m1) <- names(m2) <- c(paste("loc",a,sep=""),markernames)
return(rbind(sort(m1),sort(m2))+minloc)
}
}
}
######################################################################
# clean functions
######################################################################
clean <-
function(object, ...)
UseMethod("clean")
######################################################################
#
# clean.cross
#
# remove all of the extraneous stuff from a cross object, to get back
# to just the data
#
######################################################################
clean.cross <-
function(object, ...)
{
if(!any(class(object) == "cross"))
stop("Input should have class \"cross\".")
cross2 <- list(geno=object$geno,pheno=object$pheno)
if("cross" %in% names(object))
cross2$cross <- object$cross
if("founderGeno" %in% names(object))
cross2$founderGeno <- object$founderGeno
if(!is.null(attr(object, "alleles")))
attr(cross2, "alleles") <- attr(object, "alleles")
if(!is.null(attr(object, "scheme")))
attr(cross2, "scheme") <- attr(object, "scheme")
for(i in 1:length(object$geno)) {
cross2$geno[[i]] <- list(data=object$geno[[i]]$data,
map=object$geno[[i]]$map)
class(cross2$geno[[i]]) <- class(object$geno[[i]])
}
class(cross2) <- class(object)
cross2
}
######################################################################
#
# drop.qtlgeno
#
# remove any QTLs from the genotype data and the genetic maps
# from data simulated via sim.cross. (They all have names "QTL*")
#
######################################################################
#drop.qtlgeno <-
#function(cross)
#{
# n.chr <- nchr(cross)
# mar.names <- lapply(cross$geno, function(a) {
# m <- a$map
# if(is.matrix(m)) return(colnames(m))
# else return(names(m)) } )
#
# for(i in 1:n.chr) {
# o <- grep("^QTL[0-9]+",mar.names[[i]])
# if(length(o) != 0) {
# cross$geno[[i]]$data <- cross$geno[[i]]$data[,-o,drop=FALSE]
# if(is.matrix(cross$geno[[i]]$map))
# cross$geno[[i]]$map <- cross$geno[[i]]$map[,-o,drop=FALSE]
# else
# cross$geno[[i]]$map <- cross$geno[[i]]$map[-o]
# }
# }
# cross
#}
######################################################################
#
# drop.nullmarkers
#
# remove markers that have no genotype data from the data matrix and
# genetic maps
#
######################################################################
drop.nullmarkers <-
function(cross)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
n.chr <- nchr(cross)
keep.chr <- rep(TRUE,n.chr)
for(i in 1:n.chr) {
o <- !apply(cross$geno[[i]]$data,2,function(a) sum(!is.na(a)))
if(any(o)) { # remove from genotype data and map
mn.drop <- colnames(cross$geno[[i]]$data)[o]
if(length(mn.drop) == ncol(cross$geno[[i]]$data))
keep.chr[i] <- FALSE # removing all markers from this chromosome
cross$geno[[i]]$data <- cross$geno[[i]]$data[,!o,drop=FALSE]
if(is.matrix(cross$geno[[i]]$map))
cross$geno[[i]]$map <- cross$geno[[i]]$map[,!o,drop=FALSE]
else
cross$geno[[i]]$map <- cross$geno[[i]]$map[!o]
# results of calc.genoprob
if("prob" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$prob))
cross$geno[[i]]$prob <- cross$geno[[i]]$prob[,-o,,drop=FALSE]
}
# results of argmax.geno
if("argmax" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$argmax))
cross$geno[[i]]$argmax <- cross$geno[[i]]$argmax[,-o,drop=FALSE]
}
# results of sim.geno
if("draws" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$draws))
cross$geno[[i]]$draws <- cross$geno[[i]]$draws[,-o,,drop=FALSE]
}
# results of est.rf
if("rf" %in% names(cross)) {
attrib <- attributes(cross$rf)
o <- match(mn.drop,colnames(cross$rf))
cross$rf <- cross$rf[-o,-o]
if("onlylod" %in% names(attrib)) # save the onlylod attribute if its there
attr(cross$rf, "onlylod") <- attrib$onlylod
}
}
}
cross$geno <- cross$geno[keep.chr]
if("founderGeno" %in% names(cross))
cross$founderGeno <- cross$founderGeno[,markernames(cross)]
cross
}
######################################################################
#
# nullmarkers
#
# identify markers that have no genotype data from the data matrix and
# genetic maps
#
######################################################################
nullmarkers <-
function(cross)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
n.chr <- nchr(cross)
keep.chr <- rep(TRUE,n.chr)
all2drop <- NULL
for(i in 1:n.chr) {
o <- !apply(cross$geno[[i]]$data,2,function(a) sum(!is.na(a)))
if(any(o)) { # remove from genotype data and map
mn.drop <- colnames(cross$geno[[i]]$data)[o]
all2drop <- c(all2drop, mn.drop)
}
}
all2drop
}
######################################################################
#
# drop.markers
#
# remove a vector of markers from the data matrix and genetic maps
#
######################################################################
drop.markers <-
function(cross, markers)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
n.chr <- nchr(cross)
keep.chr <- rep(TRUE,n.chr)
found <- rep(FALSE, length(markers))
for(i in 1:n.chr) {
# find markers on this chromosome
o <- match(markers,colnames(cross$geno[[i]]$data))
found[!is.na(o)] <- TRUE
o <- o[!is.na(o)]
a <- rep(FALSE,ncol(cross$geno[[i]]$data))
a[o] <- TRUE
o <- a
if(any(o)) { # remove from genotype data and map
mn.drop <- colnames(cross$geno[[i]]$data)[o]
if(length(mn.drop) == ncol(cross$geno[[i]]$data))
keep.chr[i] <- FALSE # removing all markers from this chromosome
cross$geno[[i]]$data <- cross$geno[[i]]$data[,!o,drop=FALSE]
if(is.matrix(cross$geno[[i]]$map))
cross$geno[[i]]$map <- cross$geno[[i]]$map[,!o,drop=FALSE]
else
cross$geno[[i]]$map <- cross$geno[[i]]$map[!o]
# results of calc.genoprob
if("prob" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$prob))
cross$geno[[i]]$prob <- cross$geno[[i]]$prob[,-o,,drop=FALSE]
}
# results of argmax.geno
if("argmax" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$argmax))
cross$geno[[i]]$argmax <- cross$geno[[i]]$argmax[,-o,drop=FALSE]
}
# results of sim.geno
if("draws" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$draws))
cross$geno[[i]]$draws <- cross$geno[[i]]$draws[,-o,,drop=FALSE]
}
# results of est.rf
if("rf" %in% names(cross)) {
attrib <- attributes(cross$rf)
o <- match(mn.drop,colnames(cross$rf))
cross$rf <- cross$rf[-o,-o]
if("onlylod" %in% names(attrib))
attr(cross$rf, "onlylod") <- attrib$onlylod
}
}
}
if(any(!found))
warning("Markers not found: ", paste(markers[!found],collapse=" "))
cross$geno <- cross$geno[keep.chr]
if("founderGeno" %in% names(cross))
cross$founderGeno <- cross$founderGeno[,markernames(cross)]
cross
}
######################################################################
# pull.markers
#
# like drop.markers, but retain just those indicated
######################################################################
pull.markers <-
function(cross, markers)
{
mn <- markernames(cross)
m <- match(markers, mn)
if(any(is.na(m)))
warning("Some markers couldn't be found: ", paste(markers[is.na(m)], collapse=" "))
drop.markers(cross, mn[is.na(match(mn, markers))])
}
######################################################################
# drop.dupmarkers
#
# drop duplicate markers, retaining the consensus genotypes
######################################################################
drop.dupmarkers <-
function(cross, verbose=TRUE)
{
mn <- markernames(cross)
tab <- table(mn)
if(all(tab==1)) {
if(verbose) cat("No duplicate markers.\n")
return(cross)
}
dup <- names(tab[tab > 1])
if(verbose) cat(" ", length(dup), "duplicate markers\n")
# get consensus genotypes
g <- pull.geno(cross)[,!is.na(match(mn, dup))]
ng.omitted <- rep(NA, length(dup))
tot.omitted <- 0
nmar.omitted <- 0
for(i in seq(along=dup)) {
gg <- g[,colnames(g)==dup[i],drop=FALSE]
res <- apply(gg, 1, function(a)
{
if(all(is.na(a))) return(c(NA, 0))
a <- unique(a[!is.na(a)])
if(length(a)==1) return(c(a, 0))
return(c(NA, 1))
} )
if(verbose>1) {
cat(" ", dup[i], ":\t", sum(res[2,]), " genotypes omitted\n", sep="")
if(sum(res[2,]) > 1) cat(" ", paste(which(res[2,]>0), collapse=" "), "\n")
}
tot.omitted <- tot.omitted + sum(res[2,])
flag <- FALSE
for(j in seq(along=cross$geno)) {
mn <- colnames(cross$geno[[j]]$data)
wh <- mn==dup[i]
if(!any(wh)) next
wh <- which(wh)
if(!flag) {
flag <- TRUE
cross$geno[[j]]$data[,wh[1]] <- res[1,]
if(length(wh)==1) next
else wh <- wh[-1]
}
if(length(wh) > 0) {
nmar.omitted <- nmar.omitted + length(wh)
cross$geno[[j]]$data <- cross$geno[[j]]$data[,-wh,drop=FALSE]
if(is.matrix(cross$geno[[j]]$map))
cross$geno[[j]]$map <- cross$geno[[j]]$map[,-wh,drop=FALSE]
else
cross$geno[[j]]$map <- cross$geno[[j]]$map[-wh]
}
}
}
if(verbose) {
cat(" Total genotypes omitted:", tot.omitted, "\n")
cat(" Total markers omitted: ", nmar.omitted, "\n")
}
if("founderGeno" %in% names(cross))
cross$founderGeno <- cross$founderGeno[,markernames(cross)]
clean(cross)
}
######################################################################
#
# geno.table
#
# create table showing observed numbers of individuals with each
# genotype at each marker
#
######################################################################
geno.table <-
function(cross, chr, scanone.output=FALSE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
if(!missing(chr))
cross <- subset(cross, chr=chr)
n.chr <- nchr(cross)
type <- class(cross)[1]
is.bcs <- type == "bcsft"
cross.scheme <- attr(cross, "scheme")
if(is.bcs)
is.bcs <- (cross.scheme[2] == 0)
chrtype <- sapply(cross$geno, class)
allchrtype <- rep(chrtype, nmar(cross))
chrname <- names(cross$geno)
allchrname <- rep(chrname, nmar(cross))
## Actually plan to have our own geno.table.bcsft routine.
if(type == "f2" || (type == "bcsft" && !is.bcs)) {
n.gen <- 5
temp <- getgenonames("f2", "A", cross.attr=attributes(cross))
gen.names <- c(temp, paste("not", temp[c(3,1)]))
}
else if(type %in% c("bc", "riself", "risib", "dh", "bcsft")) {
n.gen <- 2
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
}
else if(type == "4way") {
n.gen <- 14
temp <- getgenonames("4way", "A", cross.attr=attributes(cross))
gen.names <- c(temp,
paste(temp[c(1,3)], collapse="/"),
paste(temp[c(2,4)], collapse="/"),
paste(temp[c(1,2)], collapse="/"),
paste(temp[c(3,4)], collapse="/"),
paste(temp[c(1,4)], collapse="/"),
paste(temp[c(2,3)], collapse="/"),
paste("not", temp[1]),
paste("not", temp[2]),
paste("not", temp[3]),
paste("not", temp[4]))
gen.names[5:8] <- substr(temp[c(1,2,1,3)], c(1,1,2,2), c(1,1,2,2))
}
else stop("Unknown cross type: ",type)
res <- lapply(cross$geno, function(a,ngen) {
a <- a$data; a[is.na(a)] <- 0
apply(a,2,function(b,ngen) table(factor(b,levels=0:ngen)),ngen)
},n.gen)
results <- NULL
for(i in 1:length(res))
results <- rbind(results,t(res[[i]]))
colnames(results) <- c("missing",gen.names)
rownames(results) <- unlist(lapply(cross$geno,function(a) colnames(a$data)))
pval <- rep(NA,nrow(results))
if(type %in% c("bc","risib","riself","dh") || (type=="bcsft" & is.bcs)) {
sexpgm <- getsex(cross)
if((type == "bc" || type=="bcsft") && any(chrtype == "X") && !is.null(sexpgm$sex) && any(sexpgm$sex==1)) {
for(i in which(allchrtype=="A")) {
x <- results[i,2:3]
if(sum(x) > 0)
pval[i] <- chisq.test(x,p=c(0.5,0.5))$p.value
else pval[i] <- 1
}
gn <- getgenonames("bc","X", "full", sexpgm, attributes(cross))
wh <- which(is.na(match(gn, colnames(results))))
temp <- matrix(0, nrow=nrow(results), ncol=length(wh))
colnames(temp) <- gn[wh]
results <- cbind(results, temp)
for(i in which(chrtype=="X")) {
dat <- reviseXdata("bc", "full", sexpgm, geno=cross$geno[[i]]$data,
cross.attr=attributes(cross))
dat[is.na(dat)] <- 0
tab <- t(apply(dat, 2, function(x) table(factor(x, levels=0:length(gn)))))
colnames(tab) <- c("missing", gn)
results[allchrname==chrname[i],] <- 0
results[allchrname==chrname[i],colnames(tab)] <- tab
for(j in 1:ncol(dat)) {
stat <- apply(table(sexpgm$sex, cross$geno[[i]]$data[,j]),1,
function(a) if(length(a) > 1 && sum(a)>0) return(chisq.test(a,p=c(0.5,0.5))$stat)
else return(0))
pval[allchrname==chrname[i]][j] <- 1-pchisq(sum(stat),length(stat))
}
}
results <- cbind(results, P.value=pval)
}
else {
for(i in 1:length(pval)) {
x <- results[i,2:3]
if(sum(x) > 0)
pval[i] <- chisq.test(x,p=c(0.5,0.5))$p.value
else
pval[i] <- 1
}
results <- cbind(results, P.value=pval)
}
}
else if(type == "f2" || (type == "bcsft" && !is.bcs)) {
sexpgm <- getsex(cross)
## F2 with set initial genotype probabilities.
init.geno <- c(0.25,0.5,0.25,0.75,0.75)
## BCsFt initial genotype probabilities need to be computed.
if(type == "bcsft") {
ret <- .C("genotab_em_bcsft",
as.integer(cross.scheme),
init.geno = as.double(init.geno))
init.geno <- ret$init.geno
}
for(i in which(allchrtype=="A")) {
dat <- results[i,2:6]
if(sum(dat)==0) pval[i] <- 1
else if(dat[4]==0 && dat[5]==0)
pval[i] <- chisq.test(dat[1:3], p=init.geno[1:3] )$p.value
else if(all(dat[2:4]==0))
pval[i] <- chisq.test(dat[c(1,5)],p=init.geno[c(1,5)])$p.value
else if(all(dat[c(1,2,5)]==0))
pval[i] <- chisq.test(dat[3:4], p=init.geno[3:4] )$p.value
else { # this is harder: some dominant and some not
pval[i] <- genotab.em(dat, init.geno)
}
}
for(i in which(chrtype=="X")) {
gn <- getgenonames("f2","X","full", getsex(cross), attributes(cross))
wh <- which(is.na(match(gn, colnames(results))))
temp <- matrix(0, nrow=nrow(results), ncol=length(wh))
colnames(temp) <- gn[wh]
results <- cbind(results, temp)
dat <- reviseXdata("f2", "full", sexpgm, geno=cross$geno[[i]]$data,
cross.attr=attributes(cross))
dat[is.na(dat)] <- 0
tab <- t(apply(dat, 2, function(x) table(factor(x, levels=0:length(gn)))))
colnames(tab) <- c("missing", gn)
results[allchrname==chrname[i],] <- 0
results[allchrname==chrname[i],colnames(tab)] <- tab
cn <- colnames(results)
f <- grep("f$", cn)
r <- grep("r$", cn)
if(length(f)>0 && length(r)>0) {
results <- results[,c(1:2,f,3,r,(1:ncol(results))[-c(1:3,f,r)])]
colnames(results) <- cn[c(1:2,f,3,r,(1:ncol(results))[-c(1:3,f,r)])]
}
sex <- sexpgm$sex
pgm <- sexpgm$pgm
for(j in 1:ncol(dat)) {
g <- cross$geno[[i]]$data[,j]
if(!is.null(sex)) {
if(!is.null(pgm))
g <- matrix(as.numeric(table(g,sex,pgm)),ncol=2,byrow=TRUE)
else
g <- matrix(as.numeric(table(g,sex)),ncol=2,byrow=TRUE)
}
else {
if(!is.null(pgm)) {
g <- matrix(as.numeric(table(g,pgm)),ncol=2,byrow=TRUE)
}
else g <- matrix(table(g),ncol=2,byrow=TRUE)
}
stat <- apply(g, 1,
function(a) if(sum(a)>0) return(chisq.test(a,p=c(0.5,0.5))$stat)
else return(0))
pval[allchrname==chrname[i]][j] <- 1-pchisq(sum(stat),length(stat))
}
}
results <- cbind(results, P.value=pval)
}
else if(type == "4way") {
for(i in 1:length(pval)) {
x <- results[i,2:5]
y <- results[i,-(1:5)]
if(sum(x) > 0 && sum(y)==0)
pval[i] <- chisq.test(x,p=c(0.25,0.25,0.25,0.25))$p.value
else {
if(allchrtype[i] == "A") {
res <- results[i,-1]
if(all(res[-c(1,11)]==0)) # AC/not AC
pval[i] <- chisq.test(res[c(1,11)], p=c(0.25, 0.75))$p.value
else if(all(res[-c(2,12)]==0)) # BC/not BC
pval[i] <- chisq.test(res[c(2,12)], p=c(0.25, 0.75))$p.value
else if(all(res[-c(3,13)]==0)) # AD/not AD
pval[i] <- chisq.test(res[c(3,13)], p=c(0.25, 0.75))$p.value
else if(all(res[-c(4,14)]==0)) # BD/not BD
pval[i] <- chisq.test(res[c(4,14)], p=c(0.25, 0.75))$p.value
else if(all(res[-c(5,6)]==0)) # A/B
pval[i] <- chisq.test(res[c(5,6)], p=c(0.5, 0.5))$p.value
else if(all(res[-c(7,8)]==0)) # C/D
pval[i] <- chisq.test(res[c(7,8)], p=c(0.5, 0.5))$p.value
else if(all(res[-c(9,10)]==0)) # AC/BD or AD/BC
pval[i] <- chisq.test(res[c(9,10)], p=c(0.5, 0.5))$p.value
else if(all(res[-c(2,4,5)]==0)) # BC/BD/A
pval[i] <- chisq.test(res[c(2,4,5)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(1,3,6)]==0)) # AC/AD/B
pval[i] <- chisq.test(res[c(1,3,6)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(3,4,7)]==0)) # AD/BD/C
pval[i] <- chisq.test(res[c(3,4,7)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(1,2,8)]==0)) # AC/BC/D
pval[i] <- chisq.test(res[c(1,2,8)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(2,3,9)]==0)) # AC/BD or AD or BC
pval[i] <- chisq.test(res[c(2,3,9)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(1,4,10)]==0)) # AD/BC or AC or BD
pval[i] <- chisq.test(res[c(1,4,10)], p=c(0.25, 0.25, 0.5))$p.value
}
}
}
results <- cbind(results, P.value=pval)
}
if(!scanone.output)
return(data.frame(chr=rep(names(cross$geno),nmar(cross)),results, stringsAsFactors=TRUE))
themap <- pull.map(cross)
if(is.matrix(themap[[1]]))
thepos <- unlist(lapply(themap, function(a) a[1,]))
else thepos <- unlist(themap)
temp <- results[,1:(ncol(results)-1),drop=FALSE]
res <- data.frame(chr=rep(names(cross$geno),nmar(cross)),
pos=thepos,
neglog10P=-log10(results[,ncol(results)]),
missing=temp[,1]/apply(temp, 1, sum),
temp[,-1]/apply(temp[,-1], 1, sum), stringsAsFactors=TRUE)
class(res) <- c("scanone", "data.frame")
rownames(res) <- rownames(results)
res[,1] <- factor(as.character(res[,1]), levels=unique(as.character(res[,1])))
res
}
genotab.em <-
function(dat, init.geno, tol=1e-6, maxit=10000, verbose=FALSE)
{
genotab.ll <-
function(dat, gam, init.geno)
{
p <- c(init.geno[1]*(1-gam[2]), init.geno[2]*gam[1], init.geno[3]*(1-gam[3]),
gam[2]*init.geno[4], gam[3]*init.geno[5])
if(any(p==0 & dat > 0)) return(-Inf)
return( sum((dat*log(p))[dat>0 & p>0]) )
}
n <- sum(dat)
gam <- c(sum(dat[1:3]), dat[4], dat[5])/n
curll <- genotab.ll(dat, gam, init.geno)
flag <- 0
if(verbose) cat(0, gam, curll, "\n")
for(i in 1:maxit) {
# estep
zAA <- dat[1]*gam[3]/(gam[1]+gam[3])
zBB <- dat[3]*gam[2]/(gam[1]+gam[2])
# mstep
gamnew <- c(sum(dat[1:3])-zAA-zBB, dat[4]+zBB, dat[5]+zAA)/n
newll <- genotab.ll(dat, gamnew, init.geno)
if(verbose) cat(i, gamnew, newll, "\n")
if(abs(curll-newll) < tol) {
flag <- 1
break
}
gam <- gamnew
curll <- newll
}
if(!flag) warning("Didn't converge.")
gam <- gamnew
p <- c(init.geno[1]*(1-gam[2]), init.geno[2]*gam[1], init.geno[3]*(1-gam[3]),
gam[2]*init.geno[4], gam[3]*init.geno[5])
ex <- p*n
1-pchisq(sum(((dat-ex)^2/ex)[ex>0]), 2)
}
######################################################################
# geno.crosstab
#
# Get a cross-tabulation of the genotypes at two markers,
# with the markers specified by name
######################################################################
geno.crosstab <-
function(cross, mname1, mname2, eliminate.zeros=TRUE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
if(missing(mname2) && length(mname1)>1) {
mname2 <- mname1[2]
mname1 <- mname1[1]
}
if(length(mname1) > 1 || length(mname2) > 1)
stop("mname1 and mname2 should both have lenght 1, or mname1 should have length 2 and mname1 should be missing.")
if(mname1==mname2)
stop("You must give two distinct marker names.")
pos <- find.markerpos(cross, c(mname1, mname2))
if(any(is.na(pos$chr))) {
if(all(is.na(pos$chr)))
stop("Markers ", mname1, " and ", mname2, " not found.")
else
stop("Marker ", rownames(pos)[is.na(pos$chr)], " not found.")
}
chrtype <- sapply(cross$geno[pos$chr], class)
crosstype <- class(cross)[1]
g1 <- pull.geno(cross, pos$chr[1])[,mname1, drop=FALSE]
g2 <- pull.geno(cross, pos$chr[2])[,mname2, drop=FALSE]
if(chrtype[1] == "X")
g1 <- reviseXdata(crosstype, "full", getsex(cross), geno=g1, cross.attr=attributes(cross))
if(chrtype[2] == "X")
g2 <- reviseXdata(crosstype, "full", getsex(cross), geno=g2, cross.attr=attributes(cross))
g1[is.na(g1)] <- 0
g2[is.na(g2)] <- 0
g1names <- getgenonames(crosstype, chrtype[1], "full", getsex(cross), attributes(cross))
g2names <- getgenonames(crosstype, chrtype[2], "full", getsex(cross), attributes(cross))
if(chrtype[1] != "X") {
if(crosstype == "f2")
g1names <- c(g1names, paste("not", g1names[c(3,1)]))
else if(crosstype == "bc" || crosstype == "risib" || crosstype=="riself" || crosstype=="dh") {
}
else if(crosstype == "4way") {
temp <- g1names
g1names <- c(temp,
paste(temp[c(1,3)], collapse="/"),
paste(temp[c(2,4)], collapse="/"),
paste(temp[c(1,2)], collapse="/"),
paste(temp[c(3,4)], collapse="/"),
paste(temp[c(1,4)], collapse="/"),
paste(temp[c(2,3)], collapse="/"),
paste("not", temp[1]),
paste("not", temp[2]),
paste("not", temp[3]),
paste("not", temp[4]))
g1names[5:8] <- substr(temp[c(1,2,1,3)], c(1,1,2,2), c(1,1,2,2))
}
else stop("Unknown cross type: ",crosstype)
}
if(chrtype[2] != "X") {
if(crosstype == "f2")
g2names <- c(g2names, paste("not", g2names[c(3,1)]))
else if(crosstype == "bc" || crosstype == "risib" || crosstype=="riself" || crosstype=="dh") {
}
else if(crosstype == "4way") {
temp <- g2names
g2names <- c(temp,
paste(temp[c(1,3)], collapse="/"),
paste(temp[c(2,4)], collapse="/"),
paste(temp[c(1,2)], collapse="/"),
paste(temp[c(3,4)], collapse="/"),
paste(temp[c(1,4)], collapse="/"),
paste(temp[c(2,3)], collapse="/"),
paste("not", temp[1]),
paste("not", temp[2]),
paste("not", temp[3]),
paste("not", temp[4]))
g2names[5:8] <- substr(temp[c(1,2,1,3)], c(1,1,2,2), c(1,1,2,2))
}
else stop("Unknown cross type: ",crosstype)
}
g1names <- c("-", g1names)
g2names <- c("-", g2names)
g1 <- as.character(g1)
g2 <- as.character(g2)
for(i in 1:length(g1names)) {
j <- as.character(i-1)
g1[g1==j] <- g1names[i]
}
g1 <- factor(g1, levels=g1names)
for(i in 1:length(g2names)) {
j <- as.character(i-1)
g2[g2==j] <- g2names[i]
}
g2 <- factor(g2, levels=g2names)
tab <- table(g1, g2)
names(attributes(tab)$dimnames) <- c(mname1, mname2)
if(eliminate.zeros) { # eliminate rows and columns with no data (but always include missing data row and column)
rs <- apply(tab, 1, sum); rs[1] <- 1
tab <- tab[rs>0,,drop=FALSE]
cs <- apply(tab, 2, sum); cs[1] <- 1
tab <- tab[,cs>0,drop=FALSE]
}
tab
}
# map functions
mf.k <- function(d) 0.5*tanh(d/50)
mf.h <- function(d) 0.5*(1-exp(-d/50))
imf.k <- function(r) { r[r >= 0.5] <- 0.5-1e-14; 50*atanh(2*r) }
imf.h <- function(r) { r[r >= 0.5] <- 0.5-1e-14; -50*log(1-2*r) }
mf.m <- function(d) sapply(d,function(a) min(a/100,0.5))
imf.m <- function(r) sapply(r,function(a) min(a*100,50))
# carter-falconer: mf.cf, imf.cf
imf.cf <- function(r) { r[r >= 0.5] <- 0.5-1e-14; 12.5*(log(1+2*r)-log(1-2*r))+25*atan(2*r) }
mf.cf <-
function(d)
{
d[d >= 300] <- 300
icf <- function(r,d)
imf.cf(r)-d
sapply(d,function(a) {
if(a==0) return(0)
uniroot(icf, c(0,0.5-1e-14),d=a,tol=1e-12)$root })
}
######################################################################
#
# switch.order: change the marker order on a given chromosome to some
# specified order
#
######################################################################
switch.order <-
function(cross, chr, order, error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"),
maxit=4000, tol=1e-6, sex.sp=TRUE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
map.function <- match.arg(map.function)
# check chromosome argument
if(missing(chr)) {
chr <- names(cross$geno)[1]
warning("Assuming you mean chromosome ", chr)
}
else {
if(length(chr) > 1) {
chr <- chr[1]
warning("switch.order can deal with just one chromosome at a time; assuming you want chr ", chr)
}
if(!testchr(chr, names(cross$geno)))
stop("Chr ", chr, " not found.")
}
chr <- matchchr(chr, names(cross$geno))
# check order argument
n.mar <- nmar(cross)
if(n.mar[chr] == length(order)-2 || n.mar[chr]==length(order)-1)
order <- order[1:n.mar[chr]] # useful for output from ripple()
if(n.mar[chr] != length(order))
stop("Incorrect number of markers.")
# save recombination fractions
flag <- 0
if("rf" %in% names(cross)) {
attrib <- attributes(cross$rf)
rf <- cross$rf
# determine column within rec fracs
whchr <- which(names(cross$geno)==chr)
oldcols <- cumsum(c(0,n.mar))[whchr]+seq(along=order)
newcols <- cumsum(c(0,n.mar))[whchr]+order
rf[oldcols,] <- rf[newcols,]
rf[,oldcols] <- rf[,newcols]
colnames(rf)[oldcols] <- colnames(rf)[newcols]
rownames(rf)[oldcols] <- rownames(rf)[newcols]
flag <- 1
}
# remove any intermediate calculations (except rec fracs),
# as they will no longer be meaningful
cross <- clean(cross)
if(!is.matrix(cross$geno[[chr]]$map))
first <- min(cross$geno[[chr]]$map)
else
first <- apply(cross$geno[[chr]]$map,1,min)
# re-order markers
cross$geno[[chr]]$data <- cross$geno[[chr]]$data[,order,drop=FALSE]
m <- seq(0,by=5,length=ncol(cross$geno[[chr]]$data))
names(m) <- colnames(cross$geno[[chr]]$data)
if(is.matrix(cross$geno[[chr]]$map))
cross$geno[[chr]]$map <- rbind(m,m)
else
cross$geno[[chr]]$map <- m
# re-estimate rec fracs for re-ordered chromosome
if(flag==1) {
temp <- est.rf(subset(cross, chr=chr))$rf
rf[oldcols,oldcols] <- temp
cross$rf <- rf
if("onlylod" %in% names(attrib))
attr(cross$rf, "onlylod") <- attrib$onlylod
}
# re-estimate map
newmap <- est.map(cross, chr=chr,
error.prob=error.prob, map.function=map.function,
maxit=maxit, tol=tol, sex.sp=sex.sp)
if(!is.matrix(newmap[[1]]))
cross$geno[[chr]]$map <- newmap[[1]] + first
else {
cross$geno[[chr]]$map[1,] <- newmap[[1]][1,] + first[1]
cross$geno[[chr]]$map[2,] <- newmap[[1]][2,] + first[2]
rownames(cross$geno[[chr]]$map) <- NULL
}
cross
}
######################################################################
#
# subset.cross: General subsetting function for a cross object
#
######################################################################
subset.cross <-
function(x, chr, ind, ...)
{
if(!any(class(x) == "cross"))
stop("Input should have class \"cross\".")
if(missing(chr) && missing(ind))
stop("You must specify either chr or ind.")
n.chr <- nchr(x)
n.ind <- nind(x)
# pull out relevant chromosomes
if(!missing(chr)) {
chr <- matchchr(chr, names(x$geno))
if("rf" %in% names(x)) { # pull out part of rec fracs
if(totmar(x) != ncol(x$rf))
x <- clean(x)
else {
attrib <- attributes(x$rf)
n.mar <- nmar(x)
n.chr <- n.chr
wh <- rbind(c(0,cumsum(n.mar)[-n.chr])+1,cumsum(n.mar))
dimnames(wh) <- list(NULL, names(n.mar))
wh <- wh[,chr,drop=FALSE]
wh <- unlist(apply(wh,2,function(a) a[1]:a[2]))
x$rf <- x$rf[wh,wh]
if("onlylod" %in% names(attrib)) # save the onlylod attribute if its there
attr(x$rf, "onlylod") <- attrib$onlylod
}
}
x$geno <- x$geno[chr]
if("founderGeno" %in% names(x))
x$founderGeno <- x$founderGeno[,unlist(lapply(x$geno, function(a) colnames(a$data)))]
}
if(!missing(ind)) {
theid <- getid(x)
if(is.logical(ind)) {
ind[is.na(ind)] <- FALSE
if(length(ind) != n.ind)
stop("ind argument has wrong length (", length(ind), "; should be ", n.ind, ")")
ind <- (1:n.ind)[ind]
}
else if(is.numeric(ind)) { # treat as numeric indices; don't match against individual identifiers
if(all(ind < 0)) { # drop all but these
ind <- -ind
if(any(ind > nind(x))) {
a <- -ind[ind > nind(x)]
if(length(a) > 1) a <- sample(a, 1)
stop("Invalid ind values (e.g., ", a, ")")
}
ind <- (1:n.ind)[-ind]
}
else if(all(ind > 0)) { # keep these
if(any(ind > nind(x))) {
a <- ind[ind > nind(x)]
if(length(a) > 1) a <- sample(a, 1)
stop("Invalid ind values (e.g., ", a, ")")
}
}
else {
stop("Need ind to be all > 0 or all < 0.")
}
}
else if(!is.null(theid)) { # cross has individual IDs
ind <- as.character(ind)
if(all(substr(ind, 1,1) == "-")) {
ind <- substr(ind, 2, nchar(ind))
m <- match(ind, theid)
if(all(is.na(m)))
stop("No matching individuals.")
if(any(is.na(m)))
warning("Individuals not found: ", paste(ind[is.na(m)]))
ind <- (1:n.ind)[-m[!is.na(m)]]
}
else {
m <- match(ind, theid)
if(any(is.na(m)))
warning("Individuals not found: ", paste(ind[is.na(m)], collapse=" "))
ind <- m[!is.na(m)]
}
ind <- ind[!is.na(ind)]
}
else { # no individual IDs
stop("In the absense of individual IDs, ind should be logical or numeric.")
}
# Note: ind should now be a numeric vector
if(length(ind) == 0)
stop("Must retain at least one individual.")
if(length(ind) == 1)
warning("Retained only one individual!")
x$pheno <- x$pheno[ind,,drop=FALSE]
if("cross" %in% names(x))
x$cross <- x$cross[ind,,drop=FALSE]
for(i in 1:nchr(x)) {
x$geno[[i]]$data <- x$geno[[i]]$data[ind,,drop=FALSE]
if("prob" %in% names(x$geno[[i]])) {
temp <- attributes(x$geno[[i]]$prob) # all attributes but dim and dimnames
x$geno[[i]]$prob <- x$geno[[i]]$prob[ind,,,drop=FALSE]
# put attributes back in
for(k in seq(along=temp)) {
if(names(temp)[k] != "dim" && names(temp)[k] != "dimnames")
attr(x$geno[[i]]$prob, names(temp)[k]) <- temp[[k]]
}
}
if("errorlod" %in% names(x$geno[[i]])) {
temp <- attributes(x$geno[[i]]$prob) # all attributes but dim and dimnames
x$geno[[i]]$errorlod <- x$geno[[i]]$errorlod[ind,,drop=FALSE]
# put attributes back in
for(k in seq(along=temp)) {
if(names(temp)[k] != "dim" && names(temp)[k] != "dimnames")
attr(x$geno[[i]]$errorlod, names(temp)[k]) <- temp[[k]]
}
}
if("argmax" %in% names(x$geno[[i]])) {
temp <- attributes(x$geno[[i]]$argmax) # all attributes but dim and dimnames
x$geno[[i]]$argmax <- x$geno[[i]]$argmax[ind,,drop=FALSE]
# put attributes back in
for(k in seq(along=temp)) {
if(names(temp)[k] != "dim" && names(temp)[k] != "dimnames")
attr(x$geno[[i]]$argmax, names(temp)[k]) <- temp[[k]]
}
}
if("draws" %in% names(x$geno[[i]])) {
temp <- attributes(x$geno[[i]]$draws) # all attributes but dim and dimnames
x$geno[[i]]$draws <- x$geno[[i]]$draws[ind,,,drop=FALSE]
# put attributes back in
for(k in seq(along=temp)) {
if(names(temp)[k] != "dim" && names(temp)[k] != "dimnames")
attr(x$geno[[i]]$draws, names(temp)[k]) <- temp[[k]]
}
}
}
if("qtlgeno" %in% names(x))
x$qtlgeno <- x$qtlgeno[ind,,drop=FALSE]
}
x
}
#pull.chr <-
#function(cross, chr) {
# warning("pull.chr is deprecated; use subset.cross.")
# subset.cross(cross, chr)
#}
######################################################################
#
# c.cross: Combine crosses
#
######################################################################
c.cross <-
function(...)
{
args <- list(...)
n.args <- length(args)
for(i in seq(along=args)) {
if(!any(class(args[[i]]) == "cross"))
stop("Input should have class \"cross\".")
}
# if only one cross, just return it
if(n.args==1) return(args[[1]])
if(any(sapply(args, function(a) !any(class(a) == "cross"))))
stop("All arguments must be cross objects.")
# crosses must be all the same, or must be combination of F2 and BC
classes <- sapply(args,function(a) class(a)[1])
if(length(unique(classes))==1) {
allsame <- TRUE
type <- classes[1]
}
else {
if(any(classes != "bc" & classes != "f2"))
stop("Experiments must be either the same type or be bc/f2.")
allsame <- FALSE
type <- "f2"
}
if(length(unique(sapply(args, nchr))) > 1)
stop("All arguments must have the same number of chromosomes.")
x <- args[[1]]
chr <- names(x$geno)
n.mar <- nmar(x)
marnam <- unlist(lapply(x$geno,function(b) colnames(b$data)))
marpos <- unlist(lapply(x$geno,function(b) b$map))
map.mismatch <- 0
for(i in 2:n.args) {
y <- args[[i]]
y.marnam <- unlist(lapply(y$geno, function(b) colnames(b$data)))
y.marpos <- unlist(lapply(y$geno, function(b) b$map))
if(chr != names(y$geno) || any(n.mar != nmar(y)) ||
any(marnam != y.marnam) || any(marpos != y.marpos)) {
map.mismatch <- 1
break
}
}
if(map.mismatch) { # get the maps to line up
for(i in 1:nchr(args[[1]])) {
themap <- NULL
themaps <- vector("list", n.args)
for(j in 1:n.args) {
if(is.matrix(args[[j]]$map))
stop("c.cross() won't work with sex-specific maps.")
themaps[[j]] <- args[[j]]$geno[[i]]$map
themap <- c(themap, themaps[[j]])
}
themap <- sort(themap)
mn <- unique(names(themap))
newmap <- rep(0,length(mn))
names(newmap) <- mn
for(j in 1:length(newmap))
newmap[j] <- mean(themap[names(themap) == mn[j]])
for(j in 1:n.args) {
m <- match(names(themaps[[j]]), mn)
m <- m[!is.na(m)]
if(any(diff(m)) < 0)
stop(" Markers must all be in the same order.")
if(!all(mn %in% names(themaps[[j]]))) {
temp <- matrix(ncol=length(mn), nrow=nind(args[[j]]))
colnames(temp) <- mn
temp[,names(themaps[[j]])] <- args[[j]]$geno[[i]]$data
args[[j]]$geno[[i]]$data <- temp
}
args[[j]]$geno[[i]]$map <- newmap
}
}
} # end of map mismatch fix
x <- args[[1]]
chr <- names(x$geno)
n.mar <- nmar(x)
marnam <- unlist(lapply(x$geno,function(b) colnames(b$data)))
# create genotype information
geno <- x$geno
for(j in 1:nchr(x)) { # drop extraneous stuff
geno[[j]] <- list(data=geno[[j]]$data, map=geno[[j]]$map)
class(geno[[j]]) <- class(x$geno[[j]])
}
for(i in 2:n.args)
for(j in 1:nchr(x))
geno[[j]]$data <- rbind(geno[[j]]$data,args[[i]]$geno[[j]]$data)
# get all phenotype names
phenam <- names(x$pheno)
for(i in 2:n.args)
phenam <- c(phenam, names(args[[i]]$pheno))
phenam <- unique(phenam)
# form big phenotype matrix
n.ind <- sapply(args,nind)
pheno <- matrix(nrow=sum(n.ind),ncol=length(phenam))
colnames(pheno) <- phenam
pheno <- as.data.frame(pheno, stringsAsFactors=TRUE)
if(!allsame) {
crosstype <- factor(rep(c("bc","f2")[match(classes,c("bc","f2"))],n.ind),
levels=c("bc","f2"))
pheno <- cbind(pheno,crosstype=crosstype)
}
for(i in 1:length(phenam)) {
phe <- vector("list",n.args)
for(j in 1:n.args) {
o <- match(phenam[i],names(args[[j]]$pheno))
if(is.na(o)) phe[[j]] <- rep(NA,n.ind[j])
else phe[[j]] <- args[[j]]$pheno[,o]
}
pheno[,i] <- unlist(phe)
}
# indicator of which cross
whichcross <- matrix(0,ncol=n.args,nrow=sum(n.ind))
colnames(whichcross) <- paste("cross",1:n.args,sep="")
thecross <- rep(NA, sum(n.ind))
prev <- 0
for(i in 1:n.args) {
wh <- prev + 1:n.ind[i]
prev <- prev + n.ind[i]
whichcross[wh,i] <- 1
thecross[wh] <- i
}
pheno <- cbind(pheno,cross=thecross,whichcross)
x$pheno <- pheno
if(!map.mismatch) { # keep probs and draws only if we've not re-aligned the maps
# if probs exist in each and all have the same
# set up values, keep them
wh <- sapply(args, function(a) match("prob",names(a$geno[[1]])))
step <- sapply(args,function(a) attr(a$geno[[1]]$prob,"step"))
error.prob <- sapply(args,function(a) attr(a$geno[[1]]$prob,"error.prob"))
off.end <- sapply(args,function(a) attr(a$geno[[1]]$prob,"off.end"))
map.function <- sapply(args,function(a) attr(a$geno[[1]]$prob,"map.function"))
map <- lapply(args,function(a) attr(a$geno[[1]]$prob,"map"))
if(!any(is.na(wh)) && length(unique(step))==1 &&
length(unique(error.prob))==1 && length(unique(off.end))==1 &&
length(unique(map.function))==1) {
if(allsame) { # all same cross type
for(j in 1:nchr(x)) {
geno[[j]]$prob <- array(dim=c(sum(n.ind),dim(x$geno[[j]]$prob)[-1]))
dimnames(geno[[j]]$prob) <- dimnames(x$geno[[j]]$prob)
prev <- 0
for(i in 1:n.args) {
wh <- prev + 1:n.ind[i]
prev <- prev + n.ind[i]
geno[[j]]$prob[wh,,] <- args[[i]]$geno[[j]]$prob
}
}
}
else { # mixed F2 and BC
for(j in 1:nchr(x)) {
wh <- match("f2",classes)
geno[[j]]$prob <- array(0,dim=c(sum(n.ind),dim(args[[wh]]$geno[[j]]$prob)[-1]))
dimnames(geno[[j]]$prob) <- dimnames(args[[wh]]$geno[[j]]$prob)
prev <- 0
for(i in 1:n.args) {
wh <- prev + 1:n.ind[i]
prev <- prev + n.ind[i]
if(classes[i]=="f2")
geno[[j]]$prob[wh,,] <- args[[i]]$geno[[j]]$prob
else # backcross
geno[[j]]$prob[wh,,1:2] <- args[[i]]$geno[[j]]$prob
}
}
}
for(j in 1:nchr(x)) {
wh <- sapply(args, function(a) match("prob",names(a$geno[[j]])))
step <- sapply(args,function(a) attr(a$geno[[j]]$prob,"step"))
error.prob <- sapply(args,function(a) attr(a$geno[[j]]$prob,"error.prob"))
off.end <- sapply(args,function(a) attr(a$geno[[j]]$prob,"off.end"))
map.function <- sapply(args,function(a) attr(a$geno[[j]]$prob,"map.function"))
map <- lapply(args,function(a) attr(a$geno[[j]]$prob,"map"))
attr(geno[[j]]$prob,"step") <- step[1]
attr(geno[[j]]$prob,"error.prob") <- error.prob[1]
attr(geno[[j]]$prob,"off.end") <- off.end[1]
attr(geno[[j]]$prob,"map.function") <- map.function[1]
attr(geno[[j]]$prob,"map") <- map[[1]]
}
}
# if draws exist in each and all have the same
# set up values, keep them
wh <- sapply(args, function(a) match("draws",names(a$geno[[1]])))
step <- sapply(args,function(a) attr(a$geno[[1]]$draws,"step"))
error.prob <- sapply(args,function(a) attr(a$geno[[1]]$draws,"error.prob"))
off.end <- sapply(args,function(a) attr(a$geno[[1]]$draws,"off.end"))
map.function <- sapply(args,function(a) attr(a$geno[[1]]$draws,"map.function"))
map <- lapply(args,function(a) attr(a$geno[[1]]$draws,"map"))
ndraws <- sapply(args,function(a) dim(a$geno[[1]]$draws)[3])
if(!any(is.na(wh)) && length(unique(step))==1 &&
length(unique(error.prob))==1 && length(unique(off.end))==1 &&
length(unique(map.function))==1 && length(unique(ndraws))==1) {
for(j in 1:nchr(x)) {
geno[[j]]$draws <- array(0,dim=c(sum(n.ind),dim(x$geno[[j]]$draws)[-1]))
dimnames(geno[[j]]$draws) <- dimnames(x$geno[[j]]$draws)
prev <- 0
for(i in 1:n.args) {
wh <- prev + 1:n.ind[i]
prev <- prev + n.ind[i]
geno[[j]]$draws[wh,,] <- args[[i]]$geno[[j]]$draws
}
wh <- sapply(args, function(a) match("draws",names(a$geno[[j]])))
step <- sapply(args,function(a) attr(a$geno[[j]]$draws,"step"))
error.prob <- sapply(args,function(a) attr(a$geno[[j]]$draws,"error.prob"))
off.end <- sapply(args,function(a) attr(a$geno[[j]]$draws,"off.end"))
map.function <- sapply(args,function(a) attr(a$geno[[j]]$draws,"map.function"))
map <- lapply(args,function(a) attr(a$geno[[j]]$draws,"map"))
attr(geno[[j]]$draws,"step") <- step[1]
attr(geno[[j]]$draws,"error.prob") <- error.prob[1]
attr(geno[[j]]$draws,"off.end") <- off.end[1]
attr(geno[[j]]$draws,"map.function") <- map.function[1]
attr(geno[[j]]$draws,"map") <- map[[1]]
}
}
}
x <- list(geno=geno, pheno=pheno)
class(x) <- c(type,"cross")
x
}
######################################################################
#
# fill.geno: Run argmax.geno or sim.geno and then fill in the
# genotype data with the results. This will allow
# rough genome scans by marker regression without
# holes. WE WOULD NOT PLACE ANY TRUST IN THE RESULTS!
#
# the newer method, "no_dbl_XO", fills in missing genotypes between
# markers with matching genotypes
#
######################################################################
fill.geno <-
function(cross, method=c("imp","argmax", "no_dbl_XO"), error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"))
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
method <- match.arg(method)
# don't let error.prob be exactly zero (or >1)
if(error.prob < 1e-50) error.prob <- 1e-50
if(error.prob > 1) {
error.prob <- 1-1e-50
warning("error.prob shouldn't be > 1!")
}
# remove any extraneous material
cross <- clean(cross)
n.chr <- nchr(cross)
n.mar <- nmar(cross)
if(method=="imp") {
# do one imputation
temp <- sim.geno(cross,n.draws=1,step=0,off.end=0,
error.prob=error.prob,map.function=map.function)
# replace the genotype data with the results,
# stripping off any attributes
for(i in 1:n.chr) {
nam <- colnames(cross$geno[[i]]$data)
if(n.mar[i] == 1)
cross$geno[[i]]$data <-
matrix(as.numeric(temp$geno[[i]]$draws[,2,1]),ncol=n.mar[i])
else
cross$geno[[i]]$data <-
matrix(as.numeric(temp$geno[[i]]$draws[,,1]),ncol=n.mar[i])
colnames(cross$geno[[i]]$data) <- nam
}
}
else if(method=="argmax") {
# run the Viterbi algorithm
temp <- argmax.geno(cross,step=0,off.end=0,error.prob=error.prob,
map.function=map.function)
# replace the genotype data with the results,
# stripping off any attributes
for(i in 1:n.chr) {
nam <- colnames(cross$geno[[i]]$data)
if(n.mar[i] == 1)
cross$geno[[i]]$data <-
matrix(as.numeric(temp$geno[[i]]$argmax[,2]),ncol=n.mar[i])
else
cross$geno[[i]]$data <-
matrix(as.numeric(temp$geno[[i]]$argmax),ncol=n.mar[i])
colnames(cross$geno[[i]]$data) <- nam
}
}
else {
for(i in 1:n.chr) {
dat <- cross$geno[[i]]$data
dat[is.na(dat)] <- 0
nr <- nrow(dat)
nc <- ncol(dat)
dn <- dimnames(dat)
dat <- .C("R_fill_geno_nodblXO",
as.integer(nr),
as.integer(nc),
dat=as.integer(dat),
PACKAGE="qtl")$dat
dat[dat==0] <- NA
dat <- matrix(dat, ncol=nc, nrow=nr)
dimnames(dat) <- dn
cross$geno[[i]]$data <- dat
}
}
cross
}
######################################################################
#
# checkcovar
#
# This is a utility function for scanone and scantwo. We remove
# individuals with missing phenotypes or covariates and check
# that the covariates are of the appropriate form.
#
######################################################################
checkcovar <-
function(cross, pheno.col, addcovar, intcovar, perm.strata, ind.noqtl=NULL, weights=NULL, verbose=TRUE)
{
chrtype <- sapply(cross$geno, class)
# drop individuals whose sex or pgm is missing if X chr is included
if(any(chrtype=="X")) {
sexpgm <- getsex(cross)
keep <- rep(TRUE,nind(cross))
flag <- 0
if(!is.null(sexpgm$sex)) {
if(any(is.na(sexpgm$sex))) {
keep[is.na(sexpgm$sex)] <- FALSE
flag <- 1
}
}
if(!is.null(sexpgm$pgm)) {
if(any(is.na(sexpgm$pgm))) {
keep[is.na(sexpgm$pgm)] <- FALSE
flag <- 1
}
}
if(flag) {
if(verbose) warning("Dropping ", sum(!keep), " individuals with missing sex or pgm.\n")
cross <- subset(cross, ind=keep)
if(!is.null(addcovar)) {
if(!is.matrix(addcovar)) addcovar <- addcovar[keep]
else addcovar <- addcovar[keep,]
}
if(!is.null(intcovar)) {
if(!is.matrix(intcovar)) intcovar <- intcovar[keep]
else intcovar <- intcovar[keep,]
}
}
}
# check phenotypes - we allow multiple phenotypes here
if(pheno.col < 1 || pheno.col > nphe(cross))
stop("Specified phenotype column is invalid.")
# check if all phenotypes are numeric
pheno <- cross$pheno[,pheno.col,drop=FALSE]
idx.nonnum <- which(!apply(pheno,2, is.numeric))
if(length(idx.nonnum) > 0)
stop("Following phenotypes are not numeric: Column ",
paste(idx.nonnum, collapse=","))
orig.n.ind <- nind(cross)
# drop individuals with missing phenotypes
if(any(!is.finite(unlist(pheno)))) {
keep.ind <- as.numeric(which(apply(pheno, 1, function(x) !any(!is.finite(x)))))
# keep.ind <- (1:length(pheno))[!is.na(pheno)]
n.drop <- nind(cross) - length(keep.ind)
keep.ind.boolean <- rep(FALSE, nind(cross))
keep.ind.boolean[keep.ind] <- TRUE
cross <- subset.cross(cross,ind=keep.ind.boolean)
pheno <- pheno[keep.ind,,drop=FALSE]
if(verbose) warning("Dropping ", n.drop, " individuals with missing phenotypes.\n")
}
else keep.ind <- 1:nind(cross)
n.ind <- nind(cross)
n.chr <- nchr(cross) # number of chromosomes
type <- class(cross)[1] # type of cross
n.addcovar <- n.intcovar <- 0
if(!is.null(addcovar)) { # for additive covariates
if(!is.matrix(addcovar)) {
if(!is.numeric(as.matrix(addcovar)))
stop("addcovar should be numeric")
if(is.vector(addcovar) || is.data.frame(addcovar))
addcovar <- as.matrix(addcovar)
else stop("addcovar should be a matrix")
}
if(!all(apply(addcovar,2,is.numeric)))
stop("All columns of addcovar must be numeric")
if( nrow(addcovar) != orig.n.ind ) {
# the length of additive covariates is incorrect
stop("Number of rows in additive covariates is incorrect")
}
addcovar <- addcovar[keep.ind,,drop=FALSE]
n.addcovar <- ncol(addcovar)
}
if(!is.null(intcovar)) { # interacting covariates
if(!is.matrix(intcovar)) {
if(!is.numeric(as.matrix(intcovar)))
stop("intcovar should be a numeric")
if(is.vector(intcovar) || is.data.frame(intcovar))
intcovar <- as.matrix(intcovar)
else stop("intcovar should be a matrix")
}
if(!all(apply(intcovar,2,is.numeric)))
stop("All columns of intcovar must be numeric")
if(nrow(intcovar)[1] != orig.n.ind) {
# the length of interacting covariates is incorrect
stop("The length of interacting covariates is incorrect!")
}
intcovar <- intcovar[keep.ind,,drop=FALSE]
n.intcovar <- ncol(intcovar)
}
if(!is.null(perm.strata)) perm.strata <- perm.strata[keep.ind]
if(!is.null(ind.noqtl)) ind.noqtl <- ind.noqtl[keep.ind]
if(!is.null(weights)) weights <- weights[keep.ind]
# drop individuals missing any covariates
if(!is.null(addcovar)) { # note that intcovar is contained in addcovar
wh <- apply(cbind(addcovar,intcovar),1,function(a) any(!is.finite(a)))
if(any(wh)) {
cross <- subset.cross(cross,ind=(!wh))
pheno <- pheno[!wh,,drop=FALSE]
addcovar <- addcovar[!wh,,drop=FALSE]
if(!is.null(intcovar)) intcovar <- intcovar[!wh,,drop=FALSE]
n.ind <- nind(cross)
if(!is.null(perm.strata)) perm.strata <- perm.strata[!wh]
if(!is.null(ind.noqtl)) ind.noqtl <- ind.noqtl[!wh]
if(!is.null(weights)) weights <- weights[!wh]
if(verbose) warning("Dropping ", sum(wh), " individuals with missing covariates.\n")
}
}
# make sure columns of intcovar are contained in addcovar
if(!is.null(intcovar)) {
if(is.null(addcovar)) {
addcovar <- intcovar
n.addcovar <- n.intcovar
if(verbose) warning("addcovar forced to contain all columns of intcovar\n")
}
else {
wh <- 1:n.intcovar
for(i in 1:n.intcovar) {
o <- (apply(addcovar,2,function(a,b) max(abs(a-b)),intcovar[,i])<1e-14)
if(any(o)) wh[i] <- (1:n.addcovar)[o]
else wh[i] <- NA
}
if(any(!is.finite(wh))) {
addcovar <- cbind(addcovar,intcovar[,!is.finite(wh)])
n.addcovar <- ncol(addcovar)
if(verbose) warning("addcovar forced to contain all columns of intcovar")
}
}
}
if(n.addcovar > 0) { # check rank
if(qr(cbind(1,addcovar))$rank < n.addcovar+1)
if(verbose) warning("addcovar appears to be over-specified; consider dropping columns.\n")
}
if(n.intcovar > 0) { # check rank
if(qr(cbind(1,intcovar))$rank < n.intcovar+1)
if(verbose) warning("intcovar appears to be over-specified; consider dropping columns.\n")
}
pheno <- as.matrix(pheno)
list(cross=cross, pheno=pheno, addcovar=addcovar, intcovar=intcovar,
n.addcovar=n.addcovar, n.intcovar=n.intcovar, perm.strata=perm.strata,
ind.noqtl=ind.noqtl, weights=weights)
}
# Find the nearest marker to a particular position
find.marker <-
function(cross, chr, pos, index)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
if(missing(pos) && missing(index))
stop("Give either pos or index.")
if(!missing(pos) && !missing(index))
stop("Give just one of pos or index.")
# if chr has length 1, expand if necessary
if(length(chr) == 1) {
if(!missing(pos))
chr <- rep(chr,length(pos))
else
chr <- rep(chr, length(index))
}
# otherwise, chr and pos should have same length
else {
if(!missing(pos) && length(chr) != length(pos))
stop("chr and pos must be the same length.")
if(!missing(index) && length(chr) != length(index))
stop("chr and index must be the same length.")
}
markers <- rep("",length(chr))
chrnotfound <- NULL
for(i in 1:length(chr)) {
# find chromosome
o <- match(chr[i], names(cross$geno))
if(is.na(o)) {
markers[i] <- NA # chr not matched
chrnotfound <- c(chrnotfound, chr[i])
}
else {
thismap <- cross$geno[[o]]$map # genetic map
# sex-specific map; look at female positions
if(is.matrix(thismap)) thismap <- thismap[1,]
if(!missing(pos)) {
# find closest marker
d <- abs(thismap-pos[i])
o2 <- (1:length(d))[d==min(d)]
if(length(o2)==1) markers[i] <- names(thismap)[o2]
# if multiple markers are equidistant,
# choose the one with the most data
# or choose among them at random
else {
x <- names(thismap)[o2]
n.geno <- apply(cross$geno[[o]]$data[,o2],2,function(a) sum(!is.na(a)))
o2 <- o2[n.geno==max(n.geno)]
if(length(o2) == 1)
markers[i] <- names(thismap)[o2]
else
markers[i] <- names(thismap)[sample(o2,1)]
}
}
else { # by index
if(index[i] < 1 || index[i] > length(thismap))
stop("Misspecified index ", index[i], " on chr ", chr[i])
markers[i] <- names(thismap)[index[i]]
}
}
}
if(length(chrnotfound) > 0) {
chrnotfound <- sort(unique(chrnotfound))
if(length(chrnotfound) == 1)
warning("Chromosome ", paste("\"", chrnotfound, "\"", sep=""), " not found")
else
warning("Chromosomes ", paste("\"", chrnotfound, "\"", sep="", collapse=", "), " not found")
}
markers
}
### Find the nearest pseudomarker to a particular position
find.pseudomarker <-
function(cross, chr, pos, where=c("draws","prob"), addchr=TRUE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
# if chr has length 1, expand if necessary
if(length(chr) == 1)
chr <- rep(chr,length(pos))
# otherwise, chr and pos should have same length
else if(length(chr) != length(pos))
stop("chr and pos must be the same length.")
markers <- rep("",length(chr))
where <- match.arg(where)
if(where=="draws" && !("draws" %in% names(cross$geno[[1]])))
stop("You'll need to first run sim.geno")
if(where=="prob" && !("prob" %in% names(cross$geno[[1]])))
stop("You'll need to first run calc.genoprob")
for(i in 1:length(chr)) {
# find chromosome
o <- match(chr[i], names(cross$geno))
if(is.na(o)) markers[i] <- NA # chr not matched
else {
if(where=="draws") {
if("map" %in% names(attributes(cross$geno[[o]]$draws)))
thismap <- attr(cross$geno[[o]]$draws, "map")
else {
stp <- attr(cross$geno[[o]]$draws, "step")
oe <- attr(cross$geno[[o]]$draws, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[o]]$draws)))
stpw <- attr(cross$geno[[o]]$draws, "stepwidth")
else stpw <- "fixed"
thismap <- create.map(cross$geno[[o]]$map,stp,oe,stpw)
attr(cross$geno[[o]]$draws, "map") <- thismap
}
}
else { # prob
if("map" %in% names(attributes(cross$geno[[o]]$prob)))
thismap <- attr(cross$geno[[o]]$prob, "map")
else {
stp <- attr(cross$geno[[o]]$prob, "step")
oe <- attr(cross$geno[[o]]$prob, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[o]]$prob)))
stpw <- attr(cross$geno[[o]]$prob, "stepwidth")
else stpw <- "fixed"
thismap <- create.map(cross$geno[[o]]$map,stp,oe,stpw)
attr(cross$geno[[o]]$prob, "map") <- thismap
}
}
# sex-specific map; look at female positions
if(is.matrix(thismap)) thismap <- thismap[1,]
# find closest marker
d <- abs(thismap-pos[i])
o2 <- (1:length(d))[d==min(d)]
if(length(o2)==1) themarker <- names(thismap)[o2]
else themarker <- names(thismap)[sample(o2, 1)]
if(addchr && length(grep("^loc[0-9]+\\.*[0-9]*(\\.[0-9]+)*$", themarker)) > 0)
themarker <- paste("c", chr[i], ".", themarker, sep="")
markers[i] <- themarker
}
}
markers
}
# expand recombination fractions for RI lines
adjust.rf.ri <-
function(r, type=c("self","sib"), chrtype=c("A","X"), expand=TRUE)
{
# type of RI lines
type <- match.arg(type)
chrtype <- match.arg(chrtype)
if(type=="self") {
if(expand) return(r*2/(1+2*r))
else return(r/2/(1-r))
}
else {
if(chrtype=="A") { # autosome / sib mating
if(expand) return(r*4/(1+6*r))
else return(r/(4-6*r))
}
else { # X chromosome/ sib mating
if(expand) return(8/3*r/(1+4*r))
else return(3/8*r/(1-1.5*r))
}
}
}
######################################################################
# lodint: function to get lod support interval
######################################################################
lodint <-
function(results, chr, qtl.index, drop=1.5, lodcolumn=1,
expandtomarkers=FALSE)
{
if(!("scanone" %in% class(results))) {
if(!("qtl" %in% class(results)))
stop("Input must have class \"scanone\" or \"qtl\".")
else {
if(!("lodprofile" %in% names(attributes(results))))
stop("qtl object needs to be produced by refineqtl with keeplodprofile=TRUE.")
else { # qtl object
if(lodcolumn != 1) {
warning("lod column ignored if input is a qtl object.")
lodcolumn <- 1
}
results <- attr(results, "lodprofile")
if(missing(qtl.index)) {
if(length(results)==1)
results <- results[[1]]
else
stop("You must specify qtl.index.")
}
else {
if(length(qtl.index)>1) stop("qtl.index should have length 1")
if(qtl.index < 1 || qtl.index > length(results))
stop("qtl.index misspecified.")
results <- results[[qtl.index]]
}
chr <- results[1,1]
}
}
}
else {
if(lodcolumn < 1 || lodcolumn +2 > ncol(results))
stop("Argument lodcolumn misspecified.")
if(missing(chr)) {
if(length(unique(results[,1]))>1)
stop("Give a chromosome ID.")
}
else {
if(length(chr) > 1) stop("chr should have length 1")
if(is.na(match(chr, results[,1])))
stop("Chromosome misspecified.")
results <- results[results[,1]==chr,]
}
}
if(all(is.na(results[,lodcolumn+2]))) return(NULL)
maxlod <- max(results[,lodcolumn+2],na.rm=TRUE)
w <- which(!is.na(results[,lodcolumn+2]) & results[,lodcolumn+2] == maxlod)
o <- range(which(!is.na(results[,lodcolumn+2]) & results[,lodcolumn+2] > maxlod-drop))
if(length(o)==0) o <- c(1,nrow(results))
else {
if(o[1] > 1) o[1] <- o[1]-1
if(o[2] < nrow(results)) o[2] <- o[2]+1
}
if(expandtomarkers) {
markerpos <- (1:nrow(results))[-grep("^c.+\\.loc-*[0-9]+(\\.[0-9]+)*$", rownames(results))]
if(any(markerpos <= o[1]))
o[1] <- max(markerpos[markerpos <= o[1]])
if(any(markerpos >= o[2]))
o[2] <- min(markerpos[markerpos >= o[2]])
}
rn <- rownames(results)[c(o[1],w,o[2])]
# look for duplicate rows
if(length(w)>1 && rn[length(rn)]==rn[length(rn)-1]) w <- w[-length(w)]
else if(length(w)>1 && rn[2]==rn[1]) w <- w[-1]
rn <- rownames(results)[c(o[1],w,o[2])]
# look for more duplicate rows
if(any(table(rn)> 1)) {
tab <- table(rn)
temp <- which(tab>1)
for(j in temp) {
z <- which(rn==names(tab)[j])
for(k in 2:length(z))
rn[z[k:length(z)]] <- paste(rn[z[k:length(z)]], " ", sep="")
}
}
results <- results[c(o[1],w,o[2]),]
rownames(results) <- rn
class(results) <- c("scanone","data.frame")
results
}
######################################################################
# bayesint: function to get Bayesian probability interval
######################################################################
bayesint <-
function(results, chr, qtl.index, prob=0.95, lodcolumn=1, expandtomarkers=FALSE)
{
if(!("scanone" %in% class(results))) {
if(!("qtl" %in% class(results)))
stop("Input must have class \"scanone\" or \"qtl\".")
else {
if(!("lodprofile" %in% names(attributes(results))))
stop("qtl object needs to be produced by refineqtl with keeplodprofile=TRUE.")
else { # qtl object
if(lodcolumn != 1) {
warning("lod column ignored if input is a qtl object.")
lodcolumn <- 1
}
results <- attr(results, "lodprofile")
if(missing(qtl.index)) {
if(length(results)==1)
results <- results[[1]]
else
stop("You must specify qtl.index.")
}
else {
if(length(qtl.index)>1) stop("qtl.index should have length 1")
if(qtl.index < 1 || qtl.index > length(results))
stop("qtl.index misspecified.")
results <- results[[qtl.index]]
}
chr <- results[1,1]
}
}
}
else {
if(lodcolumn < 1 || lodcolumn +2 > ncol(results))
stop("Argument lodcolumn misspecified.")
if(missing(chr)) {
if(length(unique(results[,1]))>1)
stop("Give a chromosome ID.")
}
else {
if(length(chr) > 1) stop("chr should have length 1")
if(is.na(match(chr, results[,1])))
stop("Chromosome misspecified.")
results <- results[results[,1]==chr,]
}
}
if(all(is.na(results[,lodcolumn+2]))) return(NULL)
loc <- results[,2]
width <- diff(( c(loc[1],loc) + c(loc, loc[length(loc)]) )/ 2)
area <- 10^results[,lodcolumn+2]*width
area <- area/sum(area)
o <- order(results[,lodcolumn+2], decreasing=TRUE)
cs <- cumsum(area[o])
wh <- min((1:length(loc))[cs >= prob])
int <- range(o[1:wh])
if(expandtomarkers) {
markerpos <- (1:nrow(results))[-grep("^c.+\\.loc-*[0-9]+(\\.[0-9]+)*$", rownames(results))]
if(any(markerpos <= int[1]))
int[1] <- max(markerpos[markerpos <= int[1]])
if(any(markerpos >= int[2]))
int[2] <- min(markerpos[markerpos >= int[2]])
}
rn <- rownames(results)[c(int[1],o[1],int[2])]
# look for duplicate rows
if(any(table(rn)> 1)) {
rn[2] <- paste(rn[2], "")
if(rn[1] == rn[3]) rn[3] <- paste(rn[3], " ")
}
results <- results[c(int[1],o[1],int[2]),]
rownames(results) <- rn
class(results) <- c("scanone", "data.frame")
results
}
######################################################################
# makeSSmap: convert a genetic map, or the genetic maps in a cross
# object, to be sex-specific (i.e., 2-row matrices)
######################################################################
makeSSmap <-
function(cross)
{
if(!any(class(cross) == "map")) {
# input object is a genetic map
for(i in 1:length(cross)) {
if(!is.matrix(cross[[i]]))
cross[[i]] <- rbind(cross[[i]], cross[[i]])
}
}
else { # input object is assumed to be a "cross" object
n.chr <- nchr(cross)
for(i in 1:n.chr) {
if(!is.matrix(cross$geno[[i]]$map))
cross$geno[[i]]$map <-
rbind(cross$geno[[i]]$map, cross$geno[[i]]$map)
}
}
cross
}
######################################################################
# comparecrosses: verify that two cross objects have identical
# classes, chromosomes, markers, genotypes, maps,
# and phenotypes
######################################################################
comparecrosses <-
function(cross1, cross2, tol=1e-5)
{
if(missing(cross1) || missing(cross2))
stop("Two crosses must be input.")
# both are of class "cross"
if(!any(class(cross1) == "cross") ||
!any(class(cross2) == "cross"))
stop("Input should have class \"cross\".")
# classes are the same
if(any(class(cross1) != class(cross2)))
stop("crosses are not the same type.")
if(nchr(cross1) != nchr(cross2))
stop("crosses do not have the same number of chromosomes.")
if(any(names(cross1$geno) != names(cross2$geno)))
stop("Chromosome names do not match.")
if(any(nmar(cross1) != nmar(cross2)))
stop("Number of markers per chromosome do not match.")
mnames1 <- unlist(lapply(cross1$geno, function(a) colnames(a$data)))
mnames2 <- unlist(lapply(cross2$geno, function(a) colnames(a$data)))
if(any(mnames1 != mnames2)) {
# stop("Markers names do not match.")
for(i in 1:nchr(cross1))
if(any(colnames(cross1$geno[[i]]$data) != colnames(cross2$geno[[i]]$data)))
stop("Marker names on chr ", names(cross1$geno)[i], " don't match.")
}
chrtype1 <- sapply(cross1$geno, class)
chrtype2 <- sapply(cross2$geno, class)
if(any(chrtype1 != chrtype2))
stop("Chromosome types (autosomal vs X) do not match.")
for(i in 1:nchr(cross1)) {
if(any(abs(diff(cross1$geno[[i]]$map) - diff(cross2$geno[[i]]$map)) > tol))
stop("Genetic maps for chromosome ", names(cross1$geno)[i],
" do not match.")
if(abs(cross1$geno[[i]]$map[1] - cross2$geno[[i]]$map[1]) > tol)
warning("Initial marker positions for chromosome ", names(cross1$geno)[i],
" do not match.")
}
if(nind(cross1) != nind(cross2))
stop("Number of individuals do not match.")
for(i in 1:nchr(cross1)) {
g1 <- cross1$geno[[i]]$data
g2 <- cross2$geno[[i]]$data
if(any((is.na(g1) & !is.na(g2)) | (!is.na(g1) & is.na(g2)) |
(!is.na(g1) & !is.na(g2) & g1!=g2)))
stop("Genotype data for chromosome ", names(cross1$geno)[i],
" do not match.")
}
if(nphe(cross1) != nphe(cross2))
stop("Number of phenotypes do not match.")
if(any(names(cross1$pheno) != names(cross2$pheno)))
stop("Phenotype names do not match.")
for(i in 1:nphe(cross1)) {
phe1 <- cross1$pheno[,i]
phe2 <- cross2$pheno[,i]
if(is.numeric(phe1) & is.numeric(phe2)) {
phe1[phe1 == Inf] <- max(phe1[phe1 < Inf], na.rm=TRUE)+5
phe2[phe2 == Inf] <- max(phe2[phe2 < Inf], na.rm=TRUE)+5
phe1[phe1 == -Inf] <- min(phe1[phe1 > -Inf], na.rm=TRUE)-5
phe2[phe2 == -Inf] <- min(phe2[phe2 > -Inf], na.rm=TRUE)-5
if(any((is.na(phe1) & !is.na(phe2)) | (!is.na(phe1) & is.na(phe2)) |
(!is.na(phe1) & !is.na(phe2) & abs(phe1-phe2) > tol))) {
stop("Data for phenotype ", names(cross1$pheno)[i],
" do not match.")
}
}
else {
if(any((is.na(phe1) & !is.na(phe2)) | (!is.na(phe1) & is.na(phe2)) |
(!is.na(phe1) & !is.na(phe2) &
as.character(phe1) != as.character(phe2)))) {
stop("Data for phenotype ", names(cross1$pheno)[i], " do not match.")
}
}
}
cat("\tCrosses are identical.\n")
}
######################################################################
# move marker
# Move a marker to a new chromosome...placed at the end
######################################################################
movemarker <-
function(cross, marker, newchr, newpos)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
mnames <- unlist(lapply(cross$geno,function(a) colnames(a$data)))
chr <- rep(names(cross$geno),nmar(cross))
pos <- unlist(lapply(cross$geno,function(a) 1:ncol(a$data)))
oldindex <- match(marker, mnames)
# Marker found precisely once?
if(is.na(oldindex)) stop(marker, " not found.\n")
if(length(oldindex) > 1) stop(marker, " found multiple times.\n")
chr <- chr[oldindex]
pos <- pos[oldindex]
# pull out genotype data
g <- cross$geno[[chr]]$data[,pos]
chrtype <- class(cross$geno[[chr]])
mapmatrix <- is.matrix(cross$geno[[chr]]$map)
# delete marker
if(nmar(cross)[chr] == 1) { # only marker on that chromosome, so drop the chromosome
cross$geno <- cross$geno[-match(chr,names(cross$geno))]
delchr <- TRUE
}
else {
delchr <- FALSE
cross$geno[[chr]]$data <- cross$geno[[chr]]$data[,-pos,drop=FALSE]
if(is.matrix(cross$geno[[chr]]$map))
cross$geno[[chr]]$map <- cross$geno[[chr]]$map[,-pos,drop=FALSE]
else
cross$geno[[chr]]$map <- cross$geno[[chr]]$map[-pos]
}
if(is.numeric(newchr)) newchr <- as.character(newchr)
if(!(newchr %in% names(cross$geno))) { # create a new chromosome
n <- length(cross$geno)
cross$geno[[n+1]] <- list("data"=as.matrix(g),
"map"=as.numeric(0))
names(cross$geno)[n+1] <- newchr
class(cross$geno[[n+1]]) <- chrtype
colnames(cross$geno[[n+1]]$data) <- marker
if(mapmatrix) {
if(missing(newpos)) newpos <- 0
cross$geno[[n+1]]$map <- matrix(newpos, ncol=1, nrow=2)
colnames(cross$geno[[n+1]]$map) <- marker
}
else {
if(missing(newpos)) newpos <- 0
cross$geno[[n+1]]$map <- newpos
names(cross$geno[[n+1]]$map) <- marker
}
return(cross)
}
if(missing(newpos)) {
# add marker to end of new chromosome
n.mar <- nmar(cross)[newchr]
cross$geno[[newchr]]$data <- cbind(cross$geno[[newchr]]$data,g)
colnames(cross$geno[[newchr]]$data)[n.mar+1] <- marker
if(is.matrix(cross$geno[[newchr]]$map)) {
cross$geno[[newchr]]$map <- cbind(cross$geno[[newchr]]$map,
cross$geno[[newchr]]$map[,n.mar]+10)
colnames(cross$geno[[newchr]]$map)[n.mar+1] <- marker
}
else {
cross$geno[[newchr]]$map <- c(cross$geno[[newchr]]$map,
cross$geno[[newchr]]$map[n.mar]+10)
names(cross$geno[[newchr]]$map)[n.mar+1] <- marker
}
}
else {
# add marker to the specified position
dat <- cross$geno[[newchr]]$data
map <- cross$geno[[newchr]]$map
if(length(newpos) != 1)
stop("newpos should be a single number.")
if(is.matrix(map)) { # sex-specific maps
wh <- which(map[1,] < newpos)
if(length(wh) == 0) { # place in first spot
map <- cbind(c(newpos,map[2,1]-(map[1,1]-newpos)),map)
colnames(map)[1] <- marker
}
else {
wh <- max(wh)
if(wh == ncol(map)) { # place at end of chromosome
map <- cbind(map,c(newpos,map[2,ncol(map)]+(newpos-map[1,ncol(map)])))
colnames(map)[ncol(map)] <- marker
}
else {
left <- map[,1:wh,drop=FALSE]
right <- map[,-(1:wh),drop=FALSE]
marleft <- colnames(left)
marright <- colnames(right)
left <- left[,ncol(left)]
right <- right[,1]
newpos2 <- (newpos-left[1])/(right[1]-left[1])*(right[2]-left[2])+left[2]
map <- cbind(map[,1:wh], c(newpos,newpos2), map[,-(1:wh)])
colnames(map) <- c(marleft, marker, marright)
}
}
}
else {
wh <- which(map < newpos)
if(length(wh) == 0) { # place in first position
map <- c(newpos,map)
names(map)[1] <- marker
}
else {
wh <- max(wh)
if(wh == length(map)) { # place in last position
map <- c(map,newpos)
names(map)[length(map)] <- marker
}
else {
map <- c(map[1:wh],newpos,map[-(1:wh)])
names(map)[wh+1] <- marker
}
}
}
cross$geno[[newchr]]$map <- map
if(length(wh)==0) { # place marker in first position
dat <- cbind(g, dat)
colnames(dat)[1] <- marker
}
else if(wh == ncol(dat)) { # place marker in last position
dat <- cbind(dat, g)
colnames(dat)[ncol(dat)] <- marker
}
else { # place marker in the middle
dat <- cbind(dat[,1:wh],g,dat[,-(1:wh)])
colnames(dat)[wh+1] <- marker
}
cross$geno[[newchr]]$data <- dat
# make sure the marker names for the data and the genetic map match
if(is.matrix(cross$geno[[newchr]]$map))
colnames(cross$geno[[newchr]]$data) <- colnames(cross$geno[[newchr]]$map)
else
colnames(cross$geno[[newchr]]$data) <- names(cross$geno[[newchr]]$map)
}
# update genoprob, errorlod, argmax, draws, rf
if("rf" %in% names(cross)) {
# reorder the recombination fractions
# -- a bit of pain, 'cause we need LODs in upper triangle
# and rec fracs in lower triangle
newmar <- unlist(lapply(cross$geno,function(a) colnames(a$data)))
attrib <- attributes(cross$rf)
rf <- cross$rf
lods <- rf;lods[lower.tri(rf)] <- t(rf)[lower.tri(rf)]
rf[upper.tri(rf)] <- t(rf)[upper.tri(rf)]
lods <- lods[newmar,newmar]
rf <- rf[newmar,newmar]
rf[upper.tri(rf)] <- lods[upper.tri(rf)]
cross$rf <- rf
if("onlylod" %in% names(attrib)) # save the onlylod attribute if its there
attr(cross$rf, "onlylod") <- attrib$onlylod
}
if(!delchr) thechr <- c(chr,newchr)
else thechr <- newchr
for(i in thechr) {
tempg <- cross$geno[[i]]
tempx <- subset(cross, chr=i)
if("prob" %in% names(tempg))
atp <- attributes(tempg$prob)
if("draws" %in% names(tempg)) {
at <- attributes(cross$geno[[1]]$draws)
tempg$draws <- sim.geno(tempx,
n.draws=at$dim[3],
step=at$step,
off.end=at$off.end,
map.function=at$map.function,
error.prob=at$error.prob)$geno[[1]]$draws
}
if("argmax" %in% names(tempg)) {
at <- attributes(cross$geno[[1]]$argmax)
tempg$argmax <- argmax.geno(tempx,
step=at$step,
off.end=at$off.end,
map.function=at$map.function,
error.prob=at$error.prob)$geno[[1]]$argmax
}
if("errorlod" %in% names(tempg)) {
at <- attributes(cross$geno[[1]]$errorlod)
tempg$errorlod <- argmax.geno(tempx,
map.function=at$map.function,
error.prob=at$error.prob)$geno[[1]]$errorlod
}
if("prob" %in% names(tempg))
tempg$prob <- calc.genoprob(tempx,
step=atp$step,
off.end=atp$off.end,
map.function=atp$map.function,
error.prob=atp$error.prob)$geno[[1]]$prob
cross$geno[[i]] <- tempg
}
cross
}
######################################################################
#
# summary.map
#
# Give a short summary of a genetic map object.
#
######################################################################
summaryMap <- summary.map <-
function(object, ...)
{
map <- object
if(any(class(map) == "cross")) # a cross object
map <- pull.map(map)
if(!any(class(map) == "map"))
stop("Input should have class \"cross\" or \"map\".")
n.chr <- length(map)
chrnames <- names(map)
if(is.null(chrnames)) chrnames <- 1:length(map)
if(is.matrix(map[[1]])) { # sex-specific map
sexsp <- TRUE
n.mar <- sapply(map,ncol)
tot.mar <- sum(n.mar)
fmap <- lapply(map,function(a) a[1,])
mmap <- lapply(map,function(a) a[2,])
len.f <- sapply(fmap,function(a) diff(range(a)))
len.m <- sapply(mmap,function(a) diff(range(a)))
avesp.f <- sapply(fmap,function(a) {if(length(a)<2) return(NA); mean(diff(a))})
avesp.m <- sapply(mmap,function(a) {if(length(a)<2) return(NA); mean(diff(a))})
maxsp.f <- sapply(fmap,function(a) {if(length(a)<2) return(NA); max(diff(a))})
maxsp.m <- sapply(mmap,function(a) {if(length(a)<2) return(NA); max(diff(a))})
totlen.f <- sum(len.f)
totlen.m <- sum(len.m)
tot.avesp.f <- mean(unlist(lapply(fmap,diff)), na.rm=TRUE)
tot.avesp.m <- mean(unlist(lapply(mmap,diff)), na.rm=TRUE)
tot.maxsp.f <- max(maxsp.f,na.rm=TRUE)
tot.maxsp.m <- max(maxsp.m,na.rm=TRUE)
output <- rbind(cbind(n.mar,len.f,len.m,avesp.f,avesp.m, maxsp.f, maxsp.m),
c(tot.mar,totlen.f,totlen.m,tot.avesp.f,tot.avesp.m,
tot.maxsp.f, tot.maxsp.m))
dimnames(output) <- list(c(chrnames,"overall"),
c("n.mar","length.female","length.male",
"ave.spacing.female","ave.spacing.male",
"max.spacing.female", "max.spacing.male"))
}
else {
sexsp=FALSE
n.mar <- sapply(map,length)
len <- sapply(map,function(a) diff(range(a)))
tot.mar <- sum(n.mar)
len <- sapply(map,function(a) diff(range(a)))
avesp <- sapply(map,function(a) {if(length(a)<2) return(NA); mean(diff(a))})
maxsp <- sapply(map,function(a) {if(length(a)<2) return(NA); max(diff(a))})
totlen <- sum(len)
tot.avesp <- mean(unlist(lapply(map,diff)), na.rm=TRUE)
tot.maxsp <- max(maxsp, na.rm=TRUE)
output <- rbind(cbind(n.mar,len,avesp, maxsp),
c(tot.mar,totlen,tot.avesp, tot.maxsp))
dimnames(output) <- list(c(chrnames,"overall"),
c("n.mar","length","ave.spacing", "max.spacing"))
}
output <- as.data.frame(output)
attr(output, "sexsp") <- sexsp
class(output) <- c("summary.map", "data.frame")
output
}
######################################################################
#
# print.summary.map
#
# Print out the result of summary.map()
#
######################################################################
print.summary.map <-
function(x, ...)
{
sexsp <- attr(x, "sexsp")
if(sexsp) cat("Sex-specific map\n\n")
x <- apply(x,2,round,1)
print(x)
}
######################################################################
# convert functions
######################################################################
convert <-
function(object, ...)
UseMethod("convert")
######################################################################
#
# convert.scanone
#
# Convert scanone output from the format for R/qtl ver 0.97 to
# that for R/qtl ver 0.98
# (previously, inter-maker locations named loc*.c*; now c*.loc*)
#
######################################################################
convert.scanone <-
function(object, ...)
{
if(!any(class(object) == "scanone"))
stop("Input should have class \"scanone\".")
rn <- rownames(object)
o <- grep("^loc-*[0-9]+(\\.[0-9]+)*\\.c[0-9A-Za-z]+$", rn)
if(length(o) > 0) {
temp <- rn[o]
temp <- strsplit(temp,"\\.")
temp <- sapply(temp, function(a)
paste(a[c(length(a),1:(length(a)-1))],collapse="."))
rownames(object)[o] <- temp
}
object
}
######################################################################
# convert.scantwo
#
# convert scantwo output from the format used in R/qtl version 1.03
# and earlier to that used in R/qtl version 1.04 and later.
#
# 1.03 and earlier: contained joint and interactive LOD scores
# 1.04 and later: contains joint LOD scores and LOD scores from the
# additive QTL model
######################################################################
convert.scantwo <-
function(object, ...)
{
if(!any(class(object) == "scantwo"))
stop("Input should have class \"scantwo\".")
lod <- object$lod
if(length(dim(lod)) == 2) {
u <- upper.tri(lod)
lod[u] <- t(lod)[u] - lod[u]
}
else { # multiple phenotypes
u <- upper.tri(lod[,,1])
for(i in 1:dim(lod)[3])
lod[,,i][u] <- t(lod[,,i])[u] - lod[,,i][u]
}
object$lod <- lod
object
}
######################################################################
# convert.map
#
# convert a genetic map from one map function to another
######################################################################
convert.map <-
function(object, old.map.function=c("haldane", "kosambi", "c-f", "morgan"),
new.map.function=c("haldane", "kosambi", "c-f", "morgan"), ...)
{
old.map.function <- match.arg(old.map.function)
new.map.function <- match.arg(new.map.function)
if(!("map" %in% class(object)))
stop("Input should have class \"map\".")
if(old.map.function==new.map.function) {
warning("old and new map functions are the same; no change.")
return(object)
}
mf <- switch(old.map.function,
"haldane"=mf.h,
"kosambi"=mf.k,
"c-f"=mf.cf,
"morgan"=mf.m)
imf <- switch(new.map.function,
"haldane"=imf.h,
"kosambi"=imf.k,
"c-f"=imf.cf,
"morgan"=imf.m)
if(is.matrix(object[[1]])) { # sex-specific map
for(i in seq(along=object)) {
theclass <- class(object[[i]])
thenames <- colnames(object[[i]])
for(j in 1:2)
object[[i]][j,] <- cumsum(c(object[[i]][j,1], imf(mf(diff(object[[i]][j,])))))
class(object[[i]]) <- theclass
colnames(object[[i]]) <- thenames
}
}
else {
for(i in seq(along=object)) {
theclass <- class(object[[i]])
thenames <- names(object[[i]])
object[[i]] <- cumsum(c(object[[i]][1], imf(mf(diff(object[[i]])))))
class(object[[i]]) <- theclass
names(object[[i]]) <- thenames
}
}
object
}
######################################################################
# find.pheno
#
# utility to get pheno number given pheno name
######################################################################
find.pheno <-
function( cross, pheno )
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
seq( ncol( cross$pheno ))[match(pheno,names(cross$pheno))]
}
######################################################################
# find.flanking
#
# utility to get flanking and/or closest marker to chr and pos
######################################################################
find.flanking <-
function( cross, chr, pos)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
map = pull.map(cross)
if(is.matrix(map[[1]]) && nrow(map[[1]]) > 1)
stop("This function works only for crosses with sex-averaged maps.")
if(length(chr) == 1 && length(pos) > 1) {
chr <- rep(chr,length(pos))
}
wh <- match(chr, names(cross$geno))
if(any(is.na(wh))) {
stop("Chr ", paste(chr[is.na(wh)], collapse=", "), " not found.")
wh <- wh[!is.na(wh)]
}
chr <- names(cross$geno)[wh]
marker = NULL
for (i in seq(length(chr))) {
tmp = map[[chr[i]]]-pos[i]
m = names(map[[chr[i]]])
left = sum(tmp < 0)
at = sum(tmp == 0)
right = sum(tmp > 0)
f <- if (at > 0)
left+at[c(1,length(at))]
else {
if (right > 0)
c(left,left+at+1)
else
c(left,left+at)
}
marker = rbind(marker,m[f[c(1:2,order(abs(tmp[f]))[1])]])
}
dimnames(marker) <- list(paste(chr,":",pos,sep=""),
c("left","right","close"))
as.data.frame(marker, stringsAsFactors=TRUE)
}
######################################################################
# strip.partials
#
# Removes partially informative genotypes in an intercross.
#
# Input: a cross object; if verbose=TRUE, a statement regarding the
# number of genotypes removed is printed.
######################################################################
strip.partials <-
function(cross, verbose=TRUE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
type <- class(cross)[1]
if(type != "f2")
stop("This is for intercrosses only")
n.removed <- 0
for(i in 1:nchr(cross)) {
g <- cross$geno[[i]]$data
wh <- !is.na(g) & g>3
if(any(wh)) {
g[wh] <- NA
cross$geno[[i]]$data <- g
n.removed <- n.removed + sum(wh)
}
}
if(verbose) {
if(n.removed == 0) cat(" --Didn't remove any genotypes.\n")
else cat("Removed", n.removed, "genotypes.\n")
}
cross
}
######################################################################
# comparegeno
######################################################################
comparegeno <-
function(cross, what=c("proportion","number", "both"))
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
what <- match.arg(what)
g <- pull.geno(cross)
g[is.na(g)] <- 0
n.ind <- nrow(g)
n.mar <- ncol(g)
z <- .C("R_comparegeno",
as.integer(g),
as.integer(n.ind),
as.integer(n.mar),
n.match=as.integer(rep(0,n.ind^2)),
n.missing=as.integer(rep(0,n.ind^2)),
PACKAGE="qtl")
if(what=="number") {
z <- matrix(z$n.match,n.ind,n.ind)
}
else {
if(what=="proportion") {
z <- matrix(z$n.match/(n.mar-z$n.missing),n.ind,n.ind)
diag(z) <- NA
}
else {
prop <- matrix(z$n.match/(n.mar-z$n.missing),n.ind,n.ind)
z <- matrix(z$n.match,n.ind,n.ind)
z[lower.tri(z)] <- prop[lower.tri(z)]
}
}
z
}
######################################################################
# print the installed version of R/qtl
######################################################################
qtlversion <-
function()
{
version <- unlist(packageVersion("qtl"))
# make it like #.#-#
paste(c(version,".","-")[c(1,4,2,5,3)], collapse="")
}
######################################################################
#
# locateXO
#
# Locate crossovers on a single chromosome in each individual
# Look at just the first chromosome
#
######################################################################
locateXO <-
function(cross, chr, full.info=FALSE)
{
if(!missing(chr)) {
cross <- subset(cross, chr=chr)
if(nchr(cross) != 1)
warning("locateXO works on just one chr; considering chr ", names(cross$geno)[1])
}
# individual IDs
id <- getid(cross)
if(is.null(id)) id <- as.character(1:nind(cross))
if(nmar(cross)[1] == 1) { # just one marker; don't need to do anything
warning("Just one marker.")
res <- vector("list", nind(cross))
names(res) <- id
for(i in seq(along=res)) res[[i]] <- numeric(0)
return(res)
}
geno <- cross$geno[[1]]$data
geno[is.na(geno)] <- 0
type <- class(cross)[1]
if(type != "bc" && type != "f2" && type != "riself" && type != "risib" && type!="dh")
stop("locateXO only working for backcross, intercross or RI strains.")
map <- cross$geno[[1]]$map
if(is.matrix(map)) map <- map[1,]
# map <- map - map[1] # shift first marker to 0
# bc or intercross? thetype==0 for BC and ==1 for intercross
if(type=="f2") {
if(class(geno) == "X") thetype <- 0
else thetype <- 1
}
else thetype <- 0
n.ind <- nrow(geno)
n.mar <- ncol(geno)
z <- .C("R_locate_xo",
as.integer(n.ind),
as.integer(n.mar),
as.integer(thetype),
as.integer(geno),
as.double(map),
location=as.double(rep(0,n.ind*2*(n.mar-1))),
nseen=as.integer(rep(0,n.ind)),
ileft=as.integer(rep(0,n.ind*2*(n.mar-1))),
iright=as.integer(rep(0,n.ind*2*(n.mar-1))),
left=as.double(rep(0,n.ind*2*(n.mar-1))),
right=as.double(rep(0,n.ind*2*(n.mar-1))),
gleft=as.integer(rep(0, n.ind*2*(n.mar-1))),
gright=as.integer(rep(0, n.ind*2*(n.mar-1))),
ntype=as.integer(rep(0,n.ind*2*(n.mar-1))),
as.integer(full.info),
PACKAGE="qtl")
location <- t(matrix(z$location, nrow=n.ind))
nseen <- z$nseen
if(full.info) {
ileft <- t(matrix(z$ileft, nrow=n.ind))
iright <- t(matrix(z$iright, nrow=n.ind))
left <- t(matrix(z$left, nrow=n.ind))
right <- t(matrix(z$right, nrow=n.ind))
gleft <- t(matrix(z$gleft, nrow=n.ind))
gright <- t(matrix(z$gright, nrow=n.ind))
ntype <- t(matrix(z$ntype, nrow=n.ind))
}
if(!full.info)
res <- lapply(as.data.frame(rbind(nseen, location), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
else {
location <- lapply(as.data.frame(rbind(nseen, location), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
ileft <- lapply(as.data.frame(rbind(nseen, ileft), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
iright <- lapply(as.data.frame(rbind(nseen, iright), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
left <- lapply(as.data.frame(rbind(nseen, left), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
right <- lapply(as.data.frame(rbind(nseen, right), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
gleft <- lapply(as.data.frame(rbind(nseen, gleft), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
gright <- lapply(as.data.frame(rbind(nseen, gright), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
ntype <- lapply(as.data.frame(rbind(nseen, ntype), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
res <- location
for(i in seq(along=res)) {
if(length(res[[i]])>0) {
ntype[[i]][length(ntype[[i]])] <- NA
res[[i]] <- cbind(location=location[[i]],
left=left[[i]],
right=right[[i]],
ileft=ileft[[i]],
iright=iright[[i]],
gleft=gleft[[i]],
gright=gright[[i]],
nTypedBetween=ntype[[i]])
}
}
}
names(res) <- id
res
}
# jittermap: make sure no two markers are at precisely the same position
jittermap <-
function(object, amount=1e-6) # x is either a cross object or a map
{
if(any(class(object) == "cross")) {
themap <- pull.map(object)
return.cross <- TRUE
}
else {
if(!any(class(object) == "map"))
stop("Input must be a cross or a map")
return.cross <- FALSE
themap <- object
}
for(i in 1:length(themap)) {
if(is.matrix(themap[[i]])) { # sex-specific maps
n <- ncol(themap[[i]])
if(n > 1) {
for(j in 1:2)
themap[[i]][j,] <- themap[[i]][j,] + c(0,cumsum(rep(amount, n-1)))
}
}
else {
n <- length(themap[[i]])
if(n > 1)
themap[[i]] <- themap[[i]] + c(0,cumsum(rep(amount, n-1)))
}
}
if(return.cross)
return(clean(replace.map(object, themap)))
themap
}
print.map <-
function(x, ...)
{
for(i in seq(along=x))
attr(x[[i]], "loglik") <- NULL
if(length(x) == 1)
print(unclass(x[[1]]))
else
print(unclass(lapply(x, unclass)))
}
# map function for Stahl model
mf.stahl <-
function(d, m=0, p=0)
{
if(any(d < 0))
stop("d must be >= 0\n")
if(m < 0)
stop("Must have m >= 0\n")
if(p < 0 || p > 1)
stop("Must have 0 <= p <= 1\n")
# handle missing values
if(any(!is.finite(d))) {
which.finite <- is.finite(d)
dsub <- d[which.finite]
dropped.some <- TRUE
} else {
dsub <- d
dropped.some <- FALSE
}
result <- .C("R_mf_stahl",
as.integer(length(dsub)),
as.double(dsub/100), # convert to Morgans
as.integer(m),
as.double(p),
result=as.double(rep(0,length(dsub))),
PACKAGE="qtl")$result
if(dropped.some) {
d[!which.finite] <- NA
d[which.finite] <- result
result <- d
}
result
}
# inverse map function for Stahl model
imf.stahl <-
function(r, m=0, p=0, tol=1e-12, maxit=1000)
{
if(any(r < 0 | r > 0.5))
stop("r must be >= 0 or <= 0.5\n")
if(m < 0)
stop("Must have m >= 0\n")
if(p < 0 || p > 1)
stop("Must have 0 <= p <= 1\n")
# handle missing values
if(any(!is.finite(r))) {
which.finite <- is.finite(r)
rsub <- r[which.finite]
dropped.some <- TRUE
} else {
rsub <- r
dropped.some <- FALSE
}
result <- .C("R_imf_stahl",
as.integer(length(rsub)),
as.double(rsub),
as.integer(m),
as.double(p),
result=as.double(rep(0,length(rsub))),
as.double(tol),
as.integer(maxit),
PACKAGE="qtl")$result*100 # convert to cM
if(dropped.some) {
r[!which.finite] <- NA
r[which.finite] <- result
result <- r
}
result
}
######################################################################
# getid: internal function to pull out the "ID" column from the
# phenotype data, if there is one
######################################################################
getid <-
function(cross)
{
phe <- cross$pheno
nam <- names(phe)
if("id" %in% nam) {
id <- phe$id
phenam <- "id"
}
else if("ID" %in% nam) {
id <- phe$ID
phenam <- "ID"
}
else if("Id" %in% nam) {
id <- phe$Id
phenam <- "Id"
}
else if("iD" %in% nam) {
id <- phe$iD
phenam <- "iD"
}
else {
id <- NULL
phenam <- NULL
}
if(is.factor(id))
id <- as.character(id)
attr(id, "phenam") <- phenam
id
}
######################################################################
# find the chromosome and position of a vector of markers
######################################################################
find.markerpos <-
function(cross, marker)
{
if(length(marker) != length(unique(marker))) {
warning("Dropping duplicate marker names.")
marker <- unique(marker)
}
output <- data.frame(chr=rep("", length(marker)),
pos=rep(NA, length(marker)), stringsAsFactors=TRUE)
output$chr <- as.character(output$chr)
rownames(output) <- marker
map <- pull.map(cross)
n.mar <- nmar(cross)
chr <- rep(names(map), n.mar)
if(!is.matrix(map[[1]])) {
pos <- unlist(map)
onemap <- TRUE
}
else {
pos <- unlist(lapply(map, function(a) a[1,]))
pos2 <- unlist(lapply(map, function(a) a[2,]))
onemap <- FALSE
output <- cbind(output, pos2=rep(NA, length(marker)))
colnames(output)[2:3] <- c("pos.female","pos.male")
}
mnam <- markernames(cross)
for(i in seq(along=marker)) {
wh <- match(marker[i], mnam)
if(is.na(wh)) {
output[i,1] <- NA
output[i,2] <- NA
}
else {
if(length(wh) > 1) {
warning("Marker ", marker[i], " found multiple times.")
wh <- sample(wh, 1)
}
output[i,1] <- chr[wh]
output[i,2] <- pos[wh]
if(!onemap) output[i,3] <- pos2[wh]
}
}
output
}
######################################################################
# find the chromosome and position of a vector of pseudomarkers
######################################################################
find.pseudomarkerpos <-
function(cross, marker, where=c("draws","prob"))
{
if(length(marker) != length(unique(marker))) {
warning("Dropping duplicate pseudomarker names.")
marker <- unique(marker)
}
output <- data.frame(chr=rep("", length(marker)),
pos=rep(NA, length(marker)), stringsAsFactors=TRUE)
output$chr <- as.character(output$chr)
rownames(output) <- marker
where <- match.arg(where)
if(where=="draws" && !("draws" %in% names(cross$geno[[1]])))
stop("You'll need to first run sim.geno")
if(where=="prob" && !("prob" %in% names(cross$geno[[1]])))
stop("You'll need to first run calc.genoprob")
themap <- vector("list", nchr(cross))
names(themap) <- names(cross$geno)
for(i in 1:nchr(cross)) {
if(where=="draws")
themap[[i]] <- attr(cross$geno[[i]]$draws, "map")
else
themap[[i]] <- attr(cross$geno[[i]]$prob, "map")
}
if(!is.matrix(themap[[1]])) {
pmar <- unlist(lapply(themap, names))
pos <- unlist(themap)
onemap <- TRUE
chr <- rep(names(themap), sapply(themap, length))
}
else {
pos <- unlist(lapply(themap, function(a) a[1,]))
pos2 <- unlist(lapply(themap, function(a) a[2,]))
onemap <- FALSE
pmar <- unlist(lapply(themap, colnames))
output <- cbind(output, pos2=rep(NA, length(marker)))
colnames(output)[2:3] <- c("pos.female","pos.male")
chr <- rep(names(themap), sapply(themap, ncol))
}
whnotmarker <- grep("^loc-*[0-9]*", pmar)
pmar[whnotmarker] <- paste("c", chr[whnotmarker], ".", pmar[whnotmarker], sep="")
for(i in seq(along=marker)) {
wh <- match(marker[i], pmar)
if(is.na(wh)) {
output[i,1] <- NA
output[i,2] <- NA
}
else {
if(length(wh) > 1) {
warning("Pseudomarker ", marker[i], " found multiple times.")
wh <- sample(wh, 1)
}
output[i,1] <- chr[wh]
output[i,2] <- pos[wh]
if(!onemap) output[i,3] <- pos2[wh]
}
}
output
}
######################################################################
# utility function for determining whether pheno.col (as argument
# to scanone etc) can be interpreted as a vector of phenotypes,
# versus a vector of phenotype columns
######################################################################
LikePheVector <-
function(pheno, n.ind, n.phe)
{
if(is.numeric(pheno) && length(pheno)==n.ind &&
any(pheno < 1 | pheno > n.phe | pheno!=round(pheno)))
return(TRUE)
FALSE
}
######################################################################
# for matching chromosome names
######################################################################
matchchr <-
function(selection, thechr)
{
if(is.factor(thechr)) thechr <- as.character(thechr)
if(length(thechr) > length(unique(thechr)))
stop("Duplicate chromosome names.")
if(is.logical(selection)) {
if(length(selection) != length(thechr))
stop("Logical vector to select chromosomes is the wrong length")
return(thechr[selection])
}
if(is.numeric(selection)) selection <- as.character(selection)
if(length(selection) > length(unique(selection))) {
warning("Dropping duplicate chromosomes")
selection <- unique(selection)
}
g <- grep("^-", selection)
if(length(g) > 0 && length(g) < length(selection))
stop("In selecting chromosomes, all must start with '-' or none should.")
if(length(g) > 0) {
selectomit <- TRUE
selection <- substr(selection, 2, nchar(selection))
}
else selectomit <- FALSE
wh <- match(selection, thechr)
if(any(is.na(wh))) {
warning("Chr ", paste(selection[is.na(wh)], collapse=", "), " not found")
wh <- wh[!is.na(wh)]
if(length(wh) == 0) return(thechr)
}
if(selectomit) return(thechr[-wh])
thechr[sort(wh)]
}
######################################################################
# check that chromosomes match appropriately
# TRUE = chr okay
# FALSE = problem
######################################################################
testchr <-
function(selection, thechr)
{
if(is.factor(thechr)) thechr <- as.character(thechr)
if(length(thechr) > length(unique(thechr))) {
# warning("Duplicate chromosome names.")
return(FALSE)
}
if(is.logical(selection)) {
if(length(selection) != length(thechr)) {
# warning("Logical vector to select chromosomes is the wrong length")
return(FALSE)
}
return(TRUE)
}
if(is.numeric(selection)) selection <- as.character(selection)
if(length(selection) > length(unique(selection))) {
# warning("Dropping duplicate chromosomes")
selection <- unique(selection)
}
g <- grep("^-", selection)
if(length(g) > 0 && length(g) < length(selection)) {
# stop("In selecting chromosomes, all must start with '-' or none should.")
return(FALSE)
}
if(length(g) > 0) {
selectomit <- TRUE
selection <- substr(selection, 2, nchar(selection))
}
else selectomit <- FALSE
wh <- match(selection, thechr)
if(any(is.na(wh))) return(FALSE)
TRUE
}
######################################################################
# convert2sa
#
# convert a sex-specific maps to a sex-averaged one.
# We pull out just the female map, and give a warning if the male and
# female maps are too different
######################################################################
convert2sa <-
function(map, tol=1e-4)
{
if(!("map" %in% class(map)))
stop("Input should have class \"map\".")
if(!is.matrix(map[[1]]))
stop("Input map doesn't seem to be a sex-specific map.")
theclass <- sapply(map, class)
fem <- lapply(map, function(a) a[1,])
dif <- sapply(map, function(a) { if(ncol(a)==1) return(diff(a))
a <- apply(a, 1, diff);
if(is.matrix(a)) return(max(abs(apply(a, 1, diff))))
abs(diff(a)) })
if(max(dif) > tol)
warning("Female and male inter-marker distances differ by as much as ", max(dif), ".")
for(i in seq(along=theclass)) class(fem[[i]]) <- theclass[i]
class(fem) <- "map"
fem
}
# round as character string, ensuring ending 0's are kept.
charround <-
function(x, digits=1)
{
if(digits < 1)
stop("This is intended for the case digits >= 1.")
y <- as.character(round(x, digits))
z <- strsplit(y, "\\.")
sapply(z, function(a, digits)
{
if(length(a) == 1)
b <- paste(a[1], ".", paste(rep("0", digits),collapse=""), sep="")
else {
if(nchar(a[2]) == digits)
b <- paste(a, collapse=".")
else
b <- paste(a[1], ".", a[2],
paste(rep("0", digits - nchar(a[2])), collapse=""),
sep="")
}
}, digits)
}
######################################################################
# scantwoperm2scanoneperm
#
# pull out the scanone permutations from scantwo permutation results,
# so that one may use the scantwo perms in calls to summary.scanone
######################################################################
scantwoperm2scanoneperm <-
function(scantwoperms)
{
if(!("scantwoperm" %in% class(scantwoperms)))
stop("Input must have class \"scantwoperm\".")
if(!("one" %in% names(scantwoperms)))
stop("Input doesn't contain scanone permutation results.")
scanoneperms <- scantwoperms$one
class(scanoneperms) <- c("scanoneperm")
scanoneperms
}
######################################################################
# subset.map
######################################################################
subset.map <-
function(x, ...)
{
cl <- class(x)
x <- x[...]
class(x) <- cl
x
}
`[.map` <-
function(x, ...)
{
x <- unclass(x)[...]
class(x) <- "map"
x
}
######################################################################
# subset.cross with [ ]
######################################################################
`[.cross` <-
function(x, chr, ind)
subset(x, chr, ind)
######################################################################
# findDupMarkers
#
# find markers whose genotype data is identical to another marker
# (which thus might be dropped, as extraneous)
#
# chr (Optional) set of chromosomes to consider
#
# exact.only If TRUE, look for markers with the same genotypes
# and the same pattern of missing data
# If FALSE, also identify markers whose observed
# genotypes match another marker, with no
# observed genotypes for which the other is
# missing
#
# adjacent.only If TRUE, only consider adjacent markers
######################################################################
findDupMarkers <-
function(cross, chr, exact.only=TRUE, adjacent.only=FALSE)
{
if(!missing(chr)) cross <- subset(cross, chr=chr)
g <- pull.geno(cross)
markers <- colnames(g)
markerloc <- lapply(nmar(cross), function(a) 1:a)
if(length(markerloc) > 1) {
for(j in 2:length(markerloc))
markerloc[[j]] <- markerloc[[j]] + max(markerloc[[j-1]]) + 10
}
markerloc <- unlist(markerloc)
if(exact.only) {
g[is.na(g)] <- 0
# genotype data patterns
pat <- apply(g, 2, paste, collapse="")
# table of unique values
tab <- table(pat)
# no duplicates; return
if(all(tab == 1)) return(NULL)
wh <- which(tab > 1)
theloc <- themar <- vector("list", length(wh))
for(i in seq(along=wh)) {
themar[[i]] <- names(pat)[pat==names(tab)[wh[i]]]
theloc[[i]] <- markerloc[pat==names(tab)[wh[i]]]
}
if(adjacent.only) {
extraloc <- list()
extramar <- list()
for(i in seq(along=theloc)) {
d <- diff(theloc[[i]]) # look for adjacency
if(any(d>1)) { # split into adjacent groups
temp <- which(d>1)
first <- c(1, temp+1)
last <- c(temp, length(theloc[[i]]))
for(j in 2:length(first)) {
extraloc[[length(extraloc)+1]] <- theloc[[i]][first[j]:last[j]]
extramar[[length(extramar)+1]] <- themar[[i]][first[j]:last[j]]
}
themar[[i]] <- themar[[i]][first[1]:last[1]]
theloc[[i]] <- theloc[[i]][first[1]:last[1]]
}
}
themar <- c(themar, extramar)
theloc <- c(theloc, extraloc)
nm <- sapply(themar, length)
if(all(nm==1)) return(NULL) # nothing left
themar <- themar[nm>1]
theloc <- theloc[nm>1]
}
# order by first locus
o <- order(sapply(theloc, min))
themar <- themar[o]
randompics <- sapply(themar, function(a) sample(length(a), 1))
for(i in seq(along=themar)) {
names(themar)[i] <- themar[[i]][randompics[i]]
themar[[i]] <- themar[[i]][-randompics[i]]
}
}
else {
themar <- NULL
ntyp <- ntyped(cross, "mar")
o <- order(ntyp, decreasing=TRUE)
g[is.na(g)] <- 0
z <- .C("R_findDupMarkers_notexact",
as.integer(nrow(g)),
as.integer(ncol(g)),
as.integer(g),
as.integer(o),
as.integer(markerloc),
as.integer(adjacent.only),
result=as.integer(rep(0,length(o))),
PACKAGE="qtl")
if(all(z$result==0)) return(NULL)
u <- unique(z$result[z$result != 0])
themar <- vector("list", length(u))
names(themar) <- colnames(g)[u]
for(i in seq(along=themar))
themar[[i]] <- colnames(g)[z$result==u[i]]
}
themar
}
######################################################################
# convert2riself
######################################################################
convert2riself <-
function(cross)
{
if(class(cross)[2] != "cross")
stop("input must be a cross object.")
curtype <- class(cross)[1]
chrtype <- sapply(cross$geno, class)
whX <- NULL
if(any(chrtype != "A")) { # there's an X chromosome
whX <- names(cross$geno)[chrtype != "A"]
if(length(whX) > 1)
warning("Converting chromosomes ", paste(whX, collapse=" "), " to autosomal.")
else
warning("Converting chromosome ", whX, " to autosomal.")
for(i in whX) class(cross$geno[[i]]) <- "A"
}
gtab <- table(pull.geno(cross))
usethree <- FALSE
if(!is.na(gtab["3"])) {
if(!is.na(gtab["2"]) && gtab["3"] < gtab["2"])
usethree <- FALSE
else usethree <- TRUE
}
g2omit <- g3omit <- g4omit <- 0
for(i in 1:nchr(cross)) {
dat <- cross$geno[[i]]$data
g1 <- sum(!is.na(dat) & dat==1)
g2 <- sum(!is.na(dat) & dat==2)
g3 <- sum(!is.na(dat) & dat==3)
g4 <- sum(!is.na(dat) & dat>3)
if(usethree && chrtype[i] == "A") {
dat[!is.na(dat) & dat!=1 & dat!=3] <- NA
dat[!is.na(dat) & dat==3] <- 2
g2omit <- g2omit + g2
g4omit <- g4omit + g4
}
else {
dat[!is.na(dat) & dat!=1 & dat!=2] <- NA
g3omit <- g3omit + g3
g4omit <- g4omit + g4
}
cross$geno[[i]]$data <- dat
}
if(g2omit > 0)
warning("Omitting ", g2omit, " genotypes with code==2.")
if(g3omit > 0)
warning("Omitting ", g3omit, " genotypes with code==3.")
if(g4omit > 0)
warning("Omitting ", g4omit, " genotypes with code>3.")
class(cross)[1] <- "riself"
cross
}
######################################################################
# convert2risib
######################################################################
convert2risib <-
function(cross)
{
if(class(cross)[2] != "cross")
stop("input must be a cross object.")
curtype <- class(cross)[1]
chrtype <- sapply(cross$geno, class)
gtab <- table(pull.geno(cross))
usethree <- FALSE
if(!is.na(gtab["3"])) {
if(!is.na(gtab["2"]) && gtab["3"] < gtab["2"])
usethree <- FALSE
else usethree <- TRUE
}
g2omit <- g3omit <- g4omit <- 0
for(i in 1:nchr(cross)) {
dat <- cross$geno[[i]]$data
g1 <- sum(!is.na(dat) & dat==1)
g2 <- sum(!is.na(dat) & dat==2)
g3 <- sum(!is.na(dat) & dat==3)
g4 <- sum(!is.na(dat) & dat>3)
if(usethree) {
if(chrtype[i] == "A") {
dat[!is.na(dat) & dat!=1 & dat!=3] <- NA
dat[!is.na(dat) & dat==3] <- 2
g2omit <- g2omit + g2
g4omit <- g4omit + g4
}
else { # X chromosome
if(g2 >= g3) {
dat[!is.na(dat) & dat!=1 & dat!=2] <- NA
g3omit <- g3omit + g3
g4omit <- g4omit + g4
}
else {
dat[!is.na(dat) & dat!=1 & dat!=3] <- NA
dat[!is.na(dat) & dat==3] <- 2
g2omit <- g2omit + g2
g4omit <- g4omit + g4
}
}
}
else {
dat[!is.na(dat) & dat!=1 & dat!=2] <- NA
g3omit <- g3omit + g3
g4omit <- g4omit + g4
}
cross$geno[[i]]$data <- dat
}
if(g2omit > 0)
warning("Omitting ", g2omit, " genotypes with code==2.")
if(g3omit > 0)
warning("Omitting ", g3omit, " genotypes with code==3.")
if(g4omit > 0)
warning("Omitting ", g4omit, " genotypes with code>3.")
class(cross)[1] <- "risib"
cross
}
rescalemap <-
function(object, scale=1e-6)
{
if("cross" %in% class(object)) {
for(i in 1:nchr(object)) {
object$geno[[i]]$map <-
object$geno[[i]]$map * scale
}
if(abs(scale - 1) > 1e-6)
object <- clean(object) # strip off intermediate calculations
} else if("map" %in% class(object)) {
for(i in seq(along=object)) {
object[[i]] <- object[[i]] * scale
}
} else
stop("rescalemap works only for objects of class \"cross\" or \"map\".")
object
}
shiftmap <-
function(object, offset=0)
{
if("cross" %in% class(object)) {
if(length(offset) == 1) offset <- rep(offset, nchr(object))
else if(length(offset) != nchr(object))
stop("offset must have length 1 or n.chr (", nchr(object), ")")
for(i in 1:nchr(object)) {
if(is.matrix(object$geno[[i]]$map)) {
for(j in 1:2)
object$geno[[i]]$map[j,] <- object$geno[[i]]$map[j,] - object$geno[[i]]$map[j,1] + offset[i]
} else {
object$geno[[i]]$map <- object$geno[[i]]$map - object$geno[[i]]$map[1] + offset[i]
}
}
} else if("map" %in% class(object)) {
if(length(offset) == 1) offset <- rep(offset, length(object))
else if(length(offset) != length(object))
stop("offset must have length 1 or n.chr (", length(object), ")")
for(i in seq(along=object)) {
if(is.matrix(object[[i]])) {
for(j in 1:2)
object[[i]][j,] <- object[[i]][j,] - object[[i]][j,1] + offset[i]
} else {
object[[i]] <- object[[i]] - object[[i]][1] + offset[i]
}
}
} else
stop("shiftmap works only for objects of class \"cross\" or \"map\".")
object
}
######################################################################
# switch alleles in a cross
######################################################################
switchAlleles <-
function(cross, markers, switch=c("AB","CD","ABCD", "parents"))
{
type <- class(cross)[1]
switch <- match.arg(switch)
if(type %in% c("bc", "risib", "riself", "dh")) {
if(switch != "AB")
warning("Using switch = \"AB\".")
found <- rep(FALSE, length(markers))
for(i in 1:nchr(cross)) {
cn <- colnames(cross$geno[[i]]$data)
m <- match(markers, cn)
if(any(!is.na(m))) {
found[!is.na(m)] <- TRUE
for(j in m[!is.na(m)]) {
g <- cross$geno[[i]]$data[,j]
cross$geno[[i]]$data[!is.na(g) & g==1,j] <- 2
cross$geno[[i]]$data[!is.na(g) & g==2,j] <- 1
}
}
}
}
else if(type=="f2") {
if(switch != "AB")
warning("Using switch = \"AB\".")
found <- rep(FALSE, length(markers))
for(i in 1:nchr(cross)) {
cn <- colnames(cross$geno[[i]]$data)
m <- match(markers, cn)
if(any(!is.na(m))) {
found[!is.na(m)] <- TRUE
for(j in m[!is.na(m)]) {
g <- cross$geno[[i]]$data[,j]
cross$geno[[i]]$data[!is.na(g) & g==1,j] <- 3
cross$geno[[i]]$data[!is.na(g) & g==3,j] <- 1
cross$geno[[i]]$data[!is.na(g) & g==4,j] <- 5
cross$geno[[i]]$data[!is.na(g) & g==5,j] <- 4
}
}
}
}
else if(type=="4way") {
if(switch=="AB")
newg <- c(2,1,4,3,6,5,7,8,10,9,12,11,14,13)
else if(switch=="CD")
newg <- c(3,4,1,2,5,6,8,7,10,9,13,14,11,12)
else if(switch=="ABCD")
newg <- c(4,3,2,1,6,5,8,7,9,10,14,13,12,11)
else # switch parents
newg <- c(1,3,2,4,7,8,5,6,9,10,11,13,12,14)
found <- rep(FALSE, length(markers))
for(i in 1:nchr(cross)) {
cn <- colnames(cross$geno[[i]]$data)
m <- match(markers, cn)
if(any(!is.na(m))) {
found[!is.na(m)] <- TRUE
for(j in m[!is.na(m)]) {
g <- cross$geno[[i]]$data[,j]
for(k in 1:14)
cross$geno[[i]]$data[!is.na(g) & g==k,j] <- newg[k]
}
}
}
}
if(any(!found))
warning("Some markers not found: ", paste(markers[!found], collapse=" "))
clean(cross)
}
######################################################################
#
# nqrank: Convert a set of quantitative values to the corresponding
# normal quantiles (preserving mean and SD)
#
######################################################################
nqrank <-
function(x, jitter=FALSE)
{
y <- x[!is.na(x)]
themean <- mean(y, na.rm=TRUE)
thesd <- sd(y, na.rm=TRUE)
y[y == Inf] <- max(y[y<Inf])+10
y[y == -Inf] <- min(y[y > -Inf]) - 10
if(jitter)
y <- rank(y+runif(length(y))/(sd(y)*10^8))
else y <- rank(y)
x[!is.na(x)] <- qnorm((y-0.5)/length(y))
x*thesd/sd(x, na.rm=TRUE)-mean(x,na.rm=TRUE)+themean
}
######################################################################
#
# cleanGeno: omit genotypes that are possibly in error, as indicated
# by apparent double-crossovers separated by a distance of
# no more than maxdist and having no more than maxmark
# interior typed markers
#
######################################################################
cleanGeno <-
function(cross, chr, maxdist=2.5, maxmark=2, verbose=TRUE)
{
if(class(cross)[1] != "bc")
stop("This function currently only works for a backcross.")
if(!missing(chr)) cleaned <- subset(cross, chr=chr)
else cleaned <- cross
thechr <- names(cleaned$geno)
totdrop <- 0
maxmaxdist <- max(maxdist)
for(i in thechr) {
xoloc <- locateXO(cleaned, chr=i, full.info=TRUE)
nxo <- sapply(xoloc, function(a) if(is.matrix(a)) return(nrow(a)) else return(0))
g <- pull.geno(cleaned, chr=i)
ndrop <- 0
for(j in which(nxo > 1)) {
maxd <- xoloc[[j]][-1,"right"] - xoloc[[j]][-nrow(xoloc[[j]]),"left"]
wh <- maxd <= maxmaxdist
if(any(wh)) {
for(k in which(wh)) {
nt <- sum(!is.na(g[j,(xoloc[[j]][k,"ileft"]+1):(xoloc[[j]][k+1,"iright"]-1)]))
if(nt > 0 && any(nt <= maxmark & maxd[k] < maxdist)) {
cleaned$geno[[i]]$data[j,(xoloc[[j]][k,"ileft"]+1):(xoloc[[j]][k+1,"iright"]-1)] <- NA
ndrop <- ndrop + nt
totdrop <- totdrop + nt
}
}
}
}
if(verbose && ndrop > 0) {
totgen <- sum(ntyped(subset(cross, chr=i)))
cat(" ---Dropping ", ndrop, " genotypes (out of ", totgen, ") on chr ", i, "\n", sep="")
}
}
if(verbose && nchr(cleaned)>1 && totdrop > 0) {
totgen <- sum(ntyped(subset(cross, chr=thechr)))
cat(" ---Dropped ", totdrop, " genotypes (out of ", totgen, ") in total\n", sep="")
}
for(i in names(cleaned$geno))
cross$geno[[i]] <- cleaned$geno[[i]]
cross
}
######################################################################
# typingGap: calculate gaps between typed markers
######################################################################
typingGap <-
function(cross, chr, terminal=FALSE)
{
if(!missing(chr))
cross <- subset(cross, chr)
n.ind <- nind(cross)
n.chr <- nchr(cross)
gaps <- matrix(nrow=n.ind, ncol=n.chr)
colnames(gaps) <- names(cross$geno)
for(i in 1:n.chr) {
map <- cross$geno[[i]]$map
map <- c(map[1], map, map[length(map)])
if(is.matrix(map)) stop("This function can't currently handle sex-specific maps.")
if(terminal) # just look at terminal gaps
gaps[,i] <- apply(cbind(1,cross$geno[[i]]$data,1), 1,
function(a,b) {d <- diff(b[!is.na(a)]); max(d[c(1,length(d))]) }, map)
else
gaps[,i] <- apply(cbind(1,cross$geno[[i]]$data,1), 1,
function(a,b) max(diff(b[!is.na(a)])), map)
}
if(n.chr==1) gaps <- as.numeric(gaps)
gaps
}
######################################################################
# calcPermPval
#
# calculate permutation pvalues for summary.scanone()
######################################################################
calcPermPval <-
function(peaks, perms)
{
if(!is.matrix(peaks))
peaks <- as.matrix(peaks)
if(!is.matrix(perms))
perms <- as.matrix(perms)
ncol.peaks <- ncol(peaks)
nrow.peaks <- nrow(peaks)
n.perms <- nrow(perms)
if(ncol.peaks != ncol(perms))
stop("ncol(peaks) != ncol(perms)")
pval <- .C("R_calcPermPval",
as.double(peaks),
as.integer(ncol.peaks),
as.integer(nrow.peaks),
as.double(perms),
as.integer(n.perms),
pval=as.double(rep(0,ncol.peaks*nrow.peaks)),
PACKAGE="qtl")$pval
matrix(pval, ncol=ncol.peaks, nrow=nrow.peaks)
}
######################################################################
# phenames: pull out phenotype names
######################################################################
phenames <-
function(cross)
colnames(cross$pheno)
# end of util.R
byandell/qtl documentation built on May 13, 2019, 9:28 a.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.
#####################################################################
#
# util.R
#
# copyright (c) 2001-2012, Karl W Broman
# [find.pheno, find.flanking, and a modification to create.map
# from Brian Yandell]
# last modified Nov, 2012
# first written Feb, 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: markernames, c.cross, create.map,
# clean, clean.cross, drop.nullmarkers, nullmarkers
# drop.markers, pull.markers, drop.dupmarkers
# geno.table, genotab.em
# mf.k, mf.h, imf.k, imf.h, mf.cf, imf.cf, mf.m, imf.m,
# mf.stahl, imf.stahl
# switch.order, makeSSmap,
# subset.cross, fill.geno, checkcovar, find.marker,
# find.pseudomarker,
# adjust.rf.ri, lodint, bayesint,
# comparecrosses, movemarker, summary.map (aka summaryMap),
# print.summary.map, find.pheno,
# convert, convert.scanone, convert.scantwo
# find.flanking, strip.partials, comparegeno
# qtlversion, locateXO, jittermap, getid,
# find.markerpos, geno.crosstab, LikePheVector,
# matchchr, convert2sa, charround, testchr,
# scantwoperm2scanoneperm, subset.map, [.map, [.cross,
# findDupMarkers, convert2riself, convert2risib,
# switchAlleles, nqrank, cleanGeno, typingGap,
# calcPermPval, phenames
#
######################################################################
######################################################################
#
# markernames
#
# pull out the marker names for selected chromosomes as one big vector
#
######################################################################
markernames <-
function(cross, chr)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
if(!missing(chr)) cross <- subset(cross, chr=chr)
temp <- unlist(lapply(cross$geno, function(a) colnames(a$data)))
names(temp) <- NULL
temp
}
######################################################################
#
# chrnames
#
# pull out the chrnames for a cross
#
######################################################################
chrnames <-
function(cross)
{
names(cross$geno)
}
######################################################################
#
# create.map
#
# create a new map with inserted inter-marker locations
#
# Note: map is a vector or a matrix with 2 rows
#
# stepwidth = "fixed" is the original R/qtl version;
# stepwidth="variable" is for Brian Yandell and the qtlbim package
# stepwidth="max" creates the minimal number of inserted pseudomarkers
# to have the maximum stepwidth = step
######################################################################
create.map <-
function(map, step, off.end, stepwidth = c("fixed", "variable", "max"))
{
stepwidth <- match.arg(stepwidth)
if(step<0 || off.end<0) stop("step and off.end must be > 0.")
if(!is.matrix(map)) { # sex-ave map
if(stepwidth == "variable") {
if(off.end > 0) {
tmp <- names(map)
## Append and prepend by off.end value (exact here, no roundoff).
map <- c(map[1] - off.end, map, map[length(map)] + off.end)
names(map) <- c("loc000", tmp, "loc999")
}
if(step == 0)
return(unclass(map))
## Determine differences and expansion vector.
dif <- diff(map)
expand <- pmax(1, floor(dif / step))
## Create pseudomarker map.
a <- min(map) + cumsum(c(0, rep(dif / expand, expand)))
## Names are marker names or locNNN.
namesa <- paste("loc", seq(length(a)), sep = "")
namesa[cumsum(c(1, expand))] <- names(map)
names(a) <- namesa
return(unclass(a))
}
if(stepwidth == "max") {
if(off.end > 0) {
toadd <- c(map[1] - off.end, map[length(map)]+off.end)
if(step==0) {
names(toadd) <- paste("loc", 1:2, sep="")
map <- sort(c(map, toadd))
return(unclass(map))
}
nmap <- c(map[1] - off.end, map, map[length(map)]+off.end)
}
else {
nmap <- map
toadd <- NULL
}
if(step==0 || (length(map)==1 && off.end==0)) return(unclass(map))
d <- diff(nmap)
nadd <- ceiling(d/step)-1
if(sum(nadd) > 0) {
for(j in 1:(length(nmap)-1)) {
if(nadd[j]>0)
toadd <- c(toadd, seq(nmap[j], nmap[j+1], len=nadd[j]+2)[-c(1,nadd[j]+2)])
}
}
if(length(toadd) > 0) {
names(toadd) <- paste("loc", 1:length(toadd), sep="")
map <- sort(c(map, toadd))
}
return(unclass(map))
}
if(length(map) == 1) { # just one marker!
if(off.end==0) {
if(step == 0) step <- 1
nam <- names(map)
map <- c(map,map+step)
names(map) <- c(nam,paste("loc",step,sep=""))
}
else {
if(step==0) m <- c(-off.end,off.end)
else m <- seq(-off.end,off.end,by=step)
m <- m[m!=0]
names(m) <- paste("loc",m,sep="")
map <- sort(c(m+map,map))
}
return(map)
}
minloc <- min(map)
map <- map-minloc
if(step==0 && off.end==0) return(map+minloc)
else if(step==0 && off.end > 0) {
a <- c(floor(min(map)-off.end),ceiling(max(map)+off.end))
names(a) <- paste("loc", a, sep="")
return(sort(c(a,map))+minloc)
}
else if(step>0 && off.end == 0) {
a <- seq(floor(min(map)),max(map),
by = step)
if(any(is.na(match(a, map)))) {
a <- a[is.na(match(a,map))]
names(a) <- paste("loc",a,sep="")
return(sort(c(a,map))+minloc)
}
else return(map+minloc)
}
else {
a <- seq(floor(min(map)-off.end),ceiling(max(map)+off.end+step),
by = step)
a <- a[is.na(match(a,map))]
# no more than one point above max(map)+off.end
z <- (seq(along=a))[a >= max(map)+off.end]
if(length(z) > 1) a <- a[-z[-1]]
names(a) <- paste("loc",a,sep="")
return(sort(c(a,map))+minloc)
}
} # end sex-ave map
else { # sex-specific map
if(stepwidth == "variable") {
if(off.end > 0) {
tmp <- colnames(map)
map <- cbind(map[, 1] - off.end, map, map[, ncol(map)] + off.end)
dimnames(map) <- list(NULL, c("loc000", tmp, "loc999"))
}
if(step == 0)
return(unclass(map))
## Determine differences and expansion vector.
dif <- diff(map[1, ])
expand <- pmax(1, floor(dif / step))
## Create pseudomarker map.
a <- min(map[1, ]) + cumsum(c(0, rep(dif / expand, expand)))
b <- min(map[2, ]) + cumsum(c(0, rep(diff(map[2, ]) / expand, expand)))
namesa <- paste("loc", seq(length(a)), sep = "")
namesa[cumsum(c(1, expand))] <- dimnames(map)[[2]]
map <- rbind(a,b)
dimnames(map) <- list(NULL, namesa)
return(unclass(map))
}
if(stepwidth == "max") {
if(step==0 && off.end==0) return(unclass(map))
if(step==0 && off.end>0) {
if(ncol(map)==1) { # only one marker; assume equal recomb in sexes
L1 <- L2 <- 1
}
else {
L1 <- diff(range(map[1,]))
L2 <- diff(range(map[2,]))
}
nam <- colnames(map)
nmap1 <- c(map[1,1]-off.end, map[1,], map[1,ncol(map)]+off.end)
nmap2 <- c(map[2,1]-off.end*L2/L1, map[2,], map[2,ncol(map)]+off.end*L2/L1)
map <- rbind(nmap1, nmap2)
colnames(map) <- c("loc1", nam, "loc2")
return(unclass(map))
}
if(ncol(map)==1) L1 <- L2 <- 1
else {
L1 <- diff(range(map[1,]))
L2 <- diff(range(map[2,]))
}
nam <- colnames(map)
if(off.end > 0) {
toadd1 <- c(map[1,1] - off.end, map[1,ncol(map)]+off.end)
toadd2 <- c(map[2,1] + off.end*L2/L1, map[2,ncol(map)]+off.end*L2/L1)
neword <- order(c(map[1,], toadd1))
nmap1 <- c(map[1,], toadd1)[neword]
nmap2 <- c(map[2,], toadd2)[neword]
}
else {
nmap1 <- map[1,]
nmap2 <- map[2,]
toadd1 <- toadd2 <- NULL
}
d <- diff(nmap1)
nadd <- ceiling(d/step)-1
if(sum(nadd) > 0) {
for(j in 1:(length(nmap1)-1)) {
if(nadd[j]>0) {
toadd1 <- c(toadd1, seq(nmap1[j], nmap1[j+1], len=nadd[j]+2)[-c(1,nadd[j]+2)])
toadd2 <- c(toadd2, seq(nmap2[j], nmap2[j+1], len=nadd[j]+2)[-c(1,nadd[j]+2)])
}
}
}
newnam <- paste("loc", 1:length(toadd1), sep="")
toadd1 <- sort(toadd1)
toadd2 <- sort(toadd2)
neword <- order(c(map[1,], toadd1))
nmap1 <- c(map[1,], toadd1)[neword]
nmap2 <- c(map[2,], toadd2)[neword]
map <- rbind(nmap1, nmap2)
colnames(map) <- c(nam, newnam)[neword]
return(unclass(map))
}
minloc <- c(min(map[1,]),min(map[2,]))
map <- unclass(map-minloc)
markernames <- colnames(map)
if(step==0 && off.end==0) return(map+minloc)
else if(step==0 && off.end > 0) {
map <- map+minloc
if(ncol(map)==1) { # only one marker; assume equal recomb in sexes
L1 <- L2 <- 1
}
else {
L1 <- diff(range(map[1,]))
L2 <- diff(range(map[2,]))
}
nam <- colnames(map)
nmap1 <- c(map[1,1]-off.end, map[1,], map[1,ncol(map)]+off.end)
nmap2 <- c(map[2,1]-off.end*L2/L1, map[2,], map[2,ncol(map)]+off.end*L2/L1)
map <- rbind(nmap1, nmap2)
colnames(map) <- c("loc1", nam, "loc2")
return(map)
}
else if(step>0 && off.end == 0) {
if(ncol(map)==1) return(map+minloc)
a <- seq(floor(min(map[1,])),max(map[1,]),
by = step)
a <- a[is.na(match(a,map[1,]))]
if(length(a)==0) return(map+minloc)
b <- sapply(a,function(x,y,z) {
ZZ <- min((seq(along=y))[y > x])
(x-y[ZZ-1])/(y[ZZ]-y[ZZ-1])*(z[ZZ]-z[ZZ-1])+z[ZZ-1] }, map[1,],map[2,])
m1 <- c(a,map[1,])
m2 <- c(b,map[2,])
names(m1) <- names(m2) <- c(paste("loc",a,sep=""),markernames)
return(rbind(sort(m1),sort(m2))+minloc)
}
else {
a <- seq(floor(min(map[1,])-off.end),ceiling(max(map[1,])+off.end+step),
by = step)
a <- a[is.na(match(a,map[1,]))]
# no more than one point above max(map)+off.end
z <- (seq(along=a))[a >= max(map[1,])+off.end]
if(length(z) > 1) a <- a[-z[-1]]
b <- sapply(a,function(x,y,z,ml) {
if(x < min(y)) {
return(min(z) - (min(y)-x)/diff(range(y))*diff(range(z)) - ml)
}
else if(x > max(y)) {
return(max(z) + (x - max(y))/diff(range(y))*diff(range(z)) - ml)
}
else {
ZZ <- min((seq(along=y))[y > x])
(x-y[ZZ-1])/(y[ZZ]-y[ZZ-1])*(z[ZZ]-z[ZZ-1])+z[ZZ-1]
}
}, map[1,],map[2,], minloc[2])
m1 <- c(a,map[1,])
m2 <- c(b,map[2,])
names(m1) <- names(m2) <- c(paste("loc",a,sep=""),markernames)
return(rbind(sort(m1),sort(m2))+minloc)
}
}
}
######################################################################
# clean functions
######################################################################
clean <-
function(object, ...)
UseMethod("clean")
######################################################################
#
# clean.cross
#
# remove all of the extraneous stuff from a cross object, to get back
# to just the data
#
######################################################################
clean.cross <-
function(object, ...)
{
if(!any(class(object) == "cross"))
stop("Input should have class \"cross\".")
cross2 <- list(geno=object$geno,pheno=object$pheno)
if("cross" %in% names(object))
cross2$cross <- object$cross
if("founderGeno" %in% names(object))
cross2$founderGeno <- object$founderGeno
if(!is.null(attr(object, "alleles")))
attr(cross2, "alleles") <- attr(object, "alleles")
if(!is.null(attr(object, "scheme")))
attr(cross2, "scheme") <- attr(object, "scheme")
for(i in 1:length(object$geno)) {
cross2$geno[[i]] <- list(data=object$geno[[i]]$data,
map=object$geno[[i]]$map)
class(cross2$geno[[i]]) <- class(object$geno[[i]])
}
class(cross2) <- class(object)
cross2
}
######################################################################
#
# drop.qtlgeno
#
# remove any QTLs from the genotype data and the genetic maps
# from data simulated via sim.cross. (They all have names "QTL*")
#
######################################################################
#drop.qtlgeno <-
#function(cross)
#{
# n.chr <- nchr(cross)
# mar.names <- lapply(cross$geno, function(a) {
# m <- a$map
# if(is.matrix(m)) return(colnames(m))
# else return(names(m)) } )
#
# for(i in 1:n.chr) {
# o <- grep("^QTL[0-9]+",mar.names[[i]])
# if(length(o) != 0) {
# cross$geno[[i]]$data <- cross$geno[[i]]$data[,-o,drop=FALSE]
# if(is.matrix(cross$geno[[i]]$map))
# cross$geno[[i]]$map <- cross$geno[[i]]$map[,-o,drop=FALSE]
# else
# cross$geno[[i]]$map <- cross$geno[[i]]$map[-o]
# }
# }
# cross
#}
######################################################################
#
# drop.nullmarkers
#
# remove markers that have no genotype data from the data matrix and
# genetic maps
#
######################################################################
drop.nullmarkers <-
function(cross)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
n.chr <- nchr(cross)
keep.chr <- rep(TRUE,n.chr)
for(i in 1:n.chr) {
o <- !apply(cross$geno[[i]]$data,2,function(a) sum(!is.na(a)))
if(any(o)) { # remove from genotype data and map
mn.drop <- colnames(cross$geno[[i]]$data)[o]
if(length(mn.drop) == ncol(cross$geno[[i]]$data))
keep.chr[i] <- FALSE # removing all markers from this chromosome
cross$geno[[i]]$data <- cross$geno[[i]]$data[,!o,drop=FALSE]
if(is.matrix(cross$geno[[i]]$map))
cross$geno[[i]]$map <- cross$geno[[i]]$map[,!o,drop=FALSE]
else
cross$geno[[i]]$map <- cross$geno[[i]]$map[!o]
# results of calc.genoprob
if("prob" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$prob))
cross$geno[[i]]$prob <- cross$geno[[i]]$prob[,-o,,drop=FALSE]
}
# results of argmax.geno
if("argmax" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$argmax))
cross$geno[[i]]$argmax <- cross$geno[[i]]$argmax[,-o,drop=FALSE]
}
# results of sim.geno
if("draws" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$draws))
cross$geno[[i]]$draws <- cross$geno[[i]]$draws[,-o,,drop=FALSE]
}
# results of est.rf
if("rf" %in% names(cross)) {
attrib <- attributes(cross$rf)
o <- match(mn.drop,colnames(cross$rf))
cross$rf <- cross$rf[-o,-o]
if("onlylod" %in% names(attrib)) # save the onlylod attribute if its there
attr(cross$rf, "onlylod") <- attrib$onlylod
}
}
}
cross$geno <- cross$geno[keep.chr]
if("founderGeno" %in% names(cross))
cross$founderGeno <- cross$founderGeno[,markernames(cross)]
cross
}
######################################################################
#
# nullmarkers
#
# identify markers that have no genotype data from the data matrix and
# genetic maps
#
######################################################################
nullmarkers <-
function(cross)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
n.chr <- nchr(cross)
keep.chr <- rep(TRUE,n.chr)
all2drop <- NULL
for(i in 1:n.chr) {
o <- !apply(cross$geno[[i]]$data,2,function(a) sum(!is.na(a)))
if(any(o)) { # remove from genotype data and map
mn.drop <- colnames(cross$geno[[i]]$data)[o]
all2drop <- c(all2drop, mn.drop)
}
}
all2drop
}
######################################################################
#
# drop.markers
#
# remove a vector of markers from the data matrix and genetic maps
#
######################################################################
drop.markers <-
function(cross, markers)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
n.chr <- nchr(cross)
keep.chr <- rep(TRUE,n.chr)
found <- rep(FALSE, length(markers))
for(i in 1:n.chr) {
# find markers on this chromosome
o <- match(markers,colnames(cross$geno[[i]]$data))
found[!is.na(o)] <- TRUE
o <- o[!is.na(o)]
a <- rep(FALSE,ncol(cross$geno[[i]]$data))
a[o] <- TRUE
o <- a
if(any(o)) { # remove from genotype data and map
mn.drop <- colnames(cross$geno[[i]]$data)[o]
if(length(mn.drop) == ncol(cross$geno[[i]]$data))
keep.chr[i] <- FALSE # removing all markers from this chromosome
cross$geno[[i]]$data <- cross$geno[[i]]$data[,!o,drop=FALSE]
if(is.matrix(cross$geno[[i]]$map))
cross$geno[[i]]$map <- cross$geno[[i]]$map[,!o,drop=FALSE]
else
cross$geno[[i]]$map <- cross$geno[[i]]$map[!o]
# results of calc.genoprob
if("prob" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$prob))
cross$geno[[i]]$prob <- cross$geno[[i]]$prob[,-o,,drop=FALSE]
}
# results of argmax.geno
if("argmax" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$argmax))
cross$geno[[i]]$argmax <- cross$geno[[i]]$argmax[,-o,drop=FALSE]
}
# results of sim.geno
if("draws" %in% names(cross$geno[[i]])) {
o <- match(mn.drop,colnames(cross$geno[[i]]$draws))
cross$geno[[i]]$draws <- cross$geno[[i]]$draws[,-o,,drop=FALSE]
}
# results of est.rf
if("rf" %in% names(cross)) {
attrib <- attributes(cross$rf)
o <- match(mn.drop,colnames(cross$rf))
cross$rf <- cross$rf[-o,-o]
if("onlylod" %in% names(attrib))
attr(cross$rf, "onlylod") <- attrib$onlylod
}
}
}
if(any(!found))
warning("Markers not found: ", paste(markers[!found],collapse=" "))
cross$geno <- cross$geno[keep.chr]
if("founderGeno" %in% names(cross))
cross$founderGeno <- cross$founderGeno[,markernames(cross)]
cross
}
######################################################################
# pull.markers
#
# like drop.markers, but retain just those indicated
######################################################################
pull.markers <-
function(cross, markers)
{
mn <- markernames(cross)
m <- match(markers, mn)
if(any(is.na(m)))
warning("Some markers couldn't be found: ", paste(markers[is.na(m)], collapse=" "))
drop.markers(cross, mn[is.na(match(mn, markers))])
}
######################################################################
# drop.dupmarkers
#
# drop duplicate markers, retaining the consensus genotypes
######################################################################
drop.dupmarkers <-
function(cross, verbose=TRUE)
{
mn <- markernames(cross)
tab <- table(mn)
if(all(tab==1)) {
if(verbose) cat("No duplicate markers.\n")
return(cross)
}
dup <- names(tab[tab > 1])
if(verbose) cat(" ", length(dup), "duplicate markers\n")
# get consensus genotypes
g <- pull.geno(cross)[,!is.na(match(mn, dup))]
ng.omitted <- rep(NA, length(dup))
tot.omitted <- 0
nmar.omitted <- 0
for(i in seq(along=dup)) {
gg <- g[,colnames(g)==dup[i],drop=FALSE]
res <- apply(gg, 1, function(a)
{
if(all(is.na(a))) return(c(NA, 0))
a <- unique(a[!is.na(a)])
if(length(a)==1) return(c(a, 0))
return(c(NA, 1))
} )
if(verbose>1) {
cat(" ", dup[i], ":\t", sum(res[2,]), " genotypes omitted\n", sep="")
if(sum(res[2,]) > 1) cat(" ", paste(which(res[2,]>0), collapse=" "), "\n")
}
tot.omitted <- tot.omitted + sum(res[2,])
flag <- FALSE
for(j in seq(along=cross$geno)) {
mn <- colnames(cross$geno[[j]]$data)
wh <- mn==dup[i]
if(!any(wh)) next
wh <- which(wh)
if(!flag) {
flag <- TRUE
cross$geno[[j]]$data[,wh[1]] <- res[1,]
if(length(wh)==1) next
else wh <- wh[-1]
}
if(length(wh) > 0) {
nmar.omitted <- nmar.omitted + length(wh)
cross$geno[[j]]$data <- cross$geno[[j]]$data[,-wh,drop=FALSE]
if(is.matrix(cross$geno[[j]]$map))
cross$geno[[j]]$map <- cross$geno[[j]]$map[,-wh,drop=FALSE]
else
cross$geno[[j]]$map <- cross$geno[[j]]$map[-wh]
}
}
}
if(verbose) {
cat(" Total genotypes omitted:", tot.omitted, "\n")
cat(" Total markers omitted: ", nmar.omitted, "\n")
}
if("founderGeno" %in% names(cross))
cross$founderGeno <- cross$founderGeno[,markernames(cross)]
clean(cross)
}
######################################################################
#
# geno.table
#
# create table showing observed numbers of individuals with each
# genotype at each marker
#
######################################################################
geno.table <-
function(cross, chr, scanone.output=FALSE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
if(!missing(chr))
cross <- subset(cross, chr=chr)
n.chr <- nchr(cross)
type <- class(cross)[1]
is.bcs <- type == "bcsft"
cross.scheme <- attr(cross, "scheme")
if(is.bcs)
is.bcs <- (cross.scheme[2] == 0)
chrtype <- sapply(cross$geno, class)
allchrtype <- rep(chrtype, nmar(cross))
chrname <- names(cross$geno)
allchrname <- rep(chrname, nmar(cross))
## Actually plan to have our own geno.table.bcsft routine.
if(type == "f2" || (type == "bcsft" && !is.bcs)) {
n.gen <- 5
temp <- getgenonames("f2", "A", cross.attr=attributes(cross))
gen.names <- c(temp, paste("not", temp[c(3,1)]))
}
else if(type %in% c("bc", "riself", "risib", "dh", "bcsft")) {
n.gen <- 2
gen.names <- getgenonames(type, "A", cross.attr=attributes(cross))
}
else if(type == "4way") {
n.gen <- 14
temp <- getgenonames("4way", "A", cross.attr=attributes(cross))
gen.names <- c(temp,
paste(temp[c(1,3)], collapse="/"),
paste(temp[c(2,4)], collapse="/"),
paste(temp[c(1,2)], collapse="/"),
paste(temp[c(3,4)], collapse="/"),
paste(temp[c(1,4)], collapse="/"),
paste(temp[c(2,3)], collapse="/"),
paste("not", temp[1]),
paste("not", temp[2]),
paste("not", temp[3]),
paste("not", temp[4]))
gen.names[5:8] <- substr(temp[c(1,2,1,3)], c(1,1,2,2), c(1,1,2,2))
}
else stop("Unknown cross type: ",type)
res <- lapply(cross$geno, function(a,ngen) {
a <- a$data; a[is.na(a)] <- 0
apply(a,2,function(b,ngen) table(factor(b,levels=0:ngen)),ngen)
},n.gen)
results <- NULL
for(i in 1:length(res))
results <- rbind(results,t(res[[i]]))
colnames(results) <- c("missing",gen.names)
rownames(results) <- unlist(lapply(cross$geno,function(a) colnames(a$data)))
pval <- rep(NA,nrow(results))
if(type %in% c("bc","risib","riself","dh") || (type=="bcsft" & is.bcs)) {
sexpgm <- getsex(cross)
if((type == "bc" || type=="bcsft") && any(chrtype == "X") && !is.null(sexpgm$sex) && any(sexpgm$sex==1)) {
for(i in which(allchrtype=="A")) {
x <- results[i,2:3]
if(sum(x) > 0)
pval[i] <- chisq.test(x,p=c(0.5,0.5))$p.value
else pval[i] <- 1
}
gn <- getgenonames("bc","X", "full", sexpgm, attributes(cross))
wh <- which(is.na(match(gn, colnames(results))))
temp <- matrix(0, nrow=nrow(results), ncol=length(wh))
colnames(temp) <- gn[wh]
results <- cbind(results, temp)
for(i in which(chrtype=="X")) {
dat <- reviseXdata("bc", "full", sexpgm, geno=cross$geno[[i]]$data,
cross.attr=attributes(cross))
dat[is.na(dat)] <- 0
tab <- t(apply(dat, 2, function(x) table(factor(x, levels=0:length(gn)))))
colnames(tab) <- c("missing", gn)
results[allchrname==chrname[i],] <- 0
results[allchrname==chrname[i],colnames(tab)] <- tab
for(j in 1:ncol(dat)) {
stat <- apply(table(sexpgm$sex, cross$geno[[i]]$data[,j]),1,
function(a) if(length(a) > 1 && sum(a)>0) return(chisq.test(a,p=c(0.5,0.5))$stat)
else return(0))
pval[allchrname==chrname[i]][j] <- 1-pchisq(sum(stat),length(stat))
}
}
results <- cbind(results, P.value=pval)
}
else {
for(i in 1:length(pval)) {
x <- results[i,2:3]
if(sum(x) > 0)
pval[i] <- chisq.test(x,p=c(0.5,0.5))$p.value
else
pval[i] <- 1
}
results <- cbind(results, P.value=pval)
}
}
else if(type == "f2" || (type == "bcsft" && !is.bcs)) {
sexpgm <- getsex(cross)
## F2 with set initial genotype probabilities.
init.geno <- c(0.25,0.5,0.25,0.75,0.75)
## BCsFt initial genotype probabilities need to be computed.
if(type == "bcsft") {
ret <- .C("genotab_em_bcsft",
as.integer(cross.scheme),
init.geno = as.double(init.geno))
init.geno <- ret$init.geno
}
for(i in which(allchrtype=="A")) {
dat <- results[i,2:6]
if(sum(dat)==0) pval[i] <- 1
else if(dat[4]==0 && dat[5]==0)
pval[i] <- chisq.test(dat[1:3], p=init.geno[1:3] )$p.value
else if(all(dat[2:4]==0))
pval[i] <- chisq.test(dat[c(1,5)],p=init.geno[c(1,5)])$p.value
else if(all(dat[c(1,2,5)]==0))
pval[i] <- chisq.test(dat[3:4], p=init.geno[3:4] )$p.value
else { # this is harder: some dominant and some not
pval[i] <- genotab.em(dat, init.geno)
}
}
for(i in which(chrtype=="X")) {
gn <- getgenonames("f2","X","full", getsex(cross), attributes(cross))
wh <- which(is.na(match(gn, colnames(results))))
temp <- matrix(0, nrow=nrow(results), ncol=length(wh))
colnames(temp) <- gn[wh]
results <- cbind(results, temp)
dat <- reviseXdata("f2", "full", sexpgm, geno=cross$geno[[i]]$data,
cross.attr=attributes(cross))
dat[is.na(dat)] <- 0
tab <- t(apply(dat, 2, function(x) table(factor(x, levels=0:length(gn)))))
colnames(tab) <- c("missing", gn)
results[allchrname==chrname[i],] <- 0
results[allchrname==chrname[i],colnames(tab)] <- tab
cn <- colnames(results)
f <- grep("f$", cn)
r <- grep("r$", cn)
if(length(f)>0 && length(r)>0) {
results <- results[,c(1:2,f,3,r,(1:ncol(results))[-c(1:3,f,r)])]
colnames(results) <- cn[c(1:2,f,3,r,(1:ncol(results))[-c(1:3,f,r)])]
}
sex <- sexpgm$sex
pgm <- sexpgm$pgm
for(j in 1:ncol(dat)) {
g <- cross$geno[[i]]$data[,j]
if(!is.null(sex)) {
if(!is.null(pgm))
g <- matrix(as.numeric(table(g,sex,pgm)),ncol=2,byrow=TRUE)
else
g <- matrix(as.numeric(table(g,sex)),ncol=2,byrow=TRUE)
}
else {
if(!is.null(pgm)) {
g <- matrix(as.numeric(table(g,pgm)),ncol=2,byrow=TRUE)
}
else g <- matrix(table(g),ncol=2,byrow=TRUE)
}
stat <- apply(g, 1,
function(a) if(sum(a)>0) return(chisq.test(a,p=c(0.5,0.5))$stat)
else return(0))
pval[allchrname==chrname[i]][j] <- 1-pchisq(sum(stat),length(stat))
}
}
results <- cbind(results, P.value=pval)
}
else if(type == "4way") {
for(i in 1:length(pval)) {
x <- results[i,2:5]
y <- results[i,-(1:5)]
if(sum(x) > 0 && sum(y)==0)
pval[i] <- chisq.test(x,p=c(0.25,0.25,0.25,0.25))$p.value
else {
if(allchrtype[i] == "A") {
res <- results[i,-1]
if(all(res[-c(1,11)]==0)) # AC/not AC
pval[i] <- chisq.test(res[c(1,11)], p=c(0.25, 0.75))$p.value
else if(all(res[-c(2,12)]==0)) # BC/not BC
pval[i] <- chisq.test(res[c(2,12)], p=c(0.25, 0.75))$p.value
else if(all(res[-c(3,13)]==0)) # AD/not AD
pval[i] <- chisq.test(res[c(3,13)], p=c(0.25, 0.75))$p.value
else if(all(res[-c(4,14)]==0)) # BD/not BD
pval[i] <- chisq.test(res[c(4,14)], p=c(0.25, 0.75))$p.value
else if(all(res[-c(5,6)]==0)) # A/B
pval[i] <- chisq.test(res[c(5,6)], p=c(0.5, 0.5))$p.value
else if(all(res[-c(7,8)]==0)) # C/D
pval[i] <- chisq.test(res[c(7,8)], p=c(0.5, 0.5))$p.value
else if(all(res[-c(9,10)]==0)) # AC/BD or AD/BC
pval[i] <- chisq.test(res[c(9,10)], p=c(0.5, 0.5))$p.value
else if(all(res[-c(2,4,5)]==0)) # BC/BD/A
pval[i] <- chisq.test(res[c(2,4,5)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(1,3,6)]==0)) # AC/AD/B
pval[i] <- chisq.test(res[c(1,3,6)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(3,4,7)]==0)) # AD/BD/C
pval[i] <- chisq.test(res[c(3,4,7)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(1,2,8)]==0)) # AC/BC/D
pval[i] <- chisq.test(res[c(1,2,8)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(2,3,9)]==0)) # AC/BD or AD or BC
pval[i] <- chisq.test(res[c(2,3,9)], p=c(0.25, 0.25, 0.5))$p.value
else if(all(res[-c(1,4,10)]==0)) # AD/BC or AC or BD
pval[i] <- chisq.test(res[c(1,4,10)], p=c(0.25, 0.25, 0.5))$p.value
}
}
}
results <- cbind(results, P.value=pval)
}
if(!scanone.output)
return(data.frame(chr=rep(names(cross$geno),nmar(cross)),results, stringsAsFactors=TRUE))
themap <- pull.map(cross)
if(is.matrix(themap[[1]]))
thepos <- unlist(lapply(themap, function(a) a[1,]))
else thepos <- unlist(themap)
temp <- results[,1:(ncol(results)-1),drop=FALSE]
res <- data.frame(chr=rep(names(cross$geno),nmar(cross)),
pos=thepos,
neglog10P=-log10(results[,ncol(results)]),
missing=temp[,1]/apply(temp, 1, sum),
temp[,-1]/apply(temp[,-1], 1, sum), stringsAsFactors=TRUE)
class(res) <- c("scanone", "data.frame")
rownames(res) <- rownames(results)
res[,1] <- factor(as.character(res[,1]), levels=unique(as.character(res[,1])))
res
}
genotab.em <-
function(dat, init.geno, tol=1e-6, maxit=10000, verbose=FALSE)
{
genotab.ll <-
function(dat, gam, init.geno)
{
p <- c(init.geno[1]*(1-gam[2]), init.geno[2]*gam[1], init.geno[3]*(1-gam[3]),
gam[2]*init.geno[4], gam[3]*init.geno[5])
if(any(p==0 & dat > 0)) return(-Inf)
return( sum((dat*log(p))[dat>0 & p>0]) )
}
n <- sum(dat)
gam <- c(sum(dat[1:3]), dat[4], dat[5])/n
curll <- genotab.ll(dat, gam, init.geno)
flag <- 0
if(verbose) cat(0, gam, curll, "\n")
for(i in 1:maxit) {
# estep
zAA <- dat[1]*gam[3]/(gam[1]+gam[3])
zBB <- dat[3]*gam[2]/(gam[1]+gam[2])
# mstep
gamnew <- c(sum(dat[1:3])-zAA-zBB, dat[4]+zBB, dat[5]+zAA)/n
newll <- genotab.ll(dat, gamnew, init.geno)
if(verbose) cat(i, gamnew, newll, "\n")
if(abs(curll-newll) < tol) {
flag <- 1
break
}
gam <- gamnew
curll <- newll
}
if(!flag) warning("Didn't converge.")
gam <- gamnew
p <- c(init.geno[1]*(1-gam[2]), init.geno[2]*gam[1], init.geno[3]*(1-gam[3]),
gam[2]*init.geno[4], gam[3]*init.geno[5])
ex <- p*n
1-pchisq(sum(((dat-ex)^2/ex)[ex>0]), 2)
}
######################################################################
# geno.crosstab
#
# Get a cross-tabulation of the genotypes at two markers,
# with the markers specified by name
######################################################################
geno.crosstab <-
function(cross, mname1, mname2, eliminate.zeros=TRUE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
if(missing(mname2) && length(mname1)>1) {
mname2 <- mname1[2]
mname1 <- mname1[1]
}
if(length(mname1) > 1 || length(mname2) > 1)
stop("mname1 and mname2 should both have lenght 1, or mname1 should have length 2 and mname1 should be missing.")
if(mname1==mname2)
stop("You must give two distinct marker names.")
pos <- find.markerpos(cross, c(mname1, mname2))
if(any(is.na(pos$chr))) {
if(all(is.na(pos$chr)))
stop("Markers ", mname1, " and ", mname2, " not found.")
else
stop("Marker ", rownames(pos)[is.na(pos$chr)], " not found.")
}
chrtype <- sapply(cross$geno[pos$chr], class)
crosstype <- class(cross)[1]
g1 <- pull.geno(cross, pos$chr[1])[,mname1, drop=FALSE]
g2 <- pull.geno(cross, pos$chr[2])[,mname2, drop=FALSE]
if(chrtype[1] == "X")
g1 <- reviseXdata(crosstype, "full", getsex(cross), geno=g1, cross.attr=attributes(cross))
if(chrtype[2] == "X")
g2 <- reviseXdata(crosstype, "full", getsex(cross), geno=g2, cross.attr=attributes(cross))
g1[is.na(g1)] <- 0
g2[is.na(g2)] <- 0
g1names <- getgenonames(crosstype, chrtype[1], "full", getsex(cross), attributes(cross))
g2names <- getgenonames(crosstype, chrtype[2], "full", getsex(cross), attributes(cross))
if(chrtype[1] != "X") {
if(crosstype == "f2")
g1names <- c(g1names, paste("not", g1names[c(3,1)]))
else if(crosstype == "bc" || crosstype == "risib" || crosstype=="riself" || crosstype=="dh") {
}
else if(crosstype == "4way") {
temp <- g1names
g1names <- c(temp,
paste(temp[c(1,3)], collapse="/"),
paste(temp[c(2,4)], collapse="/"),
paste(temp[c(1,2)], collapse="/"),
paste(temp[c(3,4)], collapse="/"),
paste(temp[c(1,4)], collapse="/"),
paste(temp[c(2,3)], collapse="/"),
paste("not", temp[1]),
paste("not", temp[2]),
paste("not", temp[3]),
paste("not", temp[4]))
g1names[5:8] <- substr(temp[c(1,2,1,3)], c(1,1,2,2), c(1,1,2,2))
}
else stop("Unknown cross type: ",crosstype)
}
if(chrtype[2] != "X") {
if(crosstype == "f2")
g2names <- c(g2names, paste("not", g2names[c(3,1)]))
else if(crosstype == "bc" || crosstype == "risib" || crosstype=="riself" || crosstype=="dh") {
}
else if(crosstype == "4way") {
temp <- g2names
g2names <- c(temp,
paste(temp[c(1,3)], collapse="/"),
paste(temp[c(2,4)], collapse="/"),
paste(temp[c(1,2)], collapse="/"),
paste(temp[c(3,4)], collapse="/"),
paste(temp[c(1,4)], collapse="/"),
paste(temp[c(2,3)], collapse="/"),
paste("not", temp[1]),
paste("not", temp[2]),
paste("not", temp[3]),
paste("not", temp[4]))
g2names[5:8] <- substr(temp[c(1,2,1,3)], c(1,1,2,2), c(1,1,2,2))
}
else stop("Unknown cross type: ",crosstype)
}
g1names <- c("-", g1names)
g2names <- c("-", g2names)
g1 <- as.character(g1)
g2 <- as.character(g2)
for(i in 1:length(g1names)) {
j <- as.character(i-1)
g1[g1==j] <- g1names[i]
}
g1 <- factor(g1, levels=g1names)
for(i in 1:length(g2names)) {
j <- as.character(i-1)
g2[g2==j] <- g2names[i]
}
g2 <- factor(g2, levels=g2names)
tab <- table(g1, g2)
names(attributes(tab)$dimnames) <- c(mname1, mname2)
if(eliminate.zeros) { # eliminate rows and columns with no data (but always include missing data row and column)
rs <- apply(tab, 1, sum); rs[1] <- 1
tab <- tab[rs>0,,drop=FALSE]
cs <- apply(tab, 2, sum); cs[1] <- 1
tab <- tab[,cs>0,drop=FALSE]
}
tab
}
# map functions
mf.k <- function(d) 0.5*tanh(d/50)
mf.h <- function(d) 0.5*(1-exp(-d/50))
imf.k <- function(r) { r[r >= 0.5] <- 0.5-1e-14; 50*atanh(2*r) }
imf.h <- function(r) { r[r >= 0.5] <- 0.5-1e-14; -50*log(1-2*r) }
mf.m <- function(d) sapply(d,function(a) min(a/100,0.5))
imf.m <- function(r) sapply(r,function(a) min(a*100,50))
# carter-falconer: mf.cf, imf.cf
imf.cf <- function(r) { r[r >= 0.5] <- 0.5-1e-14; 12.5*(log(1+2*r)-log(1-2*r))+25*atan(2*r) }
mf.cf <-
function(d)
{
d[d >= 300] <- 300
icf <- function(r,d)
imf.cf(r)-d
sapply(d,function(a) {
if(a==0) return(0)
uniroot(icf, c(0,0.5-1e-14),d=a,tol=1e-12)$root })
}
######################################################################
#
# switch.order: change the marker order on a given chromosome to some
# specified order
#
######################################################################
switch.order <-
function(cross, chr, order, error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"),
maxit=4000, tol=1e-6, sex.sp=TRUE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
map.function <- match.arg(map.function)
# check chromosome argument
if(missing(chr)) {
chr <- names(cross$geno)[1]
warning("Assuming you mean chromosome ", chr)
}
else {
if(length(chr) > 1) {
chr <- chr[1]
warning("switch.order can deal with just one chromosome at a time; assuming you want chr ", chr)
}
if(!testchr(chr, names(cross$geno)))
stop("Chr ", chr, " not found.")
}
chr <- matchchr(chr, names(cross$geno))
# check order argument
n.mar <- nmar(cross)
if(n.mar[chr] == length(order)-2 || n.mar[chr]==length(order)-1)
order <- order[1:n.mar[chr]] # useful for output from ripple()
if(n.mar[chr] != length(order))
stop("Incorrect number of markers.")
# save recombination fractions
flag <- 0
if("rf" %in% names(cross)) {
attrib <- attributes(cross$rf)
rf <- cross$rf
# determine column within rec fracs
whchr <- which(names(cross$geno)==chr)
oldcols <- cumsum(c(0,n.mar))[whchr]+seq(along=order)
newcols <- cumsum(c(0,n.mar))[whchr]+order
rf[oldcols,] <- rf[newcols,]
rf[,oldcols] <- rf[,newcols]
colnames(rf)[oldcols] <- colnames(rf)[newcols]
rownames(rf)[oldcols] <- rownames(rf)[newcols]
flag <- 1
}
# remove any intermediate calculations (except rec fracs),
# as they will no longer be meaningful
cross <- clean(cross)
if(!is.matrix(cross$geno[[chr]]$map))
first <- min(cross$geno[[chr]]$map)
else
first <- apply(cross$geno[[chr]]$map,1,min)
# re-order markers
cross$geno[[chr]]$data <- cross$geno[[chr]]$data[,order,drop=FALSE]
m <- seq(0,by=5,length=ncol(cross$geno[[chr]]$data))
names(m) <- colnames(cross$geno[[chr]]$data)
if(is.matrix(cross$geno[[chr]]$map))
cross$geno[[chr]]$map <- rbind(m,m)
else
cross$geno[[chr]]$map <- m
# re-estimate rec fracs for re-ordered chromosome
if(flag==1) {
temp <- est.rf(subset(cross, chr=chr))$rf
rf[oldcols,oldcols] <- temp
cross$rf <- rf
if("onlylod" %in% names(attrib))
attr(cross$rf, "onlylod") <- attrib$onlylod
}
# re-estimate map
newmap <- est.map(cross, chr=chr,
error.prob=error.prob, map.function=map.function,
maxit=maxit, tol=tol, sex.sp=sex.sp)
if(!is.matrix(newmap[[1]]))
cross$geno[[chr]]$map <- newmap[[1]] + first
else {
cross$geno[[chr]]$map[1,] <- newmap[[1]][1,] + first[1]
cross$geno[[chr]]$map[2,] <- newmap[[1]][2,] + first[2]
rownames(cross$geno[[chr]]$map) <- NULL
}
cross
}
######################################################################
#
# subset.cross: General subsetting function for a cross object
#
######################################################################
subset.cross <-
function(x, chr, ind, ...)
{
if(!any(class(x) == "cross"))
stop("Input should have class \"cross\".")
if(missing(chr) && missing(ind))
stop("You must specify either chr or ind.")
n.chr <- nchr(x)
n.ind <- nind(x)
# pull out relevant chromosomes
if(!missing(chr)) {
chr <- matchchr(chr, names(x$geno))
if("rf" %in% names(x)) { # pull out part of rec fracs
if(totmar(x) != ncol(x$rf))
x <- clean(x)
else {
attrib <- attributes(x$rf)
n.mar <- nmar(x)
n.chr <- n.chr
wh <- rbind(c(0,cumsum(n.mar)[-n.chr])+1,cumsum(n.mar))
dimnames(wh) <- list(NULL, names(n.mar))
wh <- wh[,chr,drop=FALSE]
wh <- unlist(apply(wh,2,function(a) a[1]:a[2]))
x$rf <- x$rf[wh,wh]
if("onlylod" %in% names(attrib)) # save the onlylod attribute if its there
attr(x$rf, "onlylod") <- attrib$onlylod
}
}
x$geno <- x$geno[chr]
if("founderGeno" %in% names(x))
x$founderGeno <- x$founderGeno[,unlist(lapply(x$geno, function(a) colnames(a$data)))]
}
if(!missing(ind)) {
theid <- getid(x)
if(is.logical(ind)) {
ind[is.na(ind)] <- FALSE
if(length(ind) != n.ind)
stop("ind argument has wrong length (", length(ind), "; should be ", n.ind, ")")
ind <- (1:n.ind)[ind]
}
else if(is.numeric(ind)) { # treat as numeric indices; don't match against individual identifiers
if(all(ind < 0)) { # drop all but these
ind <- -ind
if(any(ind > nind(x))) {
a <- -ind[ind > nind(x)]
if(length(a) > 1) a <- sample(a, 1)
stop("Invalid ind values (e.g., ", a, ")")
}
ind <- (1:n.ind)[-ind]
}
else if(all(ind > 0)) { # keep these
if(any(ind > nind(x))) {
a <- ind[ind > nind(x)]
if(length(a) > 1) a <- sample(a, 1)
stop("Invalid ind values (e.g., ", a, ")")
}
}
else {
stop("Need ind to be all > 0 or all < 0.")
}
}
else if(!is.null(theid)) { # cross has individual IDs
ind <- as.character(ind)
if(all(substr(ind, 1,1) == "-")) {
ind <- substr(ind, 2, nchar(ind))
m <- match(ind, theid)
if(all(is.na(m)))
stop("No matching individuals.")
if(any(is.na(m)))
warning("Individuals not found: ", paste(ind[is.na(m)]))
ind <- (1:n.ind)[-m[!is.na(m)]]
}
else {
m <- match(ind, theid)
if(any(is.na(m)))
warning("Individuals not found: ", paste(ind[is.na(m)], collapse=" "))
ind <- m[!is.na(m)]
}
ind <- ind[!is.na(ind)]
}
else { # no individual IDs
stop("In the absense of individual IDs, ind should be logical or numeric.")
}
# Note: ind should now be a numeric vector
if(length(ind) == 0)
stop("Must retain at least one individual.")
if(length(ind) == 1)
warning("Retained only one individual!")
x$pheno <- x$pheno[ind,,drop=FALSE]
if("cross" %in% names(x))
x$cross <- x$cross[ind,,drop=FALSE]
for(i in 1:nchr(x)) {
x$geno[[i]]$data <- x$geno[[i]]$data[ind,,drop=FALSE]
if("prob" %in% names(x$geno[[i]])) {
temp <- attributes(x$geno[[i]]$prob) # all attributes but dim and dimnames
x$geno[[i]]$prob <- x$geno[[i]]$prob[ind,,,drop=FALSE]
# put attributes back in
for(k in seq(along=temp)) {
if(names(temp)[k] != "dim" && names(temp)[k] != "dimnames")
attr(x$geno[[i]]$prob, names(temp)[k]) <- temp[[k]]
}
}
if("errorlod" %in% names(x$geno[[i]])) {
temp <- attributes(x$geno[[i]]$prob) # all attributes but dim and dimnames
x$geno[[i]]$errorlod <- x$geno[[i]]$errorlod[ind,,drop=FALSE]
# put attributes back in
for(k in seq(along=temp)) {
if(names(temp)[k] != "dim" && names(temp)[k] != "dimnames")
attr(x$geno[[i]]$errorlod, names(temp)[k]) <- temp[[k]]
}
}
if("argmax" %in% names(x$geno[[i]])) {
temp <- attributes(x$geno[[i]]$argmax) # all attributes but dim and dimnames
x$geno[[i]]$argmax <- x$geno[[i]]$argmax[ind,,drop=FALSE]
# put attributes back in
for(k in seq(along=temp)) {
if(names(temp)[k] != "dim" && names(temp)[k] != "dimnames")
attr(x$geno[[i]]$argmax, names(temp)[k]) <- temp[[k]]
}
}
if("draws" %in% names(x$geno[[i]])) {
temp <- attributes(x$geno[[i]]$draws) # all attributes but dim and dimnames
x$geno[[i]]$draws <- x$geno[[i]]$draws[ind,,,drop=FALSE]
# put attributes back in
for(k in seq(along=temp)) {
if(names(temp)[k] != "dim" && names(temp)[k] != "dimnames")
attr(x$geno[[i]]$draws, names(temp)[k]) <- temp[[k]]
}
}
}
if("qtlgeno" %in% names(x))
x$qtlgeno <- x$qtlgeno[ind,,drop=FALSE]
}
x
}
#pull.chr <-
#function(cross, chr) {
# warning("pull.chr is deprecated; use subset.cross.")
# subset.cross(cross, chr)
#}
######################################################################
#
# c.cross: Combine crosses
#
######################################################################
c.cross <-
function(...)
{
args <- list(...)
n.args <- length(args)
for(i in seq(along=args)) {
if(!any(class(args[[i]]) == "cross"))
stop("Input should have class \"cross\".")
}
# if only one cross, just return it
if(n.args==1) return(args[[1]])
if(any(sapply(args, function(a) !any(class(a) == "cross"))))
stop("All arguments must be cross objects.")
# crosses must be all the same, or must be combination of F2 and BC
classes <- sapply(args,function(a) class(a)[1])
if(length(unique(classes))==1) {
allsame <- TRUE
type <- classes[1]
}
else {
if(any(classes != "bc" & classes != "f2"))
stop("Experiments must be either the same type or be bc/f2.")
allsame <- FALSE
type <- "f2"
}
if(length(unique(sapply(args, nchr))) > 1)
stop("All arguments must have the same number of chromosomes.")
x <- args[[1]]
chr <- names(x$geno)
n.mar <- nmar(x)
marnam <- unlist(lapply(x$geno,function(b) colnames(b$data)))
marpos <- unlist(lapply(x$geno,function(b) b$map))
map.mismatch <- 0
for(i in 2:n.args) {
y <- args[[i]]
y.marnam <- unlist(lapply(y$geno, function(b) colnames(b$data)))
y.marpos <- unlist(lapply(y$geno, function(b) b$map))
if(chr != names(y$geno) || any(n.mar != nmar(y)) ||
any(marnam != y.marnam) || any(marpos != y.marpos)) {
map.mismatch <- 1
break
}
}
if(map.mismatch) { # get the maps to line up
for(i in 1:nchr(args[[1]])) {
themap <- NULL
themaps <- vector("list", n.args)
for(j in 1:n.args) {
if(is.matrix(args[[j]]$map))
stop("c.cross() won't work with sex-specific maps.")
themaps[[j]] <- args[[j]]$geno[[i]]$map
themap <- c(themap, themaps[[j]])
}
themap <- sort(themap)
mn <- unique(names(themap))
newmap <- rep(0,length(mn))
names(newmap) <- mn
for(j in 1:length(newmap))
newmap[j] <- mean(themap[names(themap) == mn[j]])
for(j in 1:n.args) {
m <- match(names(themaps[[j]]), mn)
m <- m[!is.na(m)]
if(any(diff(m)) < 0)
stop(" Markers must all be in the same order.")
if(!all(mn %in% names(themaps[[j]]))) {
temp <- matrix(ncol=length(mn), nrow=nind(args[[j]]))
colnames(temp) <- mn
temp[,names(themaps[[j]])] <- args[[j]]$geno[[i]]$data
args[[j]]$geno[[i]]$data <- temp
}
args[[j]]$geno[[i]]$map <- newmap
}
}
} # end of map mismatch fix
x <- args[[1]]
chr <- names(x$geno)
n.mar <- nmar(x)
marnam <- unlist(lapply(x$geno,function(b) colnames(b$data)))
# create genotype information
geno <- x$geno
for(j in 1:nchr(x)) { # drop extraneous stuff
geno[[j]] <- list(data=geno[[j]]$data, map=geno[[j]]$map)
class(geno[[j]]) <- class(x$geno[[j]])
}
for(i in 2:n.args)
for(j in 1:nchr(x))
geno[[j]]$data <- rbind(geno[[j]]$data,args[[i]]$geno[[j]]$data)
# get all phenotype names
phenam <- names(x$pheno)
for(i in 2:n.args)
phenam <- c(phenam, names(args[[i]]$pheno))
phenam <- unique(phenam)
# form big phenotype matrix
n.ind <- sapply(args,nind)
pheno <- matrix(nrow=sum(n.ind),ncol=length(phenam))
colnames(pheno) <- phenam
pheno <- as.data.frame(pheno, stringsAsFactors=TRUE)
if(!allsame) {
crosstype <- factor(rep(c("bc","f2")[match(classes,c("bc","f2"))],n.ind),
levels=c("bc","f2"))
pheno <- cbind(pheno,crosstype=crosstype)
}
for(i in 1:length(phenam)) {
phe <- vector("list",n.args)
for(j in 1:n.args) {
o <- match(phenam[i],names(args[[j]]$pheno))
if(is.na(o)) phe[[j]] <- rep(NA,n.ind[j])
else phe[[j]] <- args[[j]]$pheno[,o]
}
pheno[,i] <- unlist(phe)
}
# indicator of which cross
whichcross <- matrix(0,ncol=n.args,nrow=sum(n.ind))
colnames(whichcross) <- paste("cross",1:n.args,sep="")
thecross <- rep(NA, sum(n.ind))
prev <- 0
for(i in 1:n.args) {
wh <- prev + 1:n.ind[i]
prev <- prev + n.ind[i]
whichcross[wh,i] <- 1
thecross[wh] <- i
}
pheno <- cbind(pheno,cross=thecross,whichcross)
x$pheno <- pheno
if(!map.mismatch) { # keep probs and draws only if we've not re-aligned the maps
# if probs exist in each and all have the same
# set up values, keep them
wh <- sapply(args, function(a) match("prob",names(a$geno[[1]])))
step <- sapply(args,function(a) attr(a$geno[[1]]$prob,"step"))
error.prob <- sapply(args,function(a) attr(a$geno[[1]]$prob,"error.prob"))
off.end <- sapply(args,function(a) attr(a$geno[[1]]$prob,"off.end"))
map.function <- sapply(args,function(a) attr(a$geno[[1]]$prob,"map.function"))
map <- lapply(args,function(a) attr(a$geno[[1]]$prob,"map"))
if(!any(is.na(wh)) && length(unique(step))==1 &&
length(unique(error.prob))==1 && length(unique(off.end))==1 &&
length(unique(map.function))==1) {
if(allsame) { # all same cross type
for(j in 1:nchr(x)) {
geno[[j]]$prob <- array(dim=c(sum(n.ind),dim(x$geno[[j]]$prob)[-1]))
dimnames(geno[[j]]$prob) <- dimnames(x$geno[[j]]$prob)
prev <- 0
for(i in 1:n.args) {
wh <- prev + 1:n.ind[i]
prev <- prev + n.ind[i]
geno[[j]]$prob[wh,,] <- args[[i]]$geno[[j]]$prob
}
}
}
else { # mixed F2 and BC
for(j in 1:nchr(x)) {
wh <- match("f2",classes)
geno[[j]]$prob <- array(0,dim=c(sum(n.ind),dim(args[[wh]]$geno[[j]]$prob)[-1]))
dimnames(geno[[j]]$prob) <- dimnames(args[[wh]]$geno[[j]]$prob)
prev <- 0
for(i in 1:n.args) {
wh <- prev + 1:n.ind[i]
prev <- prev + n.ind[i]
if(classes[i]=="f2")
geno[[j]]$prob[wh,,] <- args[[i]]$geno[[j]]$prob
else # backcross
geno[[j]]$prob[wh,,1:2] <- args[[i]]$geno[[j]]$prob
}
}
}
for(j in 1:nchr(x)) {
wh <- sapply(args, function(a) match("prob",names(a$geno[[j]])))
step <- sapply(args,function(a) attr(a$geno[[j]]$prob,"step"))
error.prob <- sapply(args,function(a) attr(a$geno[[j]]$prob,"error.prob"))
off.end <- sapply(args,function(a) attr(a$geno[[j]]$prob,"off.end"))
map.function <- sapply(args,function(a) attr(a$geno[[j]]$prob,"map.function"))
map <- lapply(args,function(a) attr(a$geno[[j]]$prob,"map"))
attr(geno[[j]]$prob,"step") <- step[1]
attr(geno[[j]]$prob,"error.prob") <- error.prob[1]
attr(geno[[j]]$prob,"off.end") <- off.end[1]
attr(geno[[j]]$prob,"map.function") <- map.function[1]
attr(geno[[j]]$prob,"map") <- map[[1]]
}
}
# if draws exist in each and all have the same
# set up values, keep them
wh <- sapply(args, function(a) match("draws",names(a$geno[[1]])))
step <- sapply(args,function(a) attr(a$geno[[1]]$draws,"step"))
error.prob <- sapply(args,function(a) attr(a$geno[[1]]$draws,"error.prob"))
off.end <- sapply(args,function(a) attr(a$geno[[1]]$draws,"off.end"))
map.function <- sapply(args,function(a) attr(a$geno[[1]]$draws,"map.function"))
map <- lapply(args,function(a) attr(a$geno[[1]]$draws,"map"))
ndraws <- sapply(args,function(a) dim(a$geno[[1]]$draws)[3])
if(!any(is.na(wh)) && length(unique(step))==1 &&
length(unique(error.prob))==1 && length(unique(off.end))==1 &&
length(unique(map.function))==1 && length(unique(ndraws))==1) {
for(j in 1:nchr(x)) {
geno[[j]]$draws <- array(0,dim=c(sum(n.ind),dim(x$geno[[j]]$draws)[-1]))
dimnames(geno[[j]]$draws) <- dimnames(x$geno[[j]]$draws)
prev <- 0
for(i in 1:n.args) {
wh <- prev + 1:n.ind[i]
prev <- prev + n.ind[i]
geno[[j]]$draws[wh,,] <- args[[i]]$geno[[j]]$draws
}
wh <- sapply(args, function(a) match("draws",names(a$geno[[j]])))
step <- sapply(args,function(a) attr(a$geno[[j]]$draws,"step"))
error.prob <- sapply(args,function(a) attr(a$geno[[j]]$draws,"error.prob"))
off.end <- sapply(args,function(a) attr(a$geno[[j]]$draws,"off.end"))
map.function <- sapply(args,function(a) attr(a$geno[[j]]$draws,"map.function"))
map <- lapply(args,function(a) attr(a$geno[[j]]$draws,"map"))
attr(geno[[j]]$draws,"step") <- step[1]
attr(geno[[j]]$draws,"error.prob") <- error.prob[1]
attr(geno[[j]]$draws,"off.end") <- off.end[1]
attr(geno[[j]]$draws,"map.function") <- map.function[1]
attr(geno[[j]]$draws,"map") <- map[[1]]
}
}
}
x <- list(geno=geno, pheno=pheno)
class(x) <- c(type,"cross")
x
}
######################################################################
#
# fill.geno: Run argmax.geno or sim.geno and then fill in the
# genotype data with the results. This will allow
# rough genome scans by marker regression without
# holes. WE WOULD NOT PLACE ANY TRUST IN THE RESULTS!
#
# the newer method, "no_dbl_XO", fills in missing genotypes between
# markers with matching genotypes
#
######################################################################
fill.geno <-
function(cross, method=c("imp","argmax", "no_dbl_XO"), error.prob=0.0001,
map.function=c("haldane","kosambi","c-f","morgan"))
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
method <- match.arg(method)
# don't let error.prob be exactly zero (or >1)
if(error.prob < 1e-50) error.prob <- 1e-50
if(error.prob > 1) {
error.prob <- 1-1e-50
warning("error.prob shouldn't be > 1!")
}
# remove any extraneous material
cross <- clean(cross)
n.chr <- nchr(cross)
n.mar <- nmar(cross)
if(method=="imp") {
# do one imputation
temp <- sim.geno(cross,n.draws=1,step=0,off.end=0,
error.prob=error.prob,map.function=map.function)
# replace the genotype data with the results,
# stripping off any attributes
for(i in 1:n.chr) {
nam <- colnames(cross$geno[[i]]$data)
if(n.mar[i] == 1)
cross$geno[[i]]$data <-
matrix(as.numeric(temp$geno[[i]]$draws[,2,1]),ncol=n.mar[i])
else
cross$geno[[i]]$data <-
matrix(as.numeric(temp$geno[[i]]$draws[,,1]),ncol=n.mar[i])
colnames(cross$geno[[i]]$data) <- nam
}
}
else if(method=="argmax") {
# run the Viterbi algorithm
temp <- argmax.geno(cross,step=0,off.end=0,error.prob=error.prob,
map.function=map.function)
# replace the genotype data with the results,
# stripping off any attributes
for(i in 1:n.chr) {
nam <- colnames(cross$geno[[i]]$data)
if(n.mar[i] == 1)
cross$geno[[i]]$data <-
matrix(as.numeric(temp$geno[[i]]$argmax[,2]),ncol=n.mar[i])
else
cross$geno[[i]]$data <-
matrix(as.numeric(temp$geno[[i]]$argmax),ncol=n.mar[i])
colnames(cross$geno[[i]]$data) <- nam
}
}
else {
for(i in 1:n.chr) {
dat <- cross$geno[[i]]$data
dat[is.na(dat)] <- 0
nr <- nrow(dat)
nc <- ncol(dat)
dn <- dimnames(dat)
dat <- .C("R_fill_geno_nodblXO",
as.integer(nr),
as.integer(nc),
dat=as.integer(dat),
PACKAGE="qtl")$dat
dat[dat==0] <- NA
dat <- matrix(dat, ncol=nc, nrow=nr)
dimnames(dat) <- dn
cross$geno[[i]]$data <- dat
}
}
cross
}
######################################################################
#
# checkcovar
#
# This is a utility function for scanone and scantwo. We remove
# individuals with missing phenotypes or covariates and check
# that the covariates are of the appropriate form.
#
######################################################################
checkcovar <-
function(cross, pheno.col, addcovar, intcovar, perm.strata, ind.noqtl=NULL, weights=NULL, verbose=TRUE)
{
chrtype <- sapply(cross$geno, class)
# drop individuals whose sex or pgm is missing if X chr is included
if(any(chrtype=="X")) {
sexpgm <- getsex(cross)
keep <- rep(TRUE,nind(cross))
flag <- 0
if(!is.null(sexpgm$sex)) {
if(any(is.na(sexpgm$sex))) {
keep[is.na(sexpgm$sex)] <- FALSE
flag <- 1
}
}
if(!is.null(sexpgm$pgm)) {
if(any(is.na(sexpgm$pgm))) {
keep[is.na(sexpgm$pgm)] <- FALSE
flag <- 1
}
}
if(flag) {
if(verbose) warning("Dropping ", sum(!keep), " individuals with missing sex or pgm.\n")
cross <- subset(cross, ind=keep)
if(!is.null(addcovar)) {
if(!is.matrix(addcovar)) addcovar <- addcovar[keep]
else addcovar <- addcovar[keep,]
}
if(!is.null(intcovar)) {
if(!is.matrix(intcovar)) intcovar <- intcovar[keep]
else intcovar <- intcovar[keep,]
}
}
}
# check phenotypes - we allow multiple phenotypes here
if(pheno.col < 1 || pheno.col > nphe(cross))
stop("Specified phenotype column is invalid.")
# check if all phenotypes are numeric
pheno <- cross$pheno[,pheno.col,drop=FALSE]
idx.nonnum <- which(!apply(pheno,2, is.numeric))
if(length(idx.nonnum) > 0)
stop("Following phenotypes are not numeric: Column ",
paste(idx.nonnum, collapse=","))
orig.n.ind <- nind(cross)
# drop individuals with missing phenotypes
if(any(!is.finite(unlist(pheno)))) {
keep.ind <- as.numeric(which(apply(pheno, 1, function(x) !any(!is.finite(x)))))
# keep.ind <- (1:length(pheno))[!is.na(pheno)]
n.drop <- nind(cross) - length(keep.ind)
keep.ind.boolean <- rep(FALSE, nind(cross))
keep.ind.boolean[keep.ind] <- TRUE
cross <- subset.cross(cross,ind=keep.ind.boolean)
pheno <- pheno[keep.ind,,drop=FALSE]
if(verbose) warning("Dropping ", n.drop, " individuals with missing phenotypes.\n")
}
else keep.ind <- 1:nind(cross)
n.ind <- nind(cross)
n.chr <- nchr(cross) # number of chromosomes
type <- class(cross)[1] # type of cross
n.addcovar <- n.intcovar <- 0
if(!is.null(addcovar)) { # for additive covariates
if(!is.matrix(addcovar)) {
if(!is.numeric(as.matrix(addcovar)))
stop("addcovar should be numeric")
if(is.vector(addcovar) || is.data.frame(addcovar))
addcovar <- as.matrix(addcovar)
else stop("addcovar should be a matrix")
}
if(!all(apply(addcovar,2,is.numeric)))
stop("All columns of addcovar must be numeric")
if( nrow(addcovar) != orig.n.ind ) {
# the length of additive covariates is incorrect
stop("Number of rows in additive covariates is incorrect")
}
addcovar <- addcovar[keep.ind,,drop=FALSE]
n.addcovar <- ncol(addcovar)
}
if(!is.null(intcovar)) { # interacting covariates
if(!is.matrix(intcovar)) {
if(!is.numeric(as.matrix(intcovar)))
stop("intcovar should be a numeric")
if(is.vector(intcovar) || is.data.frame(intcovar))
intcovar <- as.matrix(intcovar)
else stop("intcovar should be a matrix")
}
if(!all(apply(intcovar,2,is.numeric)))
stop("All columns of intcovar must be numeric")
if(nrow(intcovar)[1] != orig.n.ind) {
# the length of interacting covariates is incorrect
stop("The length of interacting covariates is incorrect!")
}
intcovar <- intcovar[keep.ind,,drop=FALSE]
n.intcovar <- ncol(intcovar)
}
if(!is.null(perm.strata)) perm.strata <- perm.strata[keep.ind]
if(!is.null(ind.noqtl)) ind.noqtl <- ind.noqtl[keep.ind]
if(!is.null(weights)) weights <- weights[keep.ind]
# drop individuals missing any covariates
if(!is.null(addcovar)) { # note that intcovar is contained in addcovar
wh <- apply(cbind(addcovar,intcovar),1,function(a) any(!is.finite(a)))
if(any(wh)) {
cross <- subset.cross(cross,ind=(!wh))
pheno <- pheno[!wh,,drop=FALSE]
addcovar <- addcovar[!wh,,drop=FALSE]
if(!is.null(intcovar)) intcovar <- intcovar[!wh,,drop=FALSE]
n.ind <- nind(cross)
if(!is.null(perm.strata)) perm.strata <- perm.strata[!wh]
if(!is.null(ind.noqtl)) ind.noqtl <- ind.noqtl[!wh]
if(!is.null(weights)) weights <- weights[!wh]
if(verbose) warning("Dropping ", sum(wh), " individuals with missing covariates.\n")
}
}
# make sure columns of intcovar are contained in addcovar
if(!is.null(intcovar)) {
if(is.null(addcovar)) {
addcovar <- intcovar
n.addcovar <- n.intcovar
if(verbose) warning("addcovar forced to contain all columns of intcovar\n")
}
else {
wh <- 1:n.intcovar
for(i in 1:n.intcovar) {
o <- (apply(addcovar,2,function(a,b) max(abs(a-b)),intcovar[,i])<1e-14)
if(any(o)) wh[i] <- (1:n.addcovar)[o]
else wh[i] <- NA
}
if(any(!is.finite(wh))) {
addcovar <- cbind(addcovar,intcovar[,!is.finite(wh)])
n.addcovar <- ncol(addcovar)
if(verbose) warning("addcovar forced to contain all columns of intcovar")
}
}
}
if(n.addcovar > 0) { # check rank
if(qr(cbind(1,addcovar))$rank < n.addcovar+1)
if(verbose) warning("addcovar appears to be over-specified; consider dropping columns.\n")
}
if(n.intcovar > 0) { # check rank
if(qr(cbind(1,intcovar))$rank < n.intcovar+1)
if(verbose) warning("intcovar appears to be over-specified; consider dropping columns.\n")
}
pheno <- as.matrix(pheno)
list(cross=cross, pheno=pheno, addcovar=addcovar, intcovar=intcovar,
n.addcovar=n.addcovar, n.intcovar=n.intcovar, perm.strata=perm.strata,
ind.noqtl=ind.noqtl, weights=weights)
}
# Find the nearest marker to a particular position
find.marker <-
function(cross, chr, pos, index)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
if(missing(pos) && missing(index))
stop("Give either pos or index.")
if(!missing(pos) && !missing(index))
stop("Give just one of pos or index.")
# if chr has length 1, expand if necessary
if(length(chr) == 1) {
if(!missing(pos))
chr <- rep(chr,length(pos))
else
chr <- rep(chr, length(index))
}
# otherwise, chr and pos should have same length
else {
if(!missing(pos) && length(chr) != length(pos))
stop("chr and pos must be the same length.")
if(!missing(index) && length(chr) != length(index))
stop("chr and index must be the same length.")
}
markers <- rep("",length(chr))
chrnotfound <- NULL
for(i in 1:length(chr)) {
# find chromosome
o <- match(chr[i], names(cross$geno))
if(is.na(o)) {
markers[i] <- NA # chr not matched
chrnotfound <- c(chrnotfound, chr[i])
}
else {
thismap <- cross$geno[[o]]$map # genetic map
# sex-specific map; look at female positions
if(is.matrix(thismap)) thismap <- thismap[1,]
if(!missing(pos)) {
# find closest marker
d <- abs(thismap-pos[i])
o2 <- (1:length(d))[d==min(d)]
if(length(o2)==1) markers[i] <- names(thismap)[o2]
# if multiple markers are equidistant,
# choose the one with the most data
# or choose among them at random
else {
x <- names(thismap)[o2]
n.geno <- apply(cross$geno[[o]]$data[,o2],2,function(a) sum(!is.na(a)))
o2 <- o2[n.geno==max(n.geno)]
if(length(o2) == 1)
markers[i] <- names(thismap)[o2]
else
markers[i] <- names(thismap)[sample(o2,1)]
}
}
else { # by index
if(index[i] < 1 || index[i] > length(thismap))
stop("Misspecified index ", index[i], " on chr ", chr[i])
markers[i] <- names(thismap)[index[i]]
}
}
}
if(length(chrnotfound) > 0) {
chrnotfound <- sort(unique(chrnotfound))
if(length(chrnotfound) == 1)
warning("Chromosome ", paste("\"", chrnotfound, "\"", sep=""), " not found")
else
warning("Chromosomes ", paste("\"", chrnotfound, "\"", sep="", collapse=", "), " not found")
}
markers
}
### Find the nearest pseudomarker to a particular position
find.pseudomarker <-
function(cross, chr, pos, where=c("draws","prob"), addchr=TRUE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
# if chr has length 1, expand if necessary
if(length(chr) == 1)
chr <- rep(chr,length(pos))
# otherwise, chr and pos should have same length
else if(length(chr) != length(pos))
stop("chr and pos must be the same length.")
markers <- rep("",length(chr))
where <- match.arg(where)
if(where=="draws" && !("draws" %in% names(cross$geno[[1]])))
stop("You'll need to first run sim.geno")
if(where=="prob" && !("prob" %in% names(cross$geno[[1]])))
stop("You'll need to first run calc.genoprob")
for(i in 1:length(chr)) {
# find chromosome
o <- match(chr[i], names(cross$geno))
if(is.na(o)) markers[i] <- NA # chr not matched
else {
if(where=="draws") {
if("map" %in% names(attributes(cross$geno[[o]]$draws)))
thismap <- attr(cross$geno[[o]]$draws, "map")
else {
stp <- attr(cross$geno[[o]]$draws, "step")
oe <- attr(cross$geno[[o]]$draws, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[o]]$draws)))
stpw <- attr(cross$geno[[o]]$draws, "stepwidth")
else stpw <- "fixed"
thismap <- create.map(cross$geno[[o]]$map,stp,oe,stpw)
attr(cross$geno[[o]]$draws, "map") <- thismap
}
}
else { # prob
if("map" %in% names(attributes(cross$geno[[o]]$prob)))
thismap <- attr(cross$geno[[o]]$prob, "map")
else {
stp <- attr(cross$geno[[o]]$prob, "step")
oe <- attr(cross$geno[[o]]$prob, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[o]]$prob)))
stpw <- attr(cross$geno[[o]]$prob, "stepwidth")
else stpw <- "fixed"
thismap <- create.map(cross$geno[[o]]$map,stp,oe,stpw)
attr(cross$geno[[o]]$prob, "map") <- thismap
}
}
# sex-specific map; look at female positions
if(is.matrix(thismap)) thismap <- thismap[1,]
# find closest marker
d <- abs(thismap-pos[i])
o2 <- (1:length(d))[d==min(d)]
if(length(o2)==1) themarker <- names(thismap)[o2]
else themarker <- names(thismap)[sample(o2, 1)]
if(addchr && length(grep("^loc[0-9]+\\.*[0-9]*(\\.[0-9]+)*$", themarker)) > 0)
themarker <- paste("c", chr[i], ".", themarker, sep="")
markers[i] <- themarker
}
}
markers
}
# expand recombination fractions for RI lines
adjust.rf.ri <-
function(r, type=c("self","sib"), chrtype=c("A","X"), expand=TRUE)
{
# type of RI lines
type <- match.arg(type)
chrtype <- match.arg(chrtype)
if(type=="self") {
if(expand) return(r*2/(1+2*r))
else return(r/2/(1-r))
}
else {
if(chrtype=="A") { # autosome / sib mating
if(expand) return(r*4/(1+6*r))
else return(r/(4-6*r))
}
else { # X chromosome/ sib mating
if(expand) return(8/3*r/(1+4*r))
else return(3/8*r/(1-1.5*r))
}
}
}
######################################################################
# lodint: function to get lod support interval
######################################################################
lodint <-
function(results, chr, qtl.index, drop=1.5, lodcolumn=1,
expandtomarkers=FALSE)
{
if(!("scanone" %in% class(results))) {
if(!("qtl" %in% class(results)))
stop("Input must have class \"scanone\" or \"qtl\".")
else {
if(!("lodprofile" %in% names(attributes(results))))
stop("qtl object needs to be produced by refineqtl with keeplodprofile=TRUE.")
else { # qtl object
if(lodcolumn != 1) {
warning("lod column ignored if input is a qtl object.")
lodcolumn <- 1
}
results <- attr(results, "lodprofile")
if(missing(qtl.index)) {
if(length(results)==1)
results <- results[[1]]
else
stop("You must specify qtl.index.")
}
else {
if(length(qtl.index)>1) stop("qtl.index should have length 1")
if(qtl.index < 1 || qtl.index > length(results))
stop("qtl.index misspecified.")
results <- results[[qtl.index]]
}
chr <- results[1,1]
}
}
}
else {
if(lodcolumn < 1 || lodcolumn +2 > ncol(results))
stop("Argument lodcolumn misspecified.")
if(missing(chr)) {
if(length(unique(results[,1]))>1)
stop("Give a chromosome ID.")
}
else {
if(length(chr) > 1) stop("chr should have length 1")
if(is.na(match(chr, results[,1])))
stop("Chromosome misspecified.")
results <- results[results[,1]==chr,]
}
}
if(all(is.na(results[,lodcolumn+2]))) return(NULL)
maxlod <- max(results[,lodcolumn+2],na.rm=TRUE)
w <- which(!is.na(results[,lodcolumn+2]) & results[,lodcolumn+2] == maxlod)
o <- range(which(!is.na(results[,lodcolumn+2]) & results[,lodcolumn+2] > maxlod-drop))
if(length(o)==0) o <- c(1,nrow(results))
else {
if(o[1] > 1) o[1] <- o[1]-1
if(o[2] < nrow(results)) o[2] <- o[2]+1
}
if(expandtomarkers) {
markerpos <- (1:nrow(results))[-grep("^c.+\\.loc-*[0-9]+(\\.[0-9]+)*$", rownames(results))]
if(any(markerpos <= o[1]))
o[1] <- max(markerpos[markerpos <= o[1]])
if(any(markerpos >= o[2]))
o[2] <- min(markerpos[markerpos >= o[2]])
}
rn <- rownames(results)[c(o[1],w,o[2])]
# look for duplicate rows
if(length(w)>1 && rn[length(rn)]==rn[length(rn)-1]) w <- w[-length(w)]
else if(length(w)>1 && rn[2]==rn[1]) w <- w[-1]
rn <- rownames(results)[c(o[1],w,o[2])]
# look for more duplicate rows
if(any(table(rn)> 1)) {
tab <- table(rn)
temp <- which(tab>1)
for(j in temp) {
z <- which(rn==names(tab)[j])
for(k in 2:length(z))
rn[z[k:length(z)]] <- paste(rn[z[k:length(z)]], " ", sep="")
}
}
results <- results[c(o[1],w,o[2]),]
rownames(results) <- rn
class(results) <- c("scanone","data.frame")
results
}
######################################################################
# bayesint: function to get Bayesian probability interval
######################################################################
bayesint <-
function(results, chr, qtl.index, prob=0.95, lodcolumn=1, expandtomarkers=FALSE)
{
if(!("scanone" %in% class(results))) {
if(!("qtl" %in% class(results)))
stop("Input must have class \"scanone\" or \"qtl\".")
else {
if(!("lodprofile" %in% names(attributes(results))))
stop("qtl object needs to be produced by refineqtl with keeplodprofile=TRUE.")
else { # qtl object
if(lodcolumn != 1) {
warning("lod column ignored if input is a qtl object.")
lodcolumn <- 1
}
results <- attr(results, "lodprofile")
if(missing(qtl.index)) {
if(length(results)==1)
results <- results[[1]]
else
stop("You must specify qtl.index.")
}
else {
if(length(qtl.index)>1) stop("qtl.index should have length 1")
if(qtl.index < 1 || qtl.index > length(results))
stop("qtl.index misspecified.")
results <- results[[qtl.index]]
}
chr <- results[1,1]
}
}
}
else {
if(lodcolumn < 1 || lodcolumn +2 > ncol(results))
stop("Argument lodcolumn misspecified.")
if(missing(chr)) {
if(length(unique(results[,1]))>1)
stop("Give a chromosome ID.")
}
else {
if(length(chr) > 1) stop("chr should have length 1")
if(is.na(match(chr, results[,1])))
stop("Chromosome misspecified.")
results <- results[results[,1]==chr,]
}
}
if(all(is.na(results[,lodcolumn+2]))) return(NULL)
loc <- results[,2]
width <- diff(( c(loc[1],loc) + c(loc, loc[length(loc)]) )/ 2)
area <- 10^results[,lodcolumn+2]*width
area <- area/sum(area)
o <- order(results[,lodcolumn+2], decreasing=TRUE)
cs <- cumsum(area[o])
wh <- min((1:length(loc))[cs >= prob])
int <- range(o[1:wh])
if(expandtomarkers) {
markerpos <- (1:nrow(results))[-grep("^c.+\\.loc-*[0-9]+(\\.[0-9]+)*$", rownames(results))]
if(any(markerpos <= int[1]))
int[1] <- max(markerpos[markerpos <= int[1]])
if(any(markerpos >= int[2]))
int[2] <- min(markerpos[markerpos >= int[2]])
}
rn <- rownames(results)[c(int[1],o[1],int[2])]
# look for duplicate rows
if(any(table(rn)> 1)) {
rn[2] <- paste(rn[2], "")
if(rn[1] == rn[3]) rn[3] <- paste(rn[3], " ")
}
results <- results[c(int[1],o[1],int[2]),]
rownames(results) <- rn
class(results) <- c("scanone", "data.frame")
results
}
######################################################################
# makeSSmap: convert a genetic map, or the genetic maps in a cross
# object, to be sex-specific (i.e., 2-row matrices)
######################################################################
makeSSmap <-
function(cross)
{
if(!any(class(cross) == "map")) {
# input object is a genetic map
for(i in 1:length(cross)) {
if(!is.matrix(cross[[i]]))
cross[[i]] <- rbind(cross[[i]], cross[[i]])
}
}
else { # input object is assumed to be a "cross" object
n.chr <- nchr(cross)
for(i in 1:n.chr) {
if(!is.matrix(cross$geno[[i]]$map))
cross$geno[[i]]$map <-
rbind(cross$geno[[i]]$map, cross$geno[[i]]$map)
}
}
cross
}
######################################################################
# comparecrosses: verify that two cross objects have identical
# classes, chromosomes, markers, genotypes, maps,
# and phenotypes
######################################################################
comparecrosses <-
function(cross1, cross2, tol=1e-5)
{
if(missing(cross1) || missing(cross2))
stop("Two crosses must be input.")
# both are of class "cross"
if(!any(class(cross1) == "cross") ||
!any(class(cross2) == "cross"))
stop("Input should have class \"cross\".")
# classes are the same
if(any(class(cross1) != class(cross2)))
stop("crosses are not the same type.")
if(nchr(cross1) != nchr(cross2))
stop("crosses do not have the same number of chromosomes.")
if(any(names(cross1$geno) != names(cross2$geno)))
stop("Chromosome names do not match.")
if(any(nmar(cross1) != nmar(cross2)))
stop("Number of markers per chromosome do not match.")
mnames1 <- unlist(lapply(cross1$geno, function(a) colnames(a$data)))
mnames2 <- unlist(lapply(cross2$geno, function(a) colnames(a$data)))
if(any(mnames1 != mnames2)) {
# stop("Markers names do not match.")
for(i in 1:nchr(cross1))
if(any(colnames(cross1$geno[[i]]$data) != colnames(cross2$geno[[i]]$data)))
stop("Marker names on chr ", names(cross1$geno)[i], " don't match.")
}
chrtype1 <- sapply(cross1$geno, class)
chrtype2 <- sapply(cross2$geno, class)
if(any(chrtype1 != chrtype2))
stop("Chromosome types (autosomal vs X) do not match.")
for(i in 1:nchr(cross1)) {
if(any(abs(diff(cross1$geno[[i]]$map) - diff(cross2$geno[[i]]$map)) > tol))
stop("Genetic maps for chromosome ", names(cross1$geno)[i],
" do not match.")
if(abs(cross1$geno[[i]]$map[1] - cross2$geno[[i]]$map[1]) > tol)
warning("Initial marker positions for chromosome ", names(cross1$geno)[i],
" do not match.")
}
if(nind(cross1) != nind(cross2))
stop("Number of individuals do not match.")
for(i in 1:nchr(cross1)) {
g1 <- cross1$geno[[i]]$data
g2 <- cross2$geno[[i]]$data
if(any((is.na(g1) & !is.na(g2)) | (!is.na(g1) & is.na(g2)) |
(!is.na(g1) & !is.na(g2) & g1!=g2)))
stop("Genotype data for chromosome ", names(cross1$geno)[i],
" do not match.")
}
if(nphe(cross1) != nphe(cross2))
stop("Number of phenotypes do not match.")
if(any(names(cross1$pheno) != names(cross2$pheno)))
stop("Phenotype names do not match.")
for(i in 1:nphe(cross1)) {
phe1 <- cross1$pheno[,i]
phe2 <- cross2$pheno[,i]
if(is.numeric(phe1) & is.numeric(phe2)) {
phe1[phe1 == Inf] <- max(phe1[phe1 < Inf], na.rm=TRUE)+5
phe2[phe2 == Inf] <- max(phe2[phe2 < Inf], na.rm=TRUE)+5
phe1[phe1 == -Inf] <- min(phe1[phe1 > -Inf], na.rm=TRUE)-5
phe2[phe2 == -Inf] <- min(phe2[phe2 > -Inf], na.rm=TRUE)-5
if(any((is.na(phe1) & !is.na(phe2)) | (!is.na(phe1) & is.na(phe2)) |
(!is.na(phe1) & !is.na(phe2) & abs(phe1-phe2) > tol))) {
stop("Data for phenotype ", names(cross1$pheno)[i],
" do not match.")
}
}
else {
if(any((is.na(phe1) & !is.na(phe2)) | (!is.na(phe1) & is.na(phe2)) |
(!is.na(phe1) & !is.na(phe2) &
as.character(phe1) != as.character(phe2)))) {
stop("Data for phenotype ", names(cross1$pheno)[i], " do not match.")
}
}
}
cat("\tCrosses are identical.\n")
}
######################################################################
# move marker
# Move a marker to a new chromosome...placed at the end
######################################################################
movemarker <-
function(cross, marker, newchr, newpos)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
mnames <- unlist(lapply(cross$geno,function(a) colnames(a$data)))
chr <- rep(names(cross$geno),nmar(cross))
pos <- unlist(lapply(cross$geno,function(a) 1:ncol(a$data)))
oldindex <- match(marker, mnames)
# Marker found precisely once?
if(is.na(oldindex)) stop(marker, " not found.\n")
if(length(oldindex) > 1) stop(marker, " found multiple times.\n")
chr <- chr[oldindex]
pos <- pos[oldindex]
# pull out genotype data
g <- cross$geno[[chr]]$data[,pos]
chrtype <- class(cross$geno[[chr]])
mapmatrix <- is.matrix(cross$geno[[chr]]$map)
# delete marker
if(nmar(cross)[chr] == 1) { # only marker on that chromosome, so drop the chromosome
cross$geno <- cross$geno[-match(chr,names(cross$geno))]
delchr <- TRUE
}
else {
delchr <- FALSE
cross$geno[[chr]]$data <- cross$geno[[chr]]$data[,-pos,drop=FALSE]
if(is.matrix(cross$geno[[chr]]$map))
cross$geno[[chr]]$map <- cross$geno[[chr]]$map[,-pos,drop=FALSE]
else
cross$geno[[chr]]$map <- cross$geno[[chr]]$map[-pos]
}
if(is.numeric(newchr)) newchr <- as.character(newchr)
if(!(newchr %in% names(cross$geno))) { # create a new chromosome
n <- length(cross$geno)
cross$geno[[n+1]] <- list("data"=as.matrix(g),
"map"=as.numeric(0))
names(cross$geno)[n+1] <- newchr
class(cross$geno[[n+1]]) <- chrtype
colnames(cross$geno[[n+1]]$data) <- marker
if(mapmatrix) {
if(missing(newpos)) newpos <- 0
cross$geno[[n+1]]$map <- matrix(newpos, ncol=1, nrow=2)
colnames(cross$geno[[n+1]]$map) <- marker
}
else {
if(missing(newpos)) newpos <- 0
cross$geno[[n+1]]$map <- newpos
names(cross$geno[[n+1]]$map) <- marker
}
return(cross)
}
if(missing(newpos)) {
# add marker to end of new chromosome
n.mar <- nmar(cross)[newchr]
cross$geno[[newchr]]$data <- cbind(cross$geno[[newchr]]$data,g)
colnames(cross$geno[[newchr]]$data)[n.mar+1] <- marker
if(is.matrix(cross$geno[[newchr]]$map)) {
cross$geno[[newchr]]$map <- cbind(cross$geno[[newchr]]$map,
cross$geno[[newchr]]$map[,n.mar]+10)
colnames(cross$geno[[newchr]]$map)[n.mar+1] <- marker
}
else {
cross$geno[[newchr]]$map <- c(cross$geno[[newchr]]$map,
cross$geno[[newchr]]$map[n.mar]+10)
names(cross$geno[[newchr]]$map)[n.mar+1] <- marker
}
}
else {
# add marker to the specified position
dat <- cross$geno[[newchr]]$data
map <- cross$geno[[newchr]]$map
if(length(newpos) != 1)
stop("newpos should be a single number.")
if(is.matrix(map)) { # sex-specific maps
wh <- which(map[1,] < newpos)
if(length(wh) == 0) { # place in first spot
map <- cbind(c(newpos,map[2,1]-(map[1,1]-newpos)),map)
colnames(map)[1] <- marker
}
else {
wh <- max(wh)
if(wh == ncol(map)) { # place at end of chromosome
map <- cbind(map,c(newpos,map[2,ncol(map)]+(newpos-map[1,ncol(map)])))
colnames(map)[ncol(map)] <- marker
}
else {
left <- map[,1:wh,drop=FALSE]
right <- map[,-(1:wh),drop=FALSE]
marleft <- colnames(left)
marright <- colnames(right)
left <- left[,ncol(left)]
right <- right[,1]
newpos2 <- (newpos-left[1])/(right[1]-left[1])*(right[2]-left[2])+left[2]
map <- cbind(map[,1:wh], c(newpos,newpos2), map[,-(1:wh)])
colnames(map) <- c(marleft, marker, marright)
}
}
}
else {
wh <- which(map < newpos)
if(length(wh) == 0) { # place in first position
map <- c(newpos,map)
names(map)[1] <- marker
}
else {
wh <- max(wh)
if(wh == length(map)) { # place in last position
map <- c(map,newpos)
names(map)[length(map)] <- marker
}
else {
map <- c(map[1:wh],newpos,map[-(1:wh)])
names(map)[wh+1] <- marker
}
}
}
cross$geno[[newchr]]$map <- map
if(length(wh)==0) { # place marker in first position
dat <- cbind(g, dat)
colnames(dat)[1] <- marker
}
else if(wh == ncol(dat)) { # place marker in last position
dat <- cbind(dat, g)
colnames(dat)[ncol(dat)] <- marker
}
else { # place marker in the middle
dat <- cbind(dat[,1:wh],g,dat[,-(1:wh)])
colnames(dat)[wh+1] <- marker
}
cross$geno[[newchr]]$data <- dat
# make sure the marker names for the data and the genetic map match
if(is.matrix(cross$geno[[newchr]]$map))
colnames(cross$geno[[newchr]]$data) <- colnames(cross$geno[[newchr]]$map)
else
colnames(cross$geno[[newchr]]$data) <- names(cross$geno[[newchr]]$map)
}
# update genoprob, errorlod, argmax, draws, rf
if("rf" %in% names(cross)) {
# reorder the recombination fractions
# -- a bit of pain, 'cause we need LODs in upper triangle
# and rec fracs in lower triangle
newmar <- unlist(lapply(cross$geno,function(a) colnames(a$data)))
attrib <- attributes(cross$rf)
rf <- cross$rf
lods <- rf;lods[lower.tri(rf)] <- t(rf)[lower.tri(rf)]
rf[upper.tri(rf)] <- t(rf)[upper.tri(rf)]
lods <- lods[newmar,newmar]
rf <- rf[newmar,newmar]
rf[upper.tri(rf)] <- lods[upper.tri(rf)]
cross$rf <- rf
if("onlylod" %in% names(attrib)) # save the onlylod attribute if its there
attr(cross$rf, "onlylod") <- attrib$onlylod
}
if(!delchr) thechr <- c(chr,newchr)
else thechr <- newchr
for(i in thechr) {
tempg <- cross$geno[[i]]
tempx <- subset(cross, chr=i)
if("prob" %in% names(tempg))
atp <- attributes(tempg$prob)
if("draws" %in% names(tempg)) {
at <- attributes(cross$geno[[1]]$draws)
tempg$draws <- sim.geno(tempx,
n.draws=at$dim[3],
step=at$step,
off.end=at$off.end,
map.function=at$map.function,
error.prob=at$error.prob)$geno[[1]]$draws
}
if("argmax" %in% names(tempg)) {
at <- attributes(cross$geno[[1]]$argmax)
tempg$argmax <- argmax.geno(tempx,
step=at$step,
off.end=at$off.end,
map.function=at$map.function,
error.prob=at$error.prob)$geno[[1]]$argmax
}
if("errorlod" %in% names(tempg)) {
at <- attributes(cross$geno[[1]]$errorlod)
tempg$errorlod <- argmax.geno(tempx,
map.function=at$map.function,
error.prob=at$error.prob)$geno[[1]]$errorlod
}
if("prob" %in% names(tempg))
tempg$prob <- calc.genoprob(tempx,
step=atp$step,
off.end=atp$off.end,
map.function=atp$map.function,
error.prob=atp$error.prob)$geno[[1]]$prob
cross$geno[[i]] <- tempg
}
cross
}
######################################################################
#
# summary.map
#
# Give a short summary of a genetic map object.
#
######################################################################
summaryMap <- summary.map <-
function(object, ...)
{
map <- object
if(any(class(map) == "cross")) # a cross object
map <- pull.map(map)
if(!any(class(map) == "map"))
stop("Input should have class \"cross\" or \"map\".")
n.chr <- length(map)
chrnames <- names(map)
if(is.null(chrnames)) chrnames <- 1:length(map)
if(is.matrix(map[[1]])) { # sex-specific map
sexsp <- TRUE
n.mar <- sapply(map,ncol)
tot.mar <- sum(n.mar)
fmap <- lapply(map,function(a) a[1,])
mmap <- lapply(map,function(a) a[2,])
len.f <- sapply(fmap,function(a) diff(range(a)))
len.m <- sapply(mmap,function(a) diff(range(a)))
avesp.f <- sapply(fmap,function(a) {if(length(a)<2) return(NA); mean(diff(a))})
avesp.m <- sapply(mmap,function(a) {if(length(a)<2) return(NA); mean(diff(a))})
maxsp.f <- sapply(fmap,function(a) {if(length(a)<2) return(NA); max(diff(a))})
maxsp.m <- sapply(mmap,function(a) {if(length(a)<2) return(NA); max(diff(a))})
totlen.f <- sum(len.f)
totlen.m <- sum(len.m)
tot.avesp.f <- mean(unlist(lapply(fmap,diff)), na.rm=TRUE)
tot.avesp.m <- mean(unlist(lapply(mmap,diff)), na.rm=TRUE)
tot.maxsp.f <- max(maxsp.f,na.rm=TRUE)
tot.maxsp.m <- max(maxsp.m,na.rm=TRUE)
output <- rbind(cbind(n.mar,len.f,len.m,avesp.f,avesp.m, maxsp.f, maxsp.m),
c(tot.mar,totlen.f,totlen.m,tot.avesp.f,tot.avesp.m,
tot.maxsp.f, tot.maxsp.m))
dimnames(output) <- list(c(chrnames,"overall"),
c("n.mar","length.female","length.male",
"ave.spacing.female","ave.spacing.male",
"max.spacing.female", "max.spacing.male"))
}
else {
sexsp=FALSE
n.mar <- sapply(map,length)
len <- sapply(map,function(a) diff(range(a)))
tot.mar <- sum(n.mar)
len <- sapply(map,function(a) diff(range(a)))
avesp <- sapply(map,function(a) {if(length(a)<2) return(NA); mean(diff(a))})
maxsp <- sapply(map,function(a) {if(length(a)<2) return(NA); max(diff(a))})
totlen <- sum(len)
tot.avesp <- mean(unlist(lapply(map,diff)), na.rm=TRUE)
tot.maxsp <- max(maxsp, na.rm=TRUE)
output <- rbind(cbind(n.mar,len,avesp, maxsp),
c(tot.mar,totlen,tot.avesp, tot.maxsp))
dimnames(output) <- list(c(chrnames,"overall"),
c("n.mar","length","ave.spacing", "max.spacing"))
}
output <- as.data.frame(output)
attr(output, "sexsp") <- sexsp
class(output) <- c("summary.map", "data.frame")
output
}
######################################################################
#
# print.summary.map
#
# Print out the result of summary.map()
#
######################################################################
print.summary.map <-
function(x, ...)
{
sexsp <- attr(x, "sexsp")
if(sexsp) cat("Sex-specific map\n\n")
x <- apply(x,2,round,1)
print(x)
}
######################################################################
# convert functions
######################################################################
convert <-
function(object, ...)
UseMethod("convert")
######################################################################
#
# convert.scanone
#
# Convert scanone output from the format for R/qtl ver 0.97 to
# that for R/qtl ver 0.98
# (previously, inter-maker locations named loc*.c*; now c*.loc*)
#
######################################################################
convert.scanone <-
function(object, ...)
{
if(!any(class(object) == "scanone"))
stop("Input should have class \"scanone\".")
rn <- rownames(object)
o <- grep("^loc-*[0-9]+(\\.[0-9]+)*\\.c[0-9A-Za-z]+$", rn)
if(length(o) > 0) {
temp <- rn[o]
temp <- strsplit(temp,"\\.")
temp <- sapply(temp, function(a)
paste(a[c(length(a),1:(length(a)-1))],collapse="."))
rownames(object)[o] <- temp
}
object
}
######################################################################
# convert.scantwo
#
# convert scantwo output from the format used in R/qtl version 1.03
# and earlier to that used in R/qtl version 1.04 and later.
#
# 1.03 and earlier: contained joint and interactive LOD scores
# 1.04 and later: contains joint LOD scores and LOD scores from the
# additive QTL model
######################################################################
convert.scantwo <-
function(object, ...)
{
if(!any(class(object) == "scantwo"))
stop("Input should have class \"scantwo\".")
lod <- object$lod
if(length(dim(lod)) == 2) {
u <- upper.tri(lod)
lod[u] <- t(lod)[u] - lod[u]
}
else { # multiple phenotypes
u <- upper.tri(lod[,,1])
for(i in 1:dim(lod)[3])
lod[,,i][u] <- t(lod[,,i])[u] - lod[,,i][u]
}
object$lod <- lod
object
}
######################################################################
# convert.map
#
# convert a genetic map from one map function to another
######################################################################
convert.map <-
function(object, old.map.function=c("haldane", "kosambi", "c-f", "morgan"),
new.map.function=c("haldane", "kosambi", "c-f", "morgan"), ...)
{
old.map.function <- match.arg(old.map.function)
new.map.function <- match.arg(new.map.function)
if(!("map" %in% class(object)))
stop("Input should have class \"map\".")
if(old.map.function==new.map.function) {
warning("old and new map functions are the same; no change.")
return(object)
}
mf <- switch(old.map.function,
"haldane"=mf.h,
"kosambi"=mf.k,
"c-f"=mf.cf,
"morgan"=mf.m)
imf <- switch(new.map.function,
"haldane"=imf.h,
"kosambi"=imf.k,
"c-f"=imf.cf,
"morgan"=imf.m)
if(is.matrix(object[[1]])) { # sex-specific map
for(i in seq(along=object)) {
theclass <- class(object[[i]])
thenames <- colnames(object[[i]])
for(j in 1:2)
object[[i]][j,] <- cumsum(c(object[[i]][j,1], imf(mf(diff(object[[i]][j,])))))
class(object[[i]]) <- theclass
colnames(object[[i]]) <- thenames
}
}
else {
for(i in seq(along=object)) {
theclass <- class(object[[i]])
thenames <- names(object[[i]])
object[[i]] <- cumsum(c(object[[i]][1], imf(mf(diff(object[[i]])))))
class(object[[i]]) <- theclass
names(object[[i]]) <- thenames
}
}
object
}
######################################################################
# find.pheno
#
# utility to get pheno number given pheno name
######################################################################
find.pheno <-
function( cross, pheno )
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
seq( ncol( cross$pheno ))[match(pheno,names(cross$pheno))]
}
######################################################################
# find.flanking
#
# utility to get flanking and/or closest marker to chr and pos
######################################################################
find.flanking <-
function( cross, chr, pos)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
map = pull.map(cross)
if(is.matrix(map[[1]]) && nrow(map[[1]]) > 1)
stop("This function works only for crosses with sex-averaged maps.")
if(length(chr) == 1 && length(pos) > 1) {
chr <- rep(chr,length(pos))
}
wh <- match(chr, names(cross$geno))
if(any(is.na(wh))) {
stop("Chr ", paste(chr[is.na(wh)], collapse=", "), " not found.")
wh <- wh[!is.na(wh)]
}
chr <- names(cross$geno)[wh]
marker = NULL
for (i in seq(length(chr))) {
tmp = map[[chr[i]]]-pos[i]
m = names(map[[chr[i]]])
left = sum(tmp < 0)
at = sum(tmp == 0)
right = sum(tmp > 0)
f <- if (at > 0)
left+at[c(1,length(at))]
else {
if (right > 0)
c(left,left+at+1)
else
c(left,left+at)
}
marker = rbind(marker,m[f[c(1:2,order(abs(tmp[f]))[1])]])
}
dimnames(marker) <- list(paste(chr,":",pos,sep=""),
c("left","right","close"))
as.data.frame(marker, stringsAsFactors=TRUE)
}
######################################################################
# strip.partials
#
# Removes partially informative genotypes in an intercross.
#
# Input: a cross object; if verbose=TRUE, a statement regarding the
# number of genotypes removed is printed.
######################################################################
strip.partials <-
function(cross, verbose=TRUE)
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
type <- class(cross)[1]
if(type != "f2")
stop("This is for intercrosses only")
n.removed <- 0
for(i in 1:nchr(cross)) {
g <- cross$geno[[i]]$data
wh <- !is.na(g) & g>3
if(any(wh)) {
g[wh] <- NA
cross$geno[[i]]$data <- g
n.removed <- n.removed + sum(wh)
}
}
if(verbose) {
if(n.removed == 0) cat(" --Didn't remove any genotypes.\n")
else cat("Removed", n.removed, "genotypes.\n")
}
cross
}
######################################################################
# comparegeno
######################################################################
comparegeno <-
function(cross, what=c("proportion","number", "both"))
{
if(!any(class(cross) == "cross"))
stop("Input should have class \"cross\".")
what <- match.arg(what)
g <- pull.geno(cross)
g[is.na(g)] <- 0
n.ind <- nrow(g)
n.mar <- ncol(g)
z <- .C("R_comparegeno",
as.integer(g),
as.integer(n.ind),
as.integer(n.mar),
n.match=as.integer(rep(0,n.ind^2)),
n.missing=as.integer(rep(0,n.ind^2)),
PACKAGE="qtl")
if(what=="number") {
z <- matrix(z$n.match,n.ind,n.ind)
}
else {
if(what=="proportion") {
z <- matrix(z$n.match/(n.mar-z$n.missing),n.ind,n.ind)
diag(z) <- NA
}
else {
prop <- matrix(z$n.match/(n.mar-z$n.missing),n.ind,n.ind)
z <- matrix(z$n.match,n.ind,n.ind)
z[lower.tri(z)] <- prop[lower.tri(z)]
}
}
z
}
######################################################################
# print the installed version of R/qtl
######################################################################
qtlversion <-
function()
{
version <- unlist(packageVersion("qtl"))
# make it like #.#-#
paste(c(version,".","-")[c(1,4,2,5,3)], collapse="")
}
######################################################################
#
# locateXO
#
# Locate crossovers on a single chromosome in each individual
# Look at just the first chromosome
#
######################################################################
locateXO <-
function(cross, chr, full.info=FALSE)
{
if(!missing(chr)) {
cross <- subset(cross, chr=chr)
if(nchr(cross) != 1)
warning("locateXO works on just one chr; considering chr ", names(cross$geno)[1])
}
# individual IDs
id <- getid(cross)
if(is.null(id)) id <- as.character(1:nind(cross))
if(nmar(cross)[1] == 1) { # just one marker; don't need to do anything
warning("Just one marker.")
res <- vector("list", nind(cross))
names(res) <- id
for(i in seq(along=res)) res[[i]] <- numeric(0)
return(res)
}
geno <- cross$geno[[1]]$data
geno[is.na(geno)] <- 0
type <- class(cross)[1]
if(type != "bc" && type != "f2" && type != "riself" && type != "risib" && type!="dh")
stop("locateXO only working for backcross, intercross or RI strains.")
map <- cross$geno[[1]]$map
if(is.matrix(map)) map <- map[1,]
# map <- map - map[1] # shift first marker to 0
# bc or intercross? thetype==0 for BC and ==1 for intercross
if(type=="f2") {
if(class(geno) == "X") thetype <- 0
else thetype <- 1
}
else thetype <- 0
n.ind <- nrow(geno)
n.mar <- ncol(geno)
z <- .C("R_locate_xo",
as.integer(n.ind),
as.integer(n.mar),
as.integer(thetype),
as.integer(geno),
as.double(map),
location=as.double(rep(0,n.ind*2*(n.mar-1))),
nseen=as.integer(rep(0,n.ind)),
ileft=as.integer(rep(0,n.ind*2*(n.mar-1))),
iright=as.integer(rep(0,n.ind*2*(n.mar-1))),
left=as.double(rep(0,n.ind*2*(n.mar-1))),
right=as.double(rep(0,n.ind*2*(n.mar-1))),
gleft=as.integer(rep(0, n.ind*2*(n.mar-1))),
gright=as.integer(rep(0, n.ind*2*(n.mar-1))),
ntype=as.integer(rep(0,n.ind*2*(n.mar-1))),
as.integer(full.info),
PACKAGE="qtl")
location <- t(matrix(z$location, nrow=n.ind))
nseen <- z$nseen
if(full.info) {
ileft <- t(matrix(z$ileft, nrow=n.ind))
iright <- t(matrix(z$iright, nrow=n.ind))
left <- t(matrix(z$left, nrow=n.ind))
right <- t(matrix(z$right, nrow=n.ind))
gleft <- t(matrix(z$gleft, nrow=n.ind))
gright <- t(matrix(z$gright, nrow=n.ind))
ntype <- t(matrix(z$ntype, nrow=n.ind))
}
if(!full.info)
res <- lapply(as.data.frame(rbind(nseen, location), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
else {
location <- lapply(as.data.frame(rbind(nseen, location), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
ileft <- lapply(as.data.frame(rbind(nseen, ileft), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
iright <- lapply(as.data.frame(rbind(nseen, iright), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
left <- lapply(as.data.frame(rbind(nseen, left), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
right <- lapply(as.data.frame(rbind(nseen, right), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
gleft <- lapply(as.data.frame(rbind(nseen, gleft), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
gright <- lapply(as.data.frame(rbind(nseen, gright), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
ntype <- lapply(as.data.frame(rbind(nseen, ntype), stringsAsFactors=TRUE),
function(a) { if(a[1]==0) return(numeric(0)); a[(1:a[1])+1] })
res <- location
for(i in seq(along=res)) {
if(length(res[[i]])>0) {
ntype[[i]][length(ntype[[i]])] <- NA
res[[i]] <- cbind(location=location[[i]],
left=left[[i]],
right=right[[i]],
ileft=ileft[[i]],
iright=iright[[i]],
gleft=gleft[[i]],
gright=gright[[i]],
nTypedBetween=ntype[[i]])
}
}
}
names(res) <- id
res
}
# jittermap: make sure no two markers are at precisely the same position
jittermap <-
function(object, amount=1e-6) # x is either a cross object or a map
{
if(any(class(object) == "cross")) {
themap <- pull.map(object)
return.cross <- TRUE
}
else {
if(!any(class(object) == "map"))
stop("Input must be a cross or a map")
return.cross <- FALSE
themap <- object
}
for(i in 1:length(themap)) {
if(is.matrix(themap[[i]])) { # sex-specific maps
n <- ncol(themap[[i]])
if(n > 1) {
for(j in 1:2)
themap[[i]][j,] <- themap[[i]][j,] + c(0,cumsum(rep(amount, n-1)))
}
}
else {
n <- length(themap[[i]])
if(n > 1)
themap[[i]] <- themap[[i]] + c(0,cumsum(rep(amount, n-1)))
}
}
if(return.cross)
return(clean(replace.map(object, themap)))
themap
}
print.map <-
function(x, ...)
{
for(i in seq(along=x))
attr(x[[i]], "loglik") <- NULL
if(length(x) == 1)
print(unclass(x[[1]]))
else
print(unclass(lapply(x, unclass)))
}
# map function for Stahl model
mf.stahl <-
function(d, m=0, p=0)
{
if(any(d < 0))
stop("d must be >= 0\n")
if(m < 0)
stop("Must have m >= 0\n")
if(p < 0 || p > 1)
stop("Must have 0 <= p <= 1\n")
# handle missing values
if(any(!is.finite(d))) {
which.finite <- is.finite(d)
dsub <- d[which.finite]
dropped.some <- TRUE
} else {
dsub <- d
dropped.some <- FALSE
}
result <- .C("R_mf_stahl",
as.integer(length(dsub)),
as.double(dsub/100), # convert to Morgans
as.integer(m),
as.double(p),
result=as.double(rep(0,length(dsub))),
PACKAGE="qtl")$result
if(dropped.some) {
d[!which.finite] <- NA
d[which.finite] <- result
result <- d
}
result
}
# inverse map function for Stahl model
imf.stahl <-
function(r, m=0, p=0, tol=1e-12, maxit=1000)
{
if(any(r < 0 | r > 0.5))
stop("r must be >= 0 or <= 0.5\n")
if(m < 0)
stop("Must have m >= 0\n")
if(p < 0 || p > 1)
stop("Must have 0 <= p <= 1\n")
# handle missing values
if(any(!is.finite(r))) {
which.finite <- is.finite(r)
rsub <- r[which.finite]
dropped.some <- TRUE
} else {
rsub <- r
dropped.some <- FALSE
}
result <- .C("R_imf_stahl",
as.integer(length(rsub)),
as.double(rsub),
as.integer(m),
as.double(p),
result=as.double(rep(0,length(rsub))),
as.double(tol),
as.integer(maxit),
PACKAGE="qtl")$result*100 # convert to cM
if(dropped.some) {
r[!which.finite] <- NA
r[which.finite] <- result
result <- r
}
result
}
######################################################################
# getid: internal function to pull out the "ID" column from the
# phenotype data, if there is one
######################################################################
getid <-
function(cross)
{
phe <- cross$pheno
nam <- names(phe)
if("id" %in% nam) {
id <- phe$id
phenam <- "id"
}
else if("ID" %in% nam) {
id <- phe$ID
phenam <- "ID"
}
else if("Id" %in% nam) {
id <- phe$Id
phenam <- "Id"
}
else if("iD" %in% nam) {
id <- phe$iD
phenam <- "iD"
}
else {
id <- NULL
phenam <- NULL
}
if(is.factor(id))
id <- as.character(id)
attr(id, "phenam") <- phenam
id
}
######################################################################
# find the chromosome and position of a vector of markers
######################################################################
find.markerpos <-
function(cross, marker)
{
if(length(marker) != length(unique(marker))) {
warning("Dropping duplicate marker names.")
marker <- unique(marker)
}
output <- data.frame(chr=rep("", length(marker)),
pos=rep(NA, length(marker)), stringsAsFactors=TRUE)
output$chr <- as.character(output$chr)
rownames(output) <- marker
map <- pull.map(cross)
n.mar <- nmar(cross)
chr <- rep(names(map), n.mar)
if(!is.matrix(map[[1]])) {
pos <- unlist(map)
onemap <- TRUE
}
else {
pos <- unlist(lapply(map, function(a) a[1,]))
pos2 <- unlist(lapply(map, function(a) a[2,]))
onemap <- FALSE
output <- cbind(output, pos2=rep(NA, length(marker)))
colnames(output)[2:3] <- c("pos.female","pos.male")
}
mnam <- markernames(cross)
for(i in seq(along=marker)) {
wh <- match(marker[i], mnam)
if(is.na(wh)) {
output[i,1] <- NA
output[i,2] <- NA
}
else {
if(length(wh) > 1) {
warning("Marker ", marker[i], " found multiple times.")
wh <- sample(wh, 1)
}
output[i,1] <- chr[wh]
output[i,2] <- pos[wh]
if(!onemap) output[i,3] <- pos2[wh]
}
}
output
}
######################################################################
# find the chromosome and position of a vector of pseudomarkers
######################################################################
find.pseudomarkerpos <-
function(cross, marker, where=c("draws","prob"))
{
if(length(marker) != length(unique(marker))) {
warning("Dropping duplicate pseudomarker names.")
marker <- unique(marker)
}
output <- data.frame(chr=rep("", length(marker)),
pos=rep(NA, length(marker)), stringsAsFactors=TRUE)
output$chr <- as.character(output$chr)
rownames(output) <- marker
where <- match.arg(where)
if(where=="draws" && !("draws" %in% names(cross$geno[[1]])))
stop("You'll need to first run sim.geno")
if(where=="prob" && !("prob" %in% names(cross$geno[[1]])))
stop("You'll need to first run calc.genoprob")
themap <- vector("list", nchr(cross))
names(themap) <- names(cross$geno)
for(i in 1:nchr(cross)) {
if(where=="draws")
themap[[i]] <- attr(cross$geno[[i]]$draws, "map")
else
themap[[i]] <- attr(cross$geno[[i]]$prob, "map")
}
if(!is.matrix(themap[[1]])) {
pmar <- unlist(lapply(themap, names))
pos <- unlist(themap)
onemap <- TRUE
chr <- rep(names(themap), sapply(themap, length))
}
else {
pos <- unlist(lapply(themap, function(a) a[1,]))
pos2 <- unlist(lapply(themap, function(a) a[2,]))
onemap <- FALSE
pmar <- unlist(lapply(themap, colnames))
output <- cbind(output, pos2=rep(NA, length(marker)))
colnames(output)[2:3] <- c("pos.female","pos.male")
chr <- rep(names(themap), sapply(themap, ncol))
}
whnotmarker <- grep("^loc-*[0-9]*", pmar)
pmar[whnotmarker] <- paste("c", chr[whnotmarker], ".", pmar[whnotmarker], sep="")
for(i in seq(along=marker)) {
wh <- match(marker[i], pmar)
if(is.na(wh)) {
output[i,1] <- NA
output[i,2] <- NA
}
else {
if(length(wh) > 1) {
warning("Pseudomarker ", marker[i], " found multiple times.")
wh <- sample(wh, 1)
}
output[i,1] <- chr[wh]
output[i,2] <- pos[wh]
if(!onemap) output[i,3] <- pos2[wh]
}
}
output
}
######################################################################
# utility function for determining whether pheno.col (as argument
# to scanone etc) can be interpreted as a vector of phenotypes,
# versus a vector of phenotype columns
######################################################################
LikePheVector <-
function(pheno, n.ind, n.phe)
{
if(is.numeric(pheno) && length(pheno)==n.ind &&
any(pheno < 1 | pheno > n.phe | pheno!=round(pheno)))
return(TRUE)
FALSE
}
######################################################################
# for matching chromosome names
######################################################################
matchchr <-
function(selection, thechr)
{
if(is.factor(thechr)) thechr <- as.character(thechr)
if(length(thechr) > length(unique(thechr)))
stop("Duplicate chromosome names.")
if(is.logical(selection)) {
if(length(selection) != length(thechr))
stop("Logical vector to select chromosomes is the wrong length")
return(thechr[selection])
}
if(is.numeric(selection)) selection <- as.character(selection)
if(length(selection) > length(unique(selection))) {
warning("Dropping duplicate chromosomes")
selection <- unique(selection)
}
g <- grep("^-", selection)
if(length(g) > 0 && length(g) < length(selection))
stop("In selecting chromosomes, all must start with '-' or none should.")
if(length(g) > 0) {
selectomit <- TRUE
selection <- substr(selection, 2, nchar(selection))
}
else selectomit <- FALSE
wh <- match(selection, thechr)
if(any(is.na(wh))) {
warning("Chr ", paste(selection[is.na(wh)], collapse=", "), " not found")
wh <- wh[!is.na(wh)]
if(length(wh) == 0) return(thechr)
}
if(selectomit) return(thechr[-wh])
thechr[sort(wh)]
}
######################################################################
# check that chromosomes match appropriately
# TRUE = chr okay
# FALSE = problem
######################################################################
testchr <-
function(selection, thechr)
{
if(is.factor(thechr)) thechr <- as.character(thechr)
if(length(thechr) > length(unique(thechr))) {
# warning("Duplicate chromosome names.")
return(FALSE)
}
if(is.logical(selection)) {
if(length(selection) != length(thechr)) {
# warning("Logical vector to select chromosomes is the wrong length")
return(FALSE)
}
return(TRUE)
}
if(is.numeric(selection)) selection <- as.character(selection)
if(length(selection) > length(unique(selection))) {
# warning("Dropping duplicate chromosomes")
selection <- unique(selection)
}
g <- grep("^-", selection)
if(length(g) > 0 && length(g) < length(selection)) {
# stop("In selecting chromosomes, all must start with '-' or none should.")
return(FALSE)
}
if(length(g) > 0) {
selectomit <- TRUE
selection <- substr(selection, 2, nchar(selection))
}
else selectomit <- FALSE
wh <- match(selection, thechr)
if(any(is.na(wh))) return(FALSE)
TRUE
}
######################################################################
# convert2sa
#
# convert a sex-specific maps to a sex-averaged one.
# We pull out just the female map, and give a warning if the male and
# female maps are too different
######################################################################
convert2sa <-
function(map, tol=1e-4)
{
if(!("map" %in% class(map)))
stop("Input should have class \"map\".")
if(!is.matrix(map[[1]]))
stop("Input map doesn't seem to be a sex-specific map.")
theclass <- sapply(map, class)
fem <- lapply(map, function(a) a[1,])
dif <- sapply(map, function(a) { if(ncol(a)==1) return(diff(a))
a <- apply(a, 1, diff);
if(is.matrix(a)) return(max(abs(apply(a, 1, diff))))
abs(diff(a)) })
if(max(dif) > tol)
warning("Female and male inter-marker distances differ by as much as ", max(dif), ".")
for(i in seq(along=theclass)) class(fem[[i]]) <- theclass[i]
class(fem) <- "map"
fem
}
# round as character string, ensuring ending 0's are kept.
charround <-
function(x, digits=1)
{
if(digits < 1)
stop("This is intended for the case digits >= 1.")
y <- as.character(round(x, digits))
z <- strsplit(y, "\\.")
sapply(z, function(a, digits)
{
if(length(a) == 1)
b <- paste(a[1], ".", paste(rep("0", digits),collapse=""), sep="")
else {
if(nchar(a[2]) == digits)
b <- paste(a, collapse=".")
else
b <- paste(a[1], ".", a[2],
paste(rep("0", digits - nchar(a[2])), collapse=""),
sep="")
}
}, digits)
}
######################################################################
# scantwoperm2scanoneperm
#
# pull out the scanone permutations from scantwo permutation results,
# so that one may use the scantwo perms in calls to summary.scanone
######################################################################
scantwoperm2scanoneperm <-
function(scantwoperms)
{
if(!("scantwoperm" %in% class(scantwoperms)))
stop("Input must have class \"scantwoperm\".")
if(!("one" %in% names(scantwoperms)))
stop("Input doesn't contain scanone permutation results.")
scanoneperms <- scantwoperms$one
class(scanoneperms) <- c("scanoneperm")
scanoneperms
}
######################################################################
# subset.map
######################################################################
subset.map <-
function(x, ...)
{
cl <- class(x)
x <- x[...]
class(x) <- cl
x
}
`[.map` <-
function(x, ...)
{
x <- unclass(x)[...]
class(x) <- "map"
x
}
######################################################################
# subset.cross with [ ]
######################################################################
`[.cross` <-
function(x, chr, ind)
subset(x, chr, ind)
######################################################################
# findDupMarkers
#
# find markers whose genotype data is identical to another marker
# (which thus might be dropped, as extraneous)
#
# chr (Optional) set of chromosomes to consider
#
# exact.only If TRUE, look for markers with the same genotypes
# and the same pattern of missing data
# If FALSE, also identify markers whose observed
# genotypes match another marker, with no
# observed genotypes for which the other is
# missing
#
# adjacent.only If TRUE, only consider adjacent markers
######################################################################
findDupMarkers <-
function(cross, chr, exact.only=TRUE, adjacent.only=FALSE)
{
if(!missing(chr)) cross <- subset(cross, chr=chr)
g <- pull.geno(cross)
markers <- colnames(g)
markerloc <- lapply(nmar(cross), function(a) 1:a)
if(length(markerloc) > 1) {
for(j in 2:length(markerloc))
markerloc[[j]] <- markerloc[[j]] + max(markerloc[[j-1]]) + 10
}
markerloc <- unlist(markerloc)
if(exact.only) {
g[is.na(g)] <- 0
# genotype data patterns
pat <- apply(g, 2, paste, collapse="")
# table of unique values
tab <- table(pat)
# no duplicates; return
if(all(tab == 1)) return(NULL)
wh <- which(tab > 1)
theloc <- themar <- vector("list", length(wh))
for(i in seq(along=wh)) {
themar[[i]] <- names(pat)[pat==names(tab)[wh[i]]]
theloc[[i]] <- markerloc[pat==names(tab)[wh[i]]]
}
if(adjacent.only) {
extraloc <- list()
extramar <- list()
for(i in seq(along=theloc)) {
d <- diff(theloc[[i]]) # look for adjacency
if(any(d>1)) { # split into adjacent groups
temp <- which(d>1)
first <- c(1, temp+1)
last <- c(temp, length(theloc[[i]]))
for(j in 2:length(first)) {
extraloc[[length(extraloc)+1]] <- theloc[[i]][first[j]:last[j]]
extramar[[length(extramar)+1]] <- themar[[i]][first[j]:last[j]]
}
themar[[i]] <- themar[[i]][first[1]:last[1]]
theloc[[i]] <- theloc[[i]][first[1]:last[1]]
}
}
themar <- c(themar, extramar)
theloc <- c(theloc, extraloc)
nm <- sapply(themar, length)
if(all(nm==1)) return(NULL) # nothing left
themar <- themar[nm>1]
theloc <- theloc[nm>1]
}
# order by first locus
o <- order(sapply(theloc, min))
themar <- themar[o]
randompics <- sapply(themar, function(a) sample(length(a), 1))
for(i in seq(along=themar)) {
names(themar)[i] <- themar[[i]][randompics[i]]
themar[[i]] <- themar[[i]][-randompics[i]]
}
}
else {
themar <- NULL
ntyp <- ntyped(cross, "mar")
o <- order(ntyp, decreasing=TRUE)
g[is.na(g)] <- 0
z <- .C("R_findDupMarkers_notexact",
as.integer(nrow(g)),
as.integer(ncol(g)),
as.integer(g),
as.integer(o),
as.integer(markerloc),
as.integer(adjacent.only),
result=as.integer(rep(0,length(o))),
PACKAGE="qtl")
if(all(z$result==0)) return(NULL)
u <- unique(z$result[z$result != 0])
themar <- vector("list", length(u))
names(themar) <- colnames(g)[u]
for(i in seq(along=themar))
themar[[i]] <- colnames(g)[z$result==u[i]]
}
themar
}
######################################################################
# convert2riself
######################################################################
convert2riself <-
function(cross)
{
if(class(cross)[2] != "cross")
stop("input must be a cross object.")
curtype <- class(cross)[1]
chrtype <- sapply(cross$geno, class)
whX <- NULL
if(any(chrtype != "A")) { # there's an X chromosome
whX <- names(cross$geno)[chrtype != "A"]
if(length(whX) > 1)
warning("Converting chromosomes ", paste(whX, collapse=" "), " to autosomal.")
else
warning("Converting chromosome ", whX, " to autosomal.")
for(i in whX) class(cross$geno[[i]]) <- "A"
}
gtab <- table(pull.geno(cross))
usethree <- FALSE
if(!is.na(gtab["3"])) {
if(!is.na(gtab["2"]) && gtab["3"] < gtab["2"])
usethree <- FALSE
else usethree <- TRUE
}
g2omit <- g3omit <- g4omit <- 0
for(i in 1:nchr(cross)) {
dat <- cross$geno[[i]]$data
g1 <- sum(!is.na(dat) & dat==1)
g2 <- sum(!is.na(dat) & dat==2)
g3 <- sum(!is.na(dat) & dat==3)
g4 <- sum(!is.na(dat) & dat>3)
if(usethree && chrtype[i] == "A") {
dat[!is.na(dat) & dat!=1 & dat!=3] <- NA
dat[!is.na(dat) & dat==3] <- 2
g2omit <- g2omit + g2
g4omit <- g4omit + g4
}
else {
dat[!is.na(dat) & dat!=1 & dat!=2] <- NA
g3omit <- g3omit + g3
g4omit <- g4omit + g4
}
cross$geno[[i]]$data <- dat
}
if(g2omit > 0)
warning("Omitting ", g2omit, " genotypes with code==2.")
if(g3omit > 0)
warning("Omitting ", g3omit, " genotypes with code==3.")
if(g4omit > 0)
warning("Omitting ", g4omit, " genotypes with code>3.")
class(cross)[1] <- "riself"
cross
}
######################################################################
# convert2risib
######################################################################
convert2risib <-
function(cross)
{
if(class(cross)[2] != "cross")
stop("input must be a cross object.")
curtype <- class(cross)[1]
chrtype <- sapply(cross$geno, class)
gtab <- table(pull.geno(cross))
usethree <- FALSE
if(!is.na(gtab["3"])) {
if(!is.na(gtab["2"]) && gtab["3"] < gtab["2"])
usethree <- FALSE
else usethree <- TRUE
}
g2omit <- g3omit <- g4omit <- 0
for(i in 1:nchr(cross)) {
dat <- cross$geno[[i]]$data
g1 <- sum(!is.na(dat) & dat==1)
g2 <- sum(!is.na(dat) & dat==2)
g3 <- sum(!is.na(dat) & dat==3)
g4 <- sum(!is.na(dat) & dat>3)
if(usethree) {
if(chrtype[i] == "A") {
dat[!is.na(dat) & dat!=1 & dat!=3] <- NA
dat[!is.na(dat) & dat==3] <- 2
g2omit <- g2omit + g2
g4omit <- g4omit + g4
}
else { # X chromosome
if(g2 >= g3) {
dat[!is.na(dat) & dat!=1 & dat!=2] <- NA
g3omit <- g3omit + g3
g4omit <- g4omit + g4
}
else {
dat[!is.na(dat) & dat!=1 & dat!=3] <- NA
dat[!is.na(dat) & dat==3] <- 2
g2omit <- g2omit + g2
g4omit <- g4omit + g4
}
}
}
else {
dat[!is.na(dat) & dat!=1 & dat!=2] <- NA
g3omit <- g3omit + g3
g4omit <- g4omit + g4
}
cross$geno[[i]]$data <- dat
}
if(g2omit > 0)
warning("Omitting ", g2omit, " genotypes with code==2.")
if(g3omit > 0)
warning("Omitting ", g3omit, " genotypes with code==3.")
if(g4omit > 0)
warning("Omitting ", g4omit, " genotypes with code>3.")
class(cross)[1] <- "risib"
cross
}
rescalemap <-
function(object, scale=1e-6)
{
if("cross" %in% class(object)) {
for(i in 1:nchr(object)) {
object$geno[[i]]$map <-
object$geno[[i]]$map * scale
}
if(abs(scale - 1) > 1e-6)
object <- clean(object) # strip off intermediate calculations
} else if("map" %in% class(object)) {
for(i in seq(along=object)) {
object[[i]] <- object[[i]] * scale
}
} else
stop("rescalemap works only for objects of class \"cross\" or \"map\".")
object
}
shiftmap <-
function(object, offset=0)
{
if("cross" %in% class(object)) {
if(length(offset) == 1) offset <- rep(offset, nchr(object))
else if(length(offset) != nchr(object))
stop("offset must have length 1 or n.chr (", nchr(object), ")")
for(i in 1:nchr(object)) {
if(is.matrix(object$geno[[i]]$map)) {
for(j in 1:2)
object$geno[[i]]$map[j,] <- object$geno[[i]]$map[j,] - object$geno[[i]]$map[j,1] + offset[i]
} else {
object$geno[[i]]$map <- object$geno[[i]]$map - object$geno[[i]]$map[1] + offset[i]
}
}
} else if("map" %in% class(object)) {
if(length(offset) == 1) offset <- rep(offset, length(object))
else if(length(offset) != length(object))
stop("offset must have length 1 or n.chr (", length(object), ")")
for(i in seq(along=object)) {
if(is.matrix(object[[i]])) {
for(j in 1:2)
object[[i]][j,] <- object[[i]][j,] - object[[i]][j,1] + offset[i]
} else {
object[[i]] <- object[[i]] - object[[i]][1] + offset[i]
}
}
} else
stop("shiftmap works only for objects of class \"cross\" or \"map\".")
object
}
######################################################################
# switch alleles in a cross
######################################################################
switchAlleles <-
function(cross, markers, switch=c("AB","CD","ABCD", "parents"))
{
type <- class(cross)[1]
switch <- match.arg(switch)
if(type %in% c("bc", "risib", "riself", "dh")) {
if(switch != "AB")
warning("Using switch = \"AB\".")
found <- rep(FALSE, length(markers))
for(i in 1:nchr(cross)) {
cn <- colnames(cross$geno[[i]]$data)
m <- match(markers, cn)
if(any(!is.na(m))) {
found[!is.na(m)] <- TRUE
for(j in m[!is.na(m)]) {
g <- cross$geno[[i]]$data[,j]
cross$geno[[i]]$data[!is.na(g) & g==1,j] <- 2
cross$geno[[i]]$data[!is.na(g) & g==2,j] <- 1
}
}
}
}
else if(type=="f2") {
if(switch != "AB")
warning("Using switch = \"AB\".")
found <- rep(FALSE, length(markers))
for(i in 1:nchr(cross)) {
cn <- colnames(cross$geno[[i]]$data)
m <- match(markers, cn)
if(any(!is.na(m))) {
found[!is.na(m)] <- TRUE
for(j in m[!is.na(m)]) {
g <- cross$geno[[i]]$data[,j]
cross$geno[[i]]$data[!is.na(g) & g==1,j] <- 3
cross$geno[[i]]$data[!is.na(g) & g==3,j] <- 1
cross$geno[[i]]$data[!is.na(g) & g==4,j] <- 5
cross$geno[[i]]$data[!is.na(g) & g==5,j] <- 4
}
}
}
}
else if(type=="4way") {
if(switch=="AB")
newg <- c(2,1,4,3,6,5,7,8,10,9,12,11,14,13)
else if(switch=="CD")
newg <- c(3,4,1,2,5,6,8,7,10,9,13,14,11,12)
else if(switch=="ABCD")
newg <- c(4,3,2,1,6,5,8,7,9,10,14,13,12,11)
else # switch parents
newg <- c(1,3,2,4,7,8,5,6,9,10,11,13,12,14)
found <- rep(FALSE, length(markers))
for(i in 1:nchr(cross)) {
cn <- colnames(cross$geno[[i]]$data)
m <- match(markers, cn)
if(any(!is.na(m))) {
found[!is.na(m)] <- TRUE
for(j in m[!is.na(m)]) {
g <- cross$geno[[i]]$data[,j]
for(k in 1:14)
cross$geno[[i]]$data[!is.na(g) & g==k,j] <- newg[k]
}
}
}
}
if(any(!found))
warning("Some markers not found: ", paste(markers[!found], collapse=" "))
clean(cross)
}
######################################################################
#
# nqrank: Convert a set of quantitative values to the corresponding
# normal quantiles (preserving mean and SD)
#
######################################################################
nqrank <-
function(x, jitter=FALSE)
{
y <- x[!is.na(x)]
themean <- mean(y, na.rm=TRUE)
thesd <- sd(y, na.rm=TRUE)
y[y == Inf] <- max(y[y<Inf])+10
y[y == -Inf] <- min(y[y > -Inf]) - 10
if(jitter)
y <- rank(y+runif(length(y))/(sd(y)*10^8))
else y <- rank(y)
x[!is.na(x)] <- qnorm((y-0.5)/length(y))
x*thesd/sd(x, na.rm=TRUE)-mean(x,na.rm=TRUE)+themean
}
######################################################################
#
# cleanGeno: omit genotypes that are possibly in error, as indicated
# by apparent double-crossovers separated by a distance of
# no more than maxdist and having no more than maxmark
# interior typed markers
#
######################################################################
cleanGeno <-
function(cross, chr, maxdist=2.5, maxmark=2, verbose=TRUE)
{
if(class(cross)[1] != "bc")
stop("This function currently only works for a backcross.")
if(!missing(chr)) cleaned <- subset(cross, chr=chr)
else cleaned <- cross
thechr <- names(cleaned$geno)
totdrop <- 0
maxmaxdist <- max(maxdist)
for(i in thechr) {
xoloc <- locateXO(cleaned, chr=i, full.info=TRUE)
nxo <- sapply(xoloc, function(a) if(is.matrix(a)) return(nrow(a)) else return(0))
g <- pull.geno(cleaned, chr=i)
ndrop <- 0
for(j in which(nxo > 1)) {
maxd <- xoloc[[j]][-1,"right"] - xoloc[[j]][-nrow(xoloc[[j]]),"left"]
wh <- maxd <= maxmaxdist
if(any(wh)) {
for(k in which(wh)) {
nt <- sum(!is.na(g[j,(xoloc[[j]][k,"ileft"]+1):(xoloc[[j]][k+1,"iright"]-1)]))
if(nt > 0 && any(nt <= maxmark & maxd[k] < maxdist)) {
cleaned$geno[[i]]$data[j,(xoloc[[j]][k,"ileft"]+1):(xoloc[[j]][k+1,"iright"]-1)] <- NA
ndrop <- ndrop + nt
totdrop <- totdrop + nt
}
}
}
}
if(verbose && ndrop > 0) {
totgen <- sum(ntyped(subset(cross, chr=i)))
cat(" ---Dropping ", ndrop, " genotypes (out of ", totgen, ") on chr ", i, "\n", sep="")
}
}
if(verbose && nchr(cleaned)>1 && totdrop > 0) {
totgen <- sum(ntyped(subset(cross, chr=thechr)))
cat(" ---Dropped ", totdrop, " genotypes (out of ", totgen, ") in total\n", sep="")
}
for(i in names(cleaned$geno))
cross$geno[[i]] <- cleaned$geno[[i]]
cross
}
######################################################################
# typingGap: calculate gaps between typed markers
######################################################################
typingGap <-
function(cross, chr, terminal=FALSE)
{
if(!missing(chr))
cross <- subset(cross, chr)
n.ind <- nind(cross)
n.chr <- nchr(cross)
gaps <- matrix(nrow=n.ind, ncol=n.chr)
colnames(gaps) <- names(cross$geno)
for(i in 1:n.chr) {
map <- cross$geno[[i]]$map
map <- c(map[1], map, map[length(map)])
if(is.matrix(map)) stop("This function can't currently handle sex-specific maps.")
if(terminal) # just look at terminal gaps
gaps[,i] <- apply(cbind(1,cross$geno[[i]]$data,1), 1,
function(a,b) {d <- diff(b[!is.na(a)]); max(d[c(1,length(d))]) }, map)
else
gaps[,i] <- apply(cbind(1,cross$geno[[i]]$data,1), 1,
function(a,b) max(diff(b[!is.na(a)])), map)
}
if(n.chr==1) gaps <- as.numeric(gaps)
gaps
}
######################################################################
# calcPermPval
#
# calculate permutation pvalues for summary.scanone()
######################################################################
calcPermPval <-
function(peaks, perms)
{
if(!is.matrix(peaks))
peaks <- as.matrix(peaks)
if(!is.matrix(perms))
perms <- as.matrix(perms)
ncol.peaks <- ncol(peaks)
nrow.peaks <- nrow(peaks)
n.perms <- nrow(perms)
if(ncol.peaks != ncol(perms))
stop("ncol(peaks) != ncol(perms)")
pval <- .C("R_calcPermPval",
as.double(peaks),
as.integer(ncol.peaks),
as.integer(nrow.peaks),
as.double(perms),
as.integer(n.perms),
pval=as.double(rep(0,ncol.peaks*nrow.peaks)),
PACKAGE="qtl")$pval
matrix(pval, ncol=ncol.peaks, nrow=nrow.peaks)
}
######################################################################
# phenames: pull out phenotype names
######################################################################
phenames <-
function(cross)
colnames(cross$pheno)
# end of util.R
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.