join <- function(...) {
  paste(..., sep = "")
}

#
# Global values for sharing scales between sites (comment out for site-specific scales)
#

# maxHistCount = 7000
# maxHistRunCount = 26000
# maxRunCount = 70
# maxPepRunCount = 55

#
# Reading data
#

command_args <- commandArgs(trailingOnly = FALSE)
script_args <- commandArgs(trailingOnly = TRUE)
# print(command_args)

isCommand = length(script_args) > 0
if (isCommand) {
  reportPath <- script_args[1]
} else {
  reportPath = "C:/DIA_Course/Data/DIA_data/6600_64w_ppm10_rt10/PeakAreas.csv"
}
pdfPath = join(dirname(reportPath), "/plot/PeakAreasCvs.pdf")

cat("Reading ", reportPath, "\n")
cat("Generating", pdfPath, "\n")
peaksTable = read.csv(reportPath, stringsAsFactors = F, na.strings = "#N/A")
#sapply(peaksTable, class)
#head(peaksTable)
#summary(peaksTable)

#
# Summarize data for plotting
#

setNaAreas <- function(areasTable, acol, qcol, qcutoff = 0.01) {
  areasTable[which(areasTable[[qcol]] >= qcutoff),][[acol]] <- NA
  return (areasTable)
}

countReps = (ncol(peaksTable) - 4)/2
acols = c() # peak area columns
qcols = c() # q value columns
ucols = c() # q value columns in unique peptide table
for (n in 1:countReps) {
  # 4 descriptive columns followed by area, q value pairs
  acol = 3 + n*2
  acols = c(acols, acol)
  qcol = acol + 1
  qcols = c(qcols, qcol)
  peaksTable = setNaAreas(peaksTable, acol, qcol)
  ucols = c(ucols, 1 + n) # just best q value for each unique peptide
}
uniqueTable = aggregate(x = peaksTable[,qcols], by = list(peaksTable$ModifiedSequence), FUN = "min")

cat("queried:  precursors = ", nrow(peaksTable), "\n",
    "          unique pep = ", nrow(uniqueTable), "\n", sep = "")
peaksTable$Count = apply(peaksTable[,acols], 1, function(v) {length(which(!is.na(v)))})
#peaksTable$Count = apply(peaksTable[,qcols], 1, function(v) {length(which(v < 0.01))})
peaksTable$Min = apply(peaksTable[,qcols], 1, min)
peaksTable$Max = apply(peaksTable[,qcols], 1, max)
mqcols = c()
xqcols = c()
for (n in 1:countReps) {
  peaksTable[[join("Min", n)]] = apply(peaksTable[,qcols[1:n],drop=FALSE], 1, min)
  mqcols = c(mqcols, ncol(peaksTable))
  peaksTable[[join("Max", n)]] = apply(peaksTable[,qcols[1:n],drop=FALSE], 1, max)
  xqcols = c(xqcols, ncol(peaksTable))
}
uniqueTable$Min = apply(uniqueTable[,ucols], 1, min)
uniqueTable$Max = apply(uniqueTable[,ucols], 1, max)
mucols = c()
xucols = c()
for (n in 1:countReps) {
  uniqueTable[[join("Min", n)]] = apply(uniqueTable[,ucols[1:n],drop=FALSE], 1, min)
  mucols = c(mucols, ncol(uniqueTable))
  uniqueTable[[join("Max", n)]] = apply(uniqueTable[,ucols[1:n],drop=FALSE], 1, max)
  xucols = c(xucols, ncol(uniqueTable))
}

acceptCount <- function(colName) {length(which(peaksTable[[colName]] < 0.01))}
acceptPercent <- function(colName) {round(acceptCount(colName) / nrow(peaksTable) * 100, 1)}
acceptText <- function(colName) {join(acceptCount(colName), " - ", acceptPercent(colName), "%")}
acceptUCount <- function(colName) {length(which(uniqueTable[[colName]] < 0.01))}
acceptUPercent <- function(colName) {round(acceptUCount(colName) / nrow(uniqueTable) * 100, 1)}
acceptUText <- function(colName) {join(acceptUCount(colName), " - ", acceptUPercent(colName), "%")}
acceptSet <- function(cols, acceptFunc) {
  acceptSet = c()
  for (col in cols) {
    acceptSet = c(acceptSet, acceptFunc(col))
  }
  return (acceptSet)
}

getModifierText <- function(percentile, caps = FALSE)
{
  if (percentile < 1.0) {
    if (caps) {
      return ("at Least")
    } else {
      return ("at least")
    }
  }
  if (caps) {
    return ("All")
  }
  return ("all")
}

getRunsText <- function(runs, percentile, caps = FALSE)
{
  text = paste(getModifierText(percentile, caps), round(runs*percentile))
  if (percentile < 1.0) {
    text = paste(text, join("(of ", runs, ")"))
  }
  return (text)
}

