resultsList <- list()
i <- 0
methodN <- c()

dofabia <- TRUE
#dobeagle <- FALSE
dobeagle <- TRUE
#doplink <- FALSE
doplink <- TRUE
#dogermline <- FALSE
dogermline <- TRUE
#dodash <-  FALSE
dodash <- TRUE
dorelate <- FALSE
domcmc <- FALSE

resHeader <- c("truePSample","falseNSample","falsePSample","trueNSample","sensitivitySample","specifitySample","FDRSample","FPRSample","accuracySample","F1Sample","MatthewsSample",
                   "truePSnps","falseNSnps","falsePSnps","trueNSnps","sensitivitySnps","specifitySnps","FDRSnps","FPRSnps","accuracySnps","F1Snps","MatthewsSnps",
                   "trueP","falseN","falseP","trueN","sensitivity","specifity","FDR","FPR","accuracy","F1","Matthews")

fabiaHeader <- c("runNumber","restarts","numberIBDs",
                   "truePSample","falseNSample","falsePSample","trueNSample","sensitivitySample","specifitySample","FDRSample","FPRSample","accuracySample","F1Sample","MatthewsSample",
                   "truePSnps","falseNSnps","falsePSnps","trueNSnps","sensitivitySnps","specifitySnps","FDRSnps","FPRSnps","accuracySnps","F1Snps","MatthewsSnps",
                   "trueP","falseN","falseP","trueN","sensitivity","specifity","FDR","FPR","accuracy","F1","Matthews")

beagleHeader <- c("runNumber",
"truePSample","falseNSample","falsePSample","trueNSample","sensitivitySample","specifitySample","FDRSample","FPRSample","accuracySample","F1Sample","MatthewsSample",
"truePSnps","falseNSnps","falsePSnps","trueNSnps","sensitivitySnps","specifitySnps","FDRSnps","FPRSnps","accuracySnps","F1Snps","MatthewsSnps",
"trueP","falseN","falseP","trueN","sensitivity","specifity","FDR","FPR","accuracy","F1","Matthews")


dashHeader <- c("runNumber","numberClusters",
"truePSample","falseNSample","falsePSample","trueNSample","sensitivitySample","specifitySample","FDRSample","FPRSample","accuracySample","F1Sample","MatthewsSample",
"truePSnps","falseNSnps","falsePSnps","trueNSnps","sensitivitySnps","specifitySnps","FDRSnps","FPRSnps","accuracySnps","F1Snps","MatthewsSnps",
"trueP","falseN","falseP","trueN","sensitivity","specifity","FDR","FPR","accuracy","F1","Matthews")               


if (dofabia) {
  fabiaresM <- read.table(file="fabiares.txt",skip=1)

  colnames(fabiaresM) <- fabiaHeader
  fabiaresM <- fabiaresM[,(-1:-3)]

  i <- i+1
  resultsList[[i]] <- fabiaresM
  methodN <- c(methodN,"fabia----")
  
}

if (dobeagle) {

  beagleres1M <- read.table(file="beagleres1.txt",skip=1)

  
  colnames(beagleres1M) <- beagleHeader

  beagleres1M <- beagleres1M[,-1]

  i <- i+1
  resultsList[[i]] <- beagleres1M

  methodN <- c(methodN,"beagle1--")

  
  beagleres2M <- read.table(file="beagleres2.txt",skip=1)

  colnames(beagleres2M) <- beagleHeader

  beagleres2M <- beagleres2M[,-1]

  i <- i+1
  resultsList[[i]] <- beagleres2M

  methodN <- c(methodN,"beagle2--")
}

if (doplink) {

  plinkresM <- read.table(file="plinkres.txt",skip=1)

  plinkHeader <- beagleHeader

  colnames(plinkresM) <- plinkHeader

  plinkresM <- plinkresM[,-1]

  i <- i+1
  resultsList[[i]] <- plinkresM

  methodN <- c(methodN,"plink----")

}

