R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

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 a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("qtl-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('qtl')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "argmax.geno"
> 
> ### * argmax.geno
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: argmax.geno
> ### Title: Reconstruct underlying genotypes
> ### Aliases: argmax.geno
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2,chr=18:19)
> ## End Don't show
> fake.f2 <- argmax.geno(fake.f2, step=2, off.end=5, err=0.01)
> 
> data(fake.bc)
> ## Don't show: 
> fake.bc <- subset(fake.bc,chr=18:19)
> ## End Don't show
> fake.bc <- argmax.geno(fake.bc, step=0, off.end=0)
> 
> 
> 
> cleanEx(); ..nameEx <- "arithscan"
> 
> ### * arithscan
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: arithscan
> ### Title: Arithmetic Operators for genome scan results
> ### Aliases: +.scanone -.scanone
> ### Keywords: arith
> 
> ### ** Examples
> 
> data(fake.bc)
> fake.bc <- calc.genoprob(fake.bc)
> # genome scans without covariates
> out.nocovar <- scanone(fake.bc)
> # genome scans with covariates
> ac <- fake.bc$pheno[,c("sex","age")]
> ic <- fake.bc$pheno[,"sex"]
> out.covar <- scanone(fake.bc, pheno.col=1,
+                      addcovar=ac, intcovar=ic)
> # plot the difference of with and without additive covariate
> plot(out.covar-out.nocovar)
> 
> 
> 
> cleanEx(); ..nameEx <- "badorder"
> 
> ### * badorder
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: badorder
> ### Title: An intercross with misplaced markers
> ### Aliases: badorder
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(badorder)
> ## Don't show: 
> badorder <- subset(badorder,chr=1)
> ## End Don't show
> 
> # estimate recombination fractions
> badorder <- est.rf(badorder)
> plot.rf(badorder)
> 
> # re-estimate map
> newmap <- est.map(badorder)
> plot.map(badorder,newmap)
> 
> # assess marker order on chr 1
> rip3 <- ripple(badorder,1,window=3)
   42 total orders
> summary(rip3)
                                     obligXO
Initial  1 2 3 4 5 6 7 8  9 10 11 12     415
1        1 2 3 4 5 6 7 8 11 10  9 12     361
2        1 2 3 4 6 5 7 8  9 10 11 12     363
3        1 2 3 6 4 5 7 8  9 10 11 12     367
4        1 2 3 4 5 6 7 8 11  9 10 12     374
5        1 2 3 4 5 6 7 8  9 12 10 11     391
6        1 2 3 4 5 6 7 8 10 11  9 12     395
7        1 2 3 4 6 7 5 8  9 10 11 12     401
8        1 2 3 4 5 6 7 8  9 11 10 12     404
9        1 2 3 4 5 6 7 8  9 12 11 10     406
10       1 2 3 4 5 6 7 8  9 10 12 11     407
11       1 2 3 5 6 4 7 8  9 10 11 12     411
12       1 2 3 4 5 6 7 8  9 11 12 10     412
13       1 2 3 5 4 6 7 8  9 10 11 12     413
14       1 2 3 6 5 4 7 8  9 10 11 12     413
15       2 1 3 4 5 6 7 8  9 10 11 12     414
16       1 2 3 4 5 6 7 9 10  8 11 12     415
17       1 2 3 4 5 6 7 8 10  9 11 12     417
attr(,"class")
[1] "summary.ripple" "matrix"        
> 
> 
> 
> cleanEx(); ..nameEx <- "bayesint"
> 
> ### * bayesint
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bayesint
> ### Title: Bayesian probability interval
> ### Aliases: bayesint
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> listeria <- subset(listeria,chr=c(5,13))
> listeria <- calc.genoprob(listeria, step=0.5)
> out <- scanone(listeria, model="np")
> bayesint(out, 5)
           chr  pos      lod
c5.loc17.5   5 17.5 4.097363
c5.loc27.5   5 27.5 5.418159
c5.loc37     5 37.0 4.195485
> bayesint(out, 13)
            chr      pos      lod
c13.loc17.5  13 17.50000 5.002267
D13M147      13 26.15954 6.705877
c13.loc29    13 29.00000 5.143878
> bayesint(out, 13, 0.99)
            chr      pos      lod
c13.loc14.5  13 14.50000 4.532136
D13M147      13 26.15954 6.705877
c13.loc30.5  13 30.50000 4.545499
> 
> 
> 
> cleanEx(); ..nameEx <- "bristle3"
> 
> ### * bristle3
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bristle3
> ### Title: Data on bristle number in Drosophila
> ### Aliases: bristle3
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(bristle3)
> # Summaries
> summary(bristle3)
    Backcross

    No. individuals:  66 

    No. phenotypes:   9 
    Percent phenotyped:  100 100 100 100 100 100 100 100 100 

    No. chromosomes:  1 
    Total markers:    29 
    No. markers:      29 
    Percent genotyped:  99.9 
    Genotypes (%):      AA:51.8  AB:48.2 
> plot(bristle3)
> 
> # genome scan for each of the average phenotypes
> bristle3 <- calc.genoprob(bristle3, step=2)
> out <- vector("list", 4)
> for(i in 1:4) out[[i]] <- scanone(bristle3, pheno.col=2*i-1)
> 
> # Plot the results
>     # maximum LOD score among four phenotypes    
> ym <- max(sapply(out, function(a) max(a[,3])))
> plot(out[[1]], ylim=c(0,ym))
> for(i in 2:4)
+   plot(out[[i]], add=TRUE, col=c("blue","red","green")[i-1])
> 
> 
> 
> cleanEx(); ..nameEx <- "bristleX"
> 
> ### * bristleX
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bristleX
> ### Title: Data on bristle number in Drosophila
> ### Aliases: bristleX
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(bristleX)
> # Summaries
> summary(bristleX)
    Backcross

    No. individuals:  92 

    No. phenotypes:   9 
    Percent phenotyped:  100 100 100 100 100 100 100 100 100 

    No. chromosomes:  1 
    Total markers:    17 
    No. markers:      17 
    Percent genotyped:  97.8 
    Genotypes (%):      AA:46.2  AB:53.8 
> plot(bristleX)
> 
> # genome scan for each of the average phenotypes
> bristleX <- calc.genoprob(bristleX, step=2)
> out <- vector("list", 4)
> for(i in 1:4) out[[i]] <- scanone(bristleX, pheno.col=2*i-1)
> 
> # Plot the results
>     # maximum LOD score among four phenotypes    
> ym <- max(sapply(out, function(a) max(a[,3])))
> plot(out[[1]], ylim=c(0,ym))
> for(i in 2:4) plot(out[[i]], add=TRUE, col=c("blue","red","green")[i-1])
> 
> 
> 
> cleanEx(); ..nameEx <- "c.cross"
> 
> ### * c.cross
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: c.cross
> ### Title: Combine data for QTL experiments
> ### Aliases: c.cross
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(fake.f2)
> junk <- fake.f2
> junk <- c(fake.f2,junk)
> 
> 
> 
> cleanEx(); ..nameEx <- "calc.errorlod"
> 
> ### * calc.errorlod
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: calc.errorlod
> ### Title: Identify likely genotyping errors
> ### Aliases: calc.errorlod
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(hyper)
> ## Don't show: 
> hyper <- subset(hyper,chr=c(1,19))
> ## End Don't show
> # run calc.genoprob, then calc.errorlod; use same error.prob
> hyper <- calc.genoprob(hyper,error.prob=0.01)
> hyper <- calc.errorlod(hyper,error.prob=0.01)
> 
> # plot the error LOD scores; print those above a specified cutoff 
> plot.errorlod(hyper)
> plot.errorlod(hyper,chr=1)
> top.errorlod(hyper,cutoff=4)
   chr ind   marker errorlod
1    1 118  D1Mit14 8.372802
2    1 162  D1Mit14 8.372802
3    1 170  D1Mit14 8.372802
4    1 159  D1Mit14 8.350349
5    1  73  D1Mit14 6.165405
6    1  65  D1Mit14 6.165405
7    1  88  D1Mit14 6.165405
8    1 184  D1Mit14 6.151614
9    1 241  D1Mit14 6.151614
10   1 215 D1Mit267 5.822197
11   1 108 D1Mit267 5.822197
12   1 138 D1Mit267 5.822197
13   1 226 D1Mit267 5.822197
14   1 199 D1Mit267 5.819255
15   1  84 D1Mit267 5.808406
> 
> # plot genotype data, flagging genotypes with error LOD > cutoff
> plot.geno(hyper, chr=1, ind=160:200, cutoff=7, min.sep=2)
> 
> # Note that if markers are extremely tightly linked,
> #     any recombination event will have a high error LOD
> plot.geno(hyper, chr=1, ind=101:118, min.sep=2)
> hyper$geno[[1]]$map
D1Mit296 D1Mit123 D1Mit156 D1Mit178  D1Mit19   D1Mit7  D1Mit46 D1Mit132 
     3.3     19.7     32.8     35.0     37.2     41.5     43.7     43.7 
D1Mit334 D1Mit305  D1Mit26  D1Mit94 D1Mit218 D1Mit100 D1Mit102  D1Mit14 
    49.2     54.6     64.5     67.8     69.9     74.3     75.4     82.0 
D1Mit105 D1Mit159 D1Mit267  D1Mit15 D1Mit456 D1Mit155 
    82.0     82.0     82.0     86.3     94.0    115.8 
> 
> 
> 
> cleanEx(); ..nameEx <- "calc.genoprob"
> 
> ### * calc.genoprob
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: calc.genoprob
> ### Title: Calculate conditional genotype probabilities
> ### Aliases: calc.genoprob
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2,chr=18:19)
> ## End Don't show
> fake.f2 <- calc.genoprob(fake.f2, step=2, off.end=5)
> 
> data(fake.bc)
> ## Don't show: 
> fake.bc <- subset(fake.bc,chr=18:19)
> ## End Don't show
> fake.bc <- calc.genoprob(fake.bc, step=0, off.end=0, err=0.01)
> 
> 
> 
> cleanEx(); ..nameEx <- "clean"
> 
> ### * clean
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: clean
> ### Title: Remove derived data
> ### Aliases: clean
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(fake.f2)
> names(fake.f2$geno)
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
[16] "16" "17" "18" "19" "X" 
> fake.f2 <- calc.genoprob(fake.f2)
> names(fake.f2$geno)
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
[16] "16" "17" "18" "19" "X" 
> fake.f2 <- clean(fake.f2)
> names(fake.f2$geno)
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
[16] "16" "17" "18" "19" "X" 
> 
> 
> 
> cleanEx(); ..nameEx <- "compareGeno"
> 
> ### * compareGeno
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: comparegeno
> ### Title: Compare individuals' genotype data
> ### Aliases: comparegeno
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> 
> output <- comparegeno(listeria)
> 
> n.ind <- nind(listeria)
> image(1:n.ind,1:n.ind,output,col=gray((0:99)/99),breaks=seq(0,1,len=101))
> 
> 
> 
> cleanEx(); ..nameEx <- "comparecrosses"
> 
> ### * comparecrosses
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: comparecrosses
> ### Title: Compare two cross objects.
> ### Aliases: comparecrosses
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> comparecrosses(listeria, listeria)
	Crosses are identical.
> 
> 
> 
> cleanEx(); ..nameEx <- "convert2ss"
> 
> ### * convert2ss
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: convert2ss
> ### Title: Convert an intercross into a form to allow sex-specific maps
> ### Aliases: convert2ss
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> listeria.ss <- convert2ss(listeria)
> 
> 
> 
> cleanEx(); ..nameEx <- "drop.markers"
> 
> ### * drop.markers
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: drop.markers
> ### Title: Drop a set of markers
> ### Aliases: drop.markers
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> listeria2 <- drop.markers(listeria, c("D10M44","D1M3","D1M75"))
> 
> 
> 
> cleanEx(); ..nameEx <- "drop.nullmarkers"
> 
> ### * drop.nullmarkers
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: drop.nullmarkers
> ### Title: Drop markers without any genotype data
> ### Aliases: drop.nullmarkers
> ### Keywords: manip
> 
> ### ** Examples
> 
> # removes one marker from hyper
> data(hyper)
> hyper <- drop.nullmarkers(hyper)
> 
> # shouldn't do anything to listeria
> data(listeria)
> listeria <- drop.nullmarkers(listeria)
> 
> 
> 
> cleanEx(); ..nameEx <- "effectplot"
> 
> ### * effectplot
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: effectplot
> ### Title: Plot phenotype means against genotypes at one or two markers.
> ### Aliases: effectplot
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2, chr=c(1, 13, "X"))
> ## End Don't show
> 
> # impute gneotype data
> ## Not run: 
> ##D fake.f2.impute <- sim.geno(fake.f2, step=5, n.draws=64)
> ## End(Not run)## Don't show: 
> fake.f2.impute <- sim.geno(fake.f2, step=5, n.draws=8)
> ## End Don't show
> 
> ########################################
> # one marker plots
> ########################################
> ### plot of genotype-specific phenotype means for 1 marker
> mname <- find.marker(fake.f2, 1, 37) # marker D1M437
> effectplot(fake.f2, pheno.col=1, mname1=mname)
> 
> ### plot a pseudomarker
> effectplot(fake.f2.impute, pheno.col=1, mname1=mname)
> 
> ### plot a phenotype
> # Plot of sex-specific phenotype means,
> # note that "sex" must be a phenotype name here
> effectplot(fake.f2, mname1="sex", geno1=c("F","M"))
> # alternatively I can do it this way
> sex <- fake.f2$pheno[,2]
> effectplot(fake.f2, mname1="Sex", mark1=sex, geno1=c("F","M"))
> 
> ########################################
> # two markers plots
> ########################################
> 
> ### plot two markers
> # plot of genotype-specific phenotype means for 2 markers
> mname1 <- find.marker(fake.f2, 1, 37) # marker D1M437
> mname2 <- find.marker(fake.f2, 13, 24) # marker D13M254
> effectplot(fake.f2, mname1=mname1, mname2=mname2)
> 
> ### plot two pseudomarkers
> pmname1 <- "c1.loc35"
> pmname2 <- "c13.loc25"
> effectplot(fake.f2.impute, mname1=pmname1, mname2=pmname2)
> 
> ### plot one marker and one imputed pseudomarker
> mark1 <- fake.f2.impute$geno[[1]]$draws[,11,]
> effectplot(fake.f2, mname1="D1M437", mark1=mark1, mname2="D13M254")
> 
> ### Plot of sex- and genotype-specific phenotype means 
> mname <- find.marker(fake.f2, 13, 24) # marker D13M254
> # sex and a marker
> effectplot(fake.f2, mname1=mname, mname2="Sex",
+            mark2=fake.f2$pheno$sex, geno2=c("F","M"))
> # sex and a pseudomarker
> effectplot(fake.f2.impute, mname1=mname, mname2="Sex",
+            mark2=fake.f2$pheno$sex, geno2=c("F","M"))
> 
> # Same as above, switch role of sex and the marker
> # sex and marker
> effectplot(fake.f2, mname1="Sex", mark1=fake.f2$pheno$sex,
+            geno1=c("F","M"), mname2=mname)
> # sex and pseudomarker
> effectplot(fake.f2.impute, mname1="Sex", mark1=fake.f2$pheno$sex,
+            geno1=c("F","M"), mname2=mname)
> 
> # X chromosome marker
> mname <- find.marker(fake.f2, "X", 14) # marker DXM66
> effectplot(fake.f2, mname1=mname)
> 
> # Two markers, including one on the X
> mname1 <- find.marker(fake.f2, 13, 24) # marker D13M254
> mname2 <- find.marker(fake.f2, "X", 14) # marker DXM66
> effectplot(fake.f2, mname1=mname1, mname2=mname2)
> 
> 
> 
> cleanEx(); ..nameEx <- "effectscan"
> 
> ### * effectscan
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: effectscan
> ### Title: Plot allelic effect across the whole genome
> ### Aliases: effectscan
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(fake.f2)
> # allelic effect on whole genome
> effectscan(fake.f2, ylim=c(-6,3))
> # on first 4 chromosomes, use triangle as marker ticks
> effectscan(fake.f2, chr=1:4, ylim=c(-6,3), mtick="triangle")
> 
> 
> 
> cleanEx(); ..nameEx <- "est.map"
> 
> ### * est.map
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: est.map
> ### Title: Estimate genetic maps
> ### Aliases: est.map
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2,chr=18:19)
> ## End Don't show
> newmap <- est.map(fake.f2)
> logliks <- sapply(newmap, attr, "loglik")
> plot.map(fake.f2, newmap)
> fake.f2 <- replace.map(fake.f2, newmap)
> 
> 
> 
> cleanEx(); ..nameEx <- "est.rf"
> 
> ### * est.rf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: est.rf
> ### Title: Estimate pairwise recombination fractions
> ### Aliases: est.rf
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(badorder)
> badorder <- est.rf(badorder)
> plot.rf(badorder)
> 
> 
> 
> cleanEx(); ..nameEx <- "fake.4way"
> 
> ### * fake.4way
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fake.4way
> ### Title: Simulated data for a 4-way cross
> ### Aliases: fake.4way
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(fake.4way)
> ## Don't show: 
> fake.4way <- subset(fake.4way,chr=c(2,7))
> ## End Don't show
> plot(fake.4way)
> summary(fake.4way)
    4-way cross

    No. individuals:  250 

    No. phenotypes:   2 
    Percent phenotyped:  100 100 

    No. chromosomes:  2 
    Total markers:    16 
    No. markers:      8 8 
    Percent genotyped:  94.7 
    Genotypes (%):      1:20.6  2:22.6  3:22.5  4:24.7  7:3.1  8:3.2  10:3.3 
> 
> # estimate recombination fractions
> fake.4way <- est.rf(fake.4way)
> plot.rf(fake.4way)
> 
> # estimate genetic maps
> ssmap <- est.map(fake.4way, verbose=TRUE)
Chr 2:
      0.130 0.028 0.065 0.021 0.017 0.276 0.034 
    7 0.163 0.039 0.080 0.032 0.013 0.356 0.032 
loglik: -1216.1585

Chr 7:
      0.029 0.052 0.091 0.029 0.018 0.151 0.012 
    6 0.029 0.093 0.124 0.044 0.020 0.185 0.032 
loglik: -1031.3994

> samap <- est.map(fake.4way, sex.sp=FALSE, verbose=TRUE)
Chr 2:
      0.130 0.029 0.062 0.023 0.018 0.274 0.031 
    5 0.164 0.038 0.080 0.032 0.013 0.356 0.032 
loglik: -1222.2000

Chr 7:
      0.030 0.050 0.090 0.032 0.020 0.149 0.011 
    6 0.030 0.093 0.124 0.044 0.020 0.185 0.032 
loglik: -1036.3476

> plot(ssmap, samap)
> 
> # error lod scores
> fake.4way <- calc.genoprob(fake.4way, err=0.01)
> fake.4way <- calc.errorlod(fake.4way, err=0.01)
> top.errorlod(fake.4way, cutoff=2.5)
  chr ind marker errorlod
1   2  89   D2M4 2.512582
> 
> # genome scan
> fake.4way <- calc.genoprob(fake.4way, step=2.5)
> out.hk <- scanone(fake.4way, method="hk")
> out.em <- scanone(fake.4way, method="em")
> plot(out.em,out.hk,chr=c(2,7))
> 
> 
> 
> cleanEx(); ..nameEx <- "fake.bc"
> 
> ### * fake.bc
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fake.bc
> ### Title: Simulated data for a backcross
> ### Aliases: fake.bc
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(fake.bc)
> ## Don't show: 
> fake.bc <- subset(fake.bc,chr=c(2,5,10))
> ## End Don't show
> summary(fake.bc)
    Backcross

    No. individuals:  400 

    No. phenotypes:   4 
    Percent phenotyped:  100 100 100 100 

    No. chromosomes:  3 
    Total markers:    19 
    No. markers:      5 10 4 
    Percent genotyped:  89.5 
    Genotypes (%):      AA:50.2  AB:49.8 
> plot(fake.bc)
> 
> # genome scans without covariates
> fake.bc <- calc.genoprob(fake.bc, step=2.5)
> out1.nocovar <- scanone(fake.bc, pheno.col=1)
> out2.nocovar <- scanone(fake.bc, pheno.col=2)
> 
> # genome scans with covariates
> ac <- fake.bc$pheno[,c("sex","age")]
> ic <- fake.bc$pheno[,"sex"]
> out1.covar <- scanone(fake.bc, pheno.col=1,
+                       addcovar=ac, intcovar=ic)
> out2.covar <- scanone(fake.bc, pheno.col=2,
+                       addcovar=ac, intcovar=ic)
> 
> # summaries
> summary(out1.nocovar, 3)

           chr pos lod  AA  AB sigma
c2.loc35     2  35 3.5 5.8 6.6   1.7
c5.loc17.5   5  18 7.8 5.7 6.7   1.6

> summary(out1.covar, 3)

            chr pos  lod
c2.loc27.5    2  28 11.7
c5.loc15      5  15  9.3
c10.loc27.5  10  28  6.6

> summary(out2.nocovar, 3)

       chr pos lod  AA  AB sigma
D5M394   5 9.8 3.9 5.6 6.4     2

> summary(out2.covar, 3)

            chr  pos  lod
c2.loc35      2 35.0 12.8
D5M394        5  9.8  3.8
c10.loc27.5  10 27.5 10.6

> 
> # plots
> plot(out1.covar,out1.nocovar,chr=c(2,5,10))
> plot(out2.covar,out2.nocovar,chr=c(2,5,10))
> 
> 
> 
> cleanEx(); ..nameEx <- "fake.f2"
> 
> ### * fake.f2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fake.f2
> ### Title: Simulated data for an F2 intercross
> ### Aliases: fake.f2
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(fake.f2)
> summary(fake.f2)
    F2 intercross

    No. individuals:  200 

    No. phenotypes:   3 
    Percent phenotyped:  100 100 100 

    No. chromosomes:  20 
    Total markers:    94 
    No. markers:      7 5 3 4 10 5 4 4 6 4 4 5 11 3 3 4 4 2 3 3 
    Percent genotyped:  85 
    Genotypes (%):      AA:25.9  AB:50.2  BB:23.9  not BB:0  not AA:0 
> plot(fake.f2)
> 
> 
> 
> cleanEx(); ..nameEx <- "fill.geno"
> 
> ### * fill.geno
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fill.geno
> ### Title: Fill holes in genotype data
> ### Aliases: fill.geno
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(hyper)
> ## Don't show: 
> hyper <- subset(hyper,chr=18:19)
> ## End Don't show
> out.mr <- scantwo(fill.geno(hyper,method="argmax"), method="mr")
 --Running scanone
 --Running scantwo
 (18,18)
 (18,19)
 (19,19)
> plot(out.mr)
> 
> 
> 
> cleanEx(); ..nameEx <- "find.flanking"
> 
> ### * find.flanking
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: find.flanking
> ### Title: Find flanking markers for a specified position
> ### Aliases: find.flanking
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> find.flanking(listeria, 5, 28)
          left  right  close
chr5:28 D5M357 D5M205 D5M357
> find.flanking(listeria, c(1, 5, 13), c(81, 28, 26))
            left   right   close
chr1:81   D1M113  D1M355  D1M113
chr5:28   D5M357  D5M205  D5M357
chr13:26 D13M106 D13M147 D13M147
> 
> 
> 
> cleanEx(); ..nameEx <- "find.marker"
> 
> ### * find.marker
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: find.marker
> ### Title: Find marker closest to a specified position
> ### Aliases: find.marker
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> find.marker(listeria, 5, 28)
[1] "D5M357"
> find.marker(listeria, c(1, 5, 13), c(81, 28, 26))
[1] "D1M113"  "D5M357"  "D13M147"
> 
> 
> 
> cleanEx(); ..nameEx <- "find.pheno"
> 
> ### * find.pheno
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: find.pheno
> ### Title: Find column number for a particular phenotype.
> ### Aliases: find.pheno
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(fake.bc)
> find.pheno(fake.bc, "sex")
[1] 3
> 
> 
> 
> cleanEx(); ..nameEx <- "fitqtl"
> 
> ### * fitqtl
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fitqtl
> ### Title: Fit a multiple QTL model
> ### Aliases: fitqtl
> ### Keywords: models
> 
> ### ** Examples
> 
> data(fake.f2)
> 
> # take out several QTLs and make QTL object
> qc <- c(1, 8, 13)
> qp <- c(26, 56, 28)
> fake.f2 <- subset(fake.f2, chr=qc)
> ## Don't show: 
> fake.f2 <- subset(fake.f2, ind=1:50)
> ## End Don't show
> fake.f2 <- sim.geno(fake.f2, n.draws=8, step=2, err=0.001)
> qtl <- makeqtl(fake.f2, qc, qp)
> 
> # fit model with 3 interacting QTLs interacting
> # (performing a drop-one-term analysis)
> lod <- fitqtl(fake.f2$pheno[,1], qtl, formula=y~Q1*Q2*Q3)
> summary(lod)

		Summary for fit QTL

Method is:  imp 
Number of observations:  50 

Full model result
----------------------------------  
Model formula is:  y ~ Q1 * Q2 * Q3 

      df        SS       MS      LOD     %var Pvalue(Chi2)  Pvalue(F)
Model 26  842.1752 32.39135 14.84610 74.52255 1.142824e-05 0.01208579
Error 23  287.9192 12.51823                                          
Total 49 1130.0944                                                   


Drop one QTL at a time ANOVA table: 
----------------------------------  
                         df Type III SS       LOD      %var   F value
Chr1@26                  18   218.87462   6.13902  19.36782   0.97136
Chr8@56                  18   333.43996   8.35180  29.50550   1.47980
Chr13@28                 18   415.00802   9.69099  36.72331   1.84179
Chr1@26:Chr8@56          12     7.33882   0.27328   0.64940   0.04885
Chr1@26:Chr13@28         12    27.02590   0.97411   2.39147   0.17991
Chr8@56:Chr13@28         12    84.27680   2.78753   7.45750   0.56103
Chr1@26:Chr8@56:Chr13@28  8    -5.35831  -0.20396  -0.47415  -0.05351
                         Pvalue(Chi2) Pvalue(F)  
Chr1@26                         0.058    0.5184  
Chr8@56                         0.003    0.1861  
Chr13@28                     0.000469    0.0838 .
Chr1@26:Chr8@56                 1.000    1.0000  
Chr1@26:Chr13@28                0.973    0.9982  
Chr8@56:Chr13@28                0.381    0.8504  
Chr1@26:Chr8@56:Chr13@28        1.000    1.0000  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

> 
> ## Not run: 
> ##D # fit an additive QTL model
> ##D lod.add <- fitqtl(fake.f2$pheno[,1], qtl, formula=y~Q1+Q2+Q3)
> ##D summary(lod.add)
> ##D 
> ##D # fit the model including sex as an interacting covariate
> ##D Sex <- data.frame(Sex=fake.f2$pheno$sex)
> ##D lod.sex <- fitqtl(fake.f2$pheno[,1], qtl, formula=y~Q1*Q2*Q3*Sex, cov=Sex)
> ##D summary(lod.sex)
> ##D 
> ##D # fit the same with an additive model
> ##D lod.sex.add <- fitqtl(fake.f2$pheno[,1], qtl, formula=y~Q1+Q2+Q3+Sex, cov=Sex)
> ##D summary(lod.sex.add)
> ## End(Not run)
> 
> 
> cleanEx(); ..nameEx <- "geno.table"
> 
> ### * geno.table
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: geno.table
> ### Title: Create table of genotype distributions
> ### Aliases: geno.table
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(listeria)
> geno.table(listeria)
        chr missing AA AB BB AA.AB AB.BB      P.value
D10M44    1      19 29 46 26     0     0 0.6125657411
D1M3      1       3 26 60 31     0     0 0.7771384555
D1M75     1      24 15 52 29     0     0 0.0930144892
D1M215    1       3 22 63 32     0     0 0.3009368909
D1M309    1       2 31 52 35     0     0 0.3805637036
D1M218    1       2 30 52 36     0     0 0.3212315037
D1M451    1       0 29 57 34     0     0 0.6988400895
D1M504    1       0 31 55 34     0     0 0.6116062006
D1M113    1       0 30 57 33     0     0 0.7985162188
D1M355    1       0 30 57 33     0     0 0.7985162188
D1M291    1       5 26 59 30     0     0 0.8367241192
D1M209    1      12 26 57 25     0     0 0.8386801055
D1M155    1       5 26 59 30     0     0 0.8367241192
D2M365    2       7 24 61 28     0     0 0.6065306597
D2M37     2       1 29 60 30     0     0 0.9874740688
D2M396    2       3 27 62 28     0     0 0.8041666799
D2M493    2      30 24 44 22     0     0 0.9355069850
D2M226    2       9 33 54 24     0     0 0.4628879466
D2M148    2       0 34 62 24     0     0 0.4065696597
D3M265    3       0 33 66 21     0     0 0.1652988882
D3M51     3      33 26 46 15     0     0 0.2155671626
D3M106    3      30 28 45 17     0     0 0.2606844924
D3M257    3       2 33 69 16     0     0 0.0158582121
D3M147    3       5 29 69 17     0     0 0.0286622626
D3M19     3      30 23 50 17     0     0 0.3845984194
D4M2      4      30 18 49 23     0     0 0.5308194506
D4M178    4      44 12 46 18     0     0 0.1155681665
D4M187    4      30 13 48 29     0     0 0.0476227617
D4M251    4       0 22 62 36     0     0 0.1826835241
D5M148    5      11 27 47 35     0     0 0.1980439991
D5M232    5       7 26 53 34     0     0 0.4569479266
D5M257    5       0 32 55 33     0     0 0.6537697851
D5M83     5       0 31 56 33     0     0 0.7408182207
D5M307    5       4 32 53 31     0     0 0.6442585409
D5M357    5       0 31 58 31     0     0 0.9355069850
D5M205    5       4 27 56 33     0     0 0.6843332004
D5M398    5      57 15 30 18     0     0 0.8071177470
D5M91     5       0 27 59 34     0     0 0.6537697851
D5M338    5       0 27 61 32     0     0 0.7985162188
D5M188    5       1 28 57 34     0     0 0.6652695213
D5M29     5      58 13 25 24     0     0 0.0444716565
D5M168    5       0 30 52 38     0     0 0.2018965180
D6M223    6       1 32 59 28     0     0 0.8705279066
D6M188    6      30 21 54 15     0     0 0.1108031584
D6M284    6       0 27 76 17     0     0 0.0060967466
D6M39     6       0 24 75 21     0     0 0.0218184355
D6M254    6       0 18 77 25     0     0 0.0053803600
D6M194    6      30 13 54 23     0     0 0.0544152349
D6M290    6      53 13 35 19     0     0 0.5463597639
D6M25     6       4 20 67 29     0     0 0.1230914645
D6M339    6       0 22 68 30     0     0 0.2018965180
D6M59_    6      30 16 49 25     0     0 0.2849175167
D6M201    6       1 23 66 30     0     0 0.3256771455
D6M15     6       1 22 67 30     0     0 0.2269119375
D6M294    6       0 19 70 31     0     0 0.0568882383
D7M246    7      30 22 42 26     0     0 0.6853827910
D7M145    7       4 28 58 30     0     0 0.9661049965
D7M62     7       3 28 61 28     0     0 0.8986715993
D7M126    7      27 29 37 27     0     0 0.1375358089
D7M105    7       0 30 56 34     0     0 0.6703200460
D7M259    7      46 17 41 16     0     0 0.6402183775
D8M94     8       0 28 58 34     0     0 0.6930406201
D8M339    8       2 28 55 35     0     0 0.5033645208
D8M178    8       5 25 58 32     0     0 0.6502263262
D8M242    8      30 23 42 25     0     0 0.7831394949
D8M213    8       3 30 52 35     0     0 0.3922337029
D8M156    8      18 25 49 28     0     0 0.8464817249
D9M247    9       1 30 59 30     0     0 0.9958071340
D9M328    9      10 30 49 31     0     0 0.5149752882
D9M106    9       0 27 58 35     0     0 0.5488116361
D9M269    9      38 22 41 19     0     0 0.8960526581
D9M346    9      22 23 43 32     0     0 0.2098789192
D9M55     9      30 27 42 21     0     0 0.5488116361
D9M18     9      30 28 41 21     0     0 0.4065696597
D10M298  10       1 40 59 20     0     0 0.0345431429
D10M294  10      30 31 45 14     0     0 0.0403117975
D10M42_  10       4 30 63 23     0     0 0.4259436255
D10M10   10      11 28 61 20     0     0 0.2560492882
D10M233  10       0 29 62 29     0     0 0.9355069850
D11M78   11      30 23 47 20     0     0 0.8278784881
D11M20   11       3 23 63 31     0     0 0.4093591491
D11M242  11       9 22 57 32     0     0 0.3900651738
D11M356  11      31 18 43 28     0     0 0.3090793301
D11M327  11       0 28 54 38     0     0 0.2385125539
D11M333  11       3 27 57 33     0     0 0.7074036474
D12M105  12      30 31 35 24     0     0 0.0628712266
D12M46   12      33 34 33 20     0     0 0.0083344619
D12M34   12       2 37 62 19     0     0 0.0551165588
D12M5    12       6 38 51 25     0     0 0.1207497463
D12M99   12      30 29 43 18     0     0 0.2385125539
D12M150  12       5 33 50 32     0     0 0.3727092993
D13M59   13       0 34 14  7     0    65           NA
D13M88   13       6 33 61 20     0     0 0.1715024212
D13M21   13      17 35 47 21     0     0 0.1006489555
D13M39   13      32 33 35 20     0     0 0.0232520115
D13M167  13      31 32 35 22     0     0 0.0427799650
D13M99   13       0 48 49 23     0     0 0.0007281525
D13M233  13      14 41 43 22     0     0 0.0050294088
D13M106  13       0 45 55 20     0     0 0.0036065631
D13M147  13       0 45 55 20     0     0 0.0036065631
D13M226  13      28 27 49 16     0     0 0.2207179658
D13M290  13       1 42 58 19     0     0 0.0112972801
D13M151  13       1 38 58 23     0     0 0.1453557012
D14M14   14      32 20 44 24     0     0 0.8337529181
D14M115  14       9 27 57 27     0     0 0.9602702338
D14M265  14       4 30 61 25     0     0 0.6902581264
D14M266  14      30 19 53 18     0     0 0.2385125539
D15M226  15       8 27 58 27     0     0 0.9310627797
D15M100  15      18 27 47 28     0     0 0.7235906755
D15M209  15       1 31 55 33     0     0 0.6880116006
D15M144  15      30 24 41 25     0     0 0.6930406201
D15M68   15       3 33 51 33     0     0 0.3823042729
D15M239  15       0 32 57 31     0     0 0.8535652128
D15M241  15       5 35 54 26     0     0 0.3995600132
D15M34   15      42 23 31 24     0     0 0.1913126738
D16M154  16      30 28 35 27     0     0 0.1071705988
D16M4    16       1 33 53 33     0     0 0.4916028846
D16M139  16      31 26 39 24     0     0 0.4844606343
D16M86   16      30 24 45 21     0     0 0.9048374180
D17M260  17      30 22 46 22     0     0 0.9780228725
D17M66   17      30 19 47 24     0     0 0.6930406201
D17M88   17      15 22 57 26     0     0 0.5838593069
D17M129  17       3 27 58 32     0     0 0.8041666799
D18M94   18      30 30 40 20     0     0 0.1888756028
D18M58   18       7 32 54 27     0     0 0.7175889212
D18M106  18      40 23 40 17     0     0 0.6376281516
D18M186  18       0 35 53 32     0     0 0.4099718965
D19M68   19      30 23 45 22     0     0 0.9889503893
D19M117  19      34 26 43 17     0     0 0.3899017615
D19M65   19      32 27 43 18     0     0 0.3893869039
D19M10   19      31 26  0  0     0    63           NA
DXM186    X       1 69 50  0     0     0 0.0815562004
DXM64     X       5 52 63  0     0     0 0.3050069459
> 
> 
> 
> cleanEx(); ..nameEx <- "hyper"
> 
> ### * hyper
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hyper
> ### Title: Data on hypertension
> ### Aliases: hyper
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(hyper)
> summary(hyper)
    Backcross

    No. individuals:  250 

    No. phenotypes:   1 
    Percent phenotyped:  100 

    No. chromosomes:  20 
    Total markers:    174 
    No. markers:      22 8 6 20 14 11 7 6 5 5 14 5 5 5 11 6 12 4 4 4 
    Percent genotyped:  47.7 
    Genotypes (%):      AA:50.2  AB:49.8 
> plot(hyper)
> 
> # Note the selective genotyping
> plot.missing(hyper, reorder=TRUE)
> 
> # A marker on c14 has no data; remove it
> hyper <- drop.nullmarkers(hyper)
> 
> 
> 
> cleanEx(); ..nameEx <- "listeria"
> 
> ### * listeria
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: listeria
> ### Title: Data on Listeria monocytogenes susceptibility
> ### Aliases: listeria
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(listeria)
> ## Don't show: 
> listeria <- subset(listeria,chr=c(1,5,6,13,15))
> ## End Don't show
> # Summaries
> summary(listeria)
    F2 intercross

    No. individuals:  120 

    No. phenotypes:   1 
    Percent phenotyped:  96.7 

    No. chromosomes:  5 
    Total markers:    59 
    No. markers:      13 13 13 12 8 
    Percent genotyped:  91.5 
    Genotypes (%):      AA:25.5  AB:48.8  BB:24.7  not BB:0  not AA:1 
> plot(listeria)
> 
> # Take log of phenotype
> listeria$pheno[,2] <- log2(listeria$pheno[,1])
> names(listeria$pheno)[2] <- "logSurv"
> plot(listeria)
> 
> # Genome scan with a two-part model, using log survival
> listeria <- calc.genoprob(listeria, step=2)
> out <- scanone(listeria, model="2part", method="em",
+                upper=TRUE, pheno.col=2)
> 
> # Summary of the results
> summary(out, 4)

          chr pos lod lod.p lod.mu mu.AA mu.AB mu.BB pi.AA pi.AB pi.BB sigma
D1M355      1  81 5.4  0.57   4.88   6.4   6.7   6.9  0.21  0.30 0.394  0.35
c5.loc28    5  28 6.8  6.13   0.66   6.9   6.7   6.6  0.64  0.29 0.031  0.39
D6M223      6  10 4.1  2.58   1.52   6.7   6.7   6.4  0.13  0.28 0.535  0.38
D13M147    13  26 7.4  3.66   3.73   6.5   6.9   6.7  0.09  0.43 0.421  0.36
c15.loc16  15  16 4.6  2.06   2.48   6.9   6.6   6.5  0.22  0.44 0.137  0.37

> 
> # Plot LOD curves for interesting chromosomes
> #     (The two-part model gives three LOD scores)
> plot(out, out, out, chr=c(1,5,6,13,15), lodcolumn=3:5,
+      lty=1, col=c("black","red","blue"))
> 
> 
> 
> cleanEx(); ..nameEx <- "lodint"
> 
> ### * lodint
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lodint
> ### Title: LOD support interval
> ### Aliases: lodint
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> listeria <- subset(listeria,chr=c(5,13))
> listeria <- calc.genoprob(listeria, step=0.5)
> out <- scanone(listeria, model="np")
> lodint(out, 5, 1.5)
           chr      pos      lod
c5.loc15     5 15.00000 3.894975
c5.loc27.5   5 27.50000 5.418159
D5M338       5 38.06807 3.915563
> lodint(out, 13, 1.5)
          chr      pos      lod
c13.loc23  13 23.00000 5.172082
D13M147    13 26.15954 6.705877
c13.loc29  13 29.00000 5.143878
> 
> 
> 
> cleanEx(); ..nameEx <- "makeqtl"
> 
> ### * makeqtl
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: makeqtl
> ### Title: Make a qtl object
> ### Aliases: makeqtl
> ### Keywords: models
> 
> ### ** Examples
> 
> data(fake.f2)
> 
> # take out several QTLs and make QTL object
> qc <- c(1, 6, 13)
> qp <- c(25.8, 33.6, 18.63)
> fake.f2 <- subset(fake.f2, chr=qc)
> ## Don't show: 
> fake.f2 <- subset(fake.f2, ind=1:50)
> ## End Don't show
> fake.f2 <- sim.geno(fake.f2, n.draws=8, step=2, err=0.001)
> qtl <- makeqtl(fake.f2, qc, qp)
> 
> 
> 
> cleanEx(); ..nameEx <- "map10"
> 
> ### * map10
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: map10
> ### Title: An example genetic map
> ### Aliases: map10
> ### Keywords: datasets
> 
> ### ** Examples
> 
>   data(map10)
>   plot(map10)
> 
>   mycross <- sim.cross(map10, type="f2", n.ind=100)
> 
> 
> 
> cleanEx(); ..nameEx <- "max.scanone"
> 
> ### * max.scanone
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: max.scanone
> ### Title: Maximum peak in genome scan
> ### Aliases: max.scanone
> ### Keywords: print
> 
> ### ** Examples
> 
> data(listeria)
> ## Don't show: 
> listeria <- subset(listeria,chr=c(5,13))
> ## End Don't show
> listeria <- calc.genoprob(listeria, step=2.5)
> out <- scanone(listeria, model="2part")
> # Maximum peak for LOD(p,mu)
> max(out)

       chr pos lod lod.p lod.mu mu.AA mu.AB mu.BB pi.AA pi.AB pi.BB sigma
D5M357   5  26 6.7  0.33    6.4   207   153   106     0 0.018     0    68

> max(out,chr=13)

        chr pos lod lod.p lod.mu mu.AA mu.AB mu.BB pi.AA pi.AB pi.BB sigma
D13M147  13  26   6  0.42    5.6   110   185   171 0.023     0     0    69

> 
> # Maximum peak for LOD(p)
> max(out[,-3])

       chr pos lod.p lod.mu mu.AA mu.AB mu.BB pi.AA pi.AB pi.BB sigma
D5M148   5   0   0.6    2.0   189   152   131 0.034     0     0    74

> 
> # Maximum peak for LOD(mu)
> max(out[,-(3:4)])

           chr pos lod.mu mu.AA mu.AB mu.BB pi.AA   pi.AB pi.BB sigma
c5.loc27.5   5  28    6.6   213   153   107     0 8.4e-05 0.031    67

