#####load libraries####### library(hapFabia) library(fabia) #####define segments, overlap, filename ####### shiftSize <- 5000 segmentSize <- 10000 #fileName="KPGP1000_chr1" # without type fileName="../ALL.chr1.merged_beagle_mach.20101123.snps_indels_svs.genotypes" # without type #####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 results; parallel####### fileName="resIBD_chr1" anaRes <- analyzeIBDs(fileName=fileName,startRun=1,endRun=endRunA,shift=shiftSize,segmentSize=segmentSize) save(anaRes,file="anaRes.Rda") plot(density(anaRes$avibdPos)) hist(anaRes$avibdPos,breaks=100)