R/rho.R

In seqinr: Biological Sequences Retrieval and Analysis

rho <- function (sequence, wordsize = 2, alphabet = s2c("acgt"))
{
  wordcount <- count(sequence, wordsize, freq = FALSE, alphabet = alphabet)
  uni <- count(sequence, 1, freq = TRUE, alphabet = alphabet)

  expected_wordfreq <- function (wordsize, uni) 
  {
    if (wordsize == 1) 
        return(uni)
    else kronecker(uni, expected_wordfreq(wordsize - 1, uni))
  }

  expected_wordcount <- sum(wordcount)*expected_wordfreq(wordsize, uni)
   
  return(wordcount/expected_wordcount)
}