# efg, Stowers Institute, 24 April 2006
# Modified 11 May 2006.

#source("0-LoadData.R")

ShowReplicateHistograms <- function()
{
  counts <- table(Info$ABTID)
  duplicates <- names(counts[counts > 1])

  duplicate.indices <- which(Info$ABTID %in% duplicates)


  # By Individual Patient
  pdf("Replicates.pdf", height=8, width=8)
  par(mar=c(6,4,2,1), oma=c(1,0,1,0))

  dlog2.subset <- log2d[duplicate.indices, !ControlSet]
  d.subset     <- d[duplicate.indices, !ControlSet]
  Info.subset  <- Info[duplicate.indices, ]

  boxplot(dlog2.subset ~ row(dlog2.subset), xaxt="n",
    main=paste("Log2[Raw Data] for Replicates (", ncol(dlog2.subset), " probes)", sep=""),
    col=c("white","white", "grey90", "grey90"))
  axis(SOUTH<-1, at=1:length(duplicate.indices), labels=Info$ABTID[duplicate.indices], las=3)
  text(-0.5 + 2*(1:8), 0.5, c("L", "", "", "W", "L", "", "M", "M"))

  PlotFooter("4-Replicates 1")

  # Row Medians - Include Zeros
  row.medians <- apply(d.subset, BY.ROW<-1, median)
  d.norm1 <- sweep(dlog2.subset, BY.ROW<-1, log2(150) - log2(row.medians), FUN="+")

  boxplot(d.norm1 ~ row(d.norm1), xaxt="n",
    main=paste("Median Scaled to 150 (", ncol(dlog2.subset), " probes) [Include 0s]", sep=""),
    col=c("white","white", "grey90", "grey90"))
  axis(SOUTH<-1, at=1:length(duplicate.indices), labels=Info$ABTID[duplicate.indices], las=3)
  text(-0.5 + 2*(1:8), 0.1, c("L", "", "", "W", "L", "", "M", "M"))
  PlotFooter("4-Replicates 2a")

  # Row Medians - Exclude Zeros
  row.medians.ExcludeZeros <- numeric(nrow(d.subset))
  for (j in 1:nrow(d.subset))
  {
    PatientExpression <- d.subset[j,]
    PatientExpression <- PatientExpression[PatientExpression != 0]
    row.medians.ExcludeZeros[j] <- median(PatientExpression)
  }
  boxplot(d.norm1 ~ row(d.norm1), xaxt="n",
    main=paste("Median Scaled to 150 (", ncol(dlog2.subset), " probes) [Include 0s]", sep=""),
    col=c("white","white", "grey90", "grey90"))
  axis(SOUTH<-1, at=1:length(duplicate.indices), labels=Info$ABTID[duplicate.indices], las=3)
  text(-0.5 + 2*(1:8), 0.1, c("L", "", "", "W", "L", "", "M", "M"))

  points(1:length(row.medians),
         log2(row.medians.ExcludeZeros), pch=19, col="red")

  d.norm2 <- sweep(dlog2.subset, BY.ROW<-1,
             log2(150) - log2(row.medians.ExcludeZeros), FUN="+")
  PlotFooter("4-Replicates 2b")

  par(mfcol=c(2,1))

  # "Raw Data
  for (j in 1:nrow(dlog2.subset))
  {
    hist(dlog2.subset[j,], breaks=200, xlab="", ylim=c(0,300),
         main=paste(rownames(dlog2.subset)[j], " (", ncol(dlog2.subset), " probes)", sep=""))

    text(0.3, 300,
       paste( sprintf("%.1f",
              100 * sum(dlog2.subset[j,] == 0) / ncol(dlog2.subset)), "% 0s", sep=""),
       adj=0, col="blue")


    if (j %%  2 == 1)
    {
      mtext(paste("REPLICATES", Info.subset$CLUSTER[j], Info.subset$Onset[j],
            Info.subset$sex[j], Info.subset$race[j]),
            NORTH<-3, adj=0.025, outer=TRUE, cex=0.8, line=0, col="blue")
    } else {
      PlotFooter("4-Replicates")
    }

  }

  # Raw Median Scaled to 150 (ignore 0s)

 for (j in 1:nrow(d.norm2))
  {
    hist(d.norm2[j,], breaks=200, xlab="",
         main=rownames(d.norm2)[j], ylim=c(0,300))


    if (j %%  2 == 1)
    {
      mtext(paste("REPLICATES MedianScaledNo0s", Info.subset$CLUSTER[j], Info.subset$Onset[j],
            Info.subset$sex[j], Info.subset$race[j]),
            NORTH<-3, adj=0.025, outer=TRUE, cex=0.8, line=0, col="blue")
    } else {
      PlotFooter("4-Replicates")
    }

  }

  dev.off()
}

sink("4-Replicates.txt")
ShowReplicateHistograms()
sink()