kappaMat {BookEKM}R Documentation

Calulates IBD matrices

Description

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,

Usage

kappaMat(rho=0, relationship="UN", n1=NULL,  n2=NULL, male=TRUE)

Arguments

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.

Details

For one marker identity works generally. Other relationships can be solved in the pairwise case by simulation using twoLocusIBD or twoLocusJacquard (inbred case).

Value

A 3 by 3 matrix of IBD probabilities

Author(s)

Thore Egeland <Thore.Egeland@gmail.com>

Examples

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)


[Package BookEKM version 1.0 Index]