Skip to content

Commit

Permalink
version 0.90
Browse files Browse the repository at this point in the history
  • Loading branch information
Saunak Sen authored and gaborcsardi committed Oct 2, 2006
1 parent 2e4132e commit 2381931
Show file tree
Hide file tree
Showing 10 changed files with 302 additions and 50 deletions.
10 changes: 9 additions & 1 deletion Changelog
Original file line number Diff line number Diff line change
Expand Up @@ -26,4 +26,12 @@
now output the percent variance explained also.
- New function called "prop.var" which calculates the proportion
of variance explained by a locus given the cross type, effect,
and error variance.
and error variance.
0.81 -> 0.86
- Added the functions to calculate genome-wide thresholds and
tail probabilities for genome scans using Poisson approximations.
- Changed search interval for detectable effect calculation.
- Renamed bio.var env.var in powercalc(), and error.var().
- Changed genome length argument to G from L in thresh() and tail.prob().
- Changed the units to centiMorgans from Morgans for genome length
and marker spacing arguments.
6 changes: 3 additions & 3 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
Package: qtlDesign
Title: Design of QTL experiments
Version: 0.85
Date: 19 February 2006
Version: 0.90
Date: 02 October 2006
Author: Saunak Sen, Jaya Satagopan, Karl Broman, and Gary Churchill
Description: Tools for the design of QTL experiments
Maintainer: Saunak Sen <sen@biostat.ucsf.edu>
License: GPL version 2 or newer (http://www.gnu.org/copyleft/gpl.html)
URL: http://www.biostat.ucsf.edu/sen/
Packaged: Thu Feb 16 11:41:03 2006; sen
Packaged: Mon Oct 2 10:14:05 2006; sen
46 changes: 23 additions & 23 deletions R/powercalc.R
Original file line number Diff line number Diff line change
@@ -1,17 +1,17 @@
"powercalc" <- function(cross,n,effect,sigma2,bio.var,gen.var,
"powercalc" <- function(cross,n,effect,sigma2,env.var,gen.var,
thresh=3,sel.frac=1,theta=0,bio.reps=1)
{
# if error variance missing, calculate it
if(missing(sigma2))
{
if((missing(bio.var))|(missing(gen.var)))
stop("Need either sigma2 or both bio.var and gen.var.")
sigma2 <- error.var(cross,bio.var,gen.var,bio.reps)
if((missing(env.var))|(missing(gen.var)))
stop("Need either sigma2 or both env.var and gen.var.")
sigma2 <- error.var(cross,env.var,gen.var,bio.reps)
}
else
{
if((!missing(bio.var))|(!missing(gen.var)))
stop("Need either sigma2 or both bio.var and gen.var.")
if((!missing(env.var))|(!missing(gen.var)))
stop("Need either sigma2 or both env.var and gen.var.")
}


Expand All @@ -30,20 +30,20 @@
}

"detectable" <- function(cross,n,effect=NULL,
sigma2,bio.var,gen.var,power=0.8,thresh=3,
sigma2,env.var,gen.var,power=0.8,thresh=3,
sel.frac=1,theta=0,bio.reps=1)
{
# argument check
if(missing(sigma2))
{
if((missing(bio.var))|(missing(gen.var)))
stop("Need either sigma2 or both bio.var and gen.var.")
sigma2 <- error.var(cross,bio.var,gen.var,bio.reps)
if((missing(env.var))|(missing(gen.var)))
stop("Need either sigma2 or both env.var and gen.var.")
sigma2 <- error.var(cross,env.var,gen.var,bio.reps)
}
else
{
if((!missing(bio.var))|(!missing(gen.var)))
stop("Need either sigma2 or both bio.var and gen.var.")
if((!missing(env.var))|(!missing(gen.var)))
stop("Need either sigma2 or both env.var and gen.var.")
}

if(cross=="bc")
Expand Down Expand Up @@ -77,20 +77,20 @@
t(ans)
}

"samplesize" <- function(cross,effect,sigma2,bio.var,gen.var,power=0.8,
"samplesize" <- function(cross,effect,sigma2,env.var,gen.var,power=0.8,
thresh=3,sel.frac=1,theta=0,bio.reps=1)
{
# argument check
if(missing(sigma2))
{
if((missing(bio.var))|(missing(gen.var)))
stop("Need either sigma2 or both bio.var and gen.var.")
sigma2 <- error.var(cross,bio.var,gen.var,bio.reps)
if((missing(env.var))|(missing(gen.var)))
stop("Need either sigma2 or both env.var and gen.var.")
sigma2 <- error.var(cross,env.var,gen.var,bio.reps)
}
else
{
if((!missing(bio.var))|(!missing(gen.var)))
stop("Need either sigma2 or both bio.var and gen.var.")
if((!missing(env.var))|(!missing(gen.var)))
stop("Need either sigma2 or both env.var and gen.var.")
}

if(cross=="bc")
Expand Down Expand Up @@ -139,7 +139,7 @@
# power, threshold, selection fraction, and size of marker interval
effect <- uniroot(function(x) {
powercalc.bc(n, x, sigma2, thresh, sel.frac, theta) -
power}, interval = c(0,10*sqrt(sigma2/n)))$root
power}, interval = c(0,30*sqrt(sigma2/n)))$root
effect
}

Expand Down Expand Up @@ -229,7 +229,7 @@
# power, threshold, selection fraction, and size of marker interval
del <- uniroot(function(x) {
powercalc.f2(n, x*effect, sigma2, thresh, sel.frac, theta) - power },
interval = c(0,10*sqrt(sigma2/n)))$root
interval = c(0,30*sqrt(sigma2/n)))$root

