aroma-project.org

Author: Mark Robinson
Created on: 2009-03-28
Last updated: 2009-10-21

Christian Stratowa had been asking some questions about whether FIRMA could be used with median polish estimates (as opposed to the current standard of using affyPLM -- robust linear model -- estimates).  He has written some test code to implement FIRMA based on either and compared the results.  I have shared his code [tidied up by Henrik Bengtsson] and results below:

# - - - - - - - - - - - - - - - - - - -
firma <- function(m, method=c("rlm", "mdp")) {
## function to compute FIRMA scores (adapted from aroma.affymetrix)
## m:      dataframe containing normalized probe intensities (linear scale)
##         with 1.column containing the transcript_cluster_id
##         and 2. column containing the corresponding probeset_ids
## method: fitting model, one of "rlm" or "mdp"
  require("matrixStats") || throw("Package not loaded: matrixStats");

  # Argument 'method':
    method <- match.arg(method);

  ## convert expression levels to log2
  y <- as.matrix(log2(m[,3:ncol(m)]));

  ## dimensions
  K <- nrow(y);   # number of probes
  I <- ncol(y);   # number of arrays

  ## fit log-additive model
  if (method == "rlm") {
     fit <- .Call("R_rlm_rma_default_model", y, 0, 1.345,
                  PACKAGE="preprocessCore");
  } else if (method == "mdp") {
     mp <- medpolish(y, trace.iter=FALSE);
     fit <- list(Estimates = c(mp$overall + mp$col,  mp$row),
                 StdErrors = rep(0, length(c(mp$row, mp$col))));
  } # if

  ## extract parameters
  est <- fit$Estimates;
  se  <- fit$StdErrors;

  ## chip effects
  beta <- est[1:I];

  ## probe affinities
  if (K == 1) {
     ## if only one probe must have affinity=1 since sum constraint
     alpha <- 0;
  } else {
     ## affinities sum to zero (on log scale)
     alpha <- est[(I+1):length(est)];
     alpha[length(alpha)] <- -sum(alpha[1:(length(alpha)-1)]);
  } # if

  ## estimates on the intensity scale
  theta <- 2^beta;
  phi   <- 2^alpha;

  ## calculate residuals
  phi   <- matrix(phi,   nrow=K, ncol=I, byrow=FALSE);
  theta <- matrix(theta, nrow=K, ncol=I, byrow=TRUE);
  yhat  <- phi*theta;
  eps   <- (2^y) / yhat;  # RMA uses y/yhat

  ## estimate of standard error
  z <- log2(eps);
  z <- unlist(z, use.names=FALSE);
  u.mad <- mad(z, center=0);

  ## get probeset_ids
  id <- unique(m[,2]);

  ## FIRMA scores
  fs <- sapply(id, function(unitName) {
    rr <- which(m[,2] == unitName);
    z <- eps[rr,,drop=FALSE];
    z <- log2(z) / u.mad;
    colMedians(z);
  });
  fs <- t(fs);
  rownames(fs) <- id;

  fs;
} # firma()
# - - - - - - - - - - - - - - - - - - -

Then the following commands will give the image below.

library("preprocessCore");
score.mdp <- firma(unr, "mdp");
score.rml <- firma(unr, "rlm");

par(pty="m", mfcol=c(3,2), mar=c(5,5,4,2));
for (ii in 1:6) {
  plot(score.rml[,ii], type="l", ylim=c(-3,5),
       main=colnames(score.rml)[ii], xlab="Probeset", ylab="Score");
  abline(h=0, lty=3);
  lines(score.mdp[,ii], lty=2);
}

Figure: FIRMA scores for 30 probesets when using RLM (solid) and median polish (dashed) estimators to fit the FIRMA model.