Nothing
#####################################################################
#
# mqmpermutation.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: mqmscanfdr
# mqmpermutation
# mqmprocesspermutation
#
#
#####################################################################
mqmscanfdr <- function(cross, scanfunction=mqmscanall, thresholds=c(1,2,3,4,5,7,10,15,20), n.perm = 10, verbose=FALSE, ...){
if(verbose){cat("Calculation of FDR estimate of threshold in multitrait analysis.\n")}
results <- NULL
above.in.real.res <- NULL
res <- scanfunction(cross,...)
for(threshold in thresholds){
above.in.real <- 0
for(x in 1:nphe(cross)){
above.in.real = above.in.real + sum(res[[x]][,3] > threshold)
}
above.in.real.res <- c(above.in.real.res,above.in.real)
}
perm <- cross
if(verbose){cat("QTL's above threshold:",above.in.real,"\n")}
above.in.perm.res <- rep(0,length(thresholds))
for(x in 1:n.perm){
if(verbose){cat("Starting permutation",x,"\n")}
perm.res <- NULL
neworder <- sample(nind(cross))
for(chr in 1:nchr(cross)){
perm$geno[[chr]]$data <- perm$geno[[chr]]$data[neworder,]
}
res <- scanfunction(perm, ...)
for(threshold in thresholds){
above.in.perm <- 0
for(y in 1:nphe(cross)){
above.in.perm = above.in.perm + sum(res[[y]][,3] > threshold)
}
perm.res <- c(perm.res,above.in.perm)
#if(verbose){cat("Permutation",x,"QTL's above threshold:",above.in.perm,"\n")}
}
above.in.perm.res <- above.in.perm.res+perm.res
}
above.in.perm.res <- above.in.perm.res/n.perm
results <- cbind(above.in.real.res,above.in.perm.res,above.in.perm.res/above.in.real.res)
rownames(results) <- thresholds
results
}
######################################################################
#
# mqmpermutation: Shuffles phenotype or does parametric bootstrapping of mqmscan
#
######################################################################
mqmpermutation <- function(cross,scanfunction=scanone,pheno.col=1,multicore=TRUE,n.perm=10,file="MQM_output.txt",n.cluster=1,method=c("permutation","simulation"),cofactors=NULL,plot=FALSE,verbose=FALSE,...)
{
bootmethod <- 0
supported <- c("permutation","simulation")
bootmethod <- pmatch(method, supported)[1]-1
if(missing(cross))
stop("No cross file. Please supply a valid cross object.")
crosstype <- crosstype(cross)
if(crosstype == "f2" || crosstype == "bc" || crosstype == "riself"){
#Echo back the cross type
if(verbose) {
cat("------------------------------------------------------------------\n")
cat("Starting permutation analysis\n")
cat("Number of permutations:",n.perm,"\n")
cat("n.cluster:",n.cluster,"\n")
cat("------------------------------------------------------------------\n")
cat("INFO: Received a valid cross file type:",crosstype,".\n")
}
b <- proc.time()
if(!bootmethod){
if(verbose) cat("INFO: Shuffleling traits between individuals.\n")
}else{
if(verbose) cat("INFO: Parametric permutation\nINFO: Calculating new traits for each individual.\n")
}
#Set the Phenotype under interest as the first
cross$pheno[,1] <- cross$pheno[,pheno.col]
names(cross$pheno)[1] <- names(cross$pheno)[pheno.col]
#Scan the original
#cross <- fill.geno(cross) # <- this should be done outside of this function
res0 <- lapply(1, FUN=snowCoreALL,all.data=cross,scanfunction=scanfunction,verbose=verbose,cofactors=cofactors,...)
#Setup bootstraps by generating a list of random numbers to set as seed for each bootstrap
batchsize <- n.perm
bootstraps <- runif(n.perm)
batches <- length(bootstraps) %/% batchsize
last.batch.num <- length(bootstraps) %% batchsize
results <- NULL
if(last.batch.num > 0){
batches = batches+1
}
SUM <- 0
AVG <- 0
LEFT <- 0
if(multicore && n.cluster >1) {
updateParallelRNG(n.cluster)
if(verbose) cat("INFO: Using ",n.cluster," Cores/CPU's/PC's for calculation.\n")
for(x in 1:(batches)){
start <- proc.time()
if(verbose) {
ourline()
cat("INFO: Starting with batch",x,"/",batches,"\n")
ourline()
}
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.cluster)
on.exit(stopCluster(cl))
res <- clusterApply(cl, boots, snowCoreBOOT, all.data=cross, scanfunction=scanfunction, bootmethod=bootmethod,
cofactors=cofactors, verbose=verbose, ...)
}
else {
res <- mclapply(boots, snowCoreBOOT, all.data=cross, scanfunction=scanfunction, bootmethod=bootmethod,
cofactors=cofactors, verbose=verbose, mc.cores=n.cluster, ...)
}
results <- c(results,res)
if(plot){
temp <- c(res0,results)
class(temp) <- c(class(temp),"mqmmulti")
mqmplot.permutations(temp)
}
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) {
ourline()
cat("INFO: Starting with batch",x,"/",batches,"\n")
ourline()
}
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)]
}
res <- lapply(boots, FUN=snowCoreBOOT,all.data=cross,scanfunction=scanfunction,bootmethod=bootmethod,cofactors=cofactors,verbose=verbose,...)
results <- c(results,res)
if(plot){
temp <- c(res0,results)
class(temp) <- c(class(temp),"mqmmulti")
mqmplot.permutations(temp)
}
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 run:",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()
}
}
}
res <- c(res0,results)
#Set the class of the result to mqmmulti (so we can use our plotting routines)
class(res) <- c(class(res),"mqmmulti")
e <- proc.time()
SUM <- (e-b)[3]
AVG <- SUM/(n.perm+1)
if(verbose) {
cat("INFO: Done with MQM permutation analysis\n")
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
}else{
stop("Currently only F2, BC, and selfed RIL crosses can be analyzed by MQM.")
}
}
mqmprocesspermutation <- function(mqmpermutationresult = NULL){
if(!is.null(mqmpermutationresult) && inherits(mqmpermutationresult, "mqmmulti")){
result <- NULL
result <- sapply(mqmpermutationresult[-1], function(a) max(a[,3], na.rm=TRUE))
result <- as.matrix(result)
colnames(result) <- colnames(mqmpermutationresult[[2]])[3]
rownames(result) <- 1:(length(mqmpermutationresult)-1)
class(result) <- c("scanoneperm",class(result))
result
}else{
stop("Please supply a valid resultobject (mqmmulti).")
}
}
# end of mqmpermutation.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.