catAcceptAt <- function(colName, runs, percentile) {
  cat(getRunsText(runs, percentile), " runs: precursors = ", acceptText(colName), "\n",
      "          unique pep = ", acceptUText(colName), "\n", sep = "")
}

detectRuns = acceptSet(qcols, acceptCount)
detectRunsGrowth = acceptSet(mqcols, acceptCount)
detectRunsComplete = acceptSet(xqcols, acceptCount)
detectRunsU = acceptSet(ucols, acceptUCount)
detectRunsGrowthU = acceptSet(mucols, acceptUCount)
detectRunsCompleteU = acceptSet(xucols, acceptUCount)

percentRunDetections = round(mean(detectRuns) / acceptCount("Min") * 100, 1)
percentRunDetectionsU = round(mean(detectRunsU) / acceptUCount("Min") * 100, 1)

cat("detected: precursors = ", acceptText("Min"), " (", paste(detectRuns, collapse = ", "), ", mean = ", mean(detectRuns), ") - ", percentRunDetections, "%\n",
    "          unique pep = ", acceptUText("Min"), " (", paste(detectRunsU, collapse = ", "), ", mean = ", mean(detectRunsU), ") - ", percentRunDetectionsU, "%\n", sep = "")

catAcceptAt("Max", countReps, 1.0)

# percentiles = c(0.9, 0.8, 0.7, 0.6, 0.5)
percentiles = c(0.5)
percentileColumn <- function(p) {join("P", p*100)}
percentileText <- function(p) {join(p*100, "% complete")}
for (p in percentiles) {
  colPercent = percentileColumn(p)
  peaksTable[[colPercent]] = apply(peaksTable[,qcols], 1, function(v) {quantile(v, p, na.rm = TRUE)})
  uniqueTable[[colPercent]] = apply(uniqueTable[,ucols], 1, function(v) {quantile(v, p, na.rm = TRUE)})
  catAcceptAt(colPercent, countReps, p)
}

# Calculate CVs and plot
meanIgnoreNA <- function(v) {mean(v, na.rm = TRUE)}
stdevIgnoreNA <- function(v) {sd(v, na.rm = TRUE)}
peaksTable$Mean.Total.Area = apply(peaksTable[,acols], 1, meanIgnoreNA)
peaksTable$SD.Total.Area = apply(peaksTable[,acols], 1, stdevIgnoreNA)
peaksTable$CV.Total.Area = peaksTable$SD.Total.Area/peaksTable$Mean.Total.Area

#
# Plotting
#

# Runs bargraph plot
plotRunsBarAndLine <- function(itemType, runsCounts, runsGrowth, runsComplete, count50, isPeps = FALSE) {
  factor = 1000
  runsCounts = runsCounts/factor
  runsGrowth = runsGrowth/factor
  runsComplete = runsComplete/factor
  count50 = count50/factor
  runsMean = round(mean(runsCounts), 1)
  runsSD = round(sd(runsComplete), 1)
  
  # Plot run bars
  ymax = max(runsGrowth) * 1.05;
  if (!isPeps & exists("maxRunCount")) {
    ymax = maxRunCount
  }
  if (isPeps & exists("maxPepRunCount")) {
    ymax = maxPepRunCount
  }
  df.bar = barplot(runsCounts, names.arg = 1:length(runsCounts), ylim = c(0, ymax),
                   main = paste(itemType, "Detections by Run"), xlab = "Run Number", ylab = "Detections (thousands)", col = blues1[5])
  text(0, ymax, paste("mean:", runsMean, "\nstddev:", runsSD), srt = 0.2, adj = c(0, 1))
  
  # Plot cumulative line
  colCumulative = r2[23]
  lines(x = df.bar, y = runsGrowth, lwd = 2, col = colCumulative)
  cumulativeLast = round(runsGrowth[length(runsGrowth)], 1)
  text(countReps, cumulativeLast, cumulativeLast, srt=0.2, adj=c(0, 0), col = colCumulative)
  
  # Complete line
  colComplete = blues1[9]
  lines(x = df.bar, y = runsComplete, lwd = 2, col = colComplete)
  completeLast = round(runsComplete[length(runsComplete)], 1)
  text(countReps, completeLast, completeLast, srt=0.2, adj = c(0, 1), col = colComplete)
  
  # 50% complete line
  col50 = "red"
  abline(h = count50, lty = 2, col = col50, lwd = 2)
  
  legend(1, 1, yjust = 0, c("cumulative", getRunsText(countReps, 0.5), "all runs"),
         lty = c(1, 2, 1), lwd = c(2, 2, 2), col = c(colCumulative, col50, colComplete), bg = "white")
}
# plotRunsBarAndLine("Precursors", detectRuns, detectRunsGrowth, detectRunsComplete, acceptCount("P50"))