# decide what to return depending on delta flag
return(del*effect)
Expand Down Expand Up @@ -278,7 +278,7 @@
# power, threshold, selection fraction, and size of marker interval
effect <- uniroot(function(x) {
powercalc.ri(n, x, sigma2, thresh) - power},
interval = c(0,10*sqrt(sigma2/n)))$root
interval = c(0,30*sqrt(sigma2/n)))$root
effect
}

Expand Down Expand Up @@ -307,7 +307,7 @@



"error.var" <- function(cross,bio.var=1,gen.var=0,bio.reps=1)
"error.var" <- function(cross,env.var=1,gen.var=0,bio.reps=1)
{
# get the genetic variance multiplier
if(cross=="bc")
Expand All @@ -319,7 +319,7 @@
else
error("Cross type not recognized.")

bio.var/bio.reps + CC*gen.var
env.var/bio.reps + CC*gen.var

}

Expand Down
164 changes: 164 additions & 0 deletions R/thresh.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,164 @@
##################################################################
# function for calculating tail probabilities of genome scans
##################################################################
# t = threshold on LOD (base 10) scale
# G = genome length in centiMorgans
# cross = cross type, "bc" or "f2"
# type = type of LOD score
# "1": one-dim scan
# "2o": two-dim scan, full model
# "2i": two-dim scan, interaction only
# "2a": two-dim scan, additive only
# "3a": three-dim scan, additive only
# d = marker spacing in centiMorgans
# cov.dim = dimension of interacting covariate

tailprob <- function(t,G,cross,type="1",d=0.01,cov.dim=0)
{

# convert to Morgans
G <- G/100
d <- d/100

# convert LOD to chi-squared scale
t <- t*log(10)*2

# local drift rate
# mu <- switch(cross,bc=sqrt(2*t*d),f2=sqrt(3*t*d))
mu <- switch(cross,bc=sqrt(4*t*d),f2=sqrt(6*t*d))

# correction for terms not contributing to drift

# backcross
if((cross=="bc")&&(type=="2o"))
{
mu <- mu*sqrt(2/3)
}

if((cross=="bc")&&(type=="2x"))
{
mu <- mu
}

# intercross
if((cross=="f2")&&(type=="2o"))
{
mu <- mu*sqrt(6/8)
}

if((cross=="f2")&&(type=="2x"))
{
mu <- mu
}

# c <- 2*mu^2*nu(2*mu)
# c <- mu^2*nu(mu)/2
# c <- tau(mu)


# calculate space volume and clump volume from search dimension

# one-dimensional search
if(regexpr("1",type)==1)
{
c <- mu^2*nu(mu)/2
S <- (G/d)
}

# two-dimensional search
if(regexpr("2",type)==1)
{
if(type=="2x")
{
if(cross=="bc")
c <- (mu^2*nu(mu)/2)*((mu/2)^2*nu(mu/2)/2)
if(cross=="f2")
c <- (mu^2*nu(mu)/2)*((mu*2/3)^2*nu(mu*2/3)/2)
}
else
{
c <- mu^2*nu(mu)/2
c <- c^2
}
S <- (G/d)^2/2
}

if(regexpr("3",type)==1)
{
c <- c^3
S <- (G/d)^3/6
}




# degrees of freedom of marginal test statistic
df <- (1+cov.dim) * switch(paste(cross,type,sep="."),
bc.1=1, bc.2i=1, bc.2a=2, bc.2o=3, bc.2x=2,
f2.1=2, f2.2i=4, f2.2a=4, f2.2o=8, f2.2x=6,
bc.3a=3,f2.3a=6)

# put everything together: search space volume, (inverse of)
# expected clump size and marginal tail probability

lam <- S*c*pchisq(t,df,low=FALSE)
# lam <- 2*S*c*dchisq(t,df)

# exponentiate for tail probability
1-exp(-lam)
}


##################################################################
# function for calculating thresholds
##################################################################
# G = genome length in centiMorgans
# cross = cross type, "bc" or "f2"
# type = type of LOD score
# "1": one-dim scan
# "2o": two-dim scan, full model
# "2i": two-dim scan, interaction only
# "2a": two-dim scan, additive only
# "3a": three-dim scan, additive only
# p = vector of genome-wide significance levels
# d = marker spacing in centiMorgans
# cov.dim = dimension of interacting covariate
# interval = LOD interval to be searched for threshold


