version 1.0-2

CHANGES

@@ -0,0 +1,28 @@
+
+Changes to package bqtl - Bayesian QTL mapping toolkit 
+
+
+Version: 1.0
+Date: 2001-02-28
+
+Initial Release
+
+----------------------------------------
+
+Version: 1.0-1
+Date: 2001-03-09
+
+added keywords
+coerce make.map.frame()$marker.name to be character
+
+----------------------------------------
+
+Version: 1.0-2
+Date: 2001-05-03
+
+bugfixes - update.bqtl(), lapadj()
+
+locus( chromo=1, cM=25 ) will now work 
+
+----------------------------------------
+

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: bqtl
-Version: 1.0-1
+Version: 1.0-2
 Date: 2001-03-09
 Title: Bayesian QTL mapping toolkit 
 Author: Charles C. Berry <cberry@ucsd.edu>

R/add.s

@@ -4,5 +4,7 @@
     if (missing(scope)) stop("missing scope arg - possible bqtl syntax error")
     new.scope <- scope[grep("^add",scope)]
     if (length(new.scope)==0) stop("no 'add' terms found")
-    configs(x,...,scope=new.scope,method=method)
+    locus(x,...,
+          scope=new.scope,
+          method="fooey")    # this forces one term per locus
 }

R/bqtl.s

@@ -4,24 +4,37 @@
 {
     local.covar <-
         function (x, ..., scope = scope, method = method,
-                  prefix = NULL, bq.spec=bqtl.specials) 
+                  prefix = NULL, bq.spec=bqtl.specials,chromo=NULL,cM=NULL) 
         {
             if (missing(scope)) 
                 stop("missing scope arg - possible bqtl syntax error")
+            dots <- list(...)
+            if (is.null(chromo)&& length(names(dots))!=0 &&
+                any(cloc <- which(1==pmatch(names(dots),"chromo",0)))){
+                chromo <- dots[[cloc]]
+                if (length(dots)>1)
+                    stop("... not allowed")
+                else
+                    dots<-NULL
+            }
             x.arg <- match.call()$x
-            if  ( length(x.arg) ==1 && deparse(x.arg) == "." )  # allow . shorthand 
+            if  ( length(x.arg) ==1 && deparse(x.arg) == "all" )  # allow 'all' shorthand 
                 x <- seq(along=scope)
             if (is.call(x.arg)){
                 if (deparse(x.arg[[1]])%in%bq.spec){
                     if (is.null(x.arg$scope)) x.arg$scope <- as.name("scope")
                     if (is.null(x.arg$method)) x.arg$method <- method
+                    if (is.null(x.arg$chromo)) x.arg$chromo <- chromo
+                    if (is.null(x.arg$cM)) x.arg$cM <- cM
+
                     if (!missing(...)) stop("cannot use ... in this context")
                     paste("covar(",eval(x.arg),")",sep="")
                 }
                 else{ ## x.arg is c(1,2) or 7:8 or whatever
                     x <- eval(x.arg)
                     new.scope <- paste("covar(", scope, ")", sep = "")
-                    locus(x, ..., scope = new.scope, method = method)
+                    locus(x, ..., scope = new.scope, method = method,
+                          chromo=chromo,cM=cM)
                 }
             }
             else {
@@ -34,7 +47,8 @@
                 }
                 else {
                     new.scope <- paste("covar(", scope, ")", sep = "")
-                    locus(x, ..., scope = new.scope, method = method)
+                    locus(x, ..., scope = new.scope, method = method,
+                          chromo=chromo,cM=cM)
                 }
             }
         }
@@ -90,6 +104,8 @@
                        x$scope <- as.name("scope")
                    if ( !is.element("method",names(x)) ) #typically use default
                        x$method <- method
+                  if ( !is.element("ana.obj",names(x)) ) #typically use default
+                       x$ana.obj <- as.name("ana.obj")
                    eval(x)
                },
                scope=scope,method=method,

R/dom.s

@@ -4,5 +4,5 @@
     if (missing(scope)) stop("missing scope arg - possible bqtl syntax error")
     new.scope <- scope[grep("^dom",scope)]
     if (length(new.scope)==0) stop("no 'dom' terms found")
-    configs(x,...,scope=new.scope,method=method)
+    locus(x,...,scope=new.scope,method="fooey")
 }

