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. > ### *
> ### > 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("bqtl-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('bqtl') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "bqtl" > > ### * bqtl > > flush(stderr()); flush(stdout()) > > ### Name: bqtl > ### Title: Bayesian QTL Model Fitting > ### Aliases: bqtl > ### Keywords: regression > > ### ** Examples > > > data(little.ana.bc ) # load BC1 dataset > > loglik( bqtl( bc.phenotype ~ 1, little.ana.bc ) ) #null loglikelihood [1] -191.5687 > #on chr 1 near cM 25 > loglik(bqtl(bc.phenotype~locus(chromo=1,cM=25),little.ana.bc)) [1] -184.0409 > > little.bqtl <- # two genes with epistasis + bqtl(bc.phenotype ~ m.12 * m.24, little.ana.bc) > summary(little.bqtl) $coefficients Intercept m.12 m.24 m.12:m.24 -0.074264764 -0.017991584 -0.002598244 -0.036974458 $loglik [1] -191.4156 $std.res [1] 0.8668414 $N N N.omit N.used 150 0 150 > > several.epi <- # 20 epistatic models + bqtl( bc.phenotype ~ m.12 * locus(31:50), little.ana.bc) > several.main <- # main effects only + bqtl( bc.phenotype ~ m.12 + locus(31:50), little.ana.bc) > > max.loglik <- max( loglik(several.epi) - loglik(several.main) ) > > round( + c( Chi.Square=2*max.loglik, df=1, p.value=1-pchisq(2*max.loglik,1)) + ,2) Chi.Square df p.value 4.21 1.00 0.04 > > five.gene <- ## a five gene model + bqtl( bc.phenotype ~ locus( 12, 32, 44, 22, 76 ), little.ana.bc , return.hess=TRUE ) > > regr.coef.table <- summary(five.gene)$coefficients > > round( regr.coef.table[,"Value"] + # coefs inside 95% CI + qnorm(0.025) * regr.coef.table[,"Std.Err"] %o% + c("Lower CI"=1,"Estimate"=0,"Upper CI"=-1),3) Lower CI Estimate Upper CI Intercept -0.188 -0.056 0.075 m.5 0.145 0.283 0.421 m.14 -0.133 -0.003 0.128 m.21 -0.176 -0.045 0.086 m.9 -0.113 0.026 0.165 m.34 0.007 0.138 0.268 > > ## Don't show: > rm(five.gene,little.ana.bc,little.bqtl,max.loglik, + regr.coef.table,several.epi,several.main) > ## End Don't show > > > > cleanEx(); ..nameEx <- "linear.bayes" > > ### * linear.bayes > > flush(stderr()); flush(stdout()) > > ### Name: linear.bayes > ### Title: Bayesian QTL mapping via Linearized Likelihood > ### Aliases: linear.bayes > ### Keywords: regression > > ### ** Examples > > data( little.ana.bc ) > little.lin <- linear.bayes( bc.phenotype~locus(all), little.ana.bc, rparm=1 ) > par(mfrow=c(2,3)) > plot( little.ana.bc, little.lin$loc.posterior, type="h" ) > little.lin$odds [1] 1126.133 1257600.980 551610.008 194460.858 60312.782 > par(mfrow=c(1,1)) > plot(fitted(little.lin), residuals(little.lin)) > ## Don't show: > rm(little.lin,little.ana.bc) > ## End Don't show > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "make.analysis.obj" > > ### * make.analysis.obj > > flush(stderr()); flush(stdout()) > > ### Name: make.analysis.obj > ### Title: Set up data for QTL mapping > ### Aliases: make.analysis.obj > ### Keywords: manip > > ### ** Examples > > data( little.bc.pheno ) > data( little.mf.5 ) > data( little.bc.markers ) > names(little.bc.pheno) [1] "bc.phenotype" > little.ana.bc <- make.analysis.obj(little.bc.pheno$bc.phenotype, + little.mf.5,little.bc.markers, + method="BC1") > summary( little.ana.bc ) $map marker pseudo-marker cM.low cM.high 1 10 13 0 84.66 2 10 10 0 70.35 3 10 13 0 84.84 4 10 16 0 99.53 5 10 12 0 89.06 $N [1] 150 $method [1] "BC1" $variable.names [1] "data" $mode.mat [,1] AA 1 Aa -1 $call make.analysis.obj(data = little.bc.pheno$bc.phenotype, map.frame = little.mf.5, marker.frame = little.bc.markers, method = "BC1") > > ## Don't show: > rm( little.ana.bc, little.bc.pheno, little.mf.5, little.bc.markers ) > ## End Don't show > > > > cleanEx(); ..nameEx <- "make.loc.right" > > ### * make.loc.right > > flush(stderr()); flush(stdout()) > > ### Name: make.loc.right > ### Title: Keep track of fully informative markers or states > ### Aliases: make.loc.right > ### Keywords: manip > > ### ** Examples > > > > > > cleanEx(); ..nameEx <- "make.map.frame" > > ### * make.map.frame > > flush(stderr()); flush(stdout()) > > ### Name: make.map.frame > ### Title: Create marker map specifications > ### Aliases: make.map.frame > ### Keywords: manip > > ### ** Examples > > > data( little.map.dx ) > little.map.frame <- make.map.frame( little.map.dx ) > plot( little.map.frame ) # there is a plot method > # add 'virtual' markers to map > little.mf.5 <- make.map.frame(little.map.frame,reso=5) > print(little.mf.5[1:10,],digits=1) # show a few rows marker.name cM prior pos.type is.marker pos.plot lambda locus chr.num 1 m.1 0.0 0.003 2 TRUE 0.0 1 C.1.0 1 2 m.2 0.9 0.002 3 TRUE 0.9 1 C.1.0.86 1 3 m.3 1.9 0.006 3 TRUE 1.9 1 C.1.1.94 1 4 m.3.1 6.4 0.010 3 FALSE 6.4 1 C.1.6.361 1 5 m.3.2 10.8 0.010 3 FALSE 10.8 1 C.1.10.783 1 6 m.3.3 15.2 0.010 3 FALSE 15.2 1 C.1.15.204 1 7 m.3.4 19.6 0.010 3 FALSE 19.6 1 C.1.19.625 1 8 m.3.5 24.0 0.010 3 FALSE 24.0 1 C.1.24.046 1 9 m.3.6 28.5 0.010 3 FALSE 28.5 1 C.1.28.468 1 10 m.3.7 32.9 0.010 3 FALSE 32.9 1 C.1.32.889 1 > plot( little.mf.5 ) # notice the 'virtual' markers added > ## Don't show: > rm( little.map.dx,little.map.frame ,little.mf.5) > ## End Don't show > > > > cleanEx(); ..nameEx <- "map.index" > > ### * map.index > > flush(stderr()); flush(stdout()) > > ### Name: map.index > ### Title: Look up numerical index(es) of map locations > ### Aliases: map.index map.index.default map.index.analysis.object > ### Keywords: methods > > ### ** Examples > > data(little.ana.bc) > map.index(little.ana.bc,chromo=1,cM=25) # locus nearest 1,25 [1] 8 > index.chr.1 <- map.index(little.ana.bc,chromo=1) > fit.on.1 <- bqtl(bc.phenotype~locus(index.chr.1),little.ana.bc) > summary( loglik( fit.on.1 ) ) Min. 1st Qu. Median Mean 3rd Qu. Max. -189.9 -188.9 -185.0 -185.8 -183.0 -181.9 > > > > > cleanEx(); ..nameEx <- "map.location" > > ### * map.location > > flush(stderr()); flush(stdout()) > > ### Name: map.location > ### Title: Report map location > ### Aliases: map.loc map.location map.location.default > ### map.location.analysis.object map.location.bqtl map.location.bqtl.list > ### Keywords: regression > > ### ** Examples > > > data(little.ana.bc) > > map.loc(little.ana.bc, c(1,15,45)) chr.num cM marker.name 1 1 0.000 m.1 15 1 49.244 m.7.1 45 3 4.450 m.22 > map.loc(little.ana.bc,chromo=3,cM=22) chr.num cM marker.name 50 3 21.46667 m.25.1 > map.loc(little.ana.bc,"m.12") chr.num cM marker.name 27 2 12.22 m.12 > rm(little.ana.bc) > > > > cleanEx(); ..nameEx <- "map.names" > > ### * map.names > > flush(stderr()); flush(stdout()) > > ### Name: map.names > ### Title: Look up names of markers or loci > ### Aliases: map.names map.names.default map.names.bqtl map.names.bqtl.list > ### map.names.analysis.object map.names.map.frame > ### Keywords: utilities > > ### ** Examples > > data(little.ana.bc) > > map.names(little.ana.bc,chromo=1,cM=24) [1] "m.3.5" > > map.names(little.ana.bc,chromo=c(1,1),cM=c(40,55)) [1] "m.6" "m.7" "m.7.1" "m.7.2" > > fit <- bqtl( bc.phenotype ~ locus(23,42) , little.ana.bc ) > > map.names( fit ) [1] "m.10" "m.19.1" > > ## Don't show: > rm(fit, little.ana.bc) > ## End Don't show > > > > cleanEx(); ..nameEx <- "marker.fill" > > ### * marker.fill > > flush(stderr()); flush(stdout()) > > ### Name: marker.fill > ### Title: Map Positions Between Markers > ### Aliases: marker.fill > ### Keywords: utilities > > ### ** Examples > > > data( little.map.frame ) > little.nint <- marker.fill( little.map.frame, reso=5, TRUE ) > cbind(nint=little.nint,cM=little.map.frame$cM)[1:10,] nint cM m.1 1 0.00 m.2 1 0.86 m.3 8 1.94 m.4 1 37.31 m.5 1 38.26 m.6 1 39.95 m.7 5 44.66 m.8 3 67.58 m.9 1 80.69 m.10 1 84.66 > rm( little.map.frame, little.nint ) > > > > cleanEx(); ..nameEx <- "marker.levels" > > ### * marker.levels > > flush(stderr()); flush(stdout()) > > ### Name: marker.levels > ### Title: Define marker level codes > ### Aliases: marker.levels bc1.levels f2.levels ri.levels > ### Keywords: manip > > ### ** Examples > > > ### show the defaults: > > f2.levels() AA Aa aa not.aa not.AA miss.val "AA" "Aa" "aa" "A-" "a-" "--" > bc1.levels() AA Aa miss.val "AA" "Aa" "nil" "nil" "nil" "--" > ri.levels() AA aa miss.val "AA" "aa" "nil" "nil" "nil" "--" > > ### suppose that 1,2,3 are codes used in F2: > > f2.levels(1,2,3) AA Aa aa not.aa not.AA miss.val "1" "2" "3" "A-" "a-" "--" > > ### show what would happen changing "Aa" to "H" > > f2.levels(Aa="H") AA Aa aa not.aa not.AA miss.val "AA" "H" "aa" "A-" "a-" "--" > bc1.levels(Aa="H") AA Aa miss.val "AA" "H" "nil" "nil" "nil" "--" > > > > > cleanEx(); ..nameEx <- "plot.map.frame" > > ### * plot.map.frame > > flush(stderr()); flush(stdout()) > > ### Name: plot.map.frame > ### Title: plots by chromosome location > ### Aliases: plot.map.frame plot.analysis.object > ### Keywords: hplot > > ### ** Examples > > > data( little.ana.bc ) > null.llk <- loglik(bqtl(bc.phenotype~1,little.ana.bc)) > llk <- loglik( bqtl( bc.phenotype~locus(all), little.ana.bc) ) - null.llk > .old.par <- par(mfrow=c(2,3)) > plot.map.frame(little.ana.bc$map.frame,llk) > par(.old.par) > ## Don't show: > rm(null.llk, llk, little.ana.bc, .old.par ) > ## End Don't show > > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "predict.bqtl" > > ### * predict.bqtl > > flush(stderr()); flush(stdout()) > > ### Name: predict.bqtl > ### Title: fitted values from QTL models > ### Aliases: predict.bqtl fitted.bqtl > ### Keywords: methods > > ### ** Examples > > > data(little.ana.bc) > > fit.pheno <- bqtl(bc.phenotype~locus(15)+locus(42),little.ana.bc) > > summary(predict(fit.pheno)) fitted.value Min. :-0.47383 1st Qu.:-0.35869 Median :-0.24354 Mean :-0.08135 3rd Qu.: 0.19673 Max. : 0.42702 > > genotype.grid <- expand.grid( c(-1,1), c(-1,1) ) # set up a grid > names(genotype.grid) <- map.names( fit.pheno ) # use matching names > > fit.vals <- predict( fit.pheno, genotype.grid ) # make predictions > cbind( genotype.grid, fit.vals ) # print them! m.7.1 m.19.1 fitted.value 1 -1 -1 -0.2491288 2 1 -1 0.4332107 3 -1 1 -0.4800235 4 1 1 0.2023160 > > > > > cleanEx(); ..nameEx <- "residuals.bqtl" > > ### * residuals.bqtl > > flush(stderr()); flush(stdout()) > > ### Name: residuals.bqtl > ### Title: Residuals from QTL models > ### Aliases: residuals.bqtl > ### Keywords: methods > > ### ** Examples > > > data(little.ana.bc) > > fit.pheno <- bqtl(bc.phenotype~locus(15)+locus(42),little.ana.bc) > > summary(residuals(fit.pheno)) object Min. :-2.491460 1st Qu.:-0.512210 Median :-0.056313 Mean : 0.004279 3rd Qu.: 0.562342 Max. : 2.235045 > > plot( fitted( fit.pheno ), residuals( fit.pheno) ) > > ## Don't show: > rm(little.ana.bc,fit.pheno) > ## End Don't show > > > > > cleanEx(); ..nameEx <- "summary.bqtl" > > ### * summary.bqtl > > flush(stderr()); flush(stdout()) > > ### Name: summary.bqtl > ### Title: Summarize bqtl object > ### Aliases: summary.bqtl > ### Keywords: methods > > ### ** Examples > > data(little.ana.bc) > fit <- bqtl( bc.phenotype~locus(4)*locus(45), little.ana.bc, + return.hess=TRUE ) > summary(fit) $coefficients Value Std.Err t.statistic p.value Intercept -0.066240519 0.06965402 -0.95099344 0.341607707 m.3.1 0.184639887 0.06965402 2.65081437 0.008029796 m.22 0.003462716 0.06965402 0.04971307 0.960351038 m.3.1:m.22 -0.061002566 0.06965402 -0.87579386 0.381142107 $loglik [1] -188.5415 $std.res [1] 0.8460853 $N N N.omit N.used 150 0 150 > ## Don't show: > rm(little.ana.bc,fit) > ## End Don't show > > > > cleanEx(); ..nameEx <- "swap" > > ### * swap > > flush(stderr()); flush(stdout()) > > ### Name: swap > ### Title: MCMC sampling of multigene models > ### Aliases: swap > ### Keywords: regression > > ### ** Examples > > data( little.ana.bc ) > little.vc <- varcov( bc.phenotype~locus(all), little.ana.bc) > little.4 <- swap( little.vc, c(1,15,55,75), rparm=1, 50, little.ana.bc ) > little.4.smry <- summary( little.4 ) > print(c("Bayes Factor (3 vs 4)"=little.4.smry$ratio$mean)) Bayes Factor (3 vs 4) 3.036744 > par(mfrow=c(3,2)) > plot( little.ana.bc, little.4.smry$loc.posterior, type="h", + ylab="E(genes)" ) > rm(little.4,little.vc,little.ana.bc) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "twohk" > > ### * twohk > > flush(stderr()); flush(stdout()) > > ### Name: twohk > ### Title: One and Two Gene Models Using Linearized Posterior > ### Aliases: twohk > ### Keywords: models > > ### ** Examples > > data(little.ana.bc) > little.vc<-make.varcov(little.ana.bc$data[,little.ana.bc$reg.names], + little.ana.bc$data$bc.phenotype) > little.2<- twohk(little.vc,little.ana.bc,rparm=1) > print( c(odds.1=sum(little.2$loc.1),odds.2=sum(little.2$loc.2)) ) odds.1 odds.2 1126.133 1257600.980 > par(mfrow=c(3,2)) > little.pe <- 2 * little.2$loc.2 / sum(little.2$loc.2) #locus-wise posterior expectation > plot(little.ana.bc,little.pe,type="h",ylab="E(genes") > rm(little.2,little.vc,little.pe,little.ana.bc) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "update.bqtl" > > ### * update.bqtl > > flush(stderr()); flush(stdout()) > > ### Name: bqtl.fitter > ### Title: Get loglikelihoods for many models of a common form > ### Aliases: bqtl.fitter > ### Keywords: models > > ### ** Examples > > data( little.ana.bc ) > little.setup <- + bqtl( bc.phenotype~locus(1)*locus(2), little.ana.bc, setup=TRUE ) > combos <- t( as.matrix( expand.grid( 1:21, 44:64 ) ) ) > little.update <- bqtl.fitter(little.setup, combos, little.ana.bc) > little.res <- matrix( little.update, nr=21 ) > image( 1:21, 44:64, little.res ) > rm(little.ana.bc, little.update, little.res ) > > > > ### *