Can you suggest any way to calculate the average structure of a molecule from multiple frames in a PDB file. I want to do it in Python preferably. I have used Biopython and MDAnalysis before but I am not sure whether these can be used to do the thing.
Can you please give me an idea how I can do this?
If you did the simulation using GROMACS, then 'gmx covar' should be useful.
It used to compute an average structure from a dynamics trajectory. For example to compute the average of the last 200 ps of a 1 ns dynamics run, then run the command:
gmx dcovar –f traj.xtc –s topol.tpr –b 801 –e 1000 –av .pdb
Something else that could help this discussion. When you calculate the average structure from a trajectory, focus only on the CA, backbone atoms, or maybe all heavy atoms. This is because the hydrogens and side chain atoms are very flexible and they may explore a large volume of phase space, and so when you visualize the average structure these atoms may not be clearly visualized.
# python; assuming one frame per PDB file; for a multi-frame PDB file, store the data (line by line) into a python list, get the list positions (indices) of the first and last atoms in each frame, and then go through the list using those indices . . .
import glob
import numpy
pdbfiles = glob.glob("*.pdb")
D1 = {}
for pdbfile in pdbfiles:
for atom in open(pdbfile):
atid = line.split()[1]
if pdbfiles.index(pdbfile) == 0:
D1[atid] = {"x": [], "y": [], "z": []}
x = float(line.split()[-3])
y = float(line.split()[-2])
z = float(line.split()[-1])
D1[atid]["x"].append(x)
D1[atid]["y"].append(y)
D1[atid]["z"].append(z)
D2 = {}
for atid in D1:
D2[atid] = {}
D2[atid]["x"] = numpy.average(D1[atid]["x"])
D2[atid]["y"] = numpy.average(D1[atid]["y"])
D2[atid]["z"] = numpy.average(D1[atid]["z"])
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics