R/statisticPlate.R

In titanQC: QC Procedures for the GeneTitan platform

#' Adds a statistic to the pData of the ExpressionSet object that indicates certain types
#' of deviations in the MA plots
#' @param eset ExpressionSet as generated by the 'preprocessing' function 
#' @return object of class ExpressionSet
#' @note this function is experimental and currently only a research idea
#' @export
statisticPlate <- function(eset) {
  
  pData(eset) <- cbind(pData(eset), statistic = NA)	
  
  for (plate in sort(unique(eset$titanPlateNo))) {
    
    MedPlate <- apply(exprs(eset)[, eset$titanPlateNo == plate], 1, median)
    esettemp <- eset[, eset$titanPlateNo == plate]
    Q1 <- quantile(MedPlate, probs = 0.25)
    Q99 <- quantile(MedPlate, probs = 0.99)
    Q01 <- quantile(MedPlate, probs = 0.01)
    
    for (row in c("A","B","C","D","E","F","G","H")) {
      for (col in c(1:12)) {
        if (nrow(pData(esettemp)[esettemp$titanRow == row & esettemp$titanColumn == col, ]) != 0) {
          intensity <- exprs(esettemp)[, esettemp$titanRow == row & esettemp$titanColumn == col]
          logratio <- intensity - MedPlate
          lowessval25 <- lowess(MedPlate[MedPlate < Q1], logratio[MedPlate < Q1])$y
          centval <- mean(c(mean(logratio[MedPlate < Q1]), mean(logratio[MedPlate > Q99])))
          statistic <- round(sum(lowessval25 - centval) / length(lowessval25), 2)
          pData(eset)[eset$titanPlateNo == plate & eset$titanRow == row & eset$titanColumn == col, "statistic"] <- statistic 
        }
      } 
    }
  }	
  eset	
}