# Histogram plot
plotHistCvs <- function(qcol, percentile)
{
  cvCutoff = 0.20 # Line and percent with CV less than 20%
  detectedCvs = peaksTable[which(peaksTable[[qcol]] < 0.01 & !is.na(peaksTable$CV.Total.Area)), ]$CV.Total.Area*100
  detectedCvsCut = peaksTable[which(peaksTable[[qcol]] < 0.01 & peaksTable$CV.Total.Area <= cvCutoff & !is.na(peaksTable$CV.Total.Area)), ]$CV.Total.Area*100
  medianCV = median(detectedCvs)
  percentLessThanCuttoff = round((length(detectedCvsCut)/length(detectedCvs))*100, 1)
  # Plot the histogram of CVs
  ymax = NULL
  if (exists("maxHistCount")) {
    hist(detectedCvs, breaks = max(detectedCvs)/2, xlim = c(0, 100), ylim = c(0, maxHistCount), col = blues1[3],
         main = paste("CVs of Precursors Detected in", getRunsText(countReps, percentile, TRUE), "Runs"), xlab = "CV (%)")
  } else {
    hist(detectedCvs, breaks = max(detectedCvs)/2, xlim = c(0, 100), col = blues1[3],
         main = paste("CVs of Precursors Detected in", getRunsText(countReps, percentile, TRUE), "Runs"), xlab = "CV (%)")
  }

  # Add median line
  abline(v = medianCV, lty = 2, col = "blue", lwd = 2)
  # Add less than cut-off line with percent label
  histMax = par("yaxp")[2] # Maximum y values
  segments(cvCutoff * 100, 0, cvCutoff * 100, histMax * 0.95, lty = 5, col = "red", lwd = 2)
  text(cvCutoff * 100, histMax * 0.95, join(format(percentLessThanCuttoff), "%"), srt=0.2, pos=3, col = "red")

  cat(getRunsText(countReps, percentile), " runs cvs: median = ", round(median(detectedCvs), 1), "%, below ", cvCutoff * 100, " = ", percentLessThanCuttoff, "%, na = ",
      length(which(peaksTable[[qcol]] < 0.01 & is.na(peaksTable$CV.Total.Area))), "\n", sep = "")
}
# plotHistCvs("P50", 1.0)

# 2D histogram colors
invisible(library(gplots))
library(RColorBrewer)
rf <- colorRampPalette(rev(brewer.pal(11,'Spectral')))
r <- rf(32)
r2 = c(rgb(255,255,255,maxColorValue = 255), r[2:length(r)])
blues1 <- brewer.pal(9, 'Blues')
#display.brewer.all()

#plot(0, type="n", ylab="", xlab="",
#     axes=FALSE, ylim=c(1,0), xlim=c(1,length(r2)))

#for (i in 1:length(r2))
#{
#  rect(i-0.5,0, i+0.5,1, border="black", col=r2[i])
#}

# 2D histogram plot
plotCvsByIntensity <- function(qcol, percentile)
{
  peaksTable$CV.Total.Area.Percent = peaksTable$CV.Total.Area*100
  peaksTable$Log.Mean.Total.Area = log10(peaksTable$Mean.Total.Area)
  detectedCvs = peaksTable[which(peaksTable[[qcol]] < 0.01 & !is.na(peaksTable$CV.Total.Area) & peaksTable$CV.Total.Area < 0.5),
                           c("Log.Mean.Total.Area", "CV.Total.Area.Percent")]

  hist2d(detectedCvs, nbins=c(150, 50), col=r2, xlim = c(1, 5),
                   xlab = "Log10 Mean Area", ylab = "CV (%)", main = paste("CVs by Intensity of Precursors Detected in",
                                                                           getRunsText(countReps, percentile, TRUE), "Runs"))
  
  abline(h = 20, lty = 2, col = "red", lwd = 2)
}

pdf(pdfPath, width = 7, height = 5)

plotRunsBarAndLine("Precursors", detectRuns, detectRunsGrowth, detectRunsComplete, acceptCount("P50"))
plotRunsBarAndLine("Peptides", detectRunsU, detectRunsGrowthU, detectRunsCompleteU, acceptUCount("P50"), TRUE)

if (exists("maxHistRunCount")) {
  hist(peaksTable[which(peaksTable$Count > 0), c("Count")], ylim = c(0, maxHistRunCount),
       main = "Precursors Detected in Runs", xlab = "Run Count", col = blues1[3])
} else {
  hist(peaksTable[which(peaksTable$Count > 0), c("Count")],
       main = "Precursors Detected in Runs", xlab = "Run Count", col = blues1[3])
}

plotHistCvs("Max", 1.0)
for (p in percentiles) {
  plotHistCvs(percentileColumn(p), p)
}
plotCvsByIntensity("Max", 1.0)
for (p in percentiles) {
  plotCvsByIntensity(percentileColumn(p), p)
}

invisible(dev.off())
