rho statistic with Biostrings

rho-Biostrings.R

library(Biostrings)

setGeneric("rho", function(x, ...) standardGeneric("rho"))

setMethod("rho", c(x="XStringSet"),
function(x, ...) {
  if (!xsbasetype(x) %in% c("DNA", "RNA")) {
    stop("`x` must be of DNA or RNA base type")
  }
  di <- oligonucleotideFrequency(x, 2L, as.prob=TRUE)
  nt <- oligonucleotideFrequency(x, 1L, as.prob=TRUE)

  denom.1 <- substring(colnames(di), 1, 1)
  denom.2 <- substring(colnames(di), 2, 2)
  di / (nt[,denom.1,drop=FALSE] * nt[,denom.2,drop=FALSE])
})

setMethod("rho", c(x="XString"),
function(x, ...) {
  if (!xsbasetype(x) %in% c("DNA", "RNA")) {
    stop("`x` must be of DNA or RNA base type")
  }
  rho(as(x, 'XStringSet'))
})

setMethod("rho", c(x="character"),
function(x, base.type=c("DNA", "RNA"), ...) {
  base.type <- match.arg(base.type)
  base.type <- switch(base.type, DNA="DNAStringSet", RNA="RNAStringSet")
  rho(as(x, base.type))
})

################################################################################
## Tests
if (FALSE) {
  tests.str <- c("GATACA", "GGCCGGCC",
                 "ATGTATATGAGAAAGGAAGAGCCTAGCGGCTCAGACAAGATTATGACTTCAGTTGTTGTT")
  rho(tests.str)
}