Nothing
###########################################################################################
scan.genome <- function(cross, pheno.col, pheno.parents, addcov, intcov,
threshold, method, ...)
{
## This currently scans one phenotype at a time.
## Want to extend to all phenotypes not in names(covM.dat).
## To do that, assume covariates are same for all i.
n.pheno <- length(pheno.col)
## *** The following could be simplified.
## *** However, it involves a rethinking of myformula.
## *** Low priority for now.
## Design matrices for qtl::scanone. Must be numeric entries.
## Design matrix for parent phenotypes
covM.dat <- pull.pheno.null(cross, pheno.parents)
## Design matrix for additive covariates.
addcovM.dat <- covM.dat
tmp <- unique(c(addcov, intcov))
if(!is.null(tmp)) {
tmp <- create.cov.matrix(cross, cov.names = tmp)
if(length(tmp))
addcovM.dat <- cbind(addcovM.dat, tmp)
}
## Design matrix for interactive covariates.
intcovM.dat <- create.cov.matrix(cross,cov.names=intcov)
## Phenotype matrices of actual data. Can be mix of numeric and factor.
addcov.dat <- pull.pheno.null(cross, unique(c(addcov, intcov)))
intcov.dat <- pull.pheno.null(cross, intcov)
## Call to scanone is the big time commitment.
ss <- scanone.summary(cross, pheno.col, addcovM.dat, intcovM.dat, threshold, method)
bic <- rep(NA, n.pheno)
cross.type <- class(cross)[1]
if(nrow(ss)) {
## Model may be different for each trait, depending on QTL.
## For loop is clumsy, but calculations are pretty fast.
for(i in seq(n.pheno)) {
## Response for linear model.
y <- cross$pheno[, pheno.col[i]]
signif.lods <- (ss[ , 1 + 2 * i] >= threshold[i])
le.markers <- sum(signif.lods)
## Build data.frame and formula for linear model fit.
if(le.markers > 0){
## Need to extend this to multiple phenotypes.
qtlo <- qtl::makeqtl(cross, chr = ss[signif.lods, 1],
pos = ss[signif.lods, 2 * i], what="prob")
geno.dat <- hk.design.matrix(qtlo=qtlo, cross.type)[,-1, drop = FALSE]
tmp <- paste("add", 1:le.markers, sep = "")
if(cross.type == "f2")
tmp <- as.vector(rbind(tmp, paste("dom", 1:le.markers, sep = "")))
dimnames(geno.dat) <- list(NULL, tmp)
form <- set.dat.form(y, covM.dat, addcov.dat, intcov.dat, geno.dat, le.markers,
cross.type)
dat <- form$dat
form <- form$form
}
else{
form <- set.dat.form(y, covM.dat, addcov.dat, intcov.dat, cross.type)
dat <- form$dat
form <- form$form
}
## Fit linear model.
fm <- stats::lm(form, dat)
## Record BIC.
bic[i] <- stats::AIC(fm, k = log(length(y)))[1]
}
}
else {
## Same model for all traits (no QTL).
for(i in seq(n.pheno)) {
if(i == 1) {
y <- cross$pheno[, pheno.col[1]]
## Need only do this once if no QTL.
form <- set.dat.form(y, covM.dat, addcov.dat, intcov.dat, cross.type)
dat <- form$dat
form <- form$form
## Fit linear model.
fm <- stats::lm(form, dat)
}
else
dat$y <- cross$pheno[, pheno.col[i]]
## Record BIC.
bic[i] <- stats::AIC(stats::update(fm, data = dat), k = log(length(y)))[1]
}
}
bic
}
###########################################################################################
scanone.summary <- function(cross, pheno.col, addcov, intcov, threshold,
method)
{
## This is the big time commitment.
scan <- qtl::scanone(cross, pheno.col = pheno.col,
addcovar = addcov, intcovar = intcov, method = method)
## Mainly intersted in this summary to determine QTLs.
summary(scan, format = ifelse(length(pheno.col) == 1, "onepheno", "allpeaks"),
threshold = threshold)
}
###########################################################################################
set.dat.form <- function(y, covM.dat=NULL, addcov.dat=NULL, intcov.dat=NULL,
geno.dat=NULL, le.markers = 0, cross.type = "f2")
{
## Set up data.frame
dat <- data.frame(y = y, stringsAsFactors = TRUE)
if(!is.null(covM.dat))
dat <- cbind.data.frame(dat, covM.dat, stringsAsFactors = TRUE)
if(!is.null(addcov.dat))
dat <- cbind.data.frame(dat, addcov.dat, stringsAsFactors = TRUE)
if(!is.null(intcov.dat))
dat <- cbind.data.frame(dat, intcov.dat, stringsAsFactors = TRUE)
if(!is.null(geno.dat))
dat <- cbind.data.frame(dat, geno.dat, stringsAsFactors = TRUE)
## Set up formula.
form <- cov.formula(c(names(covM.dat),names(addcov.dat)),
names(intcov.dat),
le.markers, cross.type)
list(dat = dat, form = form)
}
######################################################################
hk.design.matrix <- function(qtlo, cross.type="f2")
{
nr <- nrow(qtlo$prob[[1]])
ng <- length(qtlo$prob)
if(cross.type == "f2"){
tmp <- unlist(lapply(qtlo$prob,
function(x) {
if(ncol(x) == 3)
cbind(x[,1]-x[,3],x[,2])
else ## Must be X chr.
x[,1]-x[,2]
}))
hkm <- matrix(tmp,nr,2*ng)
}
if(cross.type == "bc"){
tmp <- unlist(lapply(qtlo$prob, function(x) x[,1]-x[,2]))
hkm <- matrix(tmp,nr,ng)
}
cbind(rep(1,nr),hkm)
}
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.