R/robustness.R

Defines functions ppa.filter.robust.default isa.filter.robust.default robustness.ppa robustness.isa robustness.default

setMethod("robustness", signature(normed.data="list"),
          function(normed.data, ...) robustness.default(normed.data, ...))

robustness.default <- function(normed.data, ...) {
  if (length(normed.data)==2) {
    robustness.isa(normed.data, ...)
  } else if (length(normed.data)==4) {
    robustness.ppa(normed.data, ...)
  }
}

robustness.isa <- function(normed.data, row.scores, col.scores) {

  isa.status("Calculating ISA robustness", "in")

  if (ncol(row.scores) != ncol(col.scores)) {
    stop("Row and column column dimensions don't match")
  }
  
  if (ncol(row.scores)==0) return(numeric())
  
  Ec <- normed.data$Ec
  Er <- normed.data$Er
  
  row.scores <- apply(row.scores, 2, function(x) x/sqrt(sum(x^2)))
  col.scores <- apply(col.scores, 2, function(x) x/sqrt(sum(x^2)))
  if ("hasNA" %in% names(attributes(normed.data)) && !attr(normed.data, "hasNA")) {
    rob1 <- colSums(col.scores * Er %*% row.scores)
    rob2 <- colSums(row.scores * Ec %*% col.scores)
  } else {
    rob1 <- colSums(col.scores * na.multiply(Er, row.scores))
    rob2 <- colSums(row.scores * na.multiply(Ec, col.scores))
  }

  rob1[ rob1 < 0 ] <- 0
  rob2[ rob2 < 0 ] <- 0
  res <- sqrt(rob1) * sqrt(rob2)

  isa.status("DONE", "out")

  res
}

robustness.ppa <- function(normed.data, row1.scores, row2.scores,
                           col.scores) {

  isa.status("Calculating PPA robustness", "in")

  if (!is.matrix(row1.scores) || !is.matrix(row2.scores) ||
      !is.matrix(col.scores)) {
    stop("`row1.scores', `row2.scores' and `col.scores' must be matrices")
  }
  
  if (ncol(row1.scores) != ncol(row2.scores) ||
      ncol(row1.scores) != ncol(col.scores)) {
    stop("`row1.scores', `row2.scores' and `col.scores' must have the same",
         "number of columns")
  }

  if (ncol(row1.scores)==0) return (numeric())

  row1.scores <- apply(row1.scores, 2, function(x) x/sqrt(sum(x^2)))
  row2.scores <- apply(row2.scores, 2, function(x) x/sqrt(sum(x^2)))
  col.scores  <- apply(col.scores,  2, function(x) x/sqrt(sum(x^2)))

  hasNA <- c(TRUE, TRUE)
  if ("hasNA" %in% names(attributes(normed.data))) {
    hasNA <- attr(normed.data, "hasNA")
  }

  if (!hasNA[1]) {
    rob1 <- colSums( col.scores * normed.data$Egc %*% row1.scores)
    rob4 <- colSums(row1.scores * normed.data$Ecg %*%  col.scores)
  } else {
    rob1 <- colSums( col.scores * na.multiply(normed.data$Egc, row1.scores))
    rob4 <- colSums(row1.scores * na.multiply(normed.data$Ecg,  col.scores))
  }

  if (!hasNA[2]) {
    rob2 <- colSums(row2.scores * normed.data$Ecd %*%  col.scores)
    rob3 <- colSums( col.scores * normed.data$Edc %*% row2.scores)
  } else {
    rob2 <- colSums(row2.scores * na.multiply(normed.data$Ecd,  col.scores))
    rob3 <- colSums( col.scores * na.multiply(normed.data$Edc, row2.scores))
  }

  rob1[ rob1 < 0 ] <- 0
  rob2[ rob2 < 0 ] <- 0
  rob3[ rob3 < 0 ] <- 0
  rob4[ rob4 < 0 ] <- 0
  
  res <- sqrt(rob1) * sqrt(rob2) * sqrt(rob3) * sqrt(rob4)

  isa.status("DONE", "out")

  res
}

setMethod("isa.filter.robust", signature(data="matrix"),
          function(data, ...) isa.filter.robust.default(data, ...))

