How to select a Ligand/Sugar with a Glycam nomenclature?

[Noah603]

Hey guys!

I have performed an MD sim with my Protein and a trisaccharide (2'FL) and penta- and hexa-saccharides. The Ligand(s) are defined in the prmtop file by the Glycam nomenclature: for example for 2'FL: ROH 4GA 2LB 0fA.
Is there a nice way to just call them collectivly? Or do I always call them individually, like I am currently doing. I tryed to change in the prmtop file the individuall names to just LIG and use that. That kinda works but is not really pleasing.

Best greetings, Noah

[hmacdope]

Hi Noah,

Can you give me an example (in code) of how you are running your selection?

Cheers,

Hugo

[hmacdope]

I'm going to raise that as an issue @orbeckst

[orbeckst]

I think the solution rmsd_for_atomgroups2() in #4562 (reply in thread) already does the right thing: superimpose on the reference (in the example, the protein backbone) and compute RMSD for different groups. (I just would never use np.round() on any data and only round when presenting data in a table... but that's just a small criticism of the otherwise nice talktorial on analyzing MD simulations).

The question

but it would be also nice to see the overall RMSD of the ligand

could be addressed with the selection that @hmacdope suggests in #4562 (reply in thread) by also selecting the whole sugar (here written as sugar and as the selection "resname 4GA 2LB 0FA")

ligand_name ='4GA'
ligand_name2='2LB'
ligand_name3='0FA'

# resnames for the whole carbohydrate
sugar = f"{ligand_name} {ligand_name2} {ligand_name3}"

rmsd = rmsd_for_atomgroups2(u, "backbone", ["protein", f"resname {sugar}", f"resname {ligand_name}",f"resname {ligand_name2}", f"resname {ligand_name3}"])
rmsd.plot(title="RMSD of protein and ligand")
plt.ylabel("RMSD (Å)");

This should add a RMSD calculation for the whole carbohydrate selection.

[Noah603]

Thank you both for the answers and the discussion!

Btw, I couldn't access the GLYCAM

Yeha same here, the GLYCAM server seems to be down. Can not access the side at all. You can see the nomencalture as well in the amber tutorial here on page 48.

It would be an interesting project to collect all GLYCAM residue names so that we can create a selector "sugar" (or "carbohydrate") that selects anything with those residue names (similar to how "protein" selects anything that has the 20+ amino acid residue names).

Thats what I was searching for! I think this could be very nice, if there is a simple implementation like that into MDanalysis.

ligand_name ='4GA'
ligand_name2='2LB'
ligand_name3='0FA'
resnames for the whole carbohydrate
sugar = f"{ligand_name} {ligand_name2} {ligand_name3}"

But this is also very nice way of doing it, since you can tune the residues you want to see in the analysis. If you have for example a hexasaccharide in which you want to focus during the analysis on only 3 residues. I am going to use that :)

just would never use np.round() on any data and only round when presenting data in a table.

Good spot! Gonna have that in mind. I really like that talktorial 19 and 20, because it is a great starting point for someone new to the world of MD sims. I was searching for a notebook where the intention behind it was: This is how I anaylsis my MD Protein-Ligand sims, that creates the most important graphs for you, you just need to put your file into it. For example if you get a new student you can just hand it over and get good graphs out.

[Noah603]

Hey guys!
I ran into another problem which I wanted to make you aware, incase you want to implement the Glycom nomenclature.
Imagine the sugars sequence is sequence: ROH 4GA VLB 4YB 0LB 4YB 0LB.
Now the above solutions do not work, since we select for the names, but we have twice the same name. In my case this lead to the effect that only one of them gets selected.
I worked around by renaming them in the prmtop name. This is not a particular nice soultion, but works.
Best, Noah

[orbeckst]

There may be some misunderstanding of how the input for the RMSD class works. Let me explain a bit.

Selections with duplicate resnames

Duplicated resnames should not matter in a normal selection. We just look for anything that has the correct resname:

If you do

sugar_1 = u.select_atoms("resname ROH 4GA VLB 4YB 0LB 4YB 0LB")
sugar_2 = u.select_atoms("resname ROH 4GA 4YB 0LB")

then you should get the same atoms selected

sugar_1 == sugar_2

should be True.

The sugar_1 atomgroup will contain all the residues that match the name in numerical order, so e.g.

print(list(sugar_1.residues))

should show the glyco-residues in order.

Selecting specific residues with the same name

If you have two sugar residues named 0LB then you need additional information to tell them apart, for instances, their residue ids.

Lets assume you have a sugar with ROH-4GA-VLB-4YB-0LB-4YB-0LB where the resids are

>>> s = u.select_atoms("resname 0LB")
>>> print(s.residues.resids)
[1 2 3 4 5 6 7]

(i.e., ROH is resid 1 and the last 0LB is resid 7).

If you want a specific sugar residue such as the second 4YB then you can do a selection such as

s4YB_2 = u.select_atoms("resid 6 and resname 4YB")

(You could probably just use "resid 6" but in some files, the resids are not unique and then you're safer to also specify the resname.)

Selecting individual residues with RMSD

Instead of renaming the residues in your file, you should use more precise selection strings:

rms.RMSD(u, u, select="backbone", 
         groupselections=["protein", "resid 1 and resname ROH", 
                          "resid 4 and resname 4YB", "resid 6 and resname 4YB")

Here you would calculate the RMSD for protein, ROH, and both 4YB separately.

If you want to automate generating groupselections, something like the following may work

sugar = u.select_atoms("resname ROH 4GA VLB 4YB 0LB")
resnames = sugar.residues.resnames
resids = sugar.residues.resids
groupselections = [f"resid {resid} and resname {resname}" for resid, resname in zip(resids, resnames))]
groupselections.prepend("protein")

[hmacdope]

@Noah603 could you mark as resolved if your question was answered?

[orbeckst]

@Noah603 instead of closing the discussion could you please let us know if your problem was solved (then click the "Mark as answer") or say what's not working.

I am re-opening the discussion so that it remains visible. Closing pretty much removes it from sight, which makes it hard for others to learn from it. This is a bit of a shortcoming in GH Discussions.

[Noah603]

ah sorry! mb