# Analyze Control Probes in Array Data
# efg 23 March 2006.  Updated 24 April 2006.

ShowHistogramsForControlProbes <- function()
{

  pdf("All-Controls.pdf", height=8, width=8)
  par(oma=c(1,0,1,0), mfrow=c(4,1))

  # 1. Histograms by Control Subset
  N <- sum(ControlSet1)
  hist(log2d[, ControlSet1],
     breaks=200, ylim=c(0,500),
     xlab=expression(Log[2] * (Expression)),
     main=paste("Histogram of", N,
          "mwgaracontrol probes --", nrow(log2d), "microarrays"))

  N <- sum(ControlSet2)
  hist(log2d[, ControlSet2],
     breaks=200, ylim=c(0,35),
     xlab="Log2(Expression)",
     main=paste("Histogram of", N,
          "mwghumcontrol probes --", nrow(log2d), "microarrays"))

  N <- sum(ControlSet3)
  hist(log2d[, ControlSet3],
     breaks=200, ylim=c(0,160),
     xlab="Log2(Expression)",
     main=paste("Histogram of", N,
          "mwghuman probes --", nrow(log2d), "microarrays"))

  N <- sum(ControlSet4)
  hist(log2d[, ControlSet4],
     breaks=200, ylim=c(0,350),
     xlab="Log2(Expression)",
     main=paste("Histogram of", N,
          "blank probes --", nrow(log2d), "microarrays"))

  PlotFooter("3-ControlProbes.R 1")

  # 2. Histograms For All, Control Set, All Less Control Set
  par(mfrow=c(3,1))

  # 2.1 All
  hist(log2d,
     breaks=200, ylim=c(0,40000),
     xlab="Log2(Expression)",
     main=paste("Histogram of", ncol(log2d), "probes on ",
                 nrow(log2d), "microarrays"))
  text(0.3, 40000,
       paste( sprintf("%.1f",
              100 * sum(log2d == 0) / prod(dim(log2d))), "% 0s", sep=""),
       adj=0, col="blue")

  # 2.2 Control Set
  PlotSet <- log2d[, ControlSet]
  hist(PlotSet,
     breaks=200, ylim=c(0,800),
     xlab="Log2(Expression)",
     main=paste("Histogram of", ncol(PlotSet),
                "Control/Blank Probes --", nrow(PlotSet), "microarrays"))
  text(0.3, 800,
       paste( sprintf("%.1f",
              100 * sum(PlotSet == 0) / prod(dim(PlotSet))), "% 0s", sep=""),
       adj=0, col="blue")

  # 2.3a All Less Control Set
  PlotSet <- log2d[, !ControlSet]
  hist(PlotSet,
     breaks=200, ylim=c(0,40000),
     xlab="Log2(Expression)",
     main=paste("Histogram of", ncol(PlotSet),
                "Probes (less Control/Blank Probes) --", nrow(PlotSet), "microarrays"))
  text(0.3, 40000,
       paste( sprintf("%.1f",
              100 * sum(PlotSet == 0) / prod(dim(PlotSet))), "% 0s", sep=""),
       adj=0, col="blue")
  PlotFooter("3-ControlProbes.R 2a")

  par(mfrow=c(1,1))

  # 2.3b
  hist(PlotSet,
     breaks=200, ylim=c(0,40000),
     xlab="Log2(Expression)",
     main=paste("Histogram of", ncol(PlotSet),
                "Probes (less Control/Blank Probes) --", nrow(PlotSet), "microarrays"))
  text(0.3, 40000,
       paste( sprintf("%.1f",
              100 * sum(PlotSet == 0) / prod(dim(PlotSet))), "% 0s", sep=""),
       adj=0, col="blue")
  mtext("Median, including and excluding 0s", col="red", line=0.5)
  abline(v=median(PlotSet), col="red")
  NonZeros <- PlotSet[PlotSet != 0]
  abline(v=median(NonZeros), col="red")

  PlotFooter("3-ControlProbes.R 2b")


  # 3 Additional look at Control Subsets
  par(mfrow=c(1,1))

  PlotSet <- log2d[, ControlSet1]
  boxplot(PlotSet ~ col(PlotSet), main="mwgaracontrol probes")
  PlotFooter("3-ControlProbes.R 3a")

  PlotSet <- log2d[, ControlSet2]
  boxplot(PlotSet ~ col(PlotSet), main="mwghumcontrol probes")
  PlotFooter("3-ControlProbes.R 3b")

  PlotSet <- log2d[, ControlSet3]
  boxplot(PlotSet ~ col(PlotSet), main="mwghuman probes")
  PlotFooter("3-ControlProbes.R 3c")

  PlotSet <- log2d[, ControlSet4]
  boxplot(PlotSet ~ col(PlotSet), main="Blank probes probes")
  PlotFooter("3-ControlProbes.R 3d")


  dev.off()
}

sink("3-ControlProbes.txt")
ShowHistogramsForControlProbes()
sink()