kappaMat {BookEKM} | R Documentation |
A pair of markers separated by a recombination distance of rho
is considered. The 3 by 3 IBD matrix is calculated exactly.
The function twoLocusIBD
in IBDsim calculates
this for all pairwise relationship by simulation,
kappaMat(rho=0, relationship="UN", n1=NULL, n2=NULL, male=TRUE)
rho |
Real number in [0,0.5] |
relationship |
Can be "UN" (unrelated), "PO" (parent-offspring), "MZ" (monozygotic twins, "Grand" (grand parent - grand child), "Half" (half sibs), "Uncle" (uncle nephew), "Sibs", "EHS" (extended half sibs as described by n1 and n2), "EHSX" (extended half sib as described by n1, X markers), "ES" (extended sibs as described by n1 and n2), or "PC" (extended parent child as described by n1). For X, parent child is as before whereas more distant relationships are handled by "PC", with n1 equal to the number of females in the chain save the youngest. |
n1 |
Integer. If autosomal, number of meioses on left hand side. If X, number of female meioses |
n2 |
Integer. If autosomal, number of meioses on right hand side. If X, NULL |
male |
Logical. Only relevant for X and then TRUE implies a common male ancestor for extended half sib relationships. |
For one marker identity works generally. Other relationships can be solved in the pairwise case by simulation using twoLocusIBD or twoLocusJacquard (inbred case).
A 3 by 3 matrix of IBD probabilities
Thore Egeland <Thore.Egeland@gmail.com>
library(paramlink) library(IBDsim) #Example 1 rho <- 0.25 kUN <- kappaMat(relationship = "UN") kPO <- kappaMat(relationship = "PO") kMZ <- kappaMat(relationship = "MZ") kGrand <- kappaMat(rho, relationship = "Grand") kHalf <- kappaMat(rho, relationship = "Half") kUncle <- kappaMat(rho, relationship = "Uncle") kSibs <- kappaMat(rho, relationship = "Sibs") kEHS <- kappaMat(rho, relationship = "EHS", n1=1, n2=1) stopifnot(kHalf[1:2,1:2]==kEHS) kEHSX <- kappaMat(rho, relationship = "EHSX", n1=3, male=TRUE) kES <- kappaMat(rho, relationship = "ES", n1=1, n2=1) stopifnot(kSibs==kES) kPC <- kappaMat(rho, relationship="PC", n1=4) stopifnot(abs((kPC[2,2]- ((1-rho)/2)^3))<1e-15) #Example 2 y <- halfCousinPed(5) for (i in c(1,4,8,12,16,20,6,10,14,22,25)) y <- swapSex(y,i) plot(y) plot(y, id.labels= c(1:19, "20(Q)",21:24,"25(NLH)"), title="Fifth half-cousins once removed") p <- c(0.2,0.8); rho <- 0.1; chr <- 1 b1 <- kappaMat(rho,relationship="EHS",n1=5,n2=6) ## Not run: Nsim <- 10000 m1 <- marker(y, 20, c(1,1), alleles=1:2, afreq=p, chrom =chr) m2 <- marker(y, 20, c(1,1), alleles=1:2, afreq=p, chrom =chr) a1 <- twoMarkerDistribution(y, id=25, m1, m2, theta=rho, verbose=FALSE) stopifnot(abs(a1["2/2","2/2"]- b1[1,1]*p[2]^4)<1e-12) b2 <- twoLocusIBD(y, 20, 25, rho=rho, N=Nsim, Xchrom=FALSE) ## End(Not run) #Example 3 Extended half sibs, X chromosomal markers y <- halfCousinPed(2) y <- swapSex(y,c(4,5)) y <- swapSex(y,c(6,7)) y <- swapSex(y,c(8,9)) plot(y) rho <- 0.25; p <- c(0.2,0.8); chr <- 23 b1 <- kappaMat(rho, relationship="EHSX", n1=3, male=TRUE) ## Not run: m1 <- marker(y, 13, c(1,1), alleles=1:2, afreq=p, chrom =chr) m2 <- marker(y, 13, c(1,1), alleles=1:2, afreq=p, chrom =chr) a1 <- twoMarkerDistribution(y, id=12, m1, m2, theta=rho, verbose=FALSE) b1 <- kappaMat(rho, relationship="EHSX", n1=3, male=TRUE) stopifnot(a1["2","2"]== b1[1,1]*p[2]^2) b2 <- twoLocusIBD(y, 12, 13, rho=rho, N=10000, Xchrom=TRUE) ## End(Not run) # Example First cousins y <- cousinPed(1); chr <- 1; p <- c(0.2,0.8); rho <-0.1 m1 <- marker(y, 8, c(1,1), alleles=1:2, afreq=p, chrom =chr) m2 <- marker(y, 8, c(1,1), alleles=1:2, afreq=p, chrom =chr) a1 <- twoMarkerDistribution(y, id=7, m1, m2, theta=rho, verbose=FALSE) b1 <- kappaMat(rho, relationship="ES", n1=2,n2=2) stopifnot(abs(a1["2/2","2/2"]- b1[1,1]*p[2]^4 )< 1e-10) # Example Second cousins y <- cousinPed(2); chr <- 1; p <- c(0.4,0.6); rho <-0.25 m1 <- marker(y, 12, c(1,1), alleles=1:2, afreq=p, chrom =chr) m2 <- marker(y, 12, c(1,1), alleles=1:2, afreq=p, chrom =chr) a1 <- twoMarkerDistribution(y, id=11, m1, m2, theta=rho, verbose=FALSE) b1 <- kappaMat(rho, relationship="ES", n1=3,n2=3) stopifnot(abs(a1["2/2","2/2"] - b1[1,1]*p[2]^4) <1e-10)