> 
> 
> 
> cleanEx(); ..nameEx <- "max.scantwo"
> 
> ### * max.scantwo
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: max.scantwo
> ### Title: Maximum peak in two-dimensional genome scan
> ### Aliases: max.scantwo
> ### Keywords: print
> 
> ### ** Examples
> 
> data(fake.f2)
> fake.f2 <- calc.genoprob(fake.f2, step=10)
> out.2dim <- scantwo(fake.f2, method="hk", chr=c(1,6,13))
 --Running scanone
 --Running scantwo
 (1,1)
 (1,6)
 (1,13)
 (6,6)
 (6,13)
 (13,13)
> max(out.2dim)

       pos1 pos2   LODjnt -logP   LODint -logP
c1: c6  100    0     5.12   2.6     4.15   3.1
c1:c13   30   30    13.86  10.1     0.78   0.3

> 
> 
> 
> cleanEx(); ..nameEx <- "movemarker"
> 
> ### * movemarker
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: movemarker
> ### Title: Move a marker to a new chromosome
> ### Aliases: movemarker
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(badorder)
> badorder <- movemarker(badorder, "D2M937", 3, 48.15)
> badorder <- movemarker(badorder, "D3M160", 2, 28.83)
> 
> 
> 
> cleanEx(); ..nameEx <- "nmissing"
> 
> ### * nmissing
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: nmissing
> ### Title: Number of missing genotypes
> ### Aliases: nmissing
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> 
> # plot number of missing genotypes for each individual
> plot(nmissing(listeria))
> 
> # plot number of missing genotypes for each marker
> plot(nmissing(listeria, which="mar"))
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.cross"
> 
> ### * plot.cross
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.cross
> ### Title: Plot various features of a cross object
> ### Aliases: plot.cross
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(fake.bc)
> plot(fake.bc)
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.errorlod"
> 
> ### * plot.errorlod
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.errorlod
> ### Title: Plot grid of error LOD values
> ### Aliases: plot.errorlod
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(hyper)
> ## Don't show: 
> hyper <- subset(hyper,chr=c(1,19))
> ## End Don't show
> # run calc.genoprob, then calc.errorlod; use same error.prob
> hyper <- calc.genoprob(hyper,error.prob=0.01)
> hyper <- calc.errorlod(hyper,error.prob=0.01)
> 
> # plot the error LOD scores; print those above a specified cutoff 
> plot.errorlod(hyper)
> plot.errorlod(hyper,chr=1)
> top.errorlod(hyper,cutoff=4)
   chr ind   marker errorlod
1    1 118  D1Mit14 8.372802
2    1 162  D1Mit14 8.372802
3    1 170  D1Mit14 8.372802
4    1 159  D1Mit14 8.350349
5    1  73  D1Mit14 6.165405
6    1  65  D1Mit14 6.165405
7    1  88  D1Mit14 6.165405
8    1 184  D1Mit14 6.151614
9    1 241  D1Mit14 6.151614
10   1 215 D1Mit267 5.822197
11   1 108 D1Mit267 5.822197
12   1 138 D1Mit267 5.822197
13   1 226 D1Mit267 5.822197
14   1 199 D1Mit267 5.819255
15   1  84 D1Mit267 5.808406
> 
> # plot genotype data, flagging genotypes with error LOD > cutoff
> plot.geno(hyper, chr=1, ind=160:200, cutoff=7, min.sep=2)
> 
> # Note that if markers are extremely tightly linked,
> #     any recombination event will have a high error LOD
> plot.geno(hyper, chr=1, ind=101:118, min.sep=2)
> hyper$geno[[1]]$map
D1Mit296 D1Mit123 D1Mit156 D1Mit178  D1Mit19   D1Mit7  D1Mit46 D1Mit132 
     3.3     19.7     32.8     35.0     37.2     41.5     43.7     43.7 
D1Mit334 D1Mit305  D1Mit26  D1Mit94 D1Mit218 D1Mit100 D1Mit102  D1Mit14 
    49.2     54.6     64.5     67.8     69.9     74.3     75.4     82.0 
D1Mit105 D1Mit159 D1Mit267  D1Mit15 D1Mit456 D1Mit155 
    82.0     82.0     82.0     86.3     94.0    115.8 
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.geno"
> 
> ### * plot.geno
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.geno
> ### Title: Plot observed genotypes, flagging likely errors
> ### Aliases: plot.geno
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(hyper)
> ## Don't show: 
> hyper <- subset(hyper, chr=c(1,19))
> ## End Don't show
> # run calc.genoprob, then calc.errorlod; use same error.prob
> hyper <- calc.genoprob(hyper,error.prob=0.01)
> hyper <- calc.errorlod(hyper,error.prob=0.01)
> 
> # plot the error LOD scores; print those above a specified cutoff 
> plot.errorlod(hyper)
> plot.errorlod(hyper,chr=1)
> top.errorlod(hyper,cutoff=4)
   chr ind   marker errorlod
1    1 118  D1Mit14 8.372802
2    1 162  D1Mit14 8.372802
3    1 170  D1Mit14 8.372802
4    1 159  D1Mit14 8.350349
5    1  73  D1Mit14 6.165405
6    1  65  D1Mit14 6.165405
7    1  88  D1Mit14 6.165405
8    1 184  D1Mit14 6.151614
9    1 241  D1Mit14 6.151614
10   1 215 D1Mit267 5.822197
11   1 108 D1Mit267 5.822197
12   1 138 D1Mit267 5.822197
13   1 226 D1Mit267 5.822197
14   1 199 D1Mit267 5.819255
15   1  84 D1Mit267 5.808406
> 
> # plot genotype data, flagging genotypes with error LOD > cutoff
> plot.geno(hyper, chr=1, ind=160:200, cutoff=7, min.sep=2)
> 
> # Note that if markers are extremely tightly linked,
> #     any recombination event will have a high error LOD
> plot.geno(hyper, chr=1, ind=101:118, min.sep=2)
> hyper$geno[[1]]$map
D1Mit296 D1Mit123 D1Mit156 D1Mit178  D1Mit19   D1Mit7  D1Mit46 D1Mit132 
     3.3     19.7     32.8     35.0     37.2     41.5     43.7     43.7 
D1Mit334 D1Mit305  D1Mit26  D1Mit94 D1Mit218 D1Mit100 D1Mit102  D1Mit14 
    49.2     54.6     64.5     67.8     69.9     74.3     75.4     82.0 
D1Mit105 D1Mit159 D1Mit267  D1Mit15 D1Mit456 D1Mit155 
    82.0     82.0     82.0     86.3     94.0    115.8 
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.info"
> 
> ### * plot.info
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.info
> ### Title: Plot the proportion of missing genotype information
> ### Aliases: plot.info
> ### Keywords: hplot univar
> 
> ### ** Examples
> 
> data(hyper)
> ## Don't show: 
> hyper <- subset(hyper,chr=1:4)
> ## End Don't show
> hyper <- calc.genoprob(hyper, step=2.5, off.end=5)
> plot.info(hyper,chr=c(1,4))
> 
> # save the results and view maximum missing info on each chr
> info <- plot.info(hyper)
> summary(info)

           chr   pos misinfo.entropy misinfo.variance
c1.loc-5     1  -1.7            0.61             0.54
c2.loc97.5   2 105.2            0.70             0.64
c3.loc70     3  72.2            0.69             0.63
c4.loc80     4  80.0            0.55             0.46

> 
> 
> 
> cleanEx(); ..nameEx <- "plot.map"
> 
> ### * plot.map
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.map
> ### Title: Plot genetic map
> ### Aliases: plot.map
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(fake.bc)
> ## Don't show: 
> fake.bc <- subset(fake.bc, chr=18:19)
> ## End Don't show
> plot.map(fake.bc)
> plot.map(fake.bc,horizontal=TRUE)
> newmap <- est.map(fake.bc)
> plot(newmap)
> plot.map(fake.bc, newmap)
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.missing"
> 
> ### * plot.missing
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.missing
> ### Title: Plot grid of missing genotypes
> ### Aliases: plot.missing
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(fake.f2)
> plot.missing(fake.f2)
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.pxg"
> 
> ### * plot.pxg
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.pxg
> ### Title: Plot phenotypes versus marker genotypes.
> ### Aliases: plot.pxg
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(listeria)
> mname <- find.marker(listeria, 5, 28) # marker D5M357
> plot.pxg(listeria, mname)
> 
> mname2 <- find.marker(listeria, 13, 26) # marker D13Mit147
> plot.pxg(listeria, c(mname, mname2))
> 
> data(fake.f2)
> mname <- find.marker(fake.f2, 1, 37) # marker D1M437
> plot.pxg(fake.f2, mname)
> 
> mname2 <- find.marker(fake.f2, "X", 14) # marker DXM66
> plot.pxg(fake.f2, mname2)
> 
> plot.pxg(fake.f2, c(mname,mname2))
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.rf"
> 
> ### * plot.rf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.rf
> ### Title: Plot recombination fractions
> ### Aliases: plot.rf
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(badorder)
> badorder <- est.rf(badorder)
> plot.rf(badorder)
> 
> # plot just chr 1
> plot.rf(badorder,chr=1)
> 
> # plot just the recombination fractions
> plot.rf(badorder,which="rf") 
> 
> # plot just the LOD scores, and just for chr 2 and 3
> plot.rf(badorder,which="lod")
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.scanone"
> 
> ### * plot.scanone
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.scanone
> ### Title: Plot LOD curves
> ### Aliases: plot.scanone
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2, chr=c(1,13))
> ## End Don't show
> fake.f2 <- calc.genoprob(fake.f2,step=2.5)
> out.mr <- scanone(fake.f2, method="mr")
> out.em <- scanone(fake.f2, method="em")
> plot(out.mr)
> plot(out.mr, out.em, chr=c(1,13), lty=1, col=c("violetred","black"))
> out.hk <- scanone(fake.f2, method="hk")
> plot(out.hk, chr=c(1,13), add=TRUE, col="slateblue")
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.scantwo"
> 
> ### * plot.scantwo
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.scantwo
> ### Title: Plot LOD scores for a two-dimensional genome scan
> ### Aliases: plot.scantwo
> ### Keywords: hplot
> 
> ### ** Examples
> 
> data(fake.f2)
> fake.f2 <- calc.genoprob(fake.f2, step=5)
> out.2dim <- scantwo(fake.f2, method="hk", chr=c(1,8,13,14))
 --Running scanone
 --Running scantwo
 (1,1)
 (1,8)
 (1,13)
 (1,14)
 (8,8)
 (8,13)
 (8,14)
 (13,13)
 (13,14)
 (14,14)
