Merge pull request #54 from Arcadia-Science/ter/aa-with-x

Deal with false mismatches in translation around Ns/Xs in nucleotide and amino acid sequences

scripts/extract_plmutils_nucleotide_sequences.py

@@ -54,7 +54,10 @@ def extract_nucleotide_peptide_sequences(
         )
         protein_peptide_sequence = peptide_record.seq
         if not utils.verify_translation(
-            nucleotide_peptide_sequence, protein_peptide_sequence, to_stop=False
+            nucleotide_peptide_sequence,
+            protein_peptide_sequence,
+            to_stop=False,
+            allow_wildcard_x=True,
         ):
             raise ValueError(
                 f"Warning: Translation mismatch for {transcript_id}. "

scripts/utils.py

@@ -2,10 +2,37 @@
 from Bio.Seq import Seq
 
 
-def verify_translation(nucleotide_seq, amino_acid_seq, to_stop):
-    """Verify that a nucleotide sequence translates correctly to its amino acid sequence."""
-    translated_seq = Seq(str(nucleotide_seq)).translate(to_stop=to_stop)
-    return str(translated_seq) == str(amino_acid_seq)
+def verify_translation(nucleotide_seq, amino_acid_seq, to_stop, allow_wildcard_x=False):
+    """
+    Verify that a nucleotide sequence correctly translates to its corresponding amino acid sequence.
+
+    The `allow_wildcard_x` argument permits any amino acid encoded by an "X" in the input sequence
+    to be replaced with the corresponding translated nucleotide sequence. This option is needed when
+    `amino_acid_seq` was translated by a tool like Orfipy that does not translate codons that
+    contain an "N". This is because Biopython (which is used here to translate `nucleotide_seq` into
+    an amino acid sequence) is less restrictive and will translate codons containing "N"s if they
+    can only translate into a single amino acid, regardless of the "N". For example, the leucine
+    codons all begins with CT. Therefore, even if the codon is CTN, Biopython translates it to
+    leucine.
+
+    Notes:
+    - This function requires the sequences to be the same length.
+    - This function only verifies that two translations match; it does not edit the amino acid
+      sequence.
+    """
+    translated_seq = str(Seq(str(nucleotide_seq)).translate(to_stop=to_stop))
+    amino_acid_seq = str(amino_acid_seq)
+
+    if len(translated_seq) != len(amino_acid_seq):
+        return False
+
+    if allow_wildcard_x:
+        amino_acid_seq = "".join(
+            translated_seq[i] if amino_acid == "X" else amino_acid
+            for i, amino_acid in enumerate(amino_acid_seq)
+        )
+
+    return translated_seq == amino_acid_seq
 
 
 def read_fasta(fasta_file):