Identifying peptide oligomeric states using Python LOOS

[acnash]

Dear LOOS developers and users,

I thought I would put up a discussion here to avoid reinventing the wheel (in case it exists).

Objective:
With the help of LOOS I'm looking at writing a python function that will report the oligomeric state of a system.

System:
Using Gromacs, I'm running a martini CG simulation of four transmembrane (TM) helices in a POPC:DPPC bilayer system. The simulation runs for N steps and saves M frames to a .trr gromacs trajectory.

Outcome:
On a single coordinate frame by coordinate frame basis (in the .trr file) I would like to report the oligomeric state of the four TM helices. Thus, by the end of the trajectory, I should have M (number of frames) rows of oligomeric state data. For examples, a single saved coordinate frame will of four helices report one of four oligomeric states:

• 1:1:1:1 (all are monomeric)
• 1:3 (monomer and a trimer)
• 2:2 (two dimers)
• 4 (all four are together)

I plan on automating these processes to run over several hundred short CG MD simulations. I'm happy doing the python grunt work, I would like to know (a) doing this exist in any shape or form already, (b) if not, which particular LOOS classes/functions would you envisage using?

I imagine an interaction distance of e.g., 1.4 nm, between n helices is essential. If within or less than that distance, those n helices are considered to be in one of the four above states.

Many thanks
Anthony Nash
University of Oxford

[agrossfield]

HI Anthony,

There isn't a tool that does exactly what you want, but the pieces are there to make it easy. Almost everything you need is in the class AtomicGroup. The steps as I'd see them:

• Use gmxdump2pdb.pl to make a fake psf to use as the system file instead of the gro file. This has the advantage of having connectivity information as well as a segid (which is useful for selection.
• Read the system in.
• Select the protein with something like segid

all_chains = loos.selectAtoms(system, 'segid ="PRTN"')
chains = all_chains.splitByMolecule() # or, if you're using a gro file, you'll need to do something else to get the individual chains

• As far as the actual logic, I see 2 possibilities. The first would just be the centroid-centroid distance between them, which would look something like

c0 = chains[0].centroid()
c1 = chains[1].centroid()
dist_sq = c0.distance2(c1, system.periodicBox()) # compute the distance squared, respecting periodicity
if (dist_sq < threshold_sq): 
    # do something to store that they're in contact

This will work fine if you know that the helices are going to stay more or less vertical in the membrane. If you don't want to assume that, then you could use AtomicGroup.contactWith(), which does the double loop over atoms in 2 groups, and returns true if at least a certain number of atoms pairs are within a given distance. This approach will be a bit slower, but it'll handle arbitrary configurations more gracefully. Again, you'll want to explicitly feed it the periodic box information -- almost all distance calculations in loos have 2 versions, one that computes the simple distance and one that respects periodicity.

I can see a couple of ways to figure out how many aggregates there are given the identification of pairs. I'll lay them out if you'd like, but I suspect you already have similar or better ideas. There might even be something in scipy (clustering, perhaps?) that does it, but with only 4 chains I suspect that would be overkill.

Just as an aside, I'd recommend using xtc instead of trr for your trajectories. If I recall correctly, trr is a much more verbose file format, and your storage will quickly become an issue.

Hope this helps (and please ask more questions if it wasn't!),

Alan

PS I could easily imagine this being generally useful. Let me know if you're interested in adding the tool to LOOS when it's done.

[acnash]

You're welcome. I hope it helps. Disclaimer: there is the possibility that I've done something wrong, but I'm always happy to be corrected. In the meantime, I can work by residue count and divide the system up accordingly.

[acnash]

Final update. I've looked at pulling out the individual AtomicGroups that represent the BB bead for each of the four helices. As each helix is identical in length (26), I thought a for loop through the first 26 residue atomic groups (based only on the bond connectivity) would prove to be a simple workaround. It turns out that I'm looking at 34 atomic groups rather than the expected 26. Those additional atomic groups look to be empty:

<GROUP PERIODIC='(101.969,78.428,137.396)'>
   <ATOM INDEX='9531' ID='9532' NAME='BB' RESID='6944' RESNAME='TYR' COORDS='(20.94,12.994,66.617)' VELOCITIES='(0,0,0)' ALTLOC='' CHAINID='' ICODE='' SEGID='PRTA' B='1' Q='1' CHARGE='0' MASS='1' ATOMICNUMBER='-1' MASK='203'>
  <BOND>9533</BOND>
</ATOM>
</GROUP>
<GROUP PERIODIC='(101.969,78.428,137.396)'>
</GROUP>
<GROUP PERIODIC='(101.969,78.428,137.396)'>
</GROUP>
<GROUP PERIODIC='(101.969,78.428,137.396)'>
</GROUP>
<GROUP PERIODIC='(101.969,78.428,137.396)'>
   <ATOM INDEX='9536' ID='9537' NAME='BB' RESID='6945' RESNAME='LEU' COORDS='(23.037,11.067,67.843)' VELOCITIES='(0,0,0)' ALTLOC='' CHAINID='' ICODE='' SEGID='PRTA' B='1' Q='1' CHARGE='0' MASS='1' ATOMICNUMBER='-1' MASK='203'>
  <BOND>9538</BOND>
</ATOM>
</GROUP>

Using the files I linked above, the test code for this:

#create system and load in trajectory
systemAG = loos.createSystem("test.pdb")
#traj = loos.pyloos.Trajectory("../EPOC_0/0/npt_prod.trr", system)

#select "domains"
allChainsAG = loos.selectAtoms(systemAG, 'segid=="PRTA"')  #contains all 236 beads
chainsList = allChainsAG.splitByMolecule() #this splits into 136 parts (it should be 4*26 = 104)

countInt: int = 0
for residueAG in chainsList:
    bbAG = loos.selectAtoms(residueAG, 'name=="BB"')
    print(bbAG)
    countInt = countInt + 1
    if countInt > 26:  #this should be ok, but the additional empty AtomicGroups means 26 is too short
        break

Again, this may be an unexpected consequence of something I've done wrong from the start, but I thought further information can only help.

[agrossfield]

I think your problem might be related to issue #47, which we just fixed. In any event, @tromo checked and the current version works correctly, while he can reproduce the error if he rolls it back.

So, I think all you need to do is pull and (optionally) reinstall:

git pull
scons PREFIX=$CONDA_PREFIX install

(you will pick up some new stuff, so the build will take a minute)

Let me know what happens -- if it works, I'll mark this as solved.

[acnash]

That's great. I'll give this a try and let you know how it goes.

[acnash]

Just an update as I understand I've left this issue hanging. I've just finished building a new GPU server. Hopefully, over the next week, I'll find the time to try out this change.