# calculate thresholds from tail probabilities
thresh <- function(G,cross,type="1",p=c(0.10,0.05,0.01),
d=0.01,cov.dim=0,interval=c(1,40))
{

thresh <- vector(mode="numeric",length=length(p))
for( i in 1:length(p) )
{
thresh[i] <- uniroot(function(x){
tailprob(x,G,cross,type,d,cov.dim) - p[i]},interval=interval)$root
}
thresh

}

# the nu function of Siegmund

nu <- function(mu,k.lim=10000)
{
if(abs(mu)<0.1)
{
res <- -0.583*mu
}
else
{
k <- 1:k.lim
zzz <- pnorm( -0.5 * mu * sqrt( k ) ) / k
res <- -2*sum(zzz) + log(2) - 2*log(mu)
}
exp(res)
}

tau <- function(mu)
{
2*mu^2*nu(2*mu)
}
6 changes: 6 additions & 0 deletions R/version.qtlDesign.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
"version.qtlDesign" <-
function()
{
0.90
}

29 changes: 16 additions & 13 deletions man/error.var.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -7,22 +7,22 @@
calculations}
\description{
The function \code{error.var} estimates the error variance using
estimates of the biological variance and genetic variance. The effect
estimates of the environmental variance and genetic variance. The effect
segregating at a locus, can be calculated using \code{gmeans2effect}
These are key inputs for power calculations. The function
\code{prop.var} calculates the proportion of variance explained by a
locus given the effect size and error variance.
}
\usage{
error.var(cross,bio.var=1,gen.var=0,bio.reps=1)
error.var(cross,env.var=1,gen.var=0,bio.reps=1)
gmeans2effect(cross,means)
prop.var(cross,effect,sigma2)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{cross}{String indicating cross type which is "bc", for
backcross, "f2" for intercross, and "ri" for recombinant inbred lines.}
\item{bio.var}{Biological (within genotype) variance}
\item{env.var}{Environmental (within genotype) variance}
\item{gen.var}{Genetic (between genotype) variance due to all loci
segregating between the parental lines.}
\item{bio.reps}{Number of biological replicates per unique genotype.
Expand All @@ -44,17 +44,20 @@ prop.var(cross,effect,sigma2)
\item{sigma2}{Error variance.}
}
\details{The function \code{error.var} estimates the error variance
segregating in a cross using estimates of the biological (within
segregating in a cross using estimates of the environmental (within
genotype) variance, and the genetic (between genotype variance). The
biological variance is assumed to be invariant between cross types.
environmental variance is assumed to be invariant between cross types.
The genetic variance segregating in RI lines is assumed to be double
that in F2 intercross, and four times that of the backcross. This
assumption holds if all loci are additive. The error variance at a
locus of interest is aproximately \deqn{\gamma^2/c + \tau^2/m}, where
\eqn{\gamma^2} is the genetic variance (\code{gen.var}), \eqn{c} is a
locus of interest is aproximately \deqn{\sigma_G^2/c + \sigma_E^2/m,}{%
sigmaG^2/c + sigmaE^2/m,}
where \eqn{$\sigma_G^2$}{sigmaG^2} is the genetic variance
(\code{gen.var}), \eqn{c} is a
constant depending on the cross type (1, for RI lines, 1/2 for F2
intercross, and 1/4 for backross), \eqn{tau^2} is the biological
variance (\code{bio.var}), and \eqn{m} is the number of biological
intercross, and 1/4 for backross), \eqn{$\sigma_E^2$}{sigmaE^2} is the
environmental
variance (\code{env.var}), and \eqn{m} is the number of biological
replicates per unique genotype (\code{bio.reps}).

The function \code{gmeans2effect} calculates the QTL effects from
Expand Down Expand Up @@ -84,10 +87,10 @@ Genetics, 170:447-64.}
\author{Saunak Sen, Jaya Satagopan, Karl Broman, and Gary Churchill}
\seealso{\code{\link{powercalc}}}
\examples{
error.var(cross="bc",bio.var=1,gen.var=1,bio.reps=1)
error.var(cross="f2",bio.var=1,gen.var=1,bio.reps=1)
error.var(cross="ri",bio.var=1,gen.var=1,bio.reps=1)
error.var(cross="ri",bio.var=1,gen.var=1,bio.reps=10)
error.var(cross="bc",env.var=1,gen.var=1,bio.reps=1)
error.var(cross="f2",env.var=1,gen.var=1,bio.reps=1)
error.var(cross="ri",env.var=1,gen.var=1,bio.reps=1)
error.var(cross="ri",env.var=1,gen.var=1,bio.reps=10)
gmeans2effect(cross="f2",means=c(0,1,2))
gmeans2effect(cross="f2",means=c(0,1,1))
gmeans2effect(cross="bc",means=c(0,1,1))
Expand Down

0 comments on commit 2381931

Please sign in to comment.