Feature Nucleic Acid Backbone

[fjclark]

Is this pull request to fix a bug, or to introduce new functionality?

This PR allows BSS to recognise nucleic acid backbones e.g.:

protocol = BSS.Protocol.Equilibration(restraint="backbone")

I've tried to implement this in the same way as for protein backbones, copying residue names over from MDAnalysis. A couple of slight complications are:

• I've had to use regex syntax to allow searching for atom names containing '
• I've created the methods containsAminoAcids and containsNucleicAcids in System to make this section of code cleaner.

If this introduces new functionality...

Changes proposed in this pull request:

• Add nucl_res - a list of nucleic acid residues, taken from MDAnalysis.
• Add some logic when creating the Amber mask to ensure that only atoms which exist in the system are included in the mask (otherwise Amber produces the error: Error: an illegal memory access was encountered launching kernel kClearForces)
  • Add the methods containsAminoAcids and containsNucleicAcids to System to make the above logic cleaner.

• I confirm that I have merged the latest version of devel into this branch before issuing this pull request (e.g. by running git pull origin devel): y
• I confirm that I have added a test for any new functionality in this pull request: n - I couldn't see tests related to the functionality which was present before.
• I confirm that I have added documentation (e.g. a new tutorial page or detailed guide) for any new functionality in this pull request: n
• I confirm that I have permission to release this code under the GPL3 license: y

Suggested reviewers:

@lohedges, @chryswoods

[lohedges]

Many thanks, this looks great. Will take a closer look tomorrow. The CI failure was just because the required Sire packages were still building so it didn't find the dependency for that platform. I've re-triggerred.

[lohedges]

Thanks for this. Other than my review comments, do you have any simple tests for the new logic? If so, it would be good if they could re-use some of the existing test files. If you need to add anything, then I can upload it to the website as needed.

If it's not too much effort, could you mirror the new functionality in the Exscientia Sandpit too? I can always do this post-merge if needed.

[lohedges]

Could all of these be moved outside of the class? I'd rather they were hidden from the user since they are only used internally. It would probably be best to place them in a utility module, e.g. _utils.py, then import them privately in _system.py to use where required, e.g.:

from ._utils import prot_res as _prot_res
...

For simplicity, they could even be private in the _utils.py module itself, e.g. just call it _prot_res, etc.

[fjclark]

Of course - have updated and will push shortly.

[lohedges]

Is it possible to match these to the other similar search functions, e.g. have functions called nNucleicAcids and getNucleicAcids, etc. The first would give you the number of matches, the second returns the atoms matching the search, if any. (Could probably do the same for ions, too.)

[fjclark]

Sure - I've added these. I've skipped getIons for now because this seems a bit more complicated. The ions list only contains atom names, so this would have to be getMonoatomicIons, and I'd have to make sure that the atoms aren't part of molecules. Happy to try and implement something if you would like though·

[lohedges]

Thanks for catching.

[fjclark]

Thanks for this. Other than my review comments, do you have any simple tests for the new logic? If so, it would be good if they could re-use some of the existing test files. If you need to add anything, then I can upload it to the website as needed.

Thanks very much for the quick review. I've added some tests. Are any of the website systems nucleic acids? If not, could you please add compressed versions of this RNA system: rna.tar.gz?

[lohedges]

I'm not sure. I'll add that to the website now and let you know when I'm done.

[lohedges]

All done. These can now be loaded with:

s = BSS.IO.readMolecules(BSS.IO.expand(BSS.tutorialUrl(), ["rna_6e1s.rst7", "rna_6e1s.prm7"]))

In the test suite the url is a variable within the conftest.py file, so you can just use that.

[fjclark]

Perfect, thanks very much. I'll update the tests and copy over to the sandpit now.

[lohedges]

Many thanks, this looks great 👍

[lohedges]

