library(fabia)


source("/seppdata/sepp/linkage/release/funct.R")

load(file="../individuals.Rda")


lii <- length(indi[,1])
samples <- 1:lii
indiA <- indi


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

N1 <- lengthAll%/%shift
endRunAll <- (N1-over+1) # 639
startRunAll <- 0

startRun <- 600

endRun <- 639


I <- intervallAll # All SNPs gives I/10 kbp DNA length
topLs <- 1000 # number of largest L used for histogram, default 100
IBDlength <- 50 # IBD length in kbp
thresCount <- 1e-5 # p-value of random histogram hit, default 1e-5
minSNPsFactor <- 3/4 # percentage of IBD overlap; 1/2 for large to 3/4 for for small intervals
pMAF <- 0.035  # minor allele frequency
samplesN <- 1092
# estimated for 1000 genomes: chromosome 1
# after removing 0 and private SNPs and MAF > 0.05:
# median distance 9.00   mean distance: 22.29
# 1/22.29 = 0.04486317  ==> k = 14
# 1/9 = 0.11            ==> k = 30

SS <- 10 # biclusters per iteration
A <- 40 # number iterations



ps <- 0.9 # percentage of Ls to remove / qunatile
psZ <- 0.8
inteA <- IBDlength*10 # histogram length gives inteA/10 kbp DNA length

kk <- 1
while ((I/inteA)*choose(samplesN,2)*(1-pbinom(kk,inteA,pMAF*pMAF))>thresCount) { kk <- kk +1}

thresA <- kk
minSNPs <- round(minSNPsFactor*thresA)



#######################################################################
#######################################################################

for (posAll in startRun:endRun)  {

start <- posAll*shift
end <- start + intervallAll

if (end > lengthAll)  {

  end <- lengthAll
}

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



annot <- read.table(paste("../ALL.chr1.merged_beagle_mach.20101123.snps_indels_svs.genotypes",pRange,"_annot.txt",sep=""),header = FALSE, sep = " ", quote = "",as.is = TRUE,skip=2)



for (i in 3:9) {

    annot[[i]] <- gsub(",",";",annot[[i]])

}

save(annot,file=paste("annotation_chr1",pRange,".Rda",sep=""))


source("/seppdata/sepp/linkage/release/split/run_fabia_split.R")

# plotIBDs(mergedIBD,paste("ALL.chr1.merged_beagle_mach.20101123.snps_indels_svs.genotypes",pRange,"_mat",sep=""))

ibd2excel(mergedIBD,paste("ALL.chr1.merged_beagle_mach.20101123.snps_indels_svs.genotypes",pRange,".csv",sep=""))



}