Nothing
#####################################################################
#
# mqmscanall.R
#
# Copyright (c) 2009-2019, Danny Arends
#
# Modified by Pjotr Prins and Karl Broman
#
# first written Februari 2009
# last modified Dec 2019
#
# 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: mqmscanall
# scanall
#
#
#####################################################################
mqmscanall <- function(cross, multicore=TRUE, n.clusters=1, batchsize=10,cofactors=NULL, ...) {
scanall(cross=cross, multicore=multicore, n.clusters=n.clusters, batchsize=batchsize,cofactors=cofactors, ..., scanfunction=mqmscan)
}
scanall <- function(cross, scanfunction=scanone, multicore=TRUE, n.clusters=1, batchsize=10, FF=0,cofactors=NULL, ..., plot=FALSE, verbose=FALSE)
{
if(missing(cross)){
ourstop("No cross file. Please supply a valid cross object.")
}
crosstype <- crosstype(cross)
if(!(crosstype == "f2" || crosstype == "bc" || crosstype == "riself"))
stop("Currently only F2, BC, and selfed RIL crosses can be analyzed by MQM.")
cross <- omit_x_chr(cross)
start <- proc.time()
n.pheno <- nphe(cross)
if(verbose) {
ourline()
cat("Starting R/QTL multitrait analysis\n")
cat("Number of phenotypes:",n.pheno,"\n")
cat("Batchsize:",batchsize," & n.clusters:",n.clusters,"\n")
ourline()
}
result <- NULL #BATCH result variable
res <- NULL #GLOBAL result variable
all.data <- cross
bootstraps <- 1:n.pheno
batches <- length(bootstraps) %/% batchsize
last.batch.num <- length(bootstraps) %% batchsize
if(last.batch.num > 0){
batches = batches+1
}
#INIT TIME VARS
SUM <- 0
AVG <- 0
LEFT <- 0
#TEST FOR SNOW CAPABILITIES
if(multicore && n.clusters >1) {
updateParallelRNG(n.clusters)
if(verbose) cat("INFO: Using ",n.clusters," Cores/CPU's/PC's for calculation.\n")
for(x in 1:(batches)){
start <- proc.time()
if(verbose) cat("INFO: Starting with batch",x,"/",batches,"\n")
if(x==batches && last.batch.num > 0){
boots <- bootstraps[((batchsize*(x-1))+1):((batchsize*(x-1))+last.batch.num)]
}else{
boots <- bootstraps[((batchsize*(x-1))+1):(batchsize*(x-1)+batchsize)]
}
if(Sys.info()[1] == "Windows") { # Windows doesn't support mclapply, but it's faster if available
cl <- makeCluster(n.clusters)
on.exit(stopCluster(cl))
result <- clusterApply(cl, boots, snowCoreALL, all.data=all.data, scanfunction=scanfunction, cofactors=cofactors,
verbose=verbose, ...)
}
else {
result <- mclapply(boots, snowCoreALL, all.data=all.data, scanfunction=scanfunction, cofactors=cofactors,
verbose=verbose, mc.cores=n.clusters, ...)
}
if(plot){
temp <- result
class(temp) <- c(class(temp),"mqmmulti")
mqmplot.multitrait(temp)
}
res <- c(res,result)
end <- proc.time()
SUM <- SUM + (end-start)[3]
AVG <- SUM/x
LEFT <- AVG*(batches-x)
if(verbose) {
cat("INFO: Done with batch",x,"/",batches,"\n")
cat("INFO: Calculation of batch",x,"took:",round((end-start)[3], digits=3),"seconds\n")
cat("INFO: Elapsed time:",(SUM%/%3600),":",(SUM%%3600)%/%60,":",round(SUM%%60, digits=0),"(Hour:Min:Sec)\n")
cat("INFO: Average time per batch:",round((AVG), digits=3)," per trait:",round((AVG %/% batchsize), digits=3),"seconds\n")
cat("INFO: Estimated time left:",LEFT%/%3600,":",(LEFT%%3600)%/%60,":",round(LEFT%%60,digits=0),"(Hour:Min:Sec)\n")
ourline()
}
}
}else{
if(verbose) cat("INFO: Going into singlemode.\n")
for(x in 1:(batches)){
start <- proc.time()
if(verbose) cat("INFO: Starting with batch",x,"/",batches,"\n")
if(x==batches && last.batch.num > 0){
boots <- bootstraps[((batchsize*(x-1))+1):((batchsize*(x-1))+last.batch.num)]
}else{
boots <- bootstraps[((batchsize*(x-1))+1):(batchsize*(x-1)+batchsize)]
}
result <- lapply(boots, FUN=snowCoreALL,all.data=all.data,scanfunction=scanfunction,cofactors=cofactors,verbose=verbose,...)
if(plot){
temp <- result
class(temp) <- c(class(temp),"mqmmulti")
mqmplot.multitrait(temp)
}
res <- c(res,result)
end <- proc.time()
SUM <- SUM + (end-start)[3]
AVG <- SUM/x
LEFT <- AVG*(batches-x)
if(verbose) {
cat("INFO: Done with batch",x,"/",batches,"\n")
cat("INFO: Calculation of batch",x,"took:",round((end-start)[3], digits=3),"seconds\n")
cat("INFO: Elapsed time:",(SUM%/%3600),":",(SUM%%3600)%/%60,":",round(SUM%%60, digits=0),"(Hour:Min:Sec)\n")
cat("INFO: Average time per batch:",round((AVG), digits=3)," per trait:",round((AVG %/% batchsize), digits=3),"seconds\n")
cat("INFO: Estimated time left:",LEFT%/%3600,":",(LEFT%%3600)%/%60,":",round(LEFT%%60,digits=0),"(Hour:Min:Sec)\n")
ourline()
}
}
}
if(FF){
if(verbose) cat(rownames(res[[1]]),"\n",res[[1]][,1],"\n",res[[1]][,2],"\n",file="out.frank")
for(i in 1:length(res)){
if(verbose) cat("INFO: Saving trait",i,"in frankformat\n")
qtl <- res[[i]]
if(verbose) cat(colnames(qtl)[3],qtl[,3],"\n",file="out.frank",append = TRUE)
}
}
#Return the results
if(length(res) > 1){
class(res) <- c(class(res),"mqmmulti")
}else{
class(res) <- c(class(res),"scanone")
}
#All done now plot the results
end <- proc.time()
SUM <- SUM + (end-start)[3]
AVG <- SUM/n.pheno
if(verbose) {
cat("------------------------------------------------------------------\n")
cat("INFO: Elapsed time:",(SUM%/%3600),":",(SUM%%3600)%/%60,":",round(SUM%%60, digits=0),"(Hour:Min:Sec)\n")
cat("INFO: Average time per trait:",round(AVG, digits=3),"seconds\n")
cat("------------------------------------------------------------------\n")
}
res
}
# end of scanall.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.