> plot(out.2dim)
> plot(out.2dim, nodiag=FALSE)
> plot(out.2dim, lower="joint", nodiag=FALSE)
> 
> 
> 
> cleanEx(); ..nameEx <- "print.cross"
> 
> ### * print.cross
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: print.cross
> ### Title: Print a QTL experiment
> ### Aliases: print.cross
> ### Keywords: print
> 
> ### ** Examples
> 
> data(hyper)
> hyper
  This is an object of class "cross".
  It is too complex to print, so we provide just this summary.
    Backcross

    No. individuals:  250 

    No. phenotypes:   1 
    Percent phenotyped:  100 

    No. chromosomes:  20 
    Total markers:    174 
    No. markers:      22 8 6 20 14 11 7 6 5 5 14 5 5 5 11 6 12 4 4 4 
    Percent genotyped:  47.7 
    Genotypes (%):      AA:50.2  AB:49.8 
> 
> print(hyper)
  This is an object of class "cross".
  It is too complex to print, so we provide just this summary.
    Backcross

    No. individuals:  250 

    No. phenotypes:   1 
    Percent phenotyped:  100 

    No. chromosomes:  20 
    Total markers:    174 
    No. markers:      22 8 6 20 14 11 7 6 5 5 14 5 5 5 11 6 12 4 4 4 
    Percent genotyped:  47.7 
    Genotypes (%):      AA:50.2  AB:49.8 
    Backcross

    No. individuals:  250 

    No. phenotypes:   1 
    Percent phenotyped:  100 

    No. chromosomes:  20 
    Total markers:    174 
    No. markers:      22 8 6 20 14 11 7 6 5 5 14 5 5 5 11 6 12 4 4 4 
    Percent genotyped:  47.7 
    Genotypes (%):      AA:50.2  AB:49.8 
> 
> 
> 
> cleanEx(); ..nameEx <- "pull.geno"
> 
> ### * pull.geno
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pull.geno
> ### Title: Pull out the genotype data from a cross
> ### Aliases: pull.geno
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> dat <- pull.geno(listeria)
> 
> # image of the genotype data
> image(1:ncol(dat),1:nrow(dat),t(dat),ylab="Individuals",xlab="Markers",
+       col=c("red","yellow","blue","green","violet"))
> abline(v=cumsum(c(0,nmar(listeria)))+0.5)
> abline(h=nrow(dat)+0.5)
> 
> 
> 
> cleanEx(); ..nameEx <- "pull.map"
> 
> ### * pull.map
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pull.map
> ### Title: Pull out the genetic map from a cross
> ### Aliases: pull.map
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(fake.f2)
> map <- pull.map(fake.f2)
> plot(map)
> 
> 
> 
> cleanEx(); ..nameEx <- "read.cross"
> 
> ### * read.cross
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: read.cross
> ### Title: Read data for a QTL experiment
> ### Aliases: read.cross
> ### Keywords: IO
> 
> ### ** Examples
> 
> ## Not run: 
> ##D # comma-delimited format
> ##D dat1 <- read.cross("csv", dir="Mydata", file="mydata.csv")
> ##D 
> ##D # Mapmaker format
> ##D dat2 <- read.cross("mm", dir="Mydata", file="mydata.raw",
> ##D                    mapfile="mydata.map")
> ##D 
> ##D # Map Manager QTX format
> ##D dat3 <- read.cross("qtx", dir="Mydata", file="mydata.qtx")
> ##D 
> ##D # QTL Cartographer format
> ##D dat4 <- read.cross("qtlcart", dir="Mydata", file="qtlcart.cro",
> ##D                    mapfile="qtlcart.map")
> ##D 
> ##D # Gary format
> ##D dat5 <- read.cross("gary", dir="Mydata", genfile="geno.dat",
> ##D                    mapfile="markerpos.txt", phefile="pheno.dat",
> ##D                    chridfile="chrid.dat", mnamesfile="mnames.txt",
> ##D                    pnamesfile="pnames.txt")
> ##D 
> ##D # Karl format
> ##D dat6 <- read.cross("karl", dir="Mydata", genfile="gen.txt",
> ##D                    phefile="phe.txt", mapfile="map.txt")
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "replace.map"
> 
> ### * replace.map
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: replace.map
> ### Title: Replace the genetic map of a cross
> ### Aliases: replace.map
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2, chr=18:19)
> ## End Don't show
> newmap <- est.map(fake.f2)
> plot.map(fake.f2,newmap)
> fake.f2 <- replace.map(fake.f2, newmap)
> 
> 
> 
> cleanEx(); ..nameEx <- "ripple"
> 
> ### * ripple
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ripple
> ### Title: Compare marker orders
> ### Aliases: ripple
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(badorder)
> rip7c <- ripple(badorder, chr=1, window=3)
   42 total orders
> summary(rip7c)
                                     obligXO
Initial  1 2 3 4 5 6 7 8  9 10 11 12     415
1        1 2 3 4 5 6 7 8 11 10  9 12     361
2        1 2 3 4 6 5 7 8  9 10 11 12     363
3        1 2 3 6 4 5 7 8  9 10 11 12     367
4        1 2 3 4 5 6 7 8 11  9 10 12     374
5        1 2 3 4 5 6 7 8  9 12 10 11     391
6        1 2 3 4 5 6 7 8 10 11  9 12     395
7        1 2 3 4 6 7 5 8  9 10 11 12     401
8        1 2 3 4 5 6 7 8  9 11 10 12     404
9        1 2 3 4 5 6 7 8  9 12 11 10     406
10       1 2 3 4 5 6 7 8  9 10 12 11     407
11       1 2 3 5 6 4 7 8  9 10 11 12     411
12       1 2 3 4 5 6 7 8  9 11 12 10     412
13       1 2 3 5 4 6 7 8  9 10 11 12     413
14       1 2 3 6 5 4 7 8  9 10 11 12     413
15       2 1 3 4 5 6 7 8  9 10 11 12     414
16       1 2 3 4 5 6 7 9 10  8 11 12     415
17       1 2 3 4 5 6 7 8 10  9 11 12     417
attr(,"class")
[1] "summary.ripple" "matrix"        
> 
> ## Not run: 
> ##D rip2l <- ripple(badorder, chr=1, window=2, method="likelihood")
> ##D summary(rip2l)
> ## End(Not run)
> 
> badorder <- switch.order(badorder, 1, rip7c[2,])
> 
> 
> 
> cleanEx(); ..nameEx <- "scanone"
> 
> ### * scanone
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: scanone
> ### Title: Genome scan with a single QTL model
> ### Aliases: scanone
> ### Keywords: models
> 
> ### ** Examples
> 
> ###################
> # Normal Model
> ###################
> data(hyper)
> ## Don't show: 
> hyper <- subset(hyper,chr=1:4,ind=1:100)
> ## End Don't show
> # Genotype probabilities for EM and H-K
> hyper <- calc.genoprob(hyper, step=2.5)
> out.em <- scanone(hyper, method="em")
> out.hk <- scanone(hyper, method="hk")
> 
> # Summarize results: peaks above 3
> summary(out.em, 3)

           chr pos lod  AA AB sigma
c1.loc77.5   1  81 3.6 106 96    11
D4Mit164     4  30 3.8 106 97    11

> summary(out.hk, 3)

           chr pos lod  AA AB sigma
c1.loc77.5   1  81 3.6 106 96    11
D4Mit164     4  30 3.8 106 97    11

> 
> # Plot the results
> plot(out.hk, out.em)
> plot(out.hk, out.em, chr=c(1,4), lty=1, col=c("blue","black"))
> 
> # Imputation; first need to run sim.geno
> # Do just chromosomes 1 and 4, to save time
> hyper.c1n4 <- sim.geno(subset(hyper, chr=c(1,4)),
+                        step=2.5, n.draws=8)
> out.imp <- scanone(hyper.c1n4, method="imp")
> summary(out.imp, 3)

         chr pos lod
D1Mit14    1  82 3.6
D4Mit164   4  30 3.8

> 
> # Plot all three results
> plot(out.imp, out.hk, out.em, chr=c(1,4), lty=1,
+      col=c("red","blue","black"))
> 
> # Permutation tests
> ## Not run: 
> ##D permo <- scanone(hyper, method="hk", n.perm=1000)
> ## End(Not run)## Don't show: 
> permo <- scanone(hyper, method="hk", n.perm=3)
Permutation 1 
Permutation 2 
Permutation 3 
> ## End Don't showquantile(permo, 0.95)
> 
> ###################
> # Non-parametric
> ###################
> out.np <- scanone(hyper, model="np")
> summary(out.np, 3)

           chr pos lod
c1.loc77.5   1  81 3.2
D4Mit164     4  30 3.3

> 
> # Plot with previous results
> plot(out.np, chr=c(1,4), lty=1, col="green")
> plot(out.imp, out.hk, out.em, chr=c(1,4), lty=1,
+      col=c("red","blue","black"), add=TRUE)
> 
> ###################
> # Two-part Model
> ###################
> data(listeria)
> ## Don't show: 
> listeria <- subset(listeria,chr=c(1,5,13))
> ## End Don't show
> listeria <- calc.genoprob(listeria,step=2.5)
> out.2p <- scanone(listeria, model="2part", upper=TRUE)
> summary(out.2p, 5)

           chr pos lod lod.p lod.mu mu.AA mu.AB mu.BB pi.AA pi.AB pi.BB sigma
D1M355       1  81 5.2  0.57   4.64    85   107   129  0.21  0.30 0.394    29
c5.loc27.5   5  28 6.7  6.08   0.59   122   106   101  0.64  0.29 0.031    33
D13M147     13  26 7.0  3.66   3.36    94   124   104  0.09  0.43 0.421    30

> 
> # Plot all three LOD scores together
> plot(out.2p, out.2p, out.2p, lodcolumn=c(4,5,3), lty=1, chr=c(1,5,13),
+      col=c("red","blue","black"))
> 
> # Permutation test
> ## Not run: 
> ##D permo <- scanone(listeria, model="2part", upper=TRUE,
> ##D                  n.perm=1000)
> ## End(Not run)## Don't show: 
> permo <- scanone(listeria, model="2part", upper=TRUE,
+                  n.perm=3)
Permutation 1 
Permutation 2 
Permutation 3 
> ## End Don't showapply(permo, 2, quantile, 0.95)
> 
> ###################
> # Binary model
> ###################
> listeria <- subset(listeria, ind=!is.na(listeria$pheno[,1]))
> listeria$pheno[,2] <- rep(0,nind(listeria))
> listeria$pheno[listeria$pheno[,1]==264,2] <- 1
> out.bin <- scanone(listeria, pheno.col=2, model="binary")
> summary(out.bin, 3)

           chr pos lod   AA   AB    BB
c5.loc27.5   5  28 6.1 0.64 0.29 0.031
D13M147     13  26 3.7 0.09 0.43 0.421