R/lapadj.s

@@ -130,7 +130,7 @@ function(reg.formula, ana.obj,
               reg.frame <- as.data.frame(do.call("expand.grid", rep(list(mode.mat), nloc)))
           marker.distances <- switch(ana.obj$method,
                                      "RI.self"={ana.obj$map.frame$lambda/(2-ana.obj$map.frame$lambda)},
-                                     "RI.sib"=={ana.obj$map.frame$lambda/(4-3*ana.obj$map.frame$lambda)},
+                                     "RI.sib"={ana.obj$map.frame$lambda/(4-3*ana.obj$map.frame$lambda)},
                                      ana.obj$map.frame$lambda)
       }
       else

R/locus.s

@@ -1,5 +1,5 @@
 locus<-
-    function(x,...,scope=scope,method=NULL)
+    function(x,...,scope=scope,method=NULL,chromo=NULL,cM=NULL,ana.obj=NULL)
 {
 ###
 ###   
@@ -14,28 +14,44 @@ locus<-
 ### I parenthesize everything. S+3.4 needs this workaround 
 ###  ~(a):(u+v):(w) parses correctly, but not  ~a:(u+v):w
 ###
-    if (missing(scope)) stop("missing scope arg - possible bqtl syntax error")
-    dots <- list(...)
-    x.call <- match.call()$x
-    if  ( length(x.call) ==1 && deparse(x.call) == "all" ){  # all loci ? 
-        x <- seq(along=scope)
-        if (method=="F2") dim(x) <- c(2,length(x)%/%2)
-        if (length(dots) != 0) stop(". must be only arg")
-        return(configs(x,scope=scope))
-    } #else
-
-    if (length(x)>1 && length(dots)!=0)
-        stop("only one arg allowed with vector or matrix")
-    if (method == "F2") {
-        if (length(x)==1) {
-            x.1 <- 2*x-1
-            y <- 2*x
-            if ( length(dots)==0 )
-                dots <- list(y)
-            else
-                dots <- c(list(y),lapply(dots,function(x) c(2*x-1,2*x)))
-            configs.args <-  c(list(x=x.1),dots,scope=list(scope))
-        }
+  if (missing(scope)) stop("missing scope arg - possible bqtl syntax error")
+  dots <- list(...)
+  if (is.null(chromo)&& length(names(dots))!=0 &&
+      any(cloc <- which(1==pmatch(names(dots),"chromo",0)))){
+    chromo <- dots[[cloc]]
+    if (length(dots)>1)
+      stop("... not allowed")
+    else
+      dots<-NULL
+  }
+  if (!is.null(chromo)){
+    if (!missing(x)) stop("using both x and chromo args not allowed")
+    x <- if (is.null(cM))
+        map.index(ana.obj,chromo=chromo)
+    else
+        map.index(ana.obj,chromo=chromo,cM=cM)
+  }
+
+  x.call <- match.call()$x
+  if  ( length(x.call) ==1 && deparse(x.call) == "all" ){  # all loci ? 
+    x <- seq(along=scope)
+    if (method=="F2") dim(x) <- c(2,length(x)%/%2)
+    if (length(dots) != 0) stop("... not allowed")
+    return(configs(x,scope=scope))
+  } #else
+
+  if (length(x)>1 && length(dots)!=0)
+    stop("only one arg allowed with vector or matrix")
+  if (method == "F2") {
+    if (length(x)==1) {
+      x.1 <- 2*x-1
+      y <- 2*x
+      if ( length(dots)==0 )
+        dots <- list(y)
+      else
+        dots <- c(list(y),lapply(dots,function(x) c(2*x-1,2*x)))
+      configs.args <-  c(list(x=x.1),dots,scope=list(scope))
+    }
         else {
             x.1 <- 2*x-1
             y <- 2*x
@@ -48,12 +64,12 @@ locus<-
             }
             configs.args <-  c(list(x.1),scope=list(scope))
         }
-        do.call("configs",configs.args)
-    }
-    else {
-        
-        configs(x,...,scope=scope)
-    }
-    
+    do.call("configs",configs.args)
+  }
+  else {
+    configs.args <- c(list(x),dots,scope=list(scope))
+    do.call("configs",configs.args)
+  }
+
 }
 

R/make.state.matrix.s

