qkappa {BookEKM} | R Documentation |
Based on conditional distribution given IBD from q012
,
the joint probability distribution for two individuals are given for
specified IBD probabilities.
Mutations have been removed from some previous examples.
qkappa(kappa = c(0, 1, 0), q = NULL)
kappa |
Three reals summing to 1 giving IBD (0,1,2) probabilities |
q |
The joint probability distribution for two individuals |
A matrix giving the joint distribution
Thore Egeland
To appear
FamiliasLocus
to specify mutation models.
require(Familias) require(paramlink) #SNP example p <- c(0.2, 0.8) L1<- FamiliasLocus(p, maleMutationRate = 0.001, allelenames = 1:length(p)) M <- L1$maleMutationMatrix qq <- q012(p) res1 <- qkappa(kappa=c(0,1,0), q=qq) # Check, Familias persons <- c("CH", "AF") sex <- c("male", "male") ped1 <- FamiliasPedigree(id = persons, dadid = c("AF", NA), momid = c(NA, NA), sex=sex) pedigrees <- list(isFather= ped1) CH <- c(2,2) AF <- c(1,1) datamatrix <- rbind(AF,CH) res2 <-FamiliasPosterior(ped1, L1, datamatrix)$likelihoodsPerSystem #Sibs. One SNP marker with qkappa(kappa=c(0.25,0.5,0.25),q012(p=c(0.2,0.8))) # Example. Table 8.1 in book (Egeland, Kling Mostad, 2016) reproduced and checked. # The code can easily be modified to deal with general pairwise relationships # (by changing kappa), markers having more than two alleles # (by changing afreq) and also general mutation models # (by changing the parameters for FamiliasLocus) p <- 0.7; q <- 1-p;m <- 0.05 afreq <- c(p,q) M <- FamiliasLocus(afreq, allelenames = 1:length(afreq), MutationRate = m)$maleMutationMatrix qq <- q012(afreq) #Genotype probs conditioned on IBD res1 <- qkappa(kappa=c(0,1,0), q=qq) #Unconditional genotype probs res2 <- qkappa(kappa=c(1,0.0,0), q=qq) order1 <- c(1,7,4,3,9,6,2,8,5) #Sort as in Table 8.1 LR.mut <- as.vector(res1/res2)[order1] M <- FamiliasLocus(afreq, allelenames = 1:length(afreq), MutationRate = 0.00)$maleMutationMatrix qq <- q012(afreq) #Genotype probs conditioned on IBD res1.no <- qkappa(kappa=c(0,1,0), q=qq) #Unconditional genotype probs res2.no <- qkappa(kappa=c(1,0.0,0), q=qq) LR <- as.vector(res1.no/res2.no)[order1] AF <- c(rep("AA",3),rep("AB",3),rep("BB",3)) CH <- rep(c("AA","AB","BB"),3) tab <- data.frame(AF,CH,"H1.no"=as.vector(res1.no)[order1], "H1"=as.vector(res1)[order1], "H2"=as.vector(res2.no)[order1],LR,LR.mut) # Check against formulae in Table 8.1 LRE <- c(1/p,1/(2*p),0,1/(2*p),1/(4*p*q),1/(2*q),0,1/(2*q),1/q) all(round(LR,12)==round(LRE,12)) LR.mut.E <- c((1-m)/p,((1-m)*q+m*p)/(2*p*q),m/q, 1/(2*p),1/(4*p*q),1/(2*q), m/p,((1-m)*p+m*q)/(2*p*q),(1-m)/q) all(round(LR.mut,12)==round(LR.mut.E,12)) # sum(tab[,3]*tab[,6]) #E(LR|Hp)=1.25 as it should sum(tab[,5]*tab[,6]) #E(LR|Hd)=1 as it should