if (dogermline) {

  germlineres1M <- read.table(file="germlineres1.txt",skip=1)

  germlineHeader <- beagleHeader

  colnames(germlineres1M) <- germlineHeader

  germlineres1M <- germlineres1M[,-1]
  
  i <- i+1
  resultsList[[i]] <- germlineres1M

  methodN <- c(methodN,"germline1")

  germlineres2M <- read.table(file="germlineres2.txt",skip=1)

  colnames(germlineres2M) <- germlineHeader

  germlineres2M <- germlineres2M[,-1]

  i <- i+1
  resultsList[[i]] <- germlineres2M

  methodN <- c(methodN,"germline2")
}

if (dodash) {

  dashresM <- read.table(file="dashres.txt",skip=1)

  colnames(dashresM) <- dashHeader

  dashresM <- dashresM[,c(-1,-2)]

  i <- i+1
  resultsList[[i]] <- dashresM

  methodN <- c(methodN,"dash-----")
}


write(paste("Results: ",sep=""),file="resultTest.txt",ncolumns=1000)

write(paste("Wilcoxon-Rank-Sum Wilcoxon-median t-test t-test-mean mean1 (std) mean2 (std) median1 median2",sep=""),file="resultTest.txt",,append=TRUE,ncolumns=1000)

for (j in 33:1) {

write(paste("######",resHeader[j],"######",sep=""),file="resultTest.txt",append=TRUE,ncolumns=1000)

for (i1 in 1:(i-1)) {
  for (i2 in (i1+1):i) {

 
    v1 <- resultsList[[i1]][,j]
    v2 <- resultsList[[i2]][,j]
    
    median1 <- median(v1)
    median2 <- median(v2)

    mean1 <- mean(v1)
    mean2 <- mean(v2)

    std1 <- sd(v1)
    std2 <- sd(v2)

    diffs <- length(which(abs(v1-v2)<0.000001))
    l1 <- length(v1)
    ddl <- l1-diffs 
    
    if ((length(table(v1))>1)&&(length(table(v2))>1)&&(length(table(v1-v2))>1)&&(abs(median1-median2)>0.0001)&&(ddl>20)) {
      resTest <- wilcox.test(v1,v2,paired=TRUE,conf.int=TRUE)
      pvalueW <- resTest$p.value
      estimW <- resTest$estimate
    } else {
      pvalueW <- 1.0
      estimW <- 0.0
    }

    if ((length(table(v1))>1)&&(length(table(v2))>1)&&(length(table(v1-v2))>1)&&(abs(mean1-mean2)>0.0001)) {

      resTest <- t.test(v1,v2,paired=TRUE)
      pvalueT <- resTest$p.value
      estimT <- resTest$estimate
    } else {
      pvalueT <- 1.0
      estimT <- 0.0
    }

    if ((abs(median1)>1)||(abs(median2)>1)) {
      nsmallV <- 0
      digitsV <- 3
    } else {
      nsmallV <- 6
      digitsV <- 3
    }
      
     write(paste("--->",methodN[i1],"---",methodN[i2],": ",format(pvalueW,digits=digitsV,nsmall=nsmallV,width=9,scientific=TRUE)," ",format(estimW,digits=digitsV,nsmall=nsmallV,width=9)," ",format(pvalueT,digits=digitsV,nsmall=nsmallV,width=9,scientific=TRUE)," ",format(estimT,digits=digitsV,nsmall=nsmallV,width=9)," ",format(mean1,digits=digitsV,nsmall=nsmallV,width=9)," (",format(std1,digits=digitsV,nsmall=nsmallV,width=9),") ",format(mean2,digits=digitsV,nsmall=nsmallV,width=9)," (",format(std2,digits=digitsV,nsmall=nsmallV,width=9),") ",format(median1,digits=digitsV,nsmall=nsmallV,width=9)," ",format(median2,digits=digitsV,nsmall=nsmallV,width=9),sep=""),file="resultTest.txt",append=TRUE,ncolumns=1000)



  }
  
}


}