diff --git a/package/CHANGELOG b/package/CHANGELOG index bca8cc647ad..52ebe8e4696 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -15,7 +15,7 @@ The rules for this file: ------------------------------------------------------------------------------- ??/??/?? IAlibay, ianmkenney, PicoCentauri, pgbarletta, p-j-smith, - richardjgowers, lilyminium + richardjgowers, lilyminium, ALescoulie * 2.7.0 @@ -29,6 +29,8 @@ Fixes * Fixes hydrogenbonds tutorial path to point to hbonds (Issue #4285, PR #4286) Enhancements + * Add faster nucleic acid Major and Minor pair distance calculators using + AnalysisBase for updated nucleicacids module (Issue #3720, PR #3735) * Adds external sidebar links (Issue #4296) * Updated lib.qcprot.CalcRMSDRotationalMatrix to accept either float32 or float64 inputs (PR #4273, part of #3927) @@ -48,10 +50,18 @@ Changes replacing the now deprecated `xdrlib` core Python library (PR #4271) * ConverterBase class moved from coordinates/base.py to converters/base.py (Issue #3404) + * Results for WatsonCrickDist nucleic acids analysis are now stored in + `analysis.nucleicacids.WatsonCrickDist.results.distances` (Issue #3720, PR #3735) Deprecations * coordinates.base.ConverterBase has been deprecated and will be removed in 3.0.0; use converters.base.ConvertBase instead (Issue #3404) + * In `nucleicacids.WatsonCrickDist`, accepting lists of `Residue` objects was deprecated + in favor of using `ResidueGroup`: using `List[Residue]` will be removed in release + 3.0.0; instead use a `ResidueGroup` (Issue #3720, PR #3735) + * In `nucleicacids.WatsonCrickDist` the result `results.pair_distances` was + deprecated and will be removed in 3.0.0; use `results.distances` (Issue #3720, + PR #3735) 28/08/23 IAlibay, hmacdope, pillose, jaclark5, tylerjereddy diff --git a/package/MDAnalysis/analysis/nucleicacids.py b/package/MDAnalysis/analysis/nucleicacids.py index 2137e2f1e49..0eccd039ba4 100644 --- a/package/MDAnalysis/analysis/nucleicacids.py +++ b/package/MDAnalysis/analysis/nucleicacids.py @@ -26,8 +26,8 @@ ========================================================================= :Author: Alia Lescoulie -:Year: 2022 -:copyright: GNU Public Licence v3 +:Year: 2022-2023 +:copyright: LGPLv2.1 The module provides classes for analyzing nucleic acids structures. This is an updated, higher performance version of previous nucleic acid tools. @@ -46,13 +46,24 @@ .. autoclass:: WatsonCrickDist :members: + :exclude-members: select_strand_atoms + :inherited-members: + +.. autoclass:: MinorPairDist + :members: + :exclude-members: select_strand_atoms + :inherited-members: + +.. autoclass:: MajorPairDist + :members: + :exclude-members: select_strand_atoms :inherited-members: .. versionadded 2.2.0 """ -from typing import List, Dict +from typing import List, Dict, Tuple, Union import warnings import numpy as np @@ -60,11 +71,20 @@ import MDAnalysis as mda from .distances import calc_bonds from .base import AnalysisBase, Results -from MDAnalysis.core.groups import Residue +from MDAnalysis.core.groups import Residue, ResidueGroup + + +# Deprecation: In 3.0.0 change type to just +# ResidueClass = ResidueGroup +ResidueClass = Union[List[Residue], ResidueGroup] +r"""A type alias for :code:`Union[List[Residue], ResidueGroup]` + +Used as an alias for methods where either class is acceptable. +""" class NucPairDist(AnalysisBase): - r"""Atom Pair distance calculation base class. + r"""Atom pair distance calculation base class. Takes two lists of :class:`~MDAnalysis.core.groups.AtomGroup` and computes the distances between them over a trajectory. Used as a @@ -83,43 +103,77 @@ class NucPairDist(AnalysisBase): kwargs: dict Arguments for :class:`~MDAnalysis.analysis.base.AnalysisBase` + Attributes ---------- - times: numpy.ndarray - Simulation times for analysis. results.pair_distances: numpy.ndarray - 2D array of pair distances. First dimension is simulation time, second - dimension contains the pair distances for each each entry pair in - selection1 and selection2. + 2D array of pair distances. First dimension is simulation time, + second dimension contains the pair distances for each each entry + pair in selection1 and selection2. .. versionadded:: 2.4.0 + .. note:: + `results.pair_distances` is slated for deprecation in + version 3.0.0, use `results.distances` instead. + .. deprecated:: 2.7.0 + `results.pair_distances` will be removed in + version 3.0.0, use :attr:`results.distances` instead. + + results.distances: numpy.ndarray + stored in a 2d numpy array with first index selecting the + Residue pair, and the second index selecting the frame number + Distances are stored in a 2d numpy array with axis 0 (first index) + indexing the trajectory frame and axis 1 (second index) selecting the + Residue pair. + + .. versionadded:: 2.7.0 + + times: numpy.ndarray + Simulation times for analysis. + Raises ------ - ValueError If the selections given are not the same length + ValueError + An :class:`~MDAnalysis.core.groups.AtomGroup` in one of the + strands not a valid nucleic acid + ValueError + If a given residue pair from the provided strands returns an empty + :class:`~MDAnalysis.core.groups.AtomGroup` when selecting the atom + pairs used in the distance calculations + + *Version Info* .. versionchanged:: 2.5.0 - The ability to access by passing selection indices to :attr:`results` is + The ability to access by passing selection indices to :attr:`results` is now removed as of MDAnalysis version 2.5.0. Please use :attr:`results.pair_distances` instead. The :attr:`results.times` was deprecated and is now removed as of - MDAnalysis 2.5.0. Please use the class attribute :attr:`times` instead. + MDAnalysis 2.5.0. + Please use the class attribute :attr:`times` instead. + + .. versionchanged:: 2.7.0 + Added static method :attr:`select_strand_atoms` as a + helper for selecting atom pairs for distance analysis. """ _s1: mda.AtomGroup _s2: mda.AtomGroup _n_sel: int _res_dict: Dict[int, List[float]] - _times = List[float] def __init__(self, selection1: List[mda.AtomGroup], selection2: List[mda.AtomGroup], **kwargs) -> None: - super(NucPairDist, self).__init__(selection1[0].universe.trajectory, **kwargs) + super( + NucPairDist, + self).__init__( + selection1[0].universe.trajectory, + **kwargs) if len(selection1) != len(selection2): raise ValueError("Selections must be same length") @@ -133,93 +187,432 @@ def __init__(self, selection1: List[mda.AtomGroup], self._s1 += selection1[i] self._s2 += selection2[i] + @staticmethod + def select_strand_atoms( + strand1: ResidueGroup, strand2: ResidueGroup, + a1_name: str, a2_name: str, g_name: str = 'G', + a_name: str = 'A', u_name: str = 'U', + t_name: str = 'T', c_name: str = 'C' + ) -> Tuple[List[mda.AtomGroup], List[mda.AtomGroup]]: + r""" + A helper method for nucleic acid pair distance analyses. + Used for selecting specific atoms from two strands of nucleic acids. + + + Parameters + ---------- + strand1: List[Residue] + The first nucleic acid strand + strand2: List[Residue] + The second nucleic acid strand + a1_name: str + The selection for the purine base of the strand pair + a2_name: str + the selection for the pyrimidine base of the strand pair + g_name: str (optional) + Name of Guanine in topology, by default assigned to G + a_name: str (optional) + Name of Adenine in topology, by default assigned to A + u_name: str (optional) + Name of Uracil in topology, by default assigned to U + t_name: str (optional) + Name of Thymine in topology, by default assigned to T + c_name: str (optional) + Name of Cytosine in topology, by default assigned to C + + Returns + ------- + Tuple[List[AtomGroup], List[AtomGroup]] + returns a tuple containing two lists of + :class:`~MDAnalysis.core.groups.AtomGroup`\s + corresponding to the provided selections from each strand. + + Raises + ------ + ValueError: + An :class:`~MDAnalysis.core.groups.AtomGroup` + in one of the strands not a valid nucleic acid + ValueError: + An :class:`~MDAnalysis.core.groups.Residue` returns an empty + :class:`~MDAnalysis.core.groups.AtomGroup` + with the provided selection + + + .. versionadded:: 2.7.0 + """ + pyrimidines: List[str] = [c_name, t_name, u_name] + purines: List[str] = [a_name, g_name] + + sel1: List[mda.AtomGroup] = [] + sel2: List[mda.AtomGroup] = [] + + for pair in zip(strand1.residues, strand2.residues): + if pair[0].resname[0] in pyrimidines: + a1, a2 = a2_name, a1_name + elif pair[0].resname[0] in purines: + a1, a2 = a1_name, a2_name + else: + raise ValueError( + f"AtomGroup in {pair} is not a valid nucleic acid" + ) + + ag1 = pair[0].atoms.select_atoms(f'name {a1}') + ag2 = pair[1].atoms.select_atoms(f'name {a2}') + + if not all(len(ag) > 0 for ag in [ag1, ag2]): + err_info: Tuple[Residue, str] = (pair[0], a1) \ + if len(ag1) == 0 else (pair[1], a2) + + raise ValueError( + ( + f"{err_info[0]} returns an empty AtomGroup" + "with selection string \"name {a2}\"" + ) + ) + + sel1.append(ag1) + sel2.append(ag2) + + return (sel1, sel2) + def _prepare(self) -> None: - self._res_array: np.ndarray = np.zeros([self.n_frames, self._n_sel]) + self._res_array: np.ndarray = np.zeros( + [self.n_frames, self._n_sel] + ) def _single_frame(self) -> None: - dist: np.ndarray = calc_bonds(self._s1.positions, self._s2.positions) + dist: np.ndarray = calc_bonds( + self._s1.positions, self._s2.positions + ) + self._res_array[self._frame_index, :] = dist def _conclude(self) -> None: - self.results['pair_distances'] = self._res_array + self.results['distances'] = self._res_array + self.results['pair_distances'] = self.results['distances'] + # TODO: remove pair_distances in 3.0.0 class WatsonCrickDist(NucPairDist): - r"""Watson-Crick basepair distance for selected residues over a trajectory. - - Takes two lists of :class:`~MDAnalysis.core.groups.Residue` objects and calculates - the Watson-Crick distance between them over the trajectory. Bases are matched by - their index in the lists given as arguments. + r""" + Watson-Crick base pair distance for selected + residues over a trajectory. + + Takes two :class:`~MDAnalysis.core.groups.ResidueGroup` + objects or two lists of :class:`~MDAnalysis.core.groups.Residue` + and calculates the distance between the nitrogen atoms in the + Watson-Crick hydrogen bond over the trajectory. Bases are matched + either by their index in the two + :class:`~MDAnalysis.core.groups.ResidueGroup` provided as arguments, + or based on the indices of the provided lists of + :class:`~MDAnalysis.core.groups.Residue` objects depending + on which is provided. + + .. note:: + Support for :class:`~MDAnalysis.core.groups.Residue` is slated for + deprecation and will raise a warning when used. It still works but + :class:`~MDAnalysis.core.groups.ResidueGroup` is preferred. Parameters ---------- - strand1: List[Residue] + strand1: ResidueClass First list of bases - strand2: List[Residue] + + .. deprecated:: 2.7.0 + Using a list of :class:`~MDAnalysis.core.groups.Residue` will + be removed in 3.0.0. Use a + :class:`~MDAnalysis.core.groups.ResidueGroup`. + + strand2: ResidueClass Second list of bases + + .. deprecated:: 2.7.0 + Using a list of :class:`~MDAnalysis.core.groups.Residue` will + be removed in 3.0.0. Use a + :class:`~MDAnalysis.core.groups.ResidueGroup`. + n1_name: str (optional) - Name of Nitrogen 1 of nucleic acids, by default assigned to N1 + Name of Nitrogen 1 of nucleic acids, by default assigned to "N1" n3_name: str (optional) - Name of Nitrogen 3 of nucleic acids, by default assigned to N3 + Name of Nitrogen 3 of nucleic acids, by default assigned to "N3" g_name: str (optional) - Name of Guanine in topology, by default assigned to G + Name of Guanine in topology, by default assigned to "G" a_name: str (optional) - Name of Adenine in topology, by default assigned to A + Name of Adenine in topology, by default assigned to "A" u_name: str (optional) - Name of Uracil in topology, by default assigned to U + Name of Uracil in topology, by default assigned to "U" t_name: str (optional) - Name of Thymine in topology, by default assigned to T + Name of Thymine in topology, by default assigned to "T" c_name: str (optional) Name of Cytosine in topology, by default assigned to C **kwargs: dict - arguments for :class:`~MDAnalysis.analysis.base.AnalysisBase` + Key word arguments for + :class:`~MDAnalysis.analysis.base.AnalysisBase` Attributes ---------- - times: numpy.ndarray - Simulation times for analysis. + results.distances: numpy.ndarray + Distances are stored in a 2d numpy array with axis 0 (first index) + indexing the trajectory frame and axis 1 (second index) selecting the + Residue pair. + + .. versionadded:: 2.7.0 + results.pair_distances: numpy.ndarray - 2D array of Watson-Crick basepair distances. First dimension is - simulation time, second dimension contains the pair distances for - each each entry pair in strand1 and strand2. + 2D array of pair distances. First dimension is + simulation time, second dimension contains the + pair distances for each each entry pair in + selection1 and selection2. .. versionadded:: 2.4.0 + .. deprecated:: 2.7.0 + `results.pair_distances` will be removed in version 3.0.0, + use :attr:`results.distances` instead. + + times: numpy.ndarray + Simulation times for analysis. Raises ------ + TypeError + If the provided list of :class:`~MDAnalysis.core.Residue` contains + non-Residue elements + + .. deprecated:: 2.7.0 + Starting with version 3.0.0, this exception will no longer + be raised because only + :class:`~MDAnalysis.core.groups.ResidueGroup` will be allowed. + ValueError - If the residues given are not amino acids + If `strand1` and `strand2` are not the same length + ValueError: + An :class:`~MDAnalysis.core.groups.AtomGroup` + in one of the strands not a valid nucleic acid ValueError - If the selections given are not the same length + If a given residue pair from the provided strands returns an empty + :class:`~MDAnalysis.core.groups.AtomGroup` when selecting the atom + pairs used in the distance calculations + *Version Info* + .. versionchanged:: 2.5.0 - Accessing results by passing strand indices to :attr:`results` is - was deprecated and is now removed as of MDAnalysis version 2.5.0. Please - use :attr:`results.pair_distances` instead. + Accessing results by passing strand indices to :attr:`results` + was deprecated and is now removed as of MDAnalysis version 2.5.0. + Please use :attr:`results.pair_distances` instead. The :attr:`results.times` was deprecated and is now removed as of - MDAnalysis 2.5.0. Please use the class attribute :attr:`times` instead. + MDAnalysis 2.5.0. Please use the class attribute + :attr:`times` instead. + + .. versionchanged:: 2.7.0 + `strand1` and `strand2` now also accept a + :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. """ - def __init__(self, strand1: List[Residue], strand2: List[Residue], + def __init__(self, strand1: ResidueClass, strand2: ResidueClass, n1_name: str = 'N1', n3_name: str = "N3", g_name: str = 'G', a_name: str = 'A', u_name: str = 'U', t_name: str = 'T', c_name: str = 'C', **kwargs) -> None: - sel1: List[mda.AtomGroup] = [] - sel2: List[mda.AtomGroup] = [] - strand = zip(strand1, strand2) - for s in strand: - if s[0].resname[0] in [c_name, t_name, u_name]: - a1, a2 = n3_name, n1_name - elif s[0].resname[0] in [a_name, g_name]: - a1, a2 = n1_name, n3_name - else: - raise ValueError(f"{s} are not valid nucleic acids") + def verify_strand(strand: ResidueClass) -> ResidueGroup: + # Helper method to verify the strands + + if isinstance(strand, list): # Checking if a list is given + # verify list is only Residues + if not all(isinstance(resid, Residue) for resid in strand): + raise TypeError(f"{strand} contains non-Residue elements") + + warnings.warn( + DeprecationWarning( + ( + f"ResidueGroup should be used for {strand} instead" + "of giving a Residue list" + ) + ) + ) + # Convert to a ResidueGroup + strand: ResidueGroup = ResidueGroup(strand) + + return strand + + strand1: ResidueGroup = verify_strand(strand1) + strand2: ResidueGroup = verify_strand(strand2) + + strand_atomgroups: Tuple[List[mda.AtomGroup], List[mda.AtomGroup]] = \ + self.select_strand_atoms( + strand1, strand2, n1_name, n3_name, + g_name=g_name, a_name=a_name, + t_name=t_name, u_name=u_name, c_name=c_name + ) + + super(WatsonCrickDist, self).__init__( + strand_atomgroups[0], strand_atomgroups[1], **kwargs + ) + + +class MinorPairDist(NucPairDist): + r"""Minor-Pair basepair distance for selected residues over a trajectory. + + Takes two :class:`~MDAnalysis.core.groups.ResidueGroup` objects and + calculates the Minor-groove hydrogen bond length between the + nitrogen and oxygen atoms over the trajectory. Bases are + matched by their index in the two + :class:`~MDAnalysis.core.groups.ResidueGroup` provided as arguments. + + Parameters + ---------- + strand1: List[Residue] + First list of bases + strand2: List[Residue] + Second list of bases + o2_name: str (optional) + Name of Oxygen 2 of nucleic acids; + by default assigned to "O2"; + c2_name: str (optional) + Name of Carbon 2 of nucleic acids; + by default assigned to "C2"; + g_name: str (optional) + Name of Guanine in topology; + by default assigned to "G"; + a_name: str (optional) + Name of Adenine in topology + by default assigned to "A"; + u_name: str (optional) + Name of Uracil in topology; + by default assigned to "U"; + t_name: str (optional) + Name of Thymine in topology; + by default assigned to "T"; + c_name: str (optional) + Name of Cytosine in topology; + by default assigned to "C"; + **kwargs: + keyword arguments for + :class:`~MDAnalysis.analysis.base.AnalysisBase` + + Attributes + ---------- + results.distances: numpy.ndarray + stored in a 2d numpy array with first index selecting + the Residue pair, and the second index selecting the frame number + times: numpy.ndarray + Simulation times for analysis. + + Raises + ------ + ValueError + If the selections given are not the same length + A :class:`~MDAnalysis.core.Residue` in + one of the strands not a valid nucleic acid + ValueError + If a given residue pair from the provided strands returns an empty + :class:`~MDAnalysis.core.groups.AtomGroup` when selecting the atom + pairs used in the distance calculations + + + .. versionadded:: 2.7.0 + """ + + def __init__(self, strand1: ResidueGroup, strand2: ResidueGroup, + o2_name: str = 'O2', c2_name: str = "C2", + g_name: str = 'G', a_name: str = 'A', u_name: str = 'U', + t_name: str = 'T', c_name: str = 'C', + **kwargs) -> None: + + selections: Tuple[List[mda.AtomGroup], List[mda.AtomGroup]] = \ + self.select_strand_atoms( + strand1, strand2, c2_name, o2_name, + g_name=g_name, a_name=a_name, + t_name=t_name, u_name=u_name, c_name=c_name + ) + + super(MinorPairDist, self).__init__( + selections[0], selections[1], **kwargs + ) + + +class MajorPairDist(NucPairDist): + r"""Minor-Pair base pair distance for + selected residues over a trajectory. + + Takes two :class:`~MDAnalysis.core.groups.ResidueGroup` objects and + calculates the Major-groove hydrogen bond length between the nitrogen + and oxygen atoms over the trajectory. Bases are matched by their index + in the two :class:`~MDAnalysis.core.groups.ResidueGroup` + provided as arguments. + + Parameters + ---------- + strand1: List[Residue] + First list of bases + strand2: List[Residue] + Second list of bases + o6_name: str (optional) + Name of Oxygen 6 of nucleic acids; + by default assigned to "O6" + n4_name: str (optional) + Name of Nitrogen 4 of nucleic acids; + by default assigned to "N4" + g_name: str (optional) + Name of Guanine in topology; + by default assigned to "G" + a_name: str (optional) + Name of Adenine in topology; + by default assigned to "A" + u_name: str (optional) + Name of Uracil in topology; + by default assigned to "U" + t_name: str (optional) + Name of Thymine in topology; + by default assigned to "T" + c_name: str (optional) + Name of Cytosine in topology; + by default assigned to "C" + **kwargs: + arguments for :class:`~MDAnalysis.analysis.base.AnalysisBase` + + Attributes + ---------- + results.distances: numpy.ndarray + Distances are stored in a 2d numpy array with axis 0 (first index) + indexing the trajectory frame and axis 1 (second index) selecting the + Residue pair. + times: numpy.ndarray + Simulation times for analysis. + + Raises + ------ + ValueError + A :class:`~MDAnalysis.core.Residue` + in one of the strands not a valid nucleic acid + ValueError + If a given residue pair from the provided strands returns an empty + :class:`~MDAnalysis.core.groups.AtomGroup` when selecting the atom + pairs used in the distance calculations + ValueError + if the selections given are not the same length + + + .. versionadded:: 2.7.0 + """ + + def __init__(self, strand1: ResidueGroup, strand2: ResidueGroup, + n4_name: str = 'N4', o6_name: str = "O6", + g_name: str = 'G', a_name: str = 'A', u_name: str = 'U', + t_name: str = 'T', c_name: str = 'C', + **kwargs) -> None: - sel1.append(s[0].atoms.select_atoms(f'name {a1}')) - sel2.append(s[1].atoms.select_atoms(f'name {a2}')) + selections: Tuple[List[mda.AtomGroup], List[mda.AtomGroup]] = \ + self.select_strand_atoms( + strand1, strand2, o6_name, n4_name, g_name=g_name, + a_name=a_name, t_name=t_name, u_name=u_name, + c_name=c_name + ) - super(WatsonCrickDist, self).__init__(sel1, sel2, **kwargs) + super(MajorPairDist, self).__init__( + selections[0], selections[1], **kwargs + ) diff --git a/testsuite/MDAnalysisTests/analysis/test_nucleicacids.py b/testsuite/MDAnalysisTests/analysis/test_nucleicacids.py index eba429d1bb9..fb7d39374cd 100644 --- a/testsuite/MDAnalysisTests/analysis/test_nucleicacids.py +++ b/testsuite/MDAnalysisTests/analysis/test_nucleicacids.py @@ -22,10 +22,17 @@ # import MDAnalysis as mda + import pytest -from MDAnalysis.analysis.nucleicacids import WatsonCrickDist, NucPairDist + +from pytest import approx + +from MDAnalysis.analysis.nucleicacids import (NucPairDist, WatsonCrickDist, + MajorPairDist, MinorPairDist) + from MDAnalysisTests.datafiles import RNA_PSF, RNA_PDB -from numpy.testing import assert_allclose + +from MDAnalysis.core.groups import ResidueGroup @pytest.fixture(scope='module') @@ -33,11 +40,24 @@ def u(): return mda.Universe(RNA_PSF, RNA_PDB) +@pytest.fixture(scope="module") +def strand(unv=u): + unv = mda.Universe(RNA_PSF, RNA_PDB) + return unv.select_atoms("segid RNAA") + + +def test_empty_ag_error(strand): + strand1 = ResidueGroup([strand.residues[0]]) + strand2 = ResidueGroup([strand.residues[1]]) + + with pytest.raises(ValueError, match="returns an empty AtomGroup"): + NucPairDist.select_strand_atoms(strand1, strand2, 'UNK1', 'O2') + + @pytest.fixture(scope='module') -def wc_rna(u): - strand: mda.AtomGroup = u.select_atoms("segid RNAA") - strand1 = [strand.residues[0], strand.residues[21]] - strand2 = [strand.residues[1], strand.residues[22]] +def wc_rna(strand): + strand1 = ResidueGroup([strand.residues[0], strand.residues[21]]) + strand2 = ResidueGroup([strand.residues[1], strand.residues[22]]) WC = WatsonCrickDist(strand1, strand2) WC.run() @@ -45,34 +65,72 @@ def wc_rna(u): def test_wc_dist_shape(wc_rna): - assert wc_rna.results.pair_distances.shape == (1, 2) + assert wc_rna.results.distances.shape == (1, 2) def test_wc_dist_results_keys(wc_rna): - assert "pair_distances" in wc_rna.results + assert "distances" in wc_rna.results def test_wc_dist(wc_rna): - assert_allclose(wc_rna.results.pair_distances[0, 0], 4.3874702, atol=1e-3) - assert_allclose(wc_rna.results.pair_distances[0, 1], 4.1716404, atol=1e-3) + assert wc_rna.results.distances[0, 0] == approx(4.3874702, rel=1e-3) + assert wc_rna.results.distances[0, 1] == approx(4.1716404, rel=1e-3) def test_wc_dist_invalid_residue_types(u): strand = u.select_atoms("resid 1-10") - strand1 = [strand.residues[0], strand.residues[21]] - strand2 = [strand.residues[2], strand.residues[22]] - with pytest.raises(ValueError, match="are not valid nucleic acids"): + 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) -def test_selection_length_mismatch(u): - sel1 = u.select_atoms("resid 1-10") - sel2 = u.select_atoms("resid 1-5") +def test_selection_length_mismatch(strand): + sel1 = strand.residues[1:10] + sel2 = strand.residues[1:9] with pytest.raises(ValueError, match="Selections must be same length"): NucPairDist(sel1, sel2) +def test_wc_dist_deprecation_warning(strand): + strand1 = [strand.residues[0], strand.residues[21]] + strand2 = [strand.residues[1], strand.residues[22]] + + with pytest.deprecated_call(): + WatsonCrickDist(strand1, strand2) + + +def test_wc_dist_strand_verification(strand): + strand1 = [strand.residues[0], strand[0]] + strand2 = [strand.residues[1], strand.residues[22]] + + with pytest.raises(TypeError, match="contains non-Residue elements"): + WatsonCrickDist(strand1, strand2) + + @pytest.mark.parametrize("key", [0, 1, 2, "parsnips", "time", -1]) def test_wc_dis_results_keyerrs(wc_rna, key): with pytest.raises(KeyError, match=f"{key}"): wc_rna.results[key] + + +def test_minor_dist(strand): + strand1 = ResidueGroup([strand.residues[2], strand.residues[19]]) + strand2 = ResidueGroup([strand.residues[16], strand.residues[4]]) + + MI = MinorPairDist(strand1, strand2) + MI.run() + + assert MI.results.distances[0, 0] == approx(15.06506, rel=1e-3) + assert MI.results.distances[0, 1] == approx(3.219116, rel=1e-3) + + +def test_major_dist(strand): + strand1 = ResidueGroup([strand.residues[1], strand.residues[4]]) + strand2 = ResidueGroup([strand.residues[11], strand.residues[8]]) + + MA = MajorPairDist(strand1, strand2) + MA.run() + + assert MA.results.distances[0, 0] == approx(26.884272, rel=1e-3) + assert MA.results.distances[0, 1] == approx(13.578535, rel=1e-3)