isa.filter.robust.default <- function(data, normed.data, isares, perms=1,
                                      row.seeds, col.seeds) {
  
  isa.status("Filtering ISA result for robustness...", "in")
  
  if (perms <= 0) {
    stop("Number of permutations must be positive")
  }

  if (length(isares$rows) == 0) {
    return(isares)
  }
  
  if (length(unique(isares$seeddata$thr.row)) != 1) {
    warning("Different row thresholds, using only the first one.")
  }
  
  if (length(unique(isares$seeddata$thr.col)) != 1) {
    warning("Different column thresholds, using only the first one.")
  }
  
  isares$seeddata$rob <- robustness(normed.data, isares$rows, isares$columns)

  if (missing(row.seeds)) {
    row.seeds <- generate.seeds(count=isares$rundata$N,
                                length=nrow(isares$rows))
  }
  if (missing(col.seeds)) {
    col.seeds <- matrix(nrow=nrow(isares$columns), ncol=0)
  }

  rob.max <- 0
  
  for (i in seq_len(perms)) {
    data.scrambled <- data
    data.scrambled[] <- sample(data.scrambled)
    normed.data.scrambled <- isa.normalize(data.scrambled)
    
    permres <- isa.iterate(normed.data.scrambled,
                           row.seeds=row.seeds, col.seeds=col.seeds,
                           thr.row=isares$seeddata$thr.row[1],
                           thr.col=isares$seeddata$thr.col[1],
                           direction=isares$rundata$direction,
                           convergence=isares$rundata$convergence,
                           cor.limit=isares$rundata$cor.limit,
                           eps=isares$rundata$eps,
                           oscillation=isares$rundata$oscillation,
                           maxiter=isares$rundata$maxiter)

    valid <- apply(permres$rows != 0, 2, any)
    valid <- valid & apply(permres$columns !=0, 2, any)

    permres$rows <- permres$rows[,valid,drop=FALSE]
    permres$columns <- permres$columns[,valid,drop=FALSE]
    permres$seeddata <- permres$seeddata[valid,,drop=FALSE]

    rob2 <- robustness(normed.data.scrambled, permres$rows, permres$columns)
    rob.max <- max(rob2, rob.max, na.rm=TRUE)
  }

  keep <- isares$seeddata$rob > rob.max

  isares$rows <- isares$rows[, keep,drop=FALSE]
  isares$columns <- isares$columns[, keep,drop=FALSE]
  isares$seeddata <- isares$seeddata[ keep,,drop=FALSE ]
  if (nrow(isares$seeddata)>0) { isares$seeddata$rob.limit <- rob.max }
  isares$rundata$rob.perms <- perms

  isa.status("DONE.", "out")
  
  isares
}

setMethod("ppa.filter.robust", signature(data="list"),
          function(data, ...) ppa.filter.robust.default(data, ...))

ppa.filter.robust.default <- function(data, normed.data, ppares, perms=1,
                                      row1.seeds, col1.seeds,
                                      row2.seeds, col2.seeds) {

  isa.status("Filtering PPA result for robustness...", "in")

  if (perms <= 0) {
    stop("Number of permutations must be positive")
  }

  if (length(ppares$rows1)==0) {
    return(ppares)
  }

  if (length(unique(ppares$seeddata$thr.row1)) != 1) {
    warning("Different row1 thresholds, using only the first one.")
  }
  if (length(unique(ppares$seeddata$thr.row2)) != 1) {
    warning("Different row2 thresholds, using only the first one.")
  }
  if (length(unique(ppares$seeddata$thr.col)) != 1) {
    warning("Different column thresholds, using only the first one.")
  }

  ppares$seeddata$rob <- robustness(normed.data, ppares$rows1, ppares$rows2,
                                    ppares$columns)

  if (missing(row1.seeds) && missing(col1.seeds) &&
      missing(row2.seeds) && missing(col2.seeds)) {
    row1.seeds <- generate.seeds(count=ppares$rundata$N,
                                 length=nrow(ppares$rows1))
    col1.seeds <- matrix(nrow=nrow(ppares$columns), ncol=0)
    row2.seeds <- matrix(nrow=nrow(ppares$rows2), ncol=0)
    col2.seeds <- matrix(nrow=nrow(ppares$columns), ncol=0)
  }

  rob.max <- 0

  for (i in seq_len(perms)) {
    data.scrambled <- data
    data.scrambled[[1]][] <- sample(data.scrambled[[1]])
    data.scrambled[[2]][] <- sample(data.scrambled[[2]])
    normed.data.scrambled <- ppa.normalize(data.scrambled)

    permres <- ppa.iterate(normed.data.scrambled,
                           row1.seeds=row1.seeds,
#                           col1.seeds=col1.seeds,
#                           row2.seeds=row2.seeds,
#                           col2.seeds=col2.seeds,
                           thr.row1=ppares$seeddata$thr.row1[1],
                           thr.row2=ppares$seeddata$thr.row2[1],
                           thr.col =ppares$seeddata$thr.col[1],
                           direction=ppares$rundata$direction,
                           convergence=ppares$rundata$convergence,
                           cor.limit=ppares$rundata$cor.limit,
                           oscillation=ppares$rundata$oscillation,
                           maxiter=ppares$rundata$maxiter)

    valid <- apply(permres$rows1 != 0, 2, any)
    valid <- valid & apply(permres$rows2 != 0, 2, any)
    valid <- valid & apply(permres$columns != 0, 2, any)

    permres$rows1 <- permres$rows1[,valid,drop=FALSE]
    permres$rows2 <- permres$rows2[,valid,drop=FALSE]
    permres$columns <- permres$columns[,valid,drop=FALSE]
    permres$seeddata <- permres$seeddata[valid,,drop=FALSE]

    rob2 <- robustness(normed.data.scrambled, permres$rows1, permres$rows2,
                       permres$columns)
    rob.max <- max(rob2, rob.max, na.rm=TRUE)
  }

  keep <- ppares$seeddata$rob > rob.max

  ppares$rows1 <- ppares$rows1[, keep, drop=FALSE]
  ppares$rows2 <- ppares$rows2[, keep, drop=FALSE]
  ppares$columns <- ppares$columns[, keep, drop=FALSE]
  ppares$seeddata <- ppares$seeddata[keep, , drop=FALSE]
  if (nrow(ppares$seeddata)>0) { ppares$seeddata$rob.limit <- rob.max }
  ppares$rundata$rob.perms <- perms  
  
  isa.status("DONE", "out")

  ppares
}

Try the isa2 package in your browser

Any scripts or data that you put into this service are public.

isa2 documentation built on May 29, 2017, 6:44 p.m.