Write an R function that takes an amino acid sequence and stochastically returns a codon string-CodePudding

Using R, I'd like to write a "reverse translation" function that takes as input a string of amino acids and for each one returns a triplet codon that encodes that amino acid. The codon would be randomly chosen from a frequency table that I will initially hard-code the values for, but I would also like to eventually be able to adjust the frequencies so that "rarer" codons will be more likely to be chosen.

So far I have tried to define a "codon" class using the 61 codons that encode an amino acid in the standard genetic code.

# Define codon class
codon <- function(sequence, aminoacid, frequency) {
  new_codon <- list(sequence = sequence, aminoacid = aminoacid, frequency = frequency)
  class(new_codon) <- "codon"
  new_codon
}

# Define the 61 codon frequency table
# structure: nucleotide sequence, amino acid it encodes, frequency codon is used
codons <- list(
  # Alanine (Ala) codons
  codon("GCT", "A", 0.26),
  codon("GCC", "A", 0.40),
  codon("GCA", "A", 0.23),
  codon("GCG", "A", 0.11),
  # Cysteine (Cys) codons
  codon("TGT", "C", 0.45),
  codon("TGC", "C", 0.55),
  # Aspartic acid (Asp) codons
  codon("GAT", "D", 0.46),
  codon("GAC", "D", 0.54),
  # Glutamic acid (Glu) codons
  codon("GAA", "E", 0.42),
  codon("GAG", "E", 0.58),
  # Phenylalanine (Phe) codons
  codon("TTT", "F", 0.45),
  codon("TTC", "F", 0.55),
  # Glycine (Gly) codons
  codon("GGT", "G", 0.16),
  codon("GGC", "G", 0.34),
  codon("GGA", "G", 0.25),
  codon("GGG", "G", 0.25),
  # Histidine (His) codons
  codon("CAT", "H", 0.41),
  codon("CAC", "H", 0.59),
  # Isoleucine (Ile) codons
  codon("ATT", "I", 0.36),
  codon("ATC", "I", 0.48),
  codon("ATA", "I", 0.16),
  # Lysine (Lys) codons
  codon("AAA", "K", 0.42),
  codon("AAG", "K", 0.58),
  # Leucine (Leu) codons
  codon("TTA", "L", 0.07),
  codon("TTG", "L", 0.13),
  codon("CTT", "L", 0.13),
  codon("CTC", "L", 0.20),
  codon("CTA", "L", 0.07),
  codon("CTG", "L", 0.41),
  # Methionine (Met) codon
  codon("ATG", "M", 1.0),
  # Asparagine (Asn) codons
  codon("AAT", "N", 0.46),
  codon("AAC", "N", 0.54),
  # Proline (Pro) codons
  codon("CCT", "P", 0.28),
  codon("CCC", "P", 0.33),
  codon("CCA", "P", 0.27),
  codon("CCG", "P", 0.11),
  # Glutamine (Gln) codons
  codon("CAA", "Q", 0.25),
  codon("CAG", "Q", 0.75),
  # Arginine (Arg) codons
  codon("CGT", "R", 0.08),
  codon("CGC", "R", 0.19),
  codon("CGA", "R", 0.11),
  codon("CGG", "R", 0.21),
  codon("AGA", "R", 0.20),
  codon("AGG", "R", 0.20),
  # Serine (Ser) codons
  codon("TCT", "S", 0.18),
  codon("TCC", "S", 0.22),
  codon("TCA", "S", 0.15),
  codon("TCG", "S", 0.06),
  codon("AGT", "S", 0.15),
  codon("AGC", "S", 0.24),
  # Threonine (Thr) codons
  codon("ACT", "T", 0.24),
  codon("ACC", "T", 0.36),
  codon("ACA", "T", 0.28),
  codon("ACG", "T", 0.12),
  # Valine (Val) codons
  codon("GTT", "V", 0.18),
  codon("GTC", "V", 0.24),
  codon("GTA", "V", 0.11),
  codon("GTG", "V", 0.47),
  # Tryptophan (Trp) codon
  codon("TGG", "W", 1.0),
  # Tyrosine (Tyr) codons
  codon("TAT", "Y", 0.43),
  codon("TAC", "Y", 0.57)
)

