FIX: Full resname matching for nucleic acid analysis (#5365)

Author: BradyAJohnston

package/CHANGELOG

@@ -22,6 +22,13 @@ The rules for this file:
  * 2.11.0
 
 Fixes
+ * `MDAnalysis.analysis.nucleicacids.WatsonCrickDist`, `MinorPairDist`,
+   and `MajorPairDist` now match residue names against the full resname
+   instead of only the first character, fixing incorrect behaviour with
+   multi-character names such as CHARMM ``DG``/``DA`` or ``GUA``/``ADE``
+   (Issue #5360). Note: users relying on the previous truncation
+   behaviour must now pass explicit residue names via the ``g_name``,
+   ``a_name``, etc. keyword arguments.
  * `MDAnalysis.analysis.atomicdistances.AtomicDistances` results are now
    consistent with expected `analysis` documentation data type = Results
    (Issue #4819, PR #5347)
@@ -57,10 +64,10 @@ Enhancements
    calculations and supports ``backend`` selection (e.g. ``"distopia"``)
    (PR #5182)
  * Added ASV benchmark for `MDAnalysis.analysis.contacts` (PR #5291)
- * Improved performance of inverse index mapping in AtomGroup using an optimized 
+ * Improved performance of inverse index mapping in AtomGroup using an optimized
    Cython implementation in lib._cutils.inverse_int_index()
    (Issue #3387, PR #5252)
- * Added documentation for all keyword in select_atoms() and 
+ * Added documentation for all keyword in select_atoms() and
    selections.rst (Issue #5317, PR #5325)
  * Added HydrogenBondAnalysis benchmark for performance tracking (PR #5309)
  * Added `select=None` in `analysis.rms.RMSD` to perform no selection on

package/MDAnalysis/analysis/nucleicacids.py

@@ -255,6 +255,12 @@ def select_strand_atoms(
 
 
         .. versionadded:: 2.7.0
+
+        .. versionchanged:: 2.11.0
+           Residue names are now matched exactly against the full resname
+           rather than only its first character. Users with multi-character
+           residue names (e.g. CHARMM ``"GUA"``, ``"DG"``) must pass the
+           appropriate names via the keyword arguments.
         """
         pyrimidines: List[str] = [c_name, t_name, u_name]
         purines: List[str] = [a_name, g_name]
@@ -263,9 +269,9 @@ def select_strand_atoms(
         sel2: List[mda.AtomGroup] = []
 
         for pair in zip(strand1.residues, strand2.residues):
-            if pair[0].resname[0] in pyrimidines:
+            if pair[0].resname in pyrimidines:
                 a1, a2 = a2_name, a1_name
-            elif pair[0].resname[0] in purines:
+            elif pair[0].resname in purines:
                 a1, a2 = a1_name, a2_name
             else:
                 raise ValueError(
@@ -431,6 +437,12 @@ class WatsonCrickDist(NucPairDist):
         :class:`~MDAnalysis.core.groups.ResidueGroup` as input.
         The previous input type, ``List[Residue]`` is still supported,
         but it is **deprecated** and will be removed in release 3.0.0.
+
+    .. versionchanged:: 2.11.0
+       Residue names are now matched exactly against the full resname
+       rather than only its first character. Users with multi-character
+       residue names (e.g. CHARMM ``"GUA"``, ``"DG"``) must pass the
+       appropriate names via the keyword arguments.
     """
 
     def __init__(
@@ -551,6 +563,12 @@ class MinorPairDist(NucPairDist):
 
 
     .. versionadded:: 2.7.0
+
+    .. versionchanged:: 2.11.0
+       Residue names are now matched exactly against the full resname
+       rather than only its first character. Users with multi-character
+       residue names (e.g. CHARMM ``"GUA"``, ``"DG"``) must pass the
+       appropriate names via the keyword arguments.
     """
 
     def __init__(
@@ -649,6 +667,12 @@ class MajorPairDist(NucPairDist):
 
 
     .. versionadded:: 2.7.0
+
+    .. versionchanged:: 2.11.0
+       Residue names are now matched exactly against the full resname
+       rather than only its first character. Users with multi-character
+       residue names (e.g. CHARMM ``"GUA"``, ``"DG"``) must pass the
+       appropriate names via the keyword arguments.
     """
 
     def __init__(

testsuite/MDAnalysisTests/analysis/test_nucleicacids.py

@@ -22,6 +22,7 @@
 #
 
 import MDAnalysis as mda
+import numpy as np
 
 import pytest
 
@@ -55,15 +56,24 @@ def test_empty_ag_error(strand):
     strand2 = ResidueGroup([strand.residues[1]])
 
     with pytest.raises(ValueError, match="returns an empty AtomGroup"):
-        NucPairDist.select_strand_atoms(strand1, strand2, "UNK1", "O2")
+        NucPairDist.select_strand_atoms(
+            strand1, strand2, "UNK1", "O2", g_name="GUA"
+        )
 
 
 @pytest.fixture(scope="module")
 def wc_rna(strand, client_NucPairDist):
     strand1 = ResidueGroup([strand.residues[0], strand.residues[21]])
     strand2 = ResidueGroup([strand.residues[1], strand.residues[22]])
 
-    WC = WatsonCrickDist(strand1, strand2)
+    WC = WatsonCrickDist(
+        strand1,
+        strand2,
+        g_name="GUA",
+        a_name="ADE",
+        c_name="CYT",
+        u_name="URA",
+    )
     WC.run(**client_NucPairDist)
     return WC
 
@@ -83,10 +93,18 @@ def test_wc_dist(wc_rna):
 
 def test_wc_dist_invalid_residue_types(u):
     strand = u.select_atoms("resid 1-10")
+    # residues[0]=GUA, residues[21]=POT (non-nucleic), residues[22]=POT
     strand1 = ResidueGroup([strand.residues[0], strand.residues[21]])
     strand2 = ResidueGroup([strand.residues[2], strand.residues[22]])
     with pytest.raises(ValueError, match="is not a valid nucleic acid"):
-        WatsonCrickDist(strand1, strand2)
+        WatsonCrickDist(
+            strand1,
+            strand2,
+            g_name="GUA",
+            a_name="ADE",
+            c_name="CYT",
+            u_name="URA",
+        )
 
 
 def test_selection_length_mismatch(strand):
@@ -101,7 +119,14 @@ def test_wc_dist_deprecation_warning(strand):
     strand2 = [strand.residues[1], strand.residues[22]]
 
     with pytest.deprecated_call():
-        WatsonCrickDist(strand1, strand2)
+        WatsonCrickDist(
+            strand1,
+            strand2,
+            g_name="GUA",
+            a_name="ADE",
+            c_name="CYT",
+            u_name="URA",
+        )
 
 
 def test_wc_dist_strand_verification(strand):
@@ -122,7 +147,14 @@ def test_minor_dist(strand, client_NucPairDist):
     strand1 = ResidueGroup([strand.residues[2], strand.residues[19]])
     strand2 = ResidueGroup([strand.residues[16], strand.residues[4]])
 
-    MI = MinorPairDist(strand1, strand2)
+    MI = MinorPairDist(
+        strand1,
+        strand2,
+        g_name="GUA",
+        a_name="ADE",
+        c_name="CYT",
+        u_name="URA",
+    )
     MI.run(**client_NucPairDist)
 
     assert MI.results.distances[0, 0] == approx(15.06506, rel=1e-3)
@@ -133,8 +165,85 @@ def test_major_dist(strand, client_NucPairDist):
     strand1 = ResidueGroup([strand.residues[1], strand.residues[4]])
     strand2 = ResidueGroup([strand.residues[11], strand.residues[8]])
 
-    MA = MajorPairDist(strand1, strand2)
+    MA = MajorPairDist(
+        strand1,
+        strand2,
+        g_name="GUA",
+        a_name="ADE",
+        c_name="CYT",
+        u_name="URA",
+    )
     MA.run(**client_NucPairDist)
 
     assert MA.results.distances[0, 0] == approx(26.884272, rel=1e-3)
     assert MA.results.distances[0, 1] == approx(13.578535, rel=1e-3)
+
+
+@pytest.fixture(scope="module")
+def dna_u():
+    """Synthetic universe with 2-letter DNA resnames (DG, DC, DA, DT).
+
+    Used to test that WatsonCrickDist matches full resnames rather than only
+    the first character.
+
+    Strand layout (one frame, all atoms at origin):
+      residue 0 – DG (purine):    atom N1
+      residue 1 – DC (pyrimidine): atom N3
+      residue 2 – DA (purine):    atom N1
+      residue 3 – DT (pyrimidine): atom N3
+    """
+    n_atoms = 4
+    n_residues = 4
+    u = mda.Universe.empty(
+        n_atoms,
+        n_residues=n_residues,
+        atom_resindex=[0, 1, 2, 3],
+        trajectory=True,
+    )
+    u.add_TopologyAttr("name", ["N1", "N3", "N1", "N3"])
+    u.add_TopologyAttr("resname", ["DG", "DC", "DA", "DT"])
+    u.add_TopologyAttr("resid", [1, 2, 3, 4])
+    coords = np.zeros((1, n_atoms, 3), dtype=np.float32)
+    u.load_new(coords)
+    return u
+
+
+def test_wc_dist_multichar_resnames(dna_u):
+    """WatsonCrickDist must match the full resname, not just resname[0]\.
+
+    With 2-letter names like DG/DA, the first character 'D' does not appear in
+    the purine/pyrimidine lists, so the buggy resname[0] check previously
+    rasied ValueError for every residue.
+    """
+    strand1 = ResidueGroup([dna_u.residues[0], dna_u.residues[2]])  # DG, DA
+    strand2 = ResidueGroup([dna_u.residues[1], dna_u.residues[3]])  # DC, DT
+
+    WC = WatsonCrickDist(
+        strand1,
+        strand2,
+        g_name="DG",
+        a_name="DA",
+        c_name="DC",
+        t_name="DT",
+    )
+    WC.run()
+    assert WC.results.distances.shape == (1, 2)
+
+
+def test_select_strand_atoms_multichar_resnames(dna_u):
+    """select_strand_atoms must recognise multi-character resnames."""
+    strand1 = ResidueGroup([dna_u.residues[0], dna_u.residues[2]])  # DG, DA
+    strand2 = ResidueGroup([dna_u.residues[1], dna_u.residues[3]])  # DC, DT
+
+    sel1, sel2 = NucPairDist.select_strand_atoms(
+        strand1,
+        strand2,
+        a1_name="N1",
+        a2_name="N3",
+        g_name="DG",
+        a_name="DA",
+        c_name="DC",
+        t_name="DT",
+    )
+    assert len(sel1) == 2
+    assert len(sel2) == 2