Autocorrelation of protein backbone

[distribut...@googlemail.com]

Hi all, 

I've Googled and searched this Google Group and I'm surprised to find no mention (could be a result of my poor search skills) of a function in MDAnalysis that calculated the autocorrelation of a protein (or a selection). I've looked through the MDAnalysis API docs and found an autocorrelation in persistence length of a polymer, and as a part of HB survival analysis. 

I thought I could retrieve the autocorrelation results using polymer.PersistenceLength but it through an error saying that there were too many bonds in the backbone. 

Have I missed a straight-forward autocorrelation function in the MDAnalysis package, or does anyone know of a User built autocorrelation function using the MDAnalysis codebase?  

Many thanks and kind regards,

Anthony 

[richard...@gmail.com]

Hi Anthony

I think the polymer Persistence length might work for you, but you'll have to identify the "backbone" manually, I think the function you found will reorder a selection but not perform a selection.  Ways I can think of identifying the backbone include based on atom names if your residues have good names, or a shortest path algorithm between each end residue if they don't.

Richard

[distribut...@googlemail.com]

Thanks Richard,

I've thrown in C and CA just as a test, but now I'm getting:

ValueError: Could not find starting point of backbone, is the backbone cyclical?

I've Ctrl-F through the source code and performed a quick Google search, but I can't find how to treat a cyclical protein. In context, the cyclical part of the system is an idealised model-collagen of x3 polypeptide chains, each of approx 36 residues in length. To replicate the mock fibril-like environment (or 1300+ residues), I've turned the collagen cyclical by forming a bond back on itself over the unit cell boundary.

The code thus far:

from MDAnalysis.analysis import polymer

#providing the topology has demarcated each chain (MMP1, each collagen polypeptide), we can retrieve the chains.
chains = universe.atoms.fragments

backbones = [chain.select_atoms("name CA C N") for chain in chains]

sortedBackbones = [polymer.sort_backbone(bb, cyclical=True) for bb in backbones]

persistenceLength = polymer.PersistenceLength(backbones)
persistenceLength.run()
print(persistenceLength.result)

Many thanks

Anthony

[Richard Gowers]

Yeah the code won't handle cyclical things, when I wrote it a long time ago I was just dealing with regular polymers for materials reasons.  Off the top of my head I wouldn't know how to do cyclical things without getting the double counting wrong.

If your input topology is already in the correct order, you could probably get away with not sorting the input, and because it's not looking at the cyclical nature it won't include the data from the ends of chains as much.  Then you could "rotate" the AtomGroups with something like ag = ag[len(ag) // 2:] + ag[:len(ag) // 2] run it again and check if that got different results (because the cycle had different behaviour in different sections).

--
You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/iyJmjxNKVbY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/f3b65417-2ef2-4153-85ae-bc1dfd16e8f6n%40googlegroups.com.

[distribut...@googlemail.com]

Brilliant; thank you for the help. I'm giving it a try now, but I still have a very basic understanding of MDAnalysis, and it's causing a bit of struggling with deciphering return types and argument types from examples. 

Firstly,  the topology order is correct (do people make topologies out of order?! I would love to know for my education reasons why this could happen). I tried the code across the four chain fragments (enzyme, x3 collagen peptides), but MDAnalysis threw an error about each chain being a different length. 

So I'm now just trying it out on one of the collagen peptides, but my understanding of the MDAnalysis data types gets a bit sketchy. I'm getting the error: TypeError: object of type 'Atom' has no len() for the code:

from MDAnalysis.analysis import polymer
MMP1_BACKBONE_ATOMS = "index 0:5820 and backbone"

topologyFile = "C:/Users/yewro/Documents/MINES_LECTURES/test_topology.pdb"
trajectoryFile = "C:/Users/yewro/Documents/MINES_LECTURES/test_trajectory.trr"

universe = mda.Universe(topologyFile, trajectoryFile)
backbone = universe.select_atoms(MMP1_BACKBONE_ATOMS)
backbone = backbone.select_atoms("name CA C N")

persistenceLength = polymer.PersistenceLength(backbone)
persistenceLength.run()
print(persistenceLength.result)

Thanks

Anthony

[distribut...@googlemail.com]

Kind of sorted. 

.PersistenceLength() expects a List of AtomGroup objects as it's argument. So, I just added a single AtomGroup into a list, passed it in, called run and pulled out the results. 

Part of the final code:

#get all the chains from the system; the enzyme, and three polypeptide chains making up the substrate

chains = universe.atoms.fragments

#get the enzyme backbone atoms

backbone = universe.select_atoms(MMP1_BACKBONE_ATOMS) #returns an AtomGroup

#**probably** don't need to do this next step, but this was part of a test
backbone = backbone.select_atoms("name CA C N")   #returns 1 AtomGroup

#PersistenceLength expects a list of AtomGroup objects; originally for identical length polymers and not for a single polypeptide aminoacid chain. Add the single chain in a list.
backboneList = [backbone]
 
persistenceLength = polymer.PersistenceLength(backboneList)
persistenceLength.run()

#results are stored in the attribute

print(persistenceLength.results)

Thanks

Anthony