Calculates the joint density of the crossover locations on a random meiotic
product, given that there are precisely two crossovers, for the gamma model.
Usage
joint.given.two(nu, L = 103, x, y, n = 20, max.conv = 25,
integr.tol = 0.00000001, max.subd = 1000, min.subd = 10)
Arguments
nu
The interference parameter in the gamma model.
L
The length of the chromsome in cM.
x
If specified, locations of the first crossover.
y
If specified, locations of the second crossover.
n
Number of points at which to calculate the density. The points
will be evenly distributed between 0 and L. Ignored if x and
y are specified.
max.conv
Maximum limit for summation in the convolutions to get
inter-crossover distance distribution from the inter-chiasma distance
distributions. This should be greater than the maximum number of chiasmata
on the 4-strand bundle.
integr.tol
Tolerance for convergence of numerical integration.
max.subd
Maximum number of subdivisions in numerical integration.
min.subd
Minimum number of subdivisions in numerical integration.
Details
Let f(x;nu) denote the density of a gamma random variable
with parameters shape=nu and rate=2 nu, and let
f_k(x;ν) denote the density of a gamma random variable
with parameters shape=k nu and rate=2 nu.
The distribution of the distance from one crossover to the next is
f*(x;nu) = sum_(k=1 to
infty) f_k(x;ν)/2^k.
The distribution of the distance from the start of the chromosome to the
first crossover is g*(x;nu) = 1 -
F*(x;nu) where F* is the cdf of f*.
Value
A data frame with three columns: x and y are the
locations (between 0 and L, in cM) at which the density was
calculated and f is the density.
Warning
We sometimes have difficulty with the numerical
integrals. You may need to use large min.subd (e.g. 25) to get
accurate results.
# Calculate the distribution of the average of the crossover locations,
# given that there are two and that they are separated by 20 cM
# (for a chromosome of length 200 cM)
L <- 200
d <- 20
x <- seq(0, L-d, by=0.5)
y <- x+d
f <- joint.given.two(4.3, L=L, x, y)
f$f <- f$f / distance.given.two(4.3, L, d)$f
plot((f$x+f$y)/2, f$f, type="l", xlim=c(0, L), ylim=c(0,max(f$f)),
lwd=2, xlab="Average location", ylab="Density")
abline(v=c(d/2,L-d/2), h=1/(L-d), lty=2, lwd=2)