> 
> # Plot LOD for binary model with LOD(p) from 2-part model
> plot(out.bin, out.2p, lodcolumn=c(3,4), lty=1, col=c("black", "red"),
+      chr=c(1,5,13))
> 
> # Permutation test
> ## Not run: 
> ##D permo <- scanone(listeria, pheno.col=2, model="binary",
> ##D                  n.perm=1000)
> ## End(Not run)## Don't show: 
> permo <- scanone(listeria, pheno.col=2, model="binary",
+                  n.perm=3)
Permutation 1 
Permutation 2 
Permutation 3 
> ## End Don't showquantile(permo, 0.95)
> 
> ###################
> # Covariates
> ###################
> data(fake.bc)
> ## Don't show: 
> fake.bc <- subset(fake.bc, chr=c(2,5,10))
> ## End Don't show
> plot(fake.bc)
> fake.bc <- calc.genoprob(fake.bc, step=2.5)
> # genome scans without covariates
> out.nocovar <- scanone(fake.bc)
> # genome scans with covariates
> ac <- fake.bc$pheno[,c("sex","age")]
> ic <- fake.bc$pheno[,"sex"]
> out.covar <- scanone(fake.bc, pheno.col=1,
+                      addcovar=ac, intcovar=ic)
> summary(out.nocovar,3)

           chr pos lod  AA  AB sigma
c2.loc35     2  35 3.5 5.8 6.6   1.7
c5.loc17.5   5  18 7.8 5.7 6.7   1.6

> summary(out.covar,3)

            chr pos  lod
c2.loc27.5    2  28 11.7
c5.loc15      5  15  9.3
c10.loc27.5  10  28  6.6

> plot(out.covar,out.nocovar,chr=c(2,5,10))
> 
> 
> 
> cleanEx(); ..nameEx <- "scanqtl"
> 
> ### * scanqtl
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: scanqtl
> ### Title: General QTL scan
> ### Aliases: scanqtl
> ### Keywords: models
> 
> ### ** Examples
> 
> data(fake.f2)
> 
> # take out several QTLs
> qc <- c(1, 8, 13)
> fake.f2 <- subset(fake.f2, chr=qc)
> 
> # imputate genotypes
> ## Don't show: 
> fake.f2 <- subset(fake.f2, ind=1:50)
> ## End Don't show
> ## Not run: 
> ##D fake.f2 <- sim.geno(fake.f2, n.draws=64, step=2, err=0.001)
> ## End(Not run)## Don't show: 
> fake.f2 <- sim.geno(fake.f2, n.draws=8, step=5, err=0.001)
> ## End Don't show
> 
> # 2-dimensional genome scan with additive 3-QTL model
> pos <- list(c(15,35), c(45,65), 28)
> result <- scanqtl(fake.f2, pheno.col=1, chr=qc, pos=pos,
+                   formula=y~Q1+Q2+Q3)
> 
> # image of the results
> # chr locations
> chr1 <- as.numeric(matrix(unlist(strsplit(colnames(result),"@")),
+                    ncol=2,byrow=TRUE)[,2])
> chr8 <- as.numeric(matrix(unlist(strsplit(rownames(result),"@")),
+                    ncol=2,byrow=TRUE)[,2])
> # image plot
> image(chr1, chr8, t(result), las=1, col=rev(rainbow(256,start=0,end=2/3)))
> 
> # do the same, allowing the QTLs on chr 1 and 13 to interact
> result2 <- scanqtl(fake.f2, pheno.col=1, chr=qc, pos=pos,
+                    formula=y~Q1+Q2+Q3+Q1:Q3)
> # image plot
> image(chr1, chr8, t(result2), las=1, col=rev(rainbow(256,start=0,end=2/3)))
> 
> 
> 
> cleanEx(); ..nameEx <- "scantwo"
> 
> ### * scantwo
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: scantwo
> ### Title: Two-dimensional genome scan with a two-QTL model
> ### Aliases: scantwo
> ### Keywords: models
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2, chr=18:19)
> ## End Don't show
> fake.f2 <- calc.genoprob(fake.f2, step=10)
> out.2dim <- scantwo(fake.f2, method="hk")
 --Running scanone
 --Running scantwo
 (18,18)
 (18,19)
 (19,19)
> plot(out.2dim)
> 
> ## Don't show: 
> permo.2dim <- scantwo(fake.f2, method="hk", n.perm=2)
Permutation 1 
Permutation 2 
> ## End Don't show## Not run: 
> ##D permo.2dim <- scantwo(fake.f2, method="hk", n.perm=1000)
> ## End(Not run)apply(permo.2dim,2,quantile,0.95)
> 
> # covariates
> data(fake.bc)
> ## Don't show: 
> fake.bc <- subset(fake.bc, chr=16:17)
> ## End Don't show
> fake.bc <- calc.genoprob(fake.bc, step=10)
> ac <- fake.bc$pheno[,c("sex","age")]
> ic <- fake.bc$pheno[,"sex"]
> out <- scantwo(fake.bc, method="hk", pheno.col=1,
+                addcovar=ac, intcovar=ic)
 --Running scanone
 --Running scantwo
 (16,16)
 (16,17)
 (17,17)
> plot(out)
> 
> 
> 
> cleanEx(); ..nameEx <- "sim.cross"
> 
> ### * sim.cross
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.cross
> ### Title: Simulate a QTL experiment
> ### Aliases: sim.cross
> ### Keywords: datagen
> 
> ### ** Examples
> 
> # simulate a genetic map
> map <- sim.map()
> 
> # simulate 250 intercross individuals with 2 QTLs
> fake <- sim.cross(map, type="f2", n.ind=250,
+                   model = rbind(c(1,45,1,1),c(5,20,0.5,-0.5)))
> 
> 
> 
> cleanEx(); ..nameEx <- "sim.geno"
> 
> ### * sim.geno
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.geno
> ### Title: Simulate genotypes given observed marker data
> ### Aliases: sim.geno
> ### Keywords: utilities
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2, chr=18:19)
> ## End Don't show
> fake.f2 <- sim.geno(fake.f2, step=2, n.draws=8)
> 
> 
> 
> cleanEx(); ..nameEx <- "sim.map"
> 
> ### * sim.map
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sim.map
> ### Title: Simulate a genetic map
> ### Aliases: sim.map
> ### Keywords: datagen
> 
> ### ** Examples
> 
> # simulate 4 autosomes, each with 10 markers
> map <- sim.map(c(100,90,80,40), 10, include.x=FALSE)
> plot.map(map)
> 
> # equally spaced markers
> map2 <- sim.map(c(100,90,80,40), 10, include.x=FALSE, eq.spacing=TRUE)
> plot(map2)
> 
> 
> 
> cleanEx(); ..nameEx <- "strip.partials"
> 
> ### * strip.partials
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: strip.partials
> ### Title: Strip partially informative genotypes
> ### Aliases: strip.partials
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(listeria)
> sum(nmissing(listeria))
[1] 1840
> listeria <- strip.partials(listeria)
Removed 128 genotypes.
> sum(nmissing(listeria))
[1] 1968
> 
> 
> 
> cleanEx(); ..nameEx <- "subset.cross"
> 
> ### * subset.cross
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: subset.cross
> ### Title: Subsetting data for QTL experiment
> ### Aliases: subset.cross
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(fake.f2)
> fake.f2.A <- subset(fake.f2, chr=c("5","13"))
> fake.f2.B <- subset(fake.f2, ind = -c(1,5,10))
> fake.f2.C <- subset(fake.f2, chr=1:5, ind=1:50)
> 
> data(listeria)
> y <- listeria$pheno[,1]
> listeriaB <- subset(listeria, ind = (!is.na(y) & y < 264))
> 
> 
> 
> cleanEx(); ..nameEx <- "summary.cross"
> 
> ### * summary.cross
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.cross
> ### Title: Print summary of QTL experiment
> ### Aliases: summary.cross nind nmar totmar nchr nphe
> ### Keywords: print
> 
> ### ** Examples
> 
> data(fake.f2)
> summary(fake.f2)
    F2 intercross

    No. individuals:  200 

    No. phenotypes:   3 
    Percent phenotyped:  100 100 100 

    No. chromosomes:  20 
    Total markers:    94 
    No. markers:      7 5 3 4 10 5 4 4 6 4 4 5 11 3 3 4 4 2 3 3 
    Percent genotyped:  85 
    Genotypes (%):      AA:25.9  AB:50.2  BB:23.9  not BB:0  not AA:0 
> nind(fake.f2)
[1] 200
> nmar(fake.f2)
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19  X 
 7  5  3  4 10  5  4  4  6  4  4  5 11  3  3  4  4  2  3  3 
> totmar(fake.f2)
[1] 94
> nchr(fake.f2)
[1] 20
> nphe(fake.f2)
[1] 3
> 
> 
> 
> cleanEx(); ..nameEx <- "summary.fitqtl"
> 
> ### * summary.fitqtl
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.fitqtl
> ### Title: Summary of fit of qtl model
> ### Aliases: summary.fitqtl
> ### Keywords: models
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "summary.map"
> 
> ### * summary.map
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.map
> ### Title: Print summary of a genetic map
> ### Aliases: summary.map
> ### Keywords: print
> 
> ### ** Examples
> 
> data(map10)
> summary(map10)
Sex-averaged map

        n.mar length ave.spacing
1          11  102.9        10.3
2           9   88.5        11.1
3           8   78.1        11.2
4          10   90.8        10.1
5           9   83.1        10.4
6           7   67.4        11.2
7           8   72.4        10.3
8           7   68.4        11.4
9           8   70.7        10.1
10          7   64.3        10.7
11          7   69.1        11.5
12          6   59.0        11.8
13          7   67.7        11.3
14          7   63.9        10.7
15          6   56.6        11.3
16          7   62.9        10.5
17          6   50.6        10.1
18          5   41.7        10.4
19          6   53.3        10.7
X           7   59.0         9.8
overall   148 1370.4        10.7
> 
> 
> 
> cleanEx(); ..nameEx <- "summary.ripple"
> 
> ### * summary.ripple
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.ripple
> ### Title: Print summary of ripple results
> ### Aliases: summary.ripple
> ### Keywords: print
> 
> ### ** Examples
>  
> ## Not run: data(badorder)
> ##D rip7c <- ripple(badorder, 1, 7)
> ##D summary(rip7c)
> ##D 
> ##D rip2l <- ripple(badorder, 1, 2, method="likelihood")
> ##D summary(rip2l)
> ##D 
> ##D badorder <- switch.order(badorder, 1, rip7c[2,])
> ## End(Not run)
> 
> 
> cleanEx(); ..nameEx <- "summary.scanone"
> 
> ### * summary.scanone
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.scanone
> ### Title: Summarize the results of one or several one-QTL genome scans
> ### Aliases: summary.scanone
> ### Keywords: print
> 
> ### ** Examples
> 
> data(fake.f2)
> fake.f2 <- calc.genoprob(fake.f2)
> # scanones on observed data
> scan1 <- scanone(fake.f2, pheno.col=1)
> # permutation tests
> ## Not run: 
> ##D out1 <- scanone(fake.f2, pheno.col=1, n.perm=100)
> ## End(Not run)## Don't show: 
> out1 <- scanone(fake.f2, pheno.col=1, n.perm=3)
Permutation 1 
Permutation 2 
Permutation 3 
> ## End Don't show
> ###### summary one result #######
> # one sig level
> summary(scan1, threshold=quantile(out1, 0.95))

        chr pos lod AA AB ABf ABr AY BB BY sigma
