pythonargparsemdanalysis

TPRParser in MDAnanlysis has no attribute 'parser'?


I'm trying to use the TPRParser to parse the TPR file generated from GROMACs. Unfortunately it throws a error which I've never seen before:

                Traceback (most recent call last):
                    line 43, in <module>
          topo = mda.topology.TPRParser.TPRParser(TPR).parser(tpr_resid_from_one=True)
          AttributeError: 'TPRParser' object has no attribute 'parser'

Could someone please help me to remove this error?
Here is my complete code and the data files are available here

import argparse

# define and parse command line arguments
parser = argparse.ArgumentParser()
parser.add_argument('input', metavar='INPUT', help='XTC input file')
parser.add_argument('--output', metavar='FILENAME', help='H5MD output file for saving trajectory data')
parser.add_argument('-v', '--verbose', action='store_true', help='be verbose')
args = parser.parse_args()

import MDAnalysis as mda
import numpy as np
import matplotlib.pyplot as plt
import os, sys
from collections import OrderedDict

XTC = args.input

TPR = os.path.splitext(XTC)[0] + '.tpr'
if not os.path.isfile(TPR):
    print("Error: TPR file {0} does not exists".format(TPR))
    sys.exit(1)

# the following list was generated manually using
# grep "^ *1[A-Z]" -B1 lp400start.gro | grep -v "^ *1[A-Z]" | cut -c 1-5

protein_seq_length = [ 38, 38, 137, 137, 252, 252, 576, 576, 1092, 1092, 1016, 1016, 1285, 1285 ]

protein_resid_range = []
i = 0
for n in protein_seq_length:
    protein_resid_range += [ slice(i, i + n) ]
    i += n
print("Expecting {0:d} residues in {1:d} proteins".format(i, len(protein_resid_range)))

# build reduced universe to match XTC
# (ignore warnings that no coordinates are found for the TPR)
topo = mda.topology.TPRParser.TPRParser(TPR).parser(tpr_resid_from_one=True)
u0 = mda.Universe(topo)
u = mda.Merge(u0.select_atoms("not resname W and not resname WF and not resname ION"))
u.load_new(XTC)

# this selection does not work, it mixes the atoms of both EPHA proteins
# FIXME delete this output
print(u.select_atoms("segid seg_0_EPHA"))

# select proteins through range of residue IDs
for idx in protein_resid_range:
    res = u.residues[idx]
    print("Protein of type {0:s} is composed of {1:d} residues and {2:n} atoms:\n{3}\n").format(res[0].segment.segid, len(res), len(res.atoms), res)

# segments 
protein_segments = u.segments

# build the fragments
# (a dictionary with the key as the protein name -- I'm using an
# OrderedDict so that the order is the same as in the TPR)
protein_fragments = OrderedDict((seg.segid[1:], seg.atoms.fragments) for seq in protein_segments)

# this doesn't seem to be helpful either FIXME delete this output
for f in protein_fragments:
    print(f)

# analyze trajectory
timeseries = []
for ts in u.trajectory[1:4]:
    coms = []
    for name, proteins in protein_fragments.items():
        # loop over all the different proteins
        # unwrap to get the true COM under PDB (double check!!)
        coms.extend([p.center_of_mass(unwrap=True) for p in proteins])
    timeseries.append(coms)
timeseries = np.array(timeseries, dtype=float)

if args.output:
    with opne(args.output, 'w') as outfile:
        # make clear distinction in frames
        outfile.write('# Array shape: {0}\n'.format(timeseries.shape))
        outfile.write('# New frame: {0}\n'.format(u.trajectory.time))
        # Iterating through a n-dim array
        for data in timeseries:
            np.savetxt(outfile, data, fmt='%-7.2f')

Solution

  • You do not need to use the TPRParser explicitly. Instead of the lines

    ### do NOT use this code
    topo = mda.topology.TPRParser.TPRParser(TPR).parser(tpr_resid_from_one=True)
    u0 = mda.Universe(topo)
    

    use just

    u0 = mda.Universe(TPR, tpr_resid_from_one=True)
    

    The Universe() class uses the correct parser (based on the filename extension) automatically. It is almost never necessary to use a parser explicitly. Just use Universe().