Sorry, just to check... I know you have tests in place, but does this work as documented:

    def getNucleotides(self, property_map={}):
        """
        Return a list containing all of the nucleotide residues in the system.
        Parameters
        ----------
        property_map : dict
            A dictionary that maps system "properties" to their user defined
            values. This allows the user to refer to properties with their
            own naming scheme, e.g. { "charge" : "my-charge" }
        Returns
        -------
        residues : [:class:`Residue <BioSimSpace._SireWrappers.Residue>`]
            The list of all nucleotide residues in the system.
        """
        search_string = "(resname " + ",".join(_nucl_res) + ")"
        try:
            residues = list(self.search(search_string, property_map))
        except:
            residues = []
        return residues

This should return a list of residues. However, by default, searches will return atoms matching the search, i.e. the atoms within the residue. Do you want to return the atoms, or should this be converted to residues. For example:

import BioSimSpace as BSS

s = BSS.IO.readMolecules(BSS.IO.expand(BSS.tutorialUrl(), ["ala.top", "ala.crd"]))

# Search for the ALA residue. This returns the atoms in residue ALA.
s.search("resname ALA")
<BioSimSpace.SearchResult: nResults=10>

# Search for the residues matching the search. Now there is a single match.
s.search("resname ALA").residues()
<BioSimSpace.SearchResult: nResults=1

[fjclark]

Sorry, just to check... I know you have tests in place, but does this work as documented:

    def getNucleotides(self, property_map={}):
        """
        Return a list containing all of the nucleotide residues in the system.
        Parameters
        ----------
        property_map : dict
            A dictionary that maps system "properties" to their user defined
            values. This allows the user to refer to properties with their
            own naming scheme, e.g. { "charge" : "my-charge" }
        Returns
        -------
        residues : [:class:`Residue <BioSimSpace._SireWrappers.Residue>`]
            The list of all nucleotide residues in the system.
        """
        search_string = "(resname " + ",".join(_nucl_res) + ")"
        try:
            residues = list(self.search(search_string, property_map))
        except:
            residues = []
        return residues

This should return a list of residues. However, by default, searches will return atoms matching the search, i.e. the atoms within the residue. Do you want to return the atoms, or should this be converted to residues. For example:

import BioSimSpace as BSS

s = BSS.IO.readMolecules(BSS.IO.expand(BSS.tutorialUrl(), ["ala.top", "ala.crd"]))

# Search for the ALA residue. This returns the atoms in residue ALA.
s.search("resname ALA")
<BioSimSpace.SearchResult: nResults=10>

# Search for the residues matching the search. Now there is a single match.
s.search("resname ALA").residues()
<BioSimSpace.SearchResult: nResults=1

Thanks for pointing this out. I've added .residues() and a check of isintance() in the tests. These functions were returning lists of residues before, but I see now that they could return lists of atoms if only one residue was found. I didn't realise that the return type depended on the number of results, e.g.:

import BioSimSpace as BSS

s = BSS.IO.readMolecules(BSS.IO.expand(BSS.tutorialUrl(), ["ala.top", "ala.crd"]))

# Search for the ALA residue. This returns the atoms in residue ALA.
s.search("resname ALA")
<BioSimSpace.SearchResult: nResults=10>

# Search for the ALA residue and the cap ACE. This returns 2 residues, rather than atoms
s.search("resname ALA,ACE")
<BioSimSpace.SearchResult: nResults=2>

[lohedges]

No problem, I agree that it's confusing. This changed in 2023.3.0 when the Sire search API was updated. As such, it's now best to use the .atoms(), .residues(), or .molecules() functions to cast the result to whatever you want. It does work with the new API, e.g.:

s = BSS.IO.readMolecules(BSS.IO.expand(BSS.tutorialUrl(), ["ala.top", "ala.crd"]))

sire_sys = sr.system.System(s._sire_object)

sire_sys["resname ALA"]
Residue( ALA:2   num_atoms=10 )

However, using the old API allows the search to respect the property map. This should all be sorted out once we fully migrate across to the new API in future.

[fjclark]

OK, thanks for the explanation.

[lohedges]

Previously I had some logic that would automatically convert the results to the minimal object, i.e. atoms if there was only one atom, residues if there was only one residue, etc, but this broke with the API update.