# Simulate 10,000 mother-child pairs with father-daughter incest with pedisimu() function
# based on the 42 STRs in "FortytwoSTR" data.
pedi<-data.frame(Person=c("F","M","C"),Father=c("RI","F","F"),Mother=c("RI","RI","M"))
II_1=II_2=data.frame(Ngs=rep(0,10000),IIphi=rep(0,10000))
for (i in 1:42) {
Genotype1<-pedisimu(af = FortytwoSTR$afmatrix[[i]],ss = 10000,pedi = pedi)
# Calculate II for each case.
II_1<-II_1+IICAL(Parent = Genotype1[,3:4],Child = Genotype1[,5:6],af=FortytwoSTR$afmatrix[[i]],
rare=FortytwoSTR$rare[i][1,1],phi=0.25)
#Simulate 10,000 non-inbred mother-child pairs
Genotype2<-pairsimu(af = FortytwoSTR$afmatrix[[i]],ss = 10000,delta = c(0,1,0),allelename = FALSE)
II_2<-II_2+IICAL(Parent = Genotype2[,1:2],Child = Genotype2[,3:4],af=FortytwoSTR$afmatrix[[i]],
rare=FortytwoSTR$rare[i][1,1],phi=0.25)
}
# histograms of CII distributions in the two groups
xmin<-floor(min(min(II_1$IIphi),min(II_2$IIphi)))
xmax<-ceiling(max(max(II_1$IIphi),max(II_2$IIphi)))
par(mfrow = c(1,2))
hist(II_2$IIphi,xlab = expression(log[10]~CII),main = "Non-inbreed cases",
xlim = c(xmin,xmax), col = "red")
hist(II_1$IIphi,xlab = expression(log[10]~CII),main = "Inbreed cases",
xlim = c(xmin,xmax), col = "blue")
Run the code above in your browser using DataLab