R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ######################################################################
> #
> # test_qtl.R
> #
> # copyright (c) 2009, Karl W Broman, Pjotr Prins
> # first written July 2009
> #
> #     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
> #
> # Some basic regression/integration testing for some of the QTL mapping routines
> #
> # You can run it with:
> #
> #   R --no-save --no-restore --no-readline --slave < ./tests/test_qtl.R 
> 
> ######################################################################
> 
> library(qtl)
> 
> version = mqm_version()
> cat("R/qtl=",version$RQTL)
R/qtl= 1.70> cat("R-MQM=",version$RMQM)
R-MQM= 0.90-pre1> cat("MQM=",version$MQM)
MQM= 0.90-pre1> 
> 
> data(listeria)
> if (nind(listeria)!=120) stop("Number of individuals incorrect")
> 
> # ---- a quick test of standard R/qtl scanone
> mr = scanone(listeria, method='mr')
Warning message:
In checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata,  :
  Dropping 4 individuals with missing phenotypes.

> test = round(mr[15,]$lod*1000)
> cat(mr[15,]$lod,test)
0.9663805 966> if (test != 966) stop("scanone_mr gives an incorrect result")
> 
> # ---- a quick test of MQM for R/qtl
> augmentedcross <- mqmaugment(listeria, minprob=1.0, verbose=TRUE)
INFO: Received a valid cross file type: f2 .
INFO: Number of individuals: 120 .
INFO: Number of chr: 19 .
INFO: Number of markers: 131 .
INFO: Done with augmentation
# Unique individuals before augmentation:120
# Unique selected individuals:120
# Marker p individual:131
# Individuals after augmentation:120
INFO: Data augmentation succesfull
INFO: DATA-Augmentation took:  0.08  seconds
Warning message:
In mqmaugment(listeria, minprob = 1, verbose = TRUE) :
  MQM not yet available for the X chromosome; omitting chr X
> nind = nind(augmentedcross)
> if (nind!=120) stop("Number of individuals incorrect: ",nind)
> result <- mqmscan(augmentedcross, logtransform=TRUE, outputmarkers = FALSE,off.end=0)
> test1 = round(result[5,5]*1000)
> test2 = round(max(result[,5]*1000))
> cat("test1 = ",test1,"\n")
test1 =  76 
> cat("test2 = ",test2,"\n")
test2 =  5384 
> if (test1 != 76) stop("MQM gives an unexpected result (1)")
> if (test2 != 5384) stop("MQM gives an unexpected result (2)")
> 
> # ---- Test for negative markerlocations
> data(hyper)
> hyper <- fill.geno(hyper)
> #Mess up the markers by shifting
> temp <- shiftmap(hyper, offset=10^7)
> out.temp <- mqmscan(temp,verb=TRUE,off.end=10)
INFO: Received a valid cross file type: bc .
INFO: Number of individuals:  250 
INFO: Number of chromosomes:  19 
INFO: Number of markers: 170 
INFO: 0 Cofactors received to be analyzed
INFO: Number of locations per chromosome:  30 
INFO: Dominance setting ignored (setting dominance to 0)
INFO: Starting C-part of the MQM analysis
INFO: Marker 6 at chr 1 is dropped
INFO: Marker 15 at chr 1 is dropped
INFO: Marker 16 at chr 1 is dropped
INFO: Marker 17 at chr 1 is dropped
INFO: Marker 42 at chr 4 is dropped
INFO: Marker 48 at chr 4 is dropped
INFO: Marker 107 at chr 11 is dropped
INFO: Marker 111 at chr 11 is dropped
INFO: Marker 133 at chr 15 is dropped
INFO: Marker 137 at chr 15 is dropped
INFO: Marker 139 at chr 15 is dropped
INFO: Marker 148 at chr 16 is dropped
INFO: Marker 150 at chr 17 is dropped
INFO: Marker 151 at chr 17 is dropped
INFO: Marker 154 at chr 17 is dropped
INFO: recombination parameters are not re-estimated
INFO: Re-estimation of the genetic map took 0 iterations, to reach a rdelta of 1.000000
INFO: Prob=0.020 Alfa=0.020000
INFO: Prob=0.019 Alfa=0.020000
INFO: dimX: 1, nInd: 250
INFO: F(Threshold, Degrees of freedom 1, Degrees of freedom 2) = Alfa
INFO: F(5.468, 1, 249) = 0.020000
INFO: F(4.003, 2, 249) = 0.020000
INFO: Log-likelihood of full model = -14206.535
INFO: Residual variance = 70.959
INFO: Trait mean= 101.611; Trait variation = 70.959
INFO: mapQTL function called
INFO: log-likelihood of full model = -14206.535300
INFO: Number of output datapoints: 570
INFO: All done in C returning to R
INFO: Calculation time (R->C,C,C-R): ( 0.01 , 0.64 , 0.56 ) (in seconds)
Warning message:
In omit_x_chr(cross) : Omitting X chromosome (X)
> if(!(rownames(out.temp)[3]=="D1Mit296")) stop("MQM something wrong with positive shifts in location")
> #Mess up the dataset by moving 1 marker infront of the chromosome
> hyper$geno[[1]]$map[1] <- -10
> res <- mqmscan(hyper,verbose=T,off.end=100)
INFO: Received a valid cross file type: bc .
INFO: Number of individuals:  250 
INFO: Number of chromosomes:  19 
INFO: Number of markers: 170 
INFO: 0 Cofactors received to be analyzed
INFO: Number of locations per chromosome:  69 
INFO: Dominance setting ignored (setting dominance to 0)
INFO: Starting C-part of the MQM analysis
INFO: recombination parameters are not re-estimated
INFO: Re-estimation of the genetic map took 0 iterations, to reach a rdelta of 1.000000
INFO: Prob=0.020 Alfa=0.020000
INFO: Prob=0.019 Alfa=0.020000
INFO: dimX: 1, nInd: 250
INFO: F(Threshold, Degrees of freedom 1, Degrees of freedom 2) = Alfa
INFO: F(5.468, 1, 249) = 0.020000
INFO: F(4.003, 2, 249) = 0.020000
INFO: Log-likelihood of full model = -14215.619
INFO: Residual variance = 70.959
INFO: Trait mean= 101.611; Trait variation = 70.959
INFO: mapQTL function called
INFO: log-likelihood of full model = -14215.619678
INFO: Number of output datapoints: 1311
INFO: All done in C returning to R
INFO: Calculation time (R->C,C,C-R): ( 0 , 1.83 , 1.06 ) (in seconds)
Warning message:
In omit_x_chr(cross) : Omitting X chromosome (X)
> if(any(is.na(res[,3]))) stop("MQM failed to handle negative locations correctly")
> if(!(rownames(res)[2]=="c1.loc-95")) stop("MQM something wrong with negative locations")   #to -15 because off.end defaults to 10
> 
> 
> cat("Version information:\n")
Version information:
> cat("R/qtl = ",version$RQTL,"\n")
R/qtl =  1.70 
> cat("R-MQM = ",version$RMQM,"\n")
R-MQM =  0.90-pre1 
> cat("MQM   = ",version$MQM,"\n\n")
MQM   =  0.90-pre1 

> 
> cat("test_qtl.R tests succesfully run!")
test_qtl.R tests succesfully run!> 
> proc.time()
   user  system elapsed 
   4.62    0.25    4.87