D1M437    1  37 4.3 22 24  NA  NA NA 26 NA   4.6
D13M254  13  24 7.9 20 25  NA  NA NA 25 NA   4.4
DXM66     X  14 7.0 21 NA  24  25 21 26 26   4.3

> # multiple sig level
> summary(scan1, threshold=quantile(out1, c(0.63,0.95,0.99)))

        chr pos lod AA AB ABf ABr AY BB BY sigma Sig
D1M437    1  37 4.3 22 24  NA  NA NA 26 NA   4.6 ***
D8M42     8  61 2.6 22 24  NA  NA NA 25 NA   4.7   *
D13M254  13  24 7.9 20 25  NA  NA NA 25 NA   4.4 ***
DXM66     X  14 7.0 21 NA  24  25 21 26 26   4.3 ***

> 
> # scan the other phenotype
> scan2 <- scanone(fake.f2, pheno.col=2)
> ## Not run: 
> ##D out2 <- scanone(fake.f2, pheno.col=2, n.perm=100)
> ## End(Not run)## Don't show: 
> out2 <- scanone(fake.f2, pheno.col=2, n.perm=3)
Permutation 1 
Permutation 2 
Permutation 3 
> ## End Don't show
> 
> ##### summary multiple results
> # one sig level
> summary(scan1,scan2, threshold=2)

        chr pos lod1 Sig1 lod2 Sig2
D1M437    1  37  4.3    * 0.18     
D8M42     8  61  2.6    * 0.24     
D13M254  13  24  7.9    * 0.28     
DXM66     X  14  7.0    * 0.00     

> # multiple sig level, same threshold
> summary(scan1,
+         scan2,threshold=matrix(quantile(out1,c(0.37,0.95,0.99)),ncol=1))   

        chr pos lod1 Sig1 lod2 Sig2
D1M437    1  37  4.3  *** 0.18     
D8M42     8  61  2.6    * 0.24     
D13M254  13  24  7.9  *** 0.28     
DXM66     X  14  7.0  *** 0.00     

> # one sig level, different threshold
> summary(scan1, scan2,threshold=c(3,0.8))

        chr  pos  lod1 Sig1 lod2 Sig2
D1M437    1 37.1 4.305    * 0.18     
D4M310    4 52.4 0.022      0.87    *
D10M166  10  0.0 0.951      0.87    *
D12M52   12 25.1 0.235      0.88    *
D13M254  13 24.0 7.932    * 0.28     
D15M42   15 42.6 0.392      1.08    *
D16M87   16  4.4 0.116      1.21    *
D19M34   19 31.7 0.360      0.87    *
DXM66     X 14.2 6.996    * 0.00     

> # multiple sig level, different threshold
> threshold <- matrix(c(quantile(out1,probs=c(0.37,0.95,0.99)),0.5,0.8,1.2), ncol=2) 
> summary(scan1, scan2, threshold=threshold)

        chr  pos  lod1 Sig1 lod2 Sig2
D1M437    1 37.1 4.305  *** 0.18     
D2M241    2 21.8 0.878      0.73    *
D4M310    4 52.4 0.022      0.87   **
D6M194    6 48.1 0.288      0.75    *
D7M190    7  0.0 0.312      0.64    *
D8M42     8 61.2 2.554    * 0.24     
D9M18     9 48.1 0.042      0.66    *
D10M166  10  0.0 0.951      0.87   **
D11M20   11  0.0 0.668      0.68    *
D12M52   12 25.1 0.235      0.88   **
D13M254  13 24.0 7.932  *** 0.28     
D14M142  14 12.0 0.253      0.69    *
D15M42   15 42.6 0.392      1.08   **
D16M87   16  4.4 0.116      1.21  ***
D19M34   19 31.7 0.360      0.87   **
DXM66     X 14.2 6.996  *** 0.00     

> 
> 
> 
> cleanEx(); ..nameEx <- "summary.scantwo"
> 
> ### * summary.scantwo
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.scantwo
> ### Title: Summarize the results of a two-dimensional genome scan
> ### Aliases: summary.scantwo
> ### Keywords: print
> 
> ### ** Examples
> 
> data(fake.f2)
> ## Don't show: 
> fake.f2 <- subset(fake.f2, chr=c(1,6,13))
> ## End Don't show
> fake.f2 <- calc.genoprob(fake.f2, step=5)
> out.2dim <- scantwo(fake.f2, method="hk")
 --Running scanone
 --Running scantwo
 (1,1)
 (1,6)
 (1,13)
 (6,6)
 (6,13)
 (13,13)
> # All pairs of loci
> summary(out.2dim)

        pos1 pos2   LODjnt -logP   LODint -logP   LODq1 -logP   LODq2 -logP
 c1: c1    5   25     6.49   3.7     0.82   0.4    0.59   0.6    1.30   1.3
 c1: c6   25    5     8.14   5.0     2.98   2.1    5.13   5.1    0.08   0.1
 c1:c13   30   30    13.86  10.1     0.78   0.3    4.27   4.3    8.04   8.0
 c6: c6   40   45     1.96   0.5     1.34   0.7    0.27   0.3    0.15   0.2
 c6:c13   35   30     9.76   6.4     0.50   0.2    0.45   0.5    8.77   8.8
c13:c13    0   40    11.71   8.1     3.79   2.8    0.63   0.6    4.58   4.6

> # Pairs meeting specified criteria
> summary(out.2dim, c(7, 3, 3))

        pos1 pos2   LODjnt -logP   LODint -logP   LODq1 -logP   LODq2 -logP
 c1:c13   30   30    13.86  10.1     0.78   0.3    4.27   4.3    8.04   8.0
c13:c13    0   40    11.71   8.1     3.79   2.8    0.63   0.6    4.58   4.6

> # Pairs with both conditional LODs > 2
> summary(out.2dim,c(0,1000,2))

       pos1 pos2   LODjnt -logP   LODint -logP   LODq1 -logP   LODq2 -logP
c1:c13   30   30    13.86  10.1     0.78   0.3    4.27   4.3    8.04     8

> # Pairs with interaction LOD is above 3
> summary(out.2dim,c(0,3,1000))

        pos1 pos2   LODjnt -logP   LODint -logP   LODq1 -logP   LODq2 -logP
c13:c13    0   40    11.71   8.1     3.79   2.8    0.63   0.6    4.58   4.6

> 
> 
> 
> cleanEx(); ..nameEx <- "switch.order"
> 
> ### * switch.order
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: switch.order
> ### Title: Switch the order of markers on a chromosome
> ### Aliases: switch.order
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(fake.f2)
> fake.f2 <- switch.order(fake.f2, 1, c(1,3,2,4:7))
> 
> 
> 
> cleanEx(); ..nameEx <- "top.errorlod"
> 
> ### * top.errorlod
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: top.errorlod
> ### Title: List genotypes with large error LOD scores
> ### Aliases: top.errorlod
> ### Keywords: manip
> 
> ### ** Examples
> 
> data(hyper)
> ## Don't show: 
> hyper <- subset(hyper, chr=c(1,19))
> ## End Don't show
> # run calc.genoprob, then calc.errorlod; use same error.prob
> hyper <- calc.genoprob(hyper,error.prob=0.01)
> hyper <- calc.errorlod(hyper,error.prob=0.01)
> 
> # plot the error LOD scores; print those above a specified cutoff 
> plot.errorlod(hyper)
> plot.errorlod(hyper,chr=1)
> top.errorlod(hyper,cutoff=4)
   chr ind   marker errorlod
1    1 118  D1Mit14 8.372802
2    1 162  D1Mit14 8.372802
3    1 170  D1Mit14 8.372802
4    1 159  D1Mit14 8.350349
5    1  73  D1Mit14 6.165405
6    1  65  D1Mit14 6.165405
7    1  88  D1Mit14 6.165405
8    1 184  D1Mit14 6.151614
9    1 241  D1Mit14 6.151614
10   1 215 D1Mit267 5.822197
11   1 108 D1Mit267 5.822197
12   1 138 D1Mit267 5.822197
13   1 226 D1Mit267 5.822197
14   1 199 D1Mit267 5.819255
15   1  84 D1Mit267 5.808406
> 
> # plot genotype data, flagging genotypes with error LOD > cutoff
> plot.geno(hyper, chr=1, ind=160:200, cutoff=7, min.sep=2)
> 
> # Note that if markers are extremely tightly linked,
> #     any recombination event will have a high error LOD
> plot.geno(hyper, chr=1, ind=101:118, min.sep=2)
> hyper$geno[[1]]$map
D1Mit296 D1Mit123 D1Mit156 D1Mit178  D1Mit19   D1Mit7  D1Mit46 D1Mit132 
     3.3     19.7     32.8     35.0     37.2     41.5     43.7     43.7 
D1Mit334 D1Mit305  D1Mit26  D1Mit94 D1Mit218 D1Mit100 D1Mit102  D1Mit14 
    49.2     54.6     64.5     67.8     69.9     74.3     75.4     82.0 
D1Mit105 D1Mit159 D1Mit267  D1Mit15 D1Mit456 D1Mit155 
    82.0     82.0     82.0     86.3     94.0    115.8 
> 
> 
> 
> cleanEx(); ..nameEx <- "write.cross"
> 
> ### * write.cross
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: write.cross
> ### Title: Write data for a QTL experiment to a file
> ### Aliases: write.cross
> ### Keywords: IO
> 
> ### ** Examples
> 
> ## Not run: data(fake.bc)
> ##D 
> ##D # comma-delimited format
> ##D write.cross(fake.bc, "csv", "Data/fakebc", c(1,5,13))
> ##D 
> ##D # Mapmaker format
> ##D write.cross(fake.bc, "mm", "Data/fakebc", c(1,5,13))
> ##D 
> ##D # QTL Cartographer format
> ##D write.cross(fake.bc, "qtlcart", "Data/fakebc", c(1,5,13))
> ##D 
> ##D # Gary's format
> ##D write.cross(fake.bc, "gary", c(1,5,13))
> ## End(Not run) 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  29.92 4.6 37.9 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')