
#####load libraries#######
library(hapFabia)
library(fabia)

#####convert from .vcf to _mat.txt: step 2. above#######
vcftoFABIA(fileName=fileName)

#####split/ generate segments: step 3. above#######
split_sparse_matrix(fileName=fileName,segmentSize=segmentSize,shiftSize=shiftSize,annotation=TRUE)

#####compute how many segments we have#######
ina <- as.numeric(readLines(paste(fileName,"_mat.txt",sep=""),n=2))
snps <- ina[2]
over <- segmentSize%/%shiftSize
N1 <- snps%/%shiftSize
endRunA <- (N1-over+2)

#####analyze each segment#######
#####may be done by parallel runs#######
iterateSegments(startRun=1,endRun=endRunA,shift=shiftSize,segmentSize=segmentSize,fileName=fileName,samples=0,upperBP=0.05,p=10,iter=40,alpha=0.03,cyc=50,IBDlength=50,Lt = 0.1,Zt = 0.2,thresCount=1e-5,minSNPsFactor=3/4,pMAF=0.035,haplotypes=TRUE)

#####identify duplicates#######
identifyDuplicates(fileName=fileName,startRun=1,endRun=endRunA,shift=shiftSize,segmentSize=segmentSize)

#####analyze results; parallel#######
anaRes <- analyzeIBDs(fileName=fileName,startRun=1,endRun=endRunA,shift=shiftSize,segmentSize=segmentSize)
print("Number haplotypes:")
print(anaRes$noIBD)
print("Haplotye length:")
print(anaRes$avibdLengthS)
print("Number of individuals possessing a haplotye:")
print(anaRes$avnoSampS)
print("Number of tagSNVs of a haplotye:")
print(anaRes$avnoSnpsS)
print("MAF of tagSNVs of a haplotye:")
print(anaRes$avnoFreqS)
print("MAF within the group of tagSNVs of a haplotye:")
print(anaRes$avnoGroupFreqS)
print("Number of changes between major and minor allele frequency:")
print(anaRes$avnosnpChangeS)
print("Number of tagSNVs per individual:")
print(anaRes$avnosnpsPerSampleS)
print("Number of individuals having a tagSNV:")
print(anaRes$avnosamplePerSnpS)

#####load result for segment 50#######
posAll <- 50 # (50-1)*5000 = 245000: segment 245000 to 255000
start <- (posAll-1)*shiftSize
end <- start + segmentSize
pRange <- paste("_",format(start,scientific=FALSE),"_",format(end,scientific=FALSE),sep="")
load(file=paste(fileName,pRange,"_resAnno",".Rda",sep=""))
IBD <- resHapFabia$mergedIBD # $

#####plot specific haplotype (the first) in segment 50#######
plotOneIBD(IBD,filename=paste(fileName,pRange,"_mat",sep=""),ibdC=1)
   ##attention: filename without type ".txt"
