R/RelativeAdaptiveness.R

In fredysiegrist/statanacoseq: Statistical Analyses of Codon Usage in Custom Genome Sequences

#' @title Compute Relative Adaptiveness
#'
#' @description
#' \code{RelativeAdaptiveness} calculates RA for a List of Entry Numbers
#'
#' @details
#' Should compute the same RA as in Darwin and should not have an error higher than (5e-04)
#' the difference of w from seqinr and relative adaptiveness (w) from codonW (by John Peden).
#'
#' @param entires List of Coding sequences or DNA string

#' @return Named (codons) numerical vector with relative adaptiveness for the 64 codons
#'
#' @author Roth, A.; Friberg, M.; Siegrist, F. and Cannarozzi, G. M. \email{gina@@cannarozzi.com}
#' @seealso  \code{\link{seqinr}} \code{\link{ComputeCarboneRA}} \code{\link{statanacoseq}} \code{\link{SetupRA}}
#' @keywords CodonBias
#' @examples
#' RelativeAdaptiveness('ATGTGGTACTCCGACTACGGAGGATAA')
#' RelativeAdaptiveness(mylist(whatout=1))
#' @import seqinr
#'
#' @export
#' @section Original code in Darwin:
#' \subsection{Compute Relative Adaptiveness as codonW (by John Peden)}{\preformatted{
#' RelativeAdaptiveness := proc(entries:list(posint))
#'  CodonCounts := CreateArray(1..64);
#'  for i in entries do
#'    dna := SearchTag('DNA', Entry(i));
#'    for j to length(dna) by 3 do
#'      cod := CodonToCInt(dna[j..j+2]);
#'      if cod=0 then next fi;   # to avoid XXX
#'      CodonCounts[cod] := CodonCounts[cod]+1;
#'    od;
#'  od;
#'  RA := CreateArray(1..64);
#'  aa := 1;
#'  for aa to 20 do
#'    codons := IntToCInt(aa);
#'    counts := [seq(CodonCounts[i], i=codons)];
#'    freqs := counts / sum(counts);
#'    for i to length(codons) do
#'      cod := codons[i];
#'      RA[cod] := freqs[i] / max(freqs);
#'    od;
#'  od;
#'  for i to length(RA) do       # set minimum RA value to 0.01
#'    if RA[i] = 0 then
#'      RA[i] := 0.01 fi od;
#'    for i in AToCInt('$') do     # set RA value of stop codons to 1
#'      RA[i] := 1; od;
#'  RA
#' end: } }
RelativeAdaptiveness <- function(entries=mylist(whatout=1)) {
  CodonCounts <- aa_ac
  for (i in 1:length(entries)) {
    dna <- c2s(entries[[i]])
    if(!(checkCDS(dna))) stop("non valid CDS)", call.=FALSE)
    for (j in seq(from=1, to=nchar(dna), by=3)) {
      codon <- substr(dna, j, j+2)
      cod <- reversecomplement(codon)
      codn <- which(substr(names(aa_ac), 5, 7) %in% cod)
      if (codn==65) break # to avoid XXX
      CodonCounts[codn] <- CodonCounts[codn]+1
    }
  }
  RA <- aa_ac
  for (aa in 1:20) {
    a <- unique(a(substr(names(aa_ac), 1, 3)[c(-59, -60, -64, -65)]))[aa]
    # warnings
    suppressWarnings(codons <-  na.omit(substr(names(CodonCounts)[a(substr(names(aa_ac), 1, 3))==a], 5, 7)))
    # warnings
    suppressWarnings(counts <- na.omit(CodonCounts[a(substr(names(aa_ac), 1, 3))==a]))
    freqs <- counts / sum(counts)
    for (i in 1:length(codons)) {
      cd <- na.omit((1:65)[substr(names(aa_ac), 5, 7)==codons[i]])
      RA[cd] <- freqs[i] / max(freqs)
    }
  }
  RA[RA<0.01] <- 0.01       # set minimum RA value to 0.01
  # warnings
  suppressWarnings(RA[is.na(a(substr(names(aa_ac), 1, 3)))] <- 1)     # set RA value of stop codons to 1
  return(RA)
}