######################################################################
#
# 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.