# Construct pedi data.frames for two types of pedigrees
pedi1 <- data.frame(Person=c("F","M","A","B"),
Father=c("RI","RI","F","F"),
Mother=c("RI","RI","M","M"))
pedi2 <- data.frame(Person=c("M","A","B"),
Father=c("RI","RI","RI"),
Mother=c("RI","M","M"))
LR_1=LR_2=data.frame(Log10CLR=rep(0,10000))
for (i in 1:42) {
# Simulate 10000 groups of A/B/P where A is full sibiling of B
Genotype1=pedisimu(af = FortytwoSTR$afmatrix[[i]],ss = 10000,pedi = pedi1)
LR_1=LR_1+LRhsip(A=Genotype1[,5:6],B=Genotype1[,7:8],P=Genotype1[,3:4],
af = FortytwoSTR$afmatrix[[i]],rare=FortytwoSTR$rare[i][1,1])
# Simulate 10000 groups of A/B/P where A is half sibling of B, i.e., the true phi=0
Genotype2=pedisimu(af = FortytwoSTR$afmatrix[[i]],ss = 10000,pedi = pedi2)
LR_2=LR_2+LRhsip(A=Genotype2[,3:4],B=Genotype2[,5:6],P=Genotype2[,1:2],
af = FortytwoSTR$afmatrix[[i]],rare=FortytwoSTR$rare[i][1,1])
}
# histograms of CLR distributions in the two groups
xmin<-floor(min(min(LR_1$Log10CLR),min(LR_2$Log10CLR)))
xmax<-ceiling(max(max(LR_1$Log10CLR),max(LR_2$Log10CLR)))
par(mfrow = c(1, 2))
hist(LR_2$Log10CLR,xlab = expression(log[10]~CLR),main = "False pedigree",
xlim = c(xmin,xmax), col = "red")
hist(LR_1$Log10CLR,xlab = expression(log[10]~CLR),main = "True cases",
xlim = c(xmin,xmax), col = "blue")
Run the code above in your browser using DataLab