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","CHB","CHS","CLM","FIN","GBR","IBS","JPT","KOR","LWK","MXL","PUR","TSI","YRI")


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
tss <- sum(tt)
tt <- tt/sum(tt)

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



}

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

}

}


}




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

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

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

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