##########################################################################
# This program reads in four text files and output a simplied version of 
# table1 in the paper.
# 
# The four input files:
# vircotable1PI.txt
# vircotable1RT.txt
# virologictable1PI.txt
# virologictable1RT.txt
#
# Output file format:
# Refer to table1
# 
# In this program, parametric and non-parametric testing were both done
# to test whether the standard deviation and median absolute deviance
# were significantly different between the two assays for wildtype isolates.
##########################################################################

baseurl <- 'http://hivdb.stanford.edu/pages/published_analysis/phenotypeCompJAID2005/transformedData/'
geturl <- function(tailend) {
  return(paste(baseurl, tailend,sep=''))
}

.eval <- function(evaltext) {
  eval(parse(text=evaltext), envir=sys.frame())
}

assays <- c('virco', 'virologic')
assaynames <- list(virco='AV', virologic='PS')
proteins <- c('PI','RT')

tables <- list()

for (protein in proteins) {
  for (assay in assays){
    filename <- paste(assay, 'table1', protein, '.txt', sep='')
    table.name <- paste(assay, protein, '.table', sep='')
    .eval(paste(table.name, ' <- read.table(geturl(filename), header=T)', sep=''))
    .eval(paste('cur.table <- ', table.name))
    cur.table$logFold <- log10(cur.table$Fold)

    drug <- unique(cur.table$Drug)
    ndrug <- length(drug)

    mean <- numeric(ndrug)
    N <- numeric(ndrug)
    exp.mean <- numeric(ndrug)
    sd <- numeric(ndrug)
    median <- numeric(ndrug)
    MAD <- numeric(ndrug)
    cMAD <- numeric(ndrug)

    for (i in 1:ndrug) {
      log.fold <- cur.table$logFold[(cur.table$Drug == drug[i])]
      fold <- 10^log.fold
      mean[i] <- 10^mean(log.fold)
      N[i] <- length(log.fold)
      sd[i] <- sd(log.fold)
      median[i] <- median(fold)
      cMAD[i] <- mad(log.fold)
      MAD[i] <- mad(log.fold, constant=1.0)
    }
    table.name <- paste(assay, protein, '.out.table', sep='')
    .eval(paste(table.name, ' <- data.frame(drug, N, mean, sd, median, cMAD, MAD)'))
    .eval(paste('cur.table <- ', table.name, sep=''))
  }

  fligner <- numeric(ndrug)
  fligner2 <- numeric(ndrug)
  ftest <- numeric(ndrug)

  for (i in 1:ndrug) {
    virco.table <- .eval(paste('virco', protein, '.table', sep=''))
    virologic.table <- .eval(paste('virologic', protein, '.table', sep=''))
    virologic <- log10(virologic.table$Fold[(virologic.table$Drug==drug[i])])
    virco <- log10(virco.table$Fold[(virco.table$Drug==drug[i])])
    fligner[i] <- fligner.test(list(virco, virologic))$p.value
    #fligner2[i] <- fligner.test(list(abs(virco - median(virco)), abs(virologic - median(virologic))))$p.value
    ftest[i] <- var.test(virco, virologic)$p.value
  }
  table.name <- paste(protein, '.table', sep='')
  .eval(paste(table.name, ' <- data.frame(drug)', sep=''))
  .eval(paste('cur.table <- ', table.name, sep=''))
  for (stat in c('N', 'mean', 'sd')) {
    for (assay in assays) {
      .eval(paste(table.name, '$', assaynames[assay], '.', stat, ' <- ', assay, protein, '.out.table$', stat, sep=''))
    }
  }
  .eval(paste(table.name, '$ftest <- ftest', sep=''))
  for (stat in c('median', 'MAD')) {
    for (assay in assays) {
      .eval(paste(table.name, '$', assaynames[assay], '.', stat, ' <- ', assay, protein, '.out.table$', stat, sep=''))
    }
  }
  .eval(paste(table.name, '$fligner <- fligner', sep=''))
  #.eval(paste(table.name, '$fligner2 <- fligner2', sep=''))
  #.eval(paste(table.name, '$fligner.lt.ftest <- (fligner < ftest)', sep=''))
  .eval(paste('cur.table <- ', table.name))
  write.table(cur.table, file=paste('Table1', protein, '.csv', sep=''), row.names=F, col.names=T, quote=F, sep=',')
  print(cur.table, digits=2)
}