load("dups.Rda")
load("analyzeResult.Rda")
load(file="individuals.Rda")
library(hapFabia)
library(fabia)
#source("prepareAncientGenome.R")
#source("plotAncientGenomes.R")

shiftSize <- 5000
segmentSize <- 10000
fileName="KPGP1000_chr1" # without type


startRun <- 1
endRun <- 633

runIndex=""


lengthAll <- 3167651
shift <- 5000
intervallAll <- 10000
over <- intervallAll%/%shift

N1 <- lengthAll%/%shift
endRunAll <- (N1-over+2) # 639


pops <- c("ASW","CEU","CEUKOR","CHB","CHS","CLM","FIN","GBR","IBS","JPT","KOR","KORdup","LWK","MXL","PUR","TSI","YRI")

# "KPGP10","KPGP11","KPGP12","KPGP3","KPGP4","KPGP5","KPGP6","KPGP7","KPGP8","KPGP88", "KPGP90","KPGP9"
#Twins: 88,89 and 90,91
# 10 is USA Caucasian
# 11 and 12 are Caucasian and Asian mixture: children of 10
#Family: 1-12
# 1 and 2: Grandparents
# 6 is external
# 3,5,9 children of 1 and 2
# 4 is children of 3
# 7 and 8 are children of 5 and 6

allCor <-  c(1,3,12,28,29,30,31,32,33,35,36,39)
korCEU <- c(1)
korCEUKOR <- c(3,12)
korDups <- c(28,29,30,31,32,33,35,36,39)

#> indi[korDups,1]
# [1] "KPGP11" "KPGP12" "KPGP3"  "KPGP4"  "KPGP5"  "KPGP6"  "KPGP7"  "KPGP88"
# [9] "KPGP8"  "KPGP90" "KPGP9" 

indi[korDups,2] <- "KORdup"
indi[korCEU,2] <- "CEU"
indi[korCEU,3] <- "EUR"
indi[korCEUKOR,2] <- "CEUKOR"
indi[korCEUKOR,3] <- "EURASN"



allCount <- 0
allCount1 <- 0

resS <- list()


for (posAll in startRun:endRun)  {

start <- (posAll-1)*shift
end <- start + intervallAll

if (end > lengthAll)  {

  end <- lengthAll
}
 
  
pRange <- paste("_",format(start,scientific=FALSE),"_",format(end,scientific=FALSE),sep="")


load(file=paste(fileName,pRange,"_resAnno",".Rda",sep=""))
IBDannot <- resHapFabia$mergedIBD 

noIBD <- length(IBDannot)

if (noIBD>0)  {

for (ibdC in 1:noIBD)  {

allCount <- allCount + 1

if (!dups[allCount]) {

allCount1 <- allCount1 + 1

samples <-  IBDannot[[ibdC]][[8]]

tt <- table(c(pops,indi[samples,2]))
tt <- tt - 1
st <- sum(tt)
tss <- st
tt <- tt/st

if (max(tt)>0.9) {
tt <- table(pops) -1
tss <-  0
} else {
samples <- setdiff(samples,allCor)
tt <- table(c(pops,indi[samples,2]))
tt <- tt - 1
st <- sum(tt)
tss <- st
if ((st>0)&&(length(samples)>1)) {
  tt <- tt/st
} else {
tt <- table(pops) -1
tss <-  0
}

}

res1 <- c(allCount1,allCount,ibdC,posAll,tss,tt)
resS[[allCount1]] <- res1



}

cat(allCount," out of 222246\r")

}

}


}




rr <- unlist(resS)
l <- length(res1)

SubPopCorMat <- matrix(rr,nrow=allCount1,ncol=l,byrow=TRUE)

colnames(SubPopCorMat) <- c("allCount1","allCount","ibdC","posAll","#Samples",pops)
                             

write.table(SubPopCorMat,file=paste("SubPopCorMat",runIndex,".txt",sep=""))
            
save(resS,SubPopCorMat,file=paste("SubPopCorMat",runIndex,".Rda",sep=""))