The function I tried below:

# Define function to generate new sequence
random_sequence <- function(protein_sequence) {
  # Initialize empty list to store new sequence
  new_sequence <- list()
  # Iterate over each amino acid in the protein sequence
  for (aa in protein_sequence) {
    # Get all codons that encode the current amino acid
    aa_codons <- codons[sapply(codons, function(x) x$aminoacid == aa)]
    # Get the frequencies of the codons
    aa_freqs <- sapply(aa_codons, function(x) x$frequency)
    # Normalize the frequencies to sum to 1
    aa_freqs <- aa_freqs / sum(aa_freqs)
    # Select a random codon based on the frequencies
    new_codon <- aa_codons[sample(1:length(aa_codons), 1, prob = aa_freqs)]
    # Add the new codon to the new sequence
    new_sequence <- c(new_sequence, new_codon[[1]][1]$sequence)
  }
  # Return the new sequence
  new_sequence <- paste(new_sequence, collapse = "")
  return(new_sequence)
}

# Example usage
protein_sequence <- "MAVVDLKECEILHTWVFPKKKHGEARNDDCQQGKPST"
random_sequence(protein_sequence)

When I run each line of this function individually I get the results I expect. Also, if I set the protein_sequence to be one letter long, it also works. However, when I try to run the function as written, I get the error message:

Error in sum(aa_freqs) : invalid 'type' (list) of argument

I tried to fix this by removing the normalize frequencies to sum to 1 safety check (# aa_freqs <- aa_freqs / sum(aa_freqs). But now I get a different error:

Error in sample.int(length(x), size, replace, prob) : incorrect number of probabilities

I am not sure how to troubleshoot the bug given that when I run through the function line by line with a one letter test amino acid, it all seems to work fine.

Any assistance or insight would be greatly appreciated! Thank you very much for your help!

CodePudding user response：

It looks like your function is trying to parse the whole amino acid sequence at once, rather than letter-by-letter. You can use strsplit or equivalent to split the passed string into component parts and then more or less do as you have indicated:

ReverseTranscribe <- function(protein_sequence) {
  paste(sapply(unlist(strsplit(protein_sequence, "")), function(cur_amino) {
    # Extract possible codons
    potential_matches <- codons[sapply(codons, function(x) {x$aminoacid == cur_amino})]
    # Convert to a weights table
    weights_table <- data.frame(matrix(sapply(potential_matches, unlist), ncol = 3, byrow = TRUE))
    # Weighted sample
    chosen_codon <- sample(weights_table[[1]], 1, prob = weights_table[[3]])
    return(chosen_codon)
  }), collapse = "-")
}

example <- "MAVVDLKECEILHTWVFPKKKHGEARNDDCQQGKPST"

> for (i in 1:3) {
    print(ReverseTranscribe(example))
  }
[1] "ATG-GCT-GTC-GTG-GAT-CTC-AAG-GAA-TGT-GAG-ATC-CTG-CAT-ACC-TGG-GTT-TTC-CCC-AAG-AAA-AAG-CAC-GGT-GAG-GCT-CGG-AAT-GAC-GAT-TGC-CAG-CAG-GGC-AAG-CCT-AGC-ACC"
[1] "ATG-GCC-GTG-GTG-GAT-CTG-AAG-GAG-TGC-GAG-ATC-CTG-CAC-ACC-TGG-GTC-TTT-CCT-AAG-AAG-AAA-CAT-GGC-GAG-GCC-CGT-AAT-GAC-GAC-TGC-CAA-CAA-GGG-AAG-CCA-TCA-ACA"
[1] "ATG-GCC-GTG-GTG-GAT-TTG-AAA-GAG-TGC-GAG-ATC-CTG-CAC-ACC-TGG-GTT-TTC-CCC-AAG-AAG-AAG-CAC-GGA-GAG-GCT-CGC-AAC-GAT-GAT-TGT-CAG-CAG-GGC-AAG-CCC-TCT-ACC"