Nothing
######################################################################
#
# makeqtl.R
#
# copyright (c) 2002-2019, Hao Wu and Karl W. Broman
# last modified Dec, 2019
# first written Apr, 2002
#
# Modified by Danny Arends
#
# 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: makeqtl, replaceqtl, addtoqtl, dropfromqtl, locatemarker
# print.qtl, summary.qtl, print.summary.qtl, reorderqtl
# plot.qtl
# print.compactqtl, summary.compactqtl, print.summary.compactqtl
#
######################################################################
######################################################################
#
# This is the function to construct an object of class "qtl"
# The phenotype data and genotype data for a given list of
# chromosome and locations will be extracted from the input
# "cross" object
#
######################################################################
makeqtl <-
function(cross, chr, pos, qtl.name, what=c("draws", "prob"))
{
if( !inherits(cross, "cross") )
stop("The first input variable must be an object of class cross")
# cross type
type <- crosstype(cross)
chr_type <- sapply(cross$geno, chrtype)
names(chr_type) <- names(cross$geno)
sexpgm <- getsex(cross)
what <- match.arg(what)
themap <- pull.map(cross)
# try to interpret chr argument
if(!is.character(chr))
chr <- as.character(chr)
# chr, pos and qtl.name must have the same length
if(length(chr) != length(pos))
stop("Input chr and pos must have the same length.")
else if( !missing(qtl.name) )
if( length(chr) != length(qtl.name) )
stop("Input chr and qtl.name must have the same length.")
# local variables
n.ind <- nrow(cross$pheno) # number of individuals
n.pos <- length(chr) # number of selected markers
n.gen <- NULL
# initialize output object
qtl <- NULL
# take out the imputed genotypes and/or genoprobs for the
# selected markers (if there are there)
if(what == "draws") { # pull out draws
if(!("draws" %in% names(cross$geno[[1]])))
stop("You must first run sim.geno.")
# take out imputed genotype data
n.draws <- dim(cross$geno[[1]]$draws)[3] # number of draws
# initialize geno matrix for selected markers
geno <- array(rep(0, n.ind*n.pos*n.draws),
dim=c(n.ind, n.pos, n.draws))
for(i in 1:n.pos) {
# get the index for this chromosome
i.chr <- which(chr[i]==names(cross$geno))
if(length(i.chr) == 0) # no this chromosome in cross
stop("There's no chromosome number ", chr[i], " in input cross object")
i.pos <- pos[i] # marker position
# make the genetic map for this chromosome
if("map" %in% names(attributes(cross$geno[[i.chr]]$draws)))
map <- attr(cross$geno[[i.chr]]$draws,"map")
else {
stp <- attr(cross$geno[[i.chr]]$draws, "step")
oe <- attr(cross$geno[[i.chr]]$draws, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[i.chr]]$draws)))
stpw <- attr(cross$geno[[i.chr]]$draws, "stepwidth")
else stpw <- "fixed"
map <- create.map(cross$geno[[i.chr]]$map,stp,oe,stpw)
}
# pull out the female map if there are sex-specific maps
if(is.matrix(map)) map <- map[1,]
# locate this marker (given chromosome and position)
marker.idx <- locatemarker(map, i.pos, i.chr, flag="draws")
if(length(marker.idx) > 1)
stop("Multiple markers at the same position; run jittermap.")
# if everything is all right, take the genotype
geno[,i,] <- cross$geno[[i.chr]]$draws[,marker.idx,]
pos[i] <- map[marker.idx]
# no. genotypes
n.gen[i] <- length(getgenonames(type,chr_type[i.chr],"full",sexpgm, attributes(cross)))
# Fix up X chromsome here
if(chr_type[i.chr]=="X" && (type=="bc" || type=="f2"))
geno[,i,] <- reviseXdata(type,"full",sexpgm,draws=geno[,i,,drop=FALSE],
cross.attr=attributes(cross))
}
# give geno dimension names
# the 2nd dimension called "Q1", "Q2", etc.
dimnames(geno) <- list(NULL, paste("Q", 1:n.pos, sep=""), NULL)
# output
qtl$geno <- geno
}
else { # pull out probs
if(!("prob" %in% names(cross$geno[[1]])))
stop("You must first run calc.genoprob.")
# initialize prob matrix
prob <- vector("list",n.pos)
# locate the marker
for(i in 1:n.pos) {
# get the index for this chromosome
i.chr <- which(chr[i]==names(cross$geno))
if(length(i.chr) == 0) # no this chromosome in cross
stop("There's no chromosome number ", chr[i], " in input cross object")
i.pos <- pos[i] # marker position
if("map" %in% names(attributes(cross$geno[[i.chr]]$prob)))
map <- attr(cross$geno[[i.chr]]$prob,"map")
else {
stp <- attr(cross$geno[[i.chr]]$prob, "step")
oe <- attr(cross$geno[[i.chr]]$prob, "off.end")
if("stepwidth" %in% names(attributes(cross$geno[[i.chr]]$prob)))
stpw <- attr(cross$geno[[i.chr]]$prob, "stepwidth")
else stpw <- "fixed"
map <- create.map(cross$geno[[i.chr]]$map,stp,oe,stpw)
}
# pull out the female map if there are sex-specific maps
if(is.matrix(map)) map <- map[1,]
# locate this marker (given chromosome and position)
marker.idx <- locatemarker(map, i.pos, i.chr, flag="prob")
if(length(marker.idx) > 1)
stop("Multiple markers at the same position; run jittermap.")
# take genoprob
if(chr_type[i.chr]=="X" && (type=="bc" || type=="f2")) { # fix X chromosome probs
prob[[i]] <- reviseXdata(type, "full", sexpgm,
prob=cross$geno[[i.chr]]$prob[,marker.idx,,drop=FALSE],
cross.attr=attributes(cross))[,1,]
}
else
prob[[i]] <- cross$geno[[i.chr]]$prob[,marker.idx,]
pos[i] <- map[marker.idx]
# no. genotypes
n.gen[i] <- ncol(prob[[i]])
}
qtl$prob <- prob
}
if(missing(qtl.name)) { # no given qtl names
dig <- 1
if(what=="draws")
step <- attr(cross$geno[[i.chr]]$draws, "step")
else
step <- attr(cross$geno[[i.chr]]$prob, "step")
if(!is.null(step)) {
if(step > 0) dig <- max(dig, -floor(log10(step)))
}
else {
if(what=="draws")
stepw <- attr(cross$geno[[i.chr]]$draws, "stepwidth")
else
stepw <- attr(cross$geno[[i.chr]]$prob, "stepwidth")
if(!is.null(stepw) && stepw > 0) dig <- max(dig, -floor(log10(stepw)))
}
# make qtl names
qtl.name <- paste( paste(chr,sep=""), charround(pos,dig), sep="@")
}
# output object
qtl$name <- qtl.name
qtl$altname <- paste("Q", 1:n.pos, sep="")
qtl$chr <- chr
qtl$pos <- pos
qtl$n.qtl <- n.pos
qtl$n.ind <- nind(cross)
qtl$n.gen <- n.gen
qtl$chrtype <- chr_type[qtl$chr]
names(qtl$chrtype) <- NULL
class(qtl) <- "qtl"
attr(qtl, "map") <- themap
qtl
}
######################################################################
#
# This is the function to replace one QTL by another.
#
######################################################################
replaceqtl <-
function(cross, qtl, index, chr, pos, qtl.name, drop.lod.profile=TRUE)
{
if(!inherits(qtl, "qtl"))
stop("qtl should have class \"qtl\".")
if(any(index < 1 | index > qtl$n.qtl))
stop("index should be between 1 and ", qtl$n.qtl)
if(length(index) != length(chr) || length(index) != length(pos))
stop("index, chr, and pos should all have the same length.")
if(!missing(qtl.name) && length(index) != length(qtl.name))
stop("index and qtl.name should have the same length.")
if("geno" %in% names(qtl)) what <- "draws"
else what <- "prob"
if(missing(qtl.name))
newqtl <- makeqtl(cross, chr, pos, what=what)
else
newqtl <- makeqtl(cross, chr, pos, qtl.name=qtl.name, what=what)
if(what=="draws") {
qtl$geno[,index,] <- newqtl$geno
}
else {
qtl$prob[index] <- newqtl$prob
}
qtl$name[index] <- newqtl$name
qtl$chr[index] <- newqtl$chr
qtl$pos[index] <- newqtl$pos
qtl$chrtype[index] <- newqtl$chrtype
if(qtl$n.ind != newqtl$n.ind) stop("Mismatch in no. individuals")
qtl$n.gen[index] <- newqtl$n.gen
if(drop.lod.profile)
attr(qtl, "lodprofile") <- NULL
qtl
}
######################################################################
#
# This is the function to add a QTL to given qtl object
#
######################################################################
addtoqtl <-
function(cross, qtl, chr, pos, qtl.name, drop.lod.profile=TRUE)
{
if(!inherits(qtl, "qtl"))
stop("qtl should have class \"qtl\".")
if("geno" %in% names(qtl)) what <- "draws"
else what <- "prob"
if(missing(qtl.name))
newqtl <- makeqtl(cross, chr, pos, what=what)
else
newqtl <- makeqtl(cross, chr, pos, qtl.name=qtl.name, what=what)
if(what=="draws") {
do <- dim(qtl$geno)
dn <- dim(newqtl$geno)
if(do[1] != dn[1] || do[3] != dn[3])
stop("Mismatch in number of individuals or number of imputations.")
temp <- array(dim=c(do[1], do[2]+dn[2], do[3]))
temp[,1:ncol(qtl$geno),] <- qtl$geno
temp[,-(1:ncol(qtl$geno)),] <- newqtl$geno
colnames(temp) <- paste("Q", 1:ncol(temp), sep="")
qtl$geno <- temp
}
else {
qtl$prob <- c(qtl$prob, newqtl$prob)
}
qtl$name <- c(qtl$name, newqtl$name)
qtl$chr <- c(qtl$chr, newqtl$chr)
qtl$pos <- c(qtl$pos, newqtl$pos)
qtl$n.qtl <- qtl$n.qtl + newqtl$n.qtl
qtl$altname <- paste("Q", 1:qtl$n.qtl, sep="")
qtl$chrtype <- c(qtl$chrtype, newqtl$chrtype)
if(qtl$n.ind != newqtl$n.ind)
stop("Mismatch in no. individuals")
qtl$n.gen <- c(qtl$n.gen, newqtl$n.gen)
attr(qtl, "formula") <- NULL
attr(qtl, "pLOD") <- NULL
if(drop.lod.profile)
attr(qtl, "lodprofile") <- NULL
qtl
}
######################################################################
#
# This is the function to drop a QTL from a given qtl object
#
######################################################################
dropfromqtl <-
function(qtl, index, chr, pos, qtl.name, drop.lod.profile=TRUE)
{
if(!inherits(qtl, "qtl"))
stop("qtl should have class \"qtl\".")
if(!missing(chr) || !missing(pos)) {
if(missing(chr) || missing(pos))
stop("Give both chr and pos, or give name, or give a numeric index")
if(!missing(qtl.name) || !missing(index))
stop("Give chr and pos or qtl.name or numeric index, but not multiple of these.")
if(length(chr) != length(pos))
stop("chr and pos must have the same lengths.")
todrop <- NULL
for(i in seq(along=chr)) {
m <- which(qtl$chr == chr[i])
if(length(m) < 1)
stop("No QTL on chr ", chr[i], " in input qtl object.")
for(j in seq(along=m)) {
d <- abs(qtl$pos[m[j]] - pos[i])
if(min(d) > 10) stop("No qtl near position ", pos[i], " on chr ", chr[i])
wh <- m[d==min(d)]
if(length(wh) > 1)
stop("Multiple QTL matching chr ", chr[i], " at pos ", pos[i])
if(min(d) > 1)
warning("No QTL on chr ", chr[i], " exactly at ", pos[i],
"; dropping that at ", qtl$pos[wh])
todrop <- c(todrop, wh)
}
}
todrop <- unique(todrop)
}
else if(!missing(qtl.name)) {
if(!missing(index))
stop("Give chr and pos or qtl.name or numeric index, but not multiple of these.")
m <- match(qtl.name, qtl$name)
if(all(is.na(m))) # if no matches, try "altname"
m <- match(qtl.name, qtl$altname)
if(any(is.na(m)))
warning("Didn't match QTL ", qtl.name[is.na(m)])
todrop <- m[!is.na(m)]
}
else {
if(missing(index))
stop("Give chr and pos or qtl.name or numeric index, but not multiple of these.")
if(any(index < 1 | index > qtl$n.qtl))
stop("index should be between 1 and ", qtl$n.qtl)
todrop <- index
}
# input drop is an integer index
# get the index for exclusing drop QTL
idx <- setdiff(1:qtl$n.qtl, todrop)
# result object
qtl$name <- qtl$name[idx]
qtl$chr <- qtl$chr[idx]
qtl$chrtype <- qtl$chrtype[idx]
qtl$pos <- qtl$pos[idx]
qtl$n.qtl <- qtl$n.qtl - length(todrop)
qtl$altname <- paste("Q", 1:qtl$n.qtl, sep="")
qtl$n.ind <- qtl$n.ind
qtl$n.gen <- qtl$n.gen[idx]
if("geno" %in% names(qtl)) {
qtl$geno <- qtl$geno[,idx,,drop=FALSE]
colnames(qtl$geno) <- paste("Q", 1:ncol(qtl$geno), sep="")
}
if("prob" %in% names(qtl))
qtl$prob <- qtl$prob[idx]
attr(qtl, "formula") <- NULL
attr(qtl, "pLOD") <- NULL
if(drop.lod.profile)
attr(qtl, "lodprofile") <- NULL
qtl
}
##################################################################
#
# locate the marker on a genetic map. Choose the nearest
# one if there's no marker or pseudomarker one the given
# location
#
# This is the internal function and not supposed to be used by user
#
###################################################################
locatemarker <-
function(map, pos, chr, flag)
{
marker.idx <- which(map == pos)
if( length(marker.idx)==0 ) {
# there's no this marker, take the nearest marker instead
# if there's a tie, take the first nearst one
m.tmp <- abs(pos-map)
marker.idx <- which(m.tmp==min(m.tmp))[[1]]
}
if(length(marker.idx) > 1)
marker.idx <- marker.idx[sample(length(marker.idx))]
marker.idx
}
# print QTL object
print.qtl <-
function(x, ...)
{
print(summary(x))
}
# summary of QTL object
summary.qtl <-
function(object, ...)
{
if(is.null(object) || length(object) == 0 || length(object$chr)==0) {
object <- numeric(0)
class(object) <- "summary.qtl"
return(object)
}
if("geno" %in% names(object)) {
type <- "draws"
n.draws <- dim(object$geno)[3]
}
else type <- "prob"
output <- data.frame(name=object$name, chr=object$chr, pos=object$pos, n.gen=object$n.gen, stringsAsFactors=TRUE)
rownames(output) <- object$altname
attr(output, "type") <- type
if(!is.null(attr(object,"mqm"))) attr(output, "mqm") <- attr(object,"mqm")
if(type=="draws") attr(output, "n.draws") <- n.draws
class(output) <- c("summary.qtl", "data.frame")
if("formula" %in% names(attributes(object)))
attr(output, "formula") <- attr(object, "formula")
if("pLOD" %in% names(attributes(object)))
attr(output, "pLOD") <- attr(object, "pLOD")
output
}
# print summary of QTL object
print.summary.qtl <-
function(x, ...)
{
if(is.null(x) || length(x) == 0) {
cat(" Null QTL model\n")
}
else {
type <- attr(x, "type")
if(type=="draws")
thetext <- paste("imputed genotypes, with", attr(x, "n.draws"), "imputations.")
else thetext <- "genotype probabilities."
if(!is.null(attr(x,"mqm"))) thetext <- paste("model created by mqmscan")
cat(" QTL object containing", thetext, "\n\n")
print.data.frame(x, digits=5)
}
if("formula" %in% names(attributes(x))) {
form <- attr(x, "formula")
if(!is.character(form)) form <- deparseQTLformula(form)
cat("\n Formula:")
w <- options("width")[[1]]
printQTLformulanicely(form, " ", w+5, w)
}
if("pLOD" %in% names(attributes(x)))
cat("\n pLOD: ", round(attr(x, "pLOD"),3), "\n")
}
######################################################################
# plot locations of QTLs on the genetic map
######################################################################
plot.qtl <-
function(x, chr, horizontal=FALSE, shift=TRUE,
show.marker.names=FALSE, alternate.chrid=FALSE,
justdots=FALSE, col="red", ...)
{
if(!inherits(x, "qtl"))
stop("input should be a qtl object")
if(length(x) == 0)
stop(" There are no QTL to plot.")
map <- attr(x, "map")
if(is.null(map))
stop("qtl object doesn't contain a genetic map.")
if(missing(chr))
chr <- names(map)
else {
chr <- matchchr(chr, names(map))
map <- map[chr]
class(map) <- "map"
}
if(horizontal)
plotMap(map, horizontal=horizontal, shift=shift,
show.marker.names=show.marker.names, alternate.chrid=alternate.chrid,
ylim=c(length(map)+0.5, 0), ...)
else
plotMap(map, horizontal=horizontal, shift=shift,
show.marker.names=show.marker.names, alternate.chrid=alternate.chrid,
xlim=c(0.5,length(map)+1), ...)
whchr <- match(x$chr, names(map))
thepos <- x$pos
thepos[is.na(whchr)] <- NA
if(any(!is.na(thepos))) {
whchr <- whchr[!is.na(whchr)]
if(shift) thepos <- thepos - sapply(map[whchr], min)
if(is.matrix(map[[1]])) whchr <- whchr - 0.3
if(length(grep("^.+@[0-9\\.]+$", x$name)) == length(x$name))
x$name <- x$altname
if(!justdots) {
if(horizontal) {
arrows(thepos, whchr - 0.35, thepos, whchr, lwd=2, col=col, length=0.05)
text(thepos, whchr-0.4, x$name, col=col, adj=c(0.5,0))
}
else {
arrows(whchr + 0.35, thepos, whchr, thepos, lwd=2, col=col, length=0.05)
text(whchr+0.4, thepos, x$name, col=col, adj=c(0,0.5))
}
}
else {
if(horizontal)
points(thepos, whchr, pch=16, col=col)
else {
points(whchr, thepos, pch=16, col=col)
}
}
}
invisible()
}
######################################################################
#
# This is the function to reorder the QTL within a QTL object
#
######################################################################
reorderqtl <-
function(qtl, neworder)
{
if(!inherits(qtl, "qtl"))
stop("qtl should have class \"qtl\".")
if(missing(neworder)) {
if(!("map" %in% names(attributes(qtl))))
stop("No map in the qtl object; you must provide argument 'neworder'.")
chr <- names(attr(qtl, "map"))
thechr <- match(qtl$chr, chr)
if(any(is.na(thechr)))
stop("Chr ", paste(qtl$chr[is.na(thechr)], " "), " not found.")
neworder <- order(thechr, qtl$pos)
}
curorder <- seq(qtl$n.qtl)
if(length(neworder) != qtl$n.qtl ||
!all(curorder == sort(neworder)))
stop("neworder should be an ordering of the integers from 1 to ", qtl$n.qtl)
if(qtl$n.qtl == 1)
stop("Nothing to do; just one qtl.")
if("geno" %in% names(qtl))
qtl$geno <- qtl$geno[,neworder,]
else
qtl$prob <- qtl$prob[neworder]
qtl$name <- qtl$name[neworder]
qtl$chr <- qtl$chr[neworder]
qtl$pos <- qtl$pos[neworder]
qtl$n.gen <- qtl$n.gen[neworder]
qtl$chrtype <- qtl$chrtype[neworder]
attr(qtl, "formula") <- NULL
attr(qtl, "pLOD") <- NULL
if("lodprofile" %in% names(attributes(qtl))) {
lodprof <- attr(qtl, "lodprofile")
if(length(lodprof) == length(neworder))
attr(qtl, "lodprofile") <- lodprof[neworder]
}
qtl
}
# print compact version of QTL object
print.compactqtl <-
function(x, ...)
{
print(summary(x))
}
summary.compactqtl <-
function(object, ...)
{
class(object) <- c("summary.compactqtl", "list")
object
}
print.summary.compactqtl <-
function(x, ...)
{
if(is.null(x) || length(x) == 0)
cat("Null QTL model\n")
else {
temp <- as.data.frame(x, stringsAsFactors=TRUE)
rownames(temp) <- paste("Q", 1:nrow(temp), sep="")
print.data.frame(temp)
}
if("formula" %in% names(attributes(x))) {
form <- attr(x, "formula")
if(!is.character(form)) form <- deparseQTLformula(form)
cat(" Formula:")
w <- options("width")[[1]]
printQTLformulanicely(form, " ", w+5, w)
}
if("pLOD" %in% names(attributes(x)))
cat(" pLOD: ", round(attr(x, "pLOD"),3), "\n")
}
# nqtl: number of qtl
nqtl <- function(qtl) ifelse(length(qtl)==0,0,length(qtl$chr))
# end of makeqtl.R
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.