resBeaglet <- read.table(file=paste("beagleres",output_index,".txt.dataSim",run,"beagle.txt.fibd",sep="")) ############################## #####RES1 with 1.0e-50 thres1 <- 1.0e-50 nnrr <- 0 while ((nnrr<1)&&(thres1<1.0e-4)) { resBeagle <- resBeaglet[which(resBeaglet[,5]0) { for (i in 1:nrow(resBeagle)) { samp1 <- resBeagle[i,1] samp2 <- resBeagle[i,2] begin <- resBeagle[i,3] end <- resBeagle[i,4] if (begin==0) { begin <- 1 } lengthIBD <- pos[end] - pos[begin] if ((lengthIBD>shunkFact*ShunkLength)) { allSamp <- union(allSamp,samp1) allSamp <- union(allSamp,samp2) range <- begin:end allSnps <- union(allSnps,range) foundAll[samp1,range] <- 1 foundAll[samp2,range] <- 1 } } } all <- samplesN*snps detected <- sum(foundAll) Pos <- sum(implantT) truePM <- foundAll*implantT trueP <- sum(truePM) falseP <- detected-trueP falseN <- Pos-trueP trueN <- all - trueP - falseN - falseP Neg <- trueN+falseP sens <- trueP/Pos spec <- trueN/Neg if (detected>0) { fdr <- falseP/detected } else { fdr <- 0 } fpr <- falseP/Neg acc <- (trueP+trueN)/all st11 <- (2*trueP+falseP+falseN) if (st11>0) { F1 <- 2*trueP/(2*trueP+falseP+falseN) } else { F1 <- 0 } stt <- (trueP+falseP)*(trueP+falseN)*(trueN+falseP)*(trueN+falseN) if (stt>0) { mcc <- (trueP*trueN - falseP*falseN)/sqrt(stt) } else { mcc <- 0 } allSampI <- intersect(implGA,allSamp) lallSamp <- length(allSamp) allSnpsI <- intersect(interA,allSnps) lallSnps <- length(allSnps) ## samples detectedsa <- lallSamp truePsa <- length(allSampI) falsePsa <- detectedsa-truePsa falseNsa <- implantedS-truePsa trueNsa <- allsa - truePsa - falseNsa - falsePsa senssa <- truePsa/Possa specsa <- trueNsa/Negsa if (detectedsa>0) { fdrsa <- falsePsa/detectedsa } else { fdrsa <- 0 } fprsa <- falsePsa/Negsa accsa <- (truePsa+trueNsa)/allsa st11sa <- (2*truePsa+falsePsa+falseNsa) if (st11sa>0) { F1sa <- 2*truePsa/(2*truePsa+falsePsa+falseNsa) } else { F1sa <- 0 } sttsa <- (truePsa+falsePsa)*(truePsa+falseNsa)*(trueNsa+falsePsa)*(trueNsa+falseNsa) if (sttsa>0) { mccsa <- (truePsa*trueNsa - falsePsa*falseNsa)/sqrt(sttsa) } else { mccsa <- 0 } ## snps detectedsn <- lallSnps truePsn <- length(allSnpsI) falsePsn <- detectedsn - truePsn falseNsn <- length(interA) - truePsn trueNsn <- snps - truePsn - falseNsn - falsePsn senssn <- truePsn/Possn specsn <- trueNsn/Negsn if (detectedsn>0) { fdrsn <- falsePsn/detectedsn } else { fdrsn <- 0 } fprsn <- falsePsn/Negsn accsn <- (truePsn+trueNsn)/allsn st11sn <- (2*truePsn+falsePsn+falseNsn) if (st11sn>0) { F1sn <- 2*truePsn/(2*truePsn+falsePsn+falseNsn) } else { F1sn <- 0 } sttsn <- (truePsn+falsePsn)*(truePsn+falseNsn)*(trueNsn+falsePsn)*(trueNsn+falseNsn) if (sttsn>0) { mccsn <- (truePsn*trueNsn - falsePsn*falseNsn)/sqrt(sttsn) } else { mccsn <- 0 } dataR <- c(run, truePsa,falseNsa,falsePsa,trueNsa,senssa,specsa,fdrsa,fprsa,accsa,F1sa,mccsa, truePsn,falseNsn,falsePsn,trueNsn,senssn,specsn,fdrsn,fprsn,accsn,F1sn,mccsn, trueP,falseN,falseP,trueN,sens,spec,fdr,fpr,acc,F1,mcc) beagleres1[[run]] <- dataR dim(dataR) <- c(1,length(dataR)) write(dataR,file=paste("beagleres1",output_index,".txt",sep=""),append=TRUE,ncolumns=100) ################################# #####RES1 with 1.0e-80 thres1 <- 1.0e-80 nnrr <- 0 while ((nnrr<1)&&(thres1<1.0e-4)) { resBeagle <- resBeaglet[which(resBeaglet[,5]0) { for (i in 1:nrow(resBeagle)) { samp1 <- resBeagle[i,1] samp2 <- resBeagle[i,2] begin <- resBeagle[i,3] end <- resBeagle[i,4] if (begin==0) { begin <- 1 } lengthIBD <- pos[end] - pos[begin] if (lengthIBD>shunkFact*ShunkLength) { allSamp <- union(allSamp,samp1) allSamp <- union(allSamp,samp2) range <- begin:end allSnps <- union(allSnps,range) foundAll[samp1,range] <- 1 foundAll[samp2,range] <- 1 } } } all <- samplesN*snps detected <- sum(foundAll) Pos <- sum(implantT) truePM <- foundAll*implantT trueP <- sum(truePM) falseP <- detected-trueP falseN <- Pos-trueP trueN <- all - trueP - falseN - falseP Neg <- trueN+falseP sens <- trueP/Pos spec <- trueN/Neg if (detected>0) { fdr <- falseP/detected } else { fdr <- 0 } fpr <- falseP/Neg acc <- (trueP+trueN)/all st11 <- (2*trueP+falseP+falseN) if (st11>0) { F1 <- 2*trueP/(2*trueP+falseP+falseN) } else { F1 <- 0 } stt <- (trueP+falseP)*(trueP+falseN)*(trueN+falseP)*(trueN+falseN) if (stt>0) { mcc <- (trueP*trueN - falseP*falseN)/sqrt(stt) } else { mcc <- 0 } allSampI <- intersect(implGA,allSamp) lallSamp <- length(allSamp) allSnpsI <- intersect(interA,allSnps) lallSnps <- length(allSnps) ## samples detectedsa <- lallSamp truePsa <- length(allSampI) falsePsa <- detectedsa -truePsa falseNsa <- implantedS-truePsa trueNsa <- allsa - truePsa - falseNsa - falsePsa senssa <- truePsa/Possa specsa <- trueNsa/Negsa if (detectedsa>0) { fdrsa <- falsePsa/detectedsa } else { fdrsa <- 0 } fprsa <- falsePsa/Negsa accsa <- (truePsa+trueNsa)/allsa st11sa <- (2*truePsa+falsePsa+falseNsa) if (st11sa>0) { F1sa <- 2*truePsa/(2*truePsa+falsePsa+falseNsa) } else { F1sa <- 0 } sttsa <- (truePsa+falsePsa)*(truePsa+falseNsa)*(trueNsa+falsePsa)*(trueNsa+falseNsa) if (sttsa>0) { mccsa <- (truePsa*trueNsa - falsePsa*falseNsa)/sqrt(sttsa) } else { mccsa <- 0 } ## snps detectedsn <- lallSnps truePsn <- length(allSnpsI) falsePsn <- detectedsn -truePsn falseNsn <- length(interA) - truePsn trueNsn <- snps - truePsn - falseNsn - falsePsn senssn <- truePsn/Possn specsn <- trueNsn/Negsn if (detectedsn>0) { fdrsn <- falsePsn/detectedsn } else { fdrsn <- 0 } fprsn <- falsePsn/Negsn accsn <- (truePsn+trueNsn)/allsn st11sn <- (2*truePsn+falsePsn+falseNsn) if (st11sn>0) { F1sn <- 2*truePsn/(2*truePsn+falsePsn+falseNsn) } else { F1sn <- 0 } sttsn <- (truePsn+falsePsn)*(truePsn+falseNsn)*(trueNsn+falsePsn)*(trueNsn+falseNsn) if (sttsn>0) { mccsn <- (truePsn*trueNsn - falsePsn*falseNsn)/sqrt(sttsn) } else { mccsn <- 0 } dataR <- c(run, truePsa,falseNsa,falsePsa,trueNsa,senssa,specsa,fdrsa,fprsa,accsa,F1sa,mccsa, truePsn,falseNsn,falsePsn,trueNsn,senssn,specsn,fdrsn,fprsn,accsn,F1sn,mccsn, trueP,falseN,falseP,trueN,sens,spec,fdr,fpr,acc,F1,mcc) beagleres2[[run]] <- dataR dim(dataR) <- c(1,length(dataR)) write(dataR,file=paste("beagleres2",output_index,".txt",sep=""),append=TRUE,ncolumns=100)