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="")) }