#' @title stepwiseqtlShape
#' @description Runs a stepwise multiple QTL model search using penalized negative $log_{10}$ of the p-value
#' @details The function built multiple QTL models using penalized negative $log_{10}$ of the p-value (as a multivariate analogue of the penalized LOD score of the R/qtl package). The search is limited so far to a model without interaction
#' @param cross
#' @param chr
#' @param pheno.col
#' @param qtl \code{\link[qtl]{summary.scanone}} object or any data.frame with 'chr' 'pos' and
#' 'lod' values. If missing a one-dimensional scan is run first.
#' @param formula
#' @param max.qtl
#' @param covar
#' @param refine.locations Logical
#' @param additive.only Logical. (Defauft FALSE). Setting it to TRUE may be used
#' to force a diplotype model (pure additive) in an intercross (Have to bare in mind
#' that it is *NOT* as in \code{\link[qtl]{summary.scanone}}).
#' @param scan.pairs Logical. (Default FALSE) #keep it default. Not use so far
#' @param penalties Genomewide threshold
#' @param disthresh numeric (Default 0). Minimal distance, between loci already in the model
#' and a candidate locus, below which the candidate locus is discarded.
#' @param keeplodprofile Logical (Default TRUE)
#' @param keeptrace Logical (Default FALSE)
#' @param verbose Logical (Default FALSE)
#' @param test Multivariate statistics (Default "Pillai". Others are : "Hotelling.Lawley",
#' "Lik.ratio", "Goodall")
#' @references Broman and Sen 2009. A guide to QTL mapping with R/qtl
#' @export
#' @seealso \code{\link[qtl]{stepwiseqtl}}
#' @keywords analysis
#' @author Nicolas Navarro
#' @return The function returns a qtl object.
#' @examples
#' data(fake.bc)
#' fake.bc <- calc.genoprob(fake.bc, step=2.5)
#' fake.bc <- update(fake.bc, phen2keep = names(fake.bc$pheno)) #Just update the class of fake.bc
#' covar <- fake.bc$pheno[, c('sex','age')]
#' out1 <- scanoneShape(fake.bc, pheno.col = 1:2, addcovar = covar,
#' test = "Pillai")
#' Q <- max(out1)
#' outStep <- stepwiseqtlShape(fake.bc, pheno.col = 1:2, qtl = Q, max.qtl = 3, covar = covar, penalties=2.5, verbose = TRUE)
stepwiseqtlShape <- function(cross, chr, pheno.col = 1, qtl, formula, max.qtl = 10,
covar = NULL,
refine.locations = TRUE,
additive.only = FALSE,
scan.pairs = FALSE, #keep default
penalties,
disthresh = 0,
keeplodprofile = TRUE,
keeptrace = FALSE,
verbose = FALSE,
test = "Pillai") {
# This is the model search version of Broman (p.250)
pheno <- as.matrix(cross$pheno[, pheno.col])
if (missing(chr))
chr <- names(cross$geno)
if (missing(penalties))
stop("you should run first 'scanoneShape' with permutations to get penalties")
if (!is.null(covar) & is.matrix(covar))
covar <- as.data.frame(covar)
if (missing(formula)) {
if (!is.null(covar)){
formula <- as.formula(paste("pheno ~", paste(names(covar), collapse = " + ")))
warning(paste("reduced model formula created from covar names:", deparse(formula)))
} else formula <- as.formula("pheno ~ 1")
} else if (is.character(formula)) {
formula <- as.formula(formula)
}
# We need to define covar anyway because stepwise is based on formula.
if (is.null(covar)) covar <- cross$pheno
if (verbose) {
cat("Null model used is :", deparse(formula), "\n")
}
# Null model (we don't need it except fm.red and except in leave1)
fm.red <- as.formula(paste("pheno", paste(deparse(formula[-2]), collapse = " + ")))
mod.red <- lm(fm.red, data = covar)
SSCPerr.red <- crossprod(mod.red$residuals)
rank.S <- qr(SSCPerr.red)$rank
# Initial model
if (missing(qtl)) {
gscan <- scanoneShape(cross, chr = chr, pheno.col = pheno.col,
addcovar = covar, test = test, formula = fm.red)
qtl <- summary(gscan)
} else if (class(qtl)[1] != "summary.scanone") {
if (all(c("chr", "pos", "lod")%in%names(qtl))) {
qtl <- data.frame(qtl$chr, qtl$pos, qtl$lod)
names(qtl) <- c("chr", "pos", "lod")
} else {
stop("qtl object must have chr, pos as variables")
}
}
#-------------------------------------------------------
# 0. remove unnessecary lodcolumn and do some checkings on the adequation between
# additive.only arg and the type of lod scores we want do the selection on
if (additive.only & "bc" %in% class(cross))
additive.only <- FALSE
if (ncol(qtl) > 3) {
qtl <- qtl[, c(1,2, 3 + additive.only), drop = FALSE]
# so either full model if add.only is FALSE or lod.add if TRUE
}
colnames(qtl)[3] <- "lod"
# 1. Select the position with largest LOD score from a single-QTL genome scan
qtl <- qtl[which.max(qtl$lod), , drop = FALSE]
n.qtls <- 1
curplod <- qtl$lod - penalties * n.qtls
curbestplod <- curplod
curbest <- qtl
if (keeptrace) {
temp <- list(chr = qtl$chr, pos = qtl$pos)
#attr(temp, "formula") <- deparseQTLformula(formula)
attr(temp, "pLOD") <- curplod
class(temp) <- c("compactqtl", "list")
temp <- list(temp)
names(temp) <- 0
thetrace <- c("0", temp)
}
#-------------------------------------------------------
# 2. Forward model search
iter <- 0
cat("Forward selection\n")
while (n.qtls < max.qtl){
iter <- iter + 1
back.qtl <- as.matrix(getGeno(cross, Q = qtl, add.only = additive.only))
n.qtls <- n.qtls + 1
if (verbose)
cat("add qtl:", iter + 1, "\n")
# 2.1 penalized genome scan
# mvGenomScan update automatically fm.red + back.qtl if back.qtl is not null
# but the LOD is not the one of the qtl model vs null (covar only)
pLod <- mvGenomScan(cross, pheno, mod.red = fm.red, covar = covar,
back.qtl = back.qtl, test = test, chr = chr,
updateRedFormula = FALSE)
# Penalization of Broman and Sen 2009, p.251
pLod[, 3 + additive.only] <- pLod[, 3 + additive.only] - penalties * n.qtls
# Set to pLod to NA for markers closely linked to QTLs already in the model
for (q in 1:(n.qtls-1)) {
pLod[as.character(pLod$chr) == as.character(qtl$chr[q]) &
abs(pLod$pos - qtl$pos[q]) <= disthresh, 3 + additive.only] <- NA
}
# Get R/qtl max over genome (handle NA)
tmp <- max(pLod, lodcolumn = 1 + additive.only)[c(1:2, 3 + additive.only)]
colnames(tmp)[3] <- "lod" #ensure colnames consistency after refining positions for rbind()
# if empty, means no position higher than 0
# penalty drops LOD scores need updating max.qtl to n.qtls - 1
if(nrow(tmp) == 0) {
cat("\n !!!!!!\n
All pLOD negatives, max.qtl updated to ", n.qtls - 1,
"\n !!!!!!\n")
max.qtl <- n.qtls - 1
n.qtls <- n.qtls - 1
break
}
qtl <- rbind(qtl, tmp)
curplod <- tmp$lod
if (verbose) {
cat("new qtl:\n")
print(qtl[,c(1,2)])
cat("pLOD = ", curplod, "\n")
}
# 2.2 Reorder and refine qtl positions
# Re-order qtls
qtl <- qtl[order(match(qtl$chr, chr), qtl$pos), ]
if (refine.locations) {
if (verbose) cat("--- Refining positions\n")
# 2.2.1 refining qtl positions
rqtls <- refine.qtl(qtl, cross = cross, fm.red = fm.red, pheno = pheno,
covar = covar, chr = NULL,
max.iter = 10, add.only = additive.only, verbose = FALSE)
if (any(rqtls$pos != qtl$pos)) { # updated positions
if (verbose) cat(" --- Moved a bit\n")
qtl <- data.frame(rqtls$chr, rqtls$pos, rqtls$lod) #row.names=rqtls$name
colnames(qtl) <- c("chr", "pos", "lod")
if(verbose) {
cat(" --- New positions\n")
print(qtl[,c(1,2)])
}
# Updating the penalized LOD score of the model
geno <- as.matrix(getGeno(cross, Q = qtl, add.only = additive.only))
tmp.form <- as.formula(paste(paste(deparse(fm.red), collapse = ""), "+ qtl"))
tmp <- lm.shape.test(geno, pheno = pheno, covar = covar, fm.full = tmp.form,
SSCPerr.red = SSCPerr.red, mod.red.rank = mod.red$rank,
rank.E = rank.S, test = test)
curplod <- tmp - penalties * n.qtls
if (verbose) cat("new pLOD = ", curplod, "\n")
}
if (verbose) cat("--- Refining done\n")
}
attr(qtl, "pLOD") <- curplod
if(curplod > curbestplod) {
if(verbose){
cat("** new best ** (pLOD increased by ",
round(curplod - curbestplod, 4), ")\n", sep = "")
}
curbestplod <- curplod
curbest <- qtl
}
if (keeptrace) {
temp <- list(chr = qtl$chr, pos = qtl$pos)
#attr(temp, "formula") <- deparseQTLformula(formula)
attr(temp, "pLOD") <- curplod
class(temp) <- c("compactqtl", "list")
temp <- list(temp)
names(temp) <- iter
thetrace <- c(thetrace, temp)
}
if(n.qtls >= max.qtl) break
}
#-------------------------------------------------------
# 3. Backward model search
cat("Backward elimination\n")
while (n.qtls > 1){
iter <- iter + 1
if (verbose) cat("remove qtl:",n.qtls,"\n")
n.qtls <- n.qtls - 1
tmp <- leave1qtl(cross, qtls = qtl, pheno = pheno, covar = covar,
formula.red = fm.red, mod.red.rank = mod.red$rank,
SSCPerr.red = SSCPerr.red, add.only = additive.only, test = test)
# look for best model (the one where the removed qtl has minimal effect)
tmp$lod <- tmp$lod - penalties * n.qtls
dropQ <- which.max(tmp$lod)
curplod <- max(tmp$lod)
qtl <- qtl[-dropQ, ]
#-------------------------------------------------------
# 3.1 Refining qtl positions
# Re-order qtl (should not make any difference; they should be already ordered)
qtl <- qtl[order(match(qtl$chr, chr), qtl$pos), ]
if(verbose) {
cat("-- QTL model\n")
print(qtl[,c(1,2)])
cat("pLOD = ", curplod, "\n")
}
if (refine.locations & n.qtls > 1) {
if (verbose) cat("--- Refining positions\n")
rqtls <- refine.qtl(qtl, cross = cross, fm.red = fm.red,
pheno = pheno, covar = covar, chr = NULL, max.iter = 10)
if (any(rqtls$pos != qtl$pos)) { # updated positions
if (verbose) cat(" --- Moved a bit\n")
qtl <- data.frame(rqtls$chr, rqtls$pos, rqtls$lod)
colnames(qtl) <- c("chr", "pos", "lod")
if(verbose) {
cat(" --- New positions\n")
print(qtl)
}
#-------------------------------------------------------
# Updating the penalized LOD score of the model
geno <- as.matrix(getGeno(cross, Q = qtl, add.only = additive.only))
tmp.form <- as.formula(paste(paste(deparse(fm.red), collapse =""),"+ qtl"))
tmp <- lm.shape.test(geno, pheno = pheno, covar = covar, fm.full = tmp.form,
SSCPerr.red = SSCPerr.red, mod.red.rank = mod.red$rank,
rank.E = rank.S, test = test)
curplod <- tmp - penalties * n.qtls
}
}
attr(qtl, "pLOD") <- curplod
if(curplod > curbestplod) {
if(verbose){
cat("** new best ** (pLOD increased by ", round(curplod - curbestplod, 4),
")\n", sep="")
}
curbestplod <- curplod
curbest <- qtl
}
if (keeptrace) {
temp <- list(chr = qtl$chr, pos = qtl$pos)
#attr(temp, "formula") <- deparseQTLformula(formula)
attr(temp, "pLOD") <- curplod
class(temp) <- c("compactqtl", "list")
temp <- list(temp)
names(temp) <- iter
thetrace <- c(thetrace, temp)
}
}
#-------------------------------------------------------
# 4. drop1qtl the current best model to make sure lod column is SS II partial lod
qq <- drop1qtl(cross, qtls = curbest, formula.red = fm.red, pheno = pheno,
covar = covar, threshold=0, add.only = additive.only, test = test)
curbest$lod <- qq$partial.logp
rownames(curbest) <- find.pseudomarker(cross, curbest$chr, curbest$pos, where="prob")
class(curbest) <- c("summary.stepwiseqtl", "data.frame")
#-------------------------------------------------------
# 5. Refining positions
# If so then output is a qtl object instead of a summary
if (refine.locations) {
curbest <- lodprofile.qtl(curbest, cross = cross, fm.red = fm.red, pheno = pheno,
covar = covar, chr = NULL, add.only = additive.only,
test = test)
if (!keeplodprofile){
# return only summary
curbest <- data.frame(chr=curbest$chr,pos=curbest$pos, lod=curbest$lod, row.names = curbest$name)
class(curbest) <- c("summary.stepwiseqtl","data.frame")
attr(curbest, which="pLOD") <- curbestplod
}
}
if(keeptrace)
attr(curbest, "trace") <- thetrace
if (verbose) {
cat(" -- final best QTL model\n")
print(curbest)
}
return(curbest)
}
leave1qtl <- function(cross, qtls, pheno, covar, formula.red, mod.red.rank,
SSCPerr.red, add.only = FALSE, test = "Pillai", verbose = FALSE) {
# Test a model with one-QTL out versus the null model
n.qtl <- nrow(qtls)
n.ind <- nrow(pheno)
fm.red <- paste(deparse(formula.red[-2]), collapse = "")
#get genotypes for all qtls
geno <- getGeno(cross, Q = qtls, add.only = add.only)
#fit the full model
gen <- colnames(geno)
tmp.dat <- cbind(covar, geno)
#fit the full minus 1 QTL
lod <- rep(0, n.qtl)
for (q in 1:n.qtl){
idx <- grep(paste("q", q, ".", sep = ""), gen, fixed=TRUE)
fm.new <- paste(fm.red, paste(gen[-idx], collapse = "+"), sep="+")
tmp.x <- model.matrix(as.formula(fm.new), data = tmp.dat)
mod.full <- .Call(C_CdqrlsShapeQTL, x = tmp.x, y = pheno, tol = 1e-07, chk = FALSE)
dfeff <- mod.full$rank - mod.red.rank
dferr <- n.ind - mod.full$rank
SSCPerr.full <- crossprod(mod.full$residuals)
rank.E <- qr(SSCPerr.full)$rank
#partial F-test: Full model vs Reduced model
SSCPfull <- SSCPerr.red - SSCPerr.full
if (pmatch(test,"Pillai", nomatch=0))
lod[q] <- Pillai.test(SSCPfull,SSCPerr.full,dfeff,dferr,rank.E)
if (pmatch(test,"Lik.ratio", nomatch=0))
lod[q] <- LikelihoodRatio.test(SSCPerr.full,SSCPerr.red,dfeff,dferr,rank.E)
if (pmatch(test,"Hotelling.Lawley", nomatch=0))
lod[q] <- Hotelling.test(SSCPfull,SSCPerr.full,dfeff,dferr,rank.E)
if (pmatch(test,"GoodallF", nomatch=0))
lod[q] <- goodallF.test(diag(SSCPerr.full), diag(SSCPerr.red), dferr,
n.ind - mod.red.rank, rank.E)
}
qtls <- data.frame(qtls$chr, qtls$pos, lod)
colnames(qtls) <- c("chr", "pos", "lod")
return(qtls)
}
drop1qtl <- function(cross, qtls, formula.red, pheno, covar, threshold,
add.only=FALSE, test="Pillai", verbose=FALSE) {
cl.qtl <- class(qtls)
n.qtl <- nrow(qtls)
if (n.qtl < 2) {
warning("Less than 2 QTLs in your model! no drop require")
partial.logp <- qtls$lod
names(partial.logp) <- rownames(qtls)
} else {
n.ind <- nrow(pheno)
fm.red <- paste(deparse(formula.red[-2]), collapse = "")
# get genotypes for all qtls
geno <- getGeno(cross, Q = qtls, add.only = add.only)
# fit the full model
gen <- colnames(geno)
form <- as.formula(paste(fm.red, paste(gen, collapse = "+"), sep="+"))
x <- model.matrix(form, data = cbind(covar,geno))
mod.full <- .Call(C_CdqrlsShapeQTL, x = x, y = pheno, tol = 1e-07, chk = FALSE)
SSCPerr.full <- crossprod(mod.full$residuals)
rank.E <- qr(SSCPerr.full)$rank
# fit the full minus 1 QTL to get SSII p-value
partial.logp <- rep(0, n.qtl)
for (q in 1:n.qtl){
idx <- grep(paste("q",q,".",sep=""), gen, fixed=TRUE)
tmp.x <- model.matrix(as.formula(paste(fm.red, paste(gen[-idx], collapse="+"), sep="+")), data=cbind(covar,geno))
mod.red <- .Call(C_CdqrlsShapeQTL, x = tmp.x, y = pheno, tol = 1e-07, chk = FALSE)
dfeff <- mod.full$rank - mod.red$rank
dferr <- n.ind - mod.full$rank
SSCPerr.red <- crossprod(mod.red$residuals)
#partial F-test: Full model vs Reduced model
SSCPfull <- SSCPerr.red - SSCPerr.full
if (pmatch(test,"Pillai",nomatch=0))
partial.logp[q] <- Pillai.test(SSCPfull,SSCPerr.full,dfeff,dferr,rank.E)
if (pmatch(test,"Lik.ratio",nomatch=0))
partial.logp[q] <- LikelihoodRatio.test(SSCPerr.full,SSCPerr.red,dfeff,dferr,rank.E)
if (pmatch(test,"Hotelling.Lawley",nomatch=0))
partial.logp[q] <- Hotelling.test(SSCPfull,SSCPerr.full,dfeff,dferr,rank.E)
if (pmatch(test,"GoodallF",nomatch=0))
partial.logp[q] <- goodallF.test(diag(SSCPerr.full), diag(SSCPerr.red), dferr,
n.ind - mod.red$rank, rank.E)
}
}
names(partial.logp) <- rownames(qtls)
if (any(partial.logp < threshold)) {
if (verbose) {
print("removed qtl")
print(qtls[partial.logp < threshold, ])
}
qtls <- qtls[partial.logp > threshold, ]
cat("some qtls have been drop!\n updating...\n")
qtls <- drop1qtl(cross, qtls = qtls, formula.red = formula.red,
pheno = pheno, covar = covar,threshold = threshold,
add.only = add.only, test = test)
} else qtls <- cbind(qtls, partial.logp)
class(qtls) <- c("summary.drop1qtl","data.frame")
return(qtls)
}
lodprofile.qtl<- function (qtls, cross, fm.red, pheno, covar, chr = NULL, add.only = FALSE,
test = "Pillai", verbose = FALSE)
{
if(any(class(qtls)%in%"summary.stepwiseqtl")){
rownames(qtls) <- qtls$name
qtls <- qtls[,c("chr","pos","lod")] #suppress name and n.gen to match summary.scanone
}
else if(!any(class(qtls)%in%"summary.scanone"))
stop("qtls must be either summary.scanone object or summary.stepwiseqtl object")
fm.red <- as.formula(paste("pheno", paste(deparse(fm.red[-2]), collapse = "")))
mod.red <- lm(fm.red, data = covar)
if (!is.null(chr)) {
if (length(chr) > 1) {
chr <- chr[1]
cat("Analysis of chr", chr, "only\n")
}
back.qtl <- qtls[qtls$chr != chr, ]
qtls <- qtls[qtls$chr == chr, ]
if (!length(back.qtl))
break
if (!length(qtls))
stop("no qtl for this chromosome")
back.geno <- getGeno(cross, Q = back.qtl, add.only = add.only)
covar <- cbind(covar, back.geno)
# Update the formula of the reduced model
fm.red <- as.formula(paste(paste(deparse(fm.red), collapse = ""),
paste(names(back.geno), collapse = "+"), sep = "+"))
}
n.qtls <- nrow(qtls)
lastOut <- NULL
newpos <- qtls
if (n.qtls == 1) {
tmp.cross <- subset(cross, chr = qtls$chr)
tmp <- mvGenomScan(tmp.cross, pheno = pheno, mod.red = fm.red,
covar = covar, back.qtl = NULL, test = test)
if (ncol(tmp) > 3)
tmp <- tmp[, 1:3]
tmp[is.na(tmp[, 3]), 3] <- 0
lastOut[[1]] <- tmp
newpos[1, ] <- max(tmp)
} else {
for (q in 1:n.qtls) {
geno <- as.matrix(getGeno(cross, Q = qtls[-q, ], add.only = add.only))
tmp.cross <- subset(cross, chr = qtls[q, "chr"])
tmp <- mvGenomScan(tmp.cross, pheno = pheno, mod.red = fm.red,
covar = covar, back.qtl = geno, test = test)
if (ncol(tmp) > 3)
tmp <- tmp[, c(1,2, add.only + 3)]
tmp[is.na(tmp[, 3]), 3] <- 0
lastOut[[q]] <- tmp
newpos[q, ] <- max(tmp)
newpos[q, "chr"] <- qtls[q, "chr"] #in case of being overriden
}
}
qtls <- newpos
qtls.name <- find.pseudomarker(cross, qtls$chr, qtls$pos, where = "prob")
QTL <- makeqtl(cross, qtls$chr, qtls$pos, qtls.name, what = "prob")
QTL$lod <- qtls$lod
names(lastOut) <- QTL$name
attr(QTL, "lodprofile") <- lastOut
return(QTL)
}
refine.qtl <- function(qtls, cross, fm.red, pheno, covar,
chr = NULL,
max.iter = 10,
add.only = FALSE,
verbose = FALSE,
test = "Pillai") {
fm.red <- as.formula(paste("pheno", paste(deparse(fm.red[-2]), collapse = "")))
mod.red <- lm(fm.red, data = covar)
# In the process, we don't update mod.red because it corresponds initially to pheno ~ covar
# and the prg if back.qtl update it to pheno ~ covar + back.qtl
run.model <- TRUE
iter <- 0
if (!is.null(chr)) {
if (length(chr) > 1) {
chr <- chr[1]
cat("Analysis of chr", chr, "only\n")
}
back.qtl <- qtls[qtls$chr != chr,]
qtls <- qtls[qtls$chr == chr, ]
if (!length(back.qtl)) break
if (!length(qtls)) stop("no qtl for this chromosome")
back.geno <- getGeno(cross, Q = back.qtl, add.only = add.only)
covar <- cbind(covar, back.geno)
# Update the formula of the reduced model
fm.red <- as.formula(paste(paste(deparse(fm.red), collapse = ""),
paste(names(back.geno), collapse = "+"), sep = "+"))
}
qtls.name <- find.pseudomarker(cross, chr = qtls$chr, pos = qtls$pos, where="prob")
QTL <- makeqtl(cross, chr = qtls$chr, pos = qtls$pos, qtl.name = qtls.name, what="prob")
n.qtls <- QTL$n.qtl
if (is.null(chr) & n.qtls == 1) {
warning("Only one qtl without background qtl: No refining done")
return(QTL)
}
map <- attr(QTL, "map")
if (is.null(map))
stop("Input qtl object should contain the genetic map.")
mind <- min(sapply(map, function(a) {
if (is.matrix(a)) a <- a[1, ]
min(diff(a))
}))/2
if (mind <= 0) mind <- 1e-06
target.qtl <- NULL
while (run.model){
iter <- iter + 1
old.qtls <- qtls
# Randomly pick a qtl. But make sure the 1st was not the first before
o <- sample(n.qtls)
if(!is.null(target.qtl))
while(o[1] == target.qtl[1]) o <- sample(n.qtls)
target.qtl <- o
outStep <- NULL
if(verbose) cat("------ scan ...")
for (q in 1:n.qtls) {
geno <- as.matrix(getGeno(cross, Q = qtls[-target.qtl[q], ], add.only))
tmp.cross <- subset(cross, chr = qtls[target.qtl[q], 'chr'])
tmp <- mvGenomScan(tmp.cross, pheno, mod.red = fm.red, covar = covar,
back.qtl = geno, test = test)
# selection on full or additive model depending on add.only arg
if (ncol(tmp) > 3) {
tmp <- tmp[, c(1,2, 3 + add.only), drop = FALSE]
}
tmp[is.na(tmp[, 3]), 3] <- 0
mx.tmp <- max(tmp) #get R/qtl max over the chromosome
if (verbose) {
print("back qtl")
print(qtls[, 1:3])
print("new qtl")
print(mx.tmp)
}
qtls[target.qtl[q], ] <- mx.tmp
qtls[target.qtl[q], 'chr'] <- as.character(mx.tmp$chr) #this is because of chr is a factor
outStep[[q]] <- tmp
}
if(verbose) cat(" done\n")
#check if positions have changed.
if (max(abs(old.qtls[target.qtl, 'pos']-qtls[target.qtl, 'pos'])) < mind){
converged <- TRUE
break
}
if (iter>=max.iter) break
}
idx <- order(as.numeric(qtls$chr) * max(qtls$pos) + qtls$pos)
qtls <- qtls[idx, ]
lastOut <- NULL
newpos <- qtls
for (q in 1:n.qtls){
geno <- as.matrix(getGeno(cross, qtls[-q, ], add.only))
tmp.cross <- subset(cross, chr = qtls[q, 'chr'])
tmp <- mvGenomScan(tmp.cross, pheno, mod.red = fm.red, covar = covar,
back.qtl = geno, test = test)
if (ncol(tmp) > 3) {
tmp <- tmp[, c(1,2, 3 + add.only), drop = FALSE]
}
tmp[is.na(tmp[, 3]), 3] <- 0
lastOut[[q]] <- tmp
mx.tmp <- max(tmp)
qtls[q, ] <- mx.tmp
qtls[q, 'chr'] <- as.character(mx.tmp$chr) #!! because chr factor
}
qtls.name <- find.pseudomarker(cross, chr = qtls$chr, pos = qtls$pos, where = "prob")
QTL <- makeqtl(cross, chr = qtls$chr, pos = qtls$pos, qtl.name = qtls.name, what = "prob")
# add the current partial lod scores
QTL$lod <- qtls$lod
# add the partial lod score profiles
names(lastOut) <- QTL$name
attr(QTL, "lodprofile") <- lastOut
return(QTL)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.