@@ -6,6 +6,9 @@
 ### for BC1 method as.numeric(marker.states)%in%c( 1,2,5,6 )
 ### or NA's
 ###  if NA's are present, they are recoded to 6
+
+                                        # no this is wrong NA cause error
+
 ###  the assumed setup is as follows:
 ###    strains are A and a
 ###
@@ -33,6 +36,17 @@
     mf.num <- as.matrix(marker.frame)
     if (any(is.na(mf.num)))
         stop("marker.frame must be all numeric")
+
+    if (n.ind==0){            # only one marker
+      if (method=="F2")
+        return(
+               array(rbind(diag(3),c(1,2,0)/3,c(0,2,1)/3,c(1,2,1)/4)[mf.num,],
+                     c(dim(mf.num),3)))
+      else
+        return(array(rbind(diag(2),NA,NA,NA,c(1,1)/2)[mf.num,],
+                     c(dim(mf.num),2)))
+    }
+
     switch(method,
            F2 = {
                tmat <- matrix(c(1, 1, 1, -1, 0, 1, 1, -1, 1), 3, byrow

R/map.index.s

@@ -2,9 +2,9 @@
     function(x,...)
     UseMethod("map.index")
 "map.index.default"<-
-    function(x,chromo,cM)
+    function(x,chromo,cM=NULL)
 {
-    if (length(chromo)==1){
+    if (length(chromo)==1 && missing(cM) || length(cM)!=2){
         subset <-    x$chr.num==chromo
         if (!missing(cM)){
             this.cM <- x[ subset ,"cM" ]
@@ -19,6 +19,7 @@
     }
     else
     {# two locations find a loci in span
+      if (length(chromo)==1) chromo[2] <- chromo
         if (length(chromo)>2) stop("> 2 values for chromo not allowed")
         if (chromo[1]>chromo[2]) stop("chromo[1]>chromo[2]")
         subset <-

R/update.bqtl.s

@@ -15,7 +15,10 @@
     if (any( loc.mat > ncol(ana.obj$loc.right))||any(loc.mat<1) )
         stop("invalid locus number")
     subset <- attr(setup,"subset")
-
+    lambda <- switch(ana.obj$method,
+                   RI.self = {ana.obj$map.frame$lambda/(2 - ana.obj$map.frame$lambda)},
+                   RI.sib = {ana.obj$map.frame$lambda/(4 - 3 * ana.obj$map.frame$lambda)},
+                   ana.obj$map.frame$lambda)
     setup <-
         c(list("upbqtl"),
           setup[-c(1,22)],
@@ -25,7 +28,7 @@
           list( res = double(n.alt)),
           list( orig.x = as.double(setup$x)) ,
           list( loc.right = as.integer(ana.obj$loc.right[subset,]-1)) ,
-          list( map.lambda = as.double(ana.obj$map.frame$lambda)) ,
+          list( map.lambda = as.double(lambda)) ,
           list( state.matrix = as.double(ana.obj$state.matrix[subset,,])) ,
           list( n.state.loc = as.integer(dim(ana.obj$state.matrix)[2])))
 

man/bqtl.Rd

@@ -57,6 +57,8 @@ and Marginal Densities. \emph{JASA}, \bold{81},82--86.
 data(little.ana.bc )                        # load  BC1 dataset
 
 loglik( bqtl( bc.phenotype ~ 1, little.ana.bc ) ) #null loglikelihood
+                                                  #on chr 1 near cM 25
+loglik(bqtl(bc.phenotype~locus(chromo=1,cM=25),little.ana.bc))
 
 little.bqtl <-                              # two genes with epistasis
    bqtl(bc.phenotype ~ m.12 * m.24, little.ana.bc)

man/locus.Rd

@@ -11,7 +11,7 @@ function.  It is used on the right hand side of a formula in the
 
 }
 \usage{
-locus(x, ..., scope=<see below>,method=<see below>)
+locus(x, ...,chromo=NULL,cM=NULL scope=<see below>,method=<see below>)
 
 add(x,...)
 
@@ -40,6 +40,15 @@ bqtl( y ~ add(14) +dom(27), my.f2.object )
   Optional arguments (usually integers) to be used when
   \code{is.atomic(x)} is TRUE. 
 }
+\item{chromo }{A chromosome number or 2 ordered numbers. The loci on the
+  chromosome or in the range of chromosome numbers are used. If
+\code{chromo} is used, \code{x} must not be used.}
+\item{cM }{(Optional) map distance or two giving a location near a locus
+  or range of locations from which loci will be included. If the one
+  chromosome number is specified in \code{chromo}, \code{cM} must be 
+  ordered. If \code{cM} is omitted, all loci on the chromosome(s) will
+  be included.)
+
 \item{scope}{ (Optional and)
   Usually not supplied by the user.  Rather \code{bqtl} fills this in
   automatically.  A vector of regressor names, like the \code{reg.names}
@@ -48,17 +57,20 @@ bqtl( y ~ add(14) +dom(27), my.f2.object )
 \item{method}{ (Optional and) 
   Usually not supplied by the user. Like \code{scope}, \code{bqtl} takes
   care of filling this in with "BC1", "F2", et cetera as appropriate.}
-}
-\details{ \code{locus} is used in the model formula notation of
-\code{bqtl}, possibly more than once, and possibly with regressors named
-in the usual manner.  \code{locus} is intended to speed up the
-specification and examination of genetic models by allowing many models
-to be specified in a shorthand notation in a single model formula.  The
-names of genetic loci can consist of marker names, names that encode
-chromosome number and location, or other shorthand notations.  The names
-of terms in genetic models will typically include the names of the locus
-and may prepend "add." or "dom." or similar abbreviations for the
-'additive' and 'dominance' terms associated with the locus. 
+\item{ana.obj}{Usually not specified by the user. This is the
+  \code{analysis.object} to be used to lookup loci if a \code{chromo}
+  argument is used.}  }
+\details{ \code{locus} is used in the model
+formula notation of \code{bqtl}, possibly more than once, and possibly
+with regressors named in the usual manner.  \code{locus} is intended to
+speed up the specification and examination of genetic models by allowing
+many models to be specified in a shorthand notation in a single model
+formula.  The names of genetic loci can consist of marker names, names
+that encode chromosome number and location, or other shorthand
+notations.  The names of terms in genetic models will typically include
+the names of the locus and may prepend "add." or "dom." or similar
+abbreviations for the 'additive' and 'dominance' terms associated with
+the locus.
 
 When used as in \code{bqtl( y ~ locus(34), my.analysis.obj )}, it will
 look up the term or terms corresponding to the 34$^{th}$ locus. When
@@ -80,6 +92,11 @@ models.  And more generally, whenever \code{is.array(x)} is TRUE, the
 columns (or slices) specify \code{ dim(x)[1]/length(x)} different
 models. 
 
+The \code{chromo} argument performs a lookup of loci on the chromosome
+via the function \code{\link{map.index}}. If \code{cM} is also given,
+the locus nearest that location is used. If two values are given for
+\code{cM} all loci in the range are used.
+
 \code{add(x)} and \code{dom(x)} are alternatives that specify that only
 the \emph{additive} or \emph{dominance} terms in an F2 intercross. 
 }

src/upbqtl.c

@@ -163,7 +163,7 @@ void upbqtl(longint *crsType, longint *nparm,
       for (i=0; i<N; i++)
 	for (j=0;j< NLOC-1; j++)
 	  needProd[i+N*j] = 
-	    (longint) (loc_right[i+N*cur_loci[j]] < cur_loci[j+1]);
+	    (longint) (loc_right[i+N*cur_loci[j]] > cur_loci[j+1]);
       /* calc lambda */
       for (j=0;j< (NLOC-1); j++){
 	tmp = 1.0;
@@ -174,11 +174,11 @@ void upbqtl(longint *crsType, longint *nparm,
       /* get perm vector for x */
 
       for (i=0;i<NLOC;i++) {
-	newRadixProd[loc_order[i]]=radixProd[i];
+	newRadixProd[loc_order[i+NLOC*i_alt]]=radixProd[i];
       }
       for (i=0;i<NRX;i++){
 	k=0;
-	for (j=0;j<NLOC;j++) k+= ((i/radixProd[j]) % radix) * newRadixProd[j];
+	for (j=0;j<NLOC;j++) k+= ((i/newRadixProd[j]) % radix) * radixProd[j];
 	perm_indx[i]=k;
       }
       /* revise x */