MDAnalysis.analysis.rms.RMSD "select" group appears to overwride prior sub-selection

[Chris Neale]

Hello,

with python 2.7.12 and MDAnalysis 0.16.2 I find that I can make a selection with .select_atoms() and that selection is honored during print, but then the RMSD() command seems to ignore it and jump back to a section from the entire structure. The issue is that I have an overall structure with some variable number of proteins and I want to compute the RMSD of each of those sampled proteins to a single reference structure.

here is the relevant code that I used:
<pre>
thisProtein = pdbPatch.select_atoms("resid %d:%d" % (proteinStartResidue, proteinStartResidue + numResPerRas - 1))
print ("sample has %d residues" % len(thisProtein.residues))
print ("reference has %d residues" % len(pdbRef.residues))
R = RMSD(thisProtein, pdbRef, select="name BB", filename="meow.cat", weights=None, tol_mass=999.99)
R.run()
</pre>

and the output for the printing of number of residues in the initial selection is:
<pre>
sample has 186 residues
reference has 186 residues
</pre>

but then the RMSD() function gives the error:
<pre>
MDAnalysis.exceptions.SelectionError: Reference and trajectory atom selections do not contain the same number of atoms: N_ref=184, N_traj=920
</pre>

Here is the actual complete code snippit:
<pre>
            import MDAnalysis as mda
            from MDAnalysis.analysis import align
            from MDAnalysis.analysis.rms import RMSD

            # read in the patch to MDAnalysis
            pdbPatch = mda.Universe(fnameTopo, fnameTraj)

            fnameRef = "/tmp/reference.pdb"
            pdbRef = mda.Universe(fnameRef, fnameRef)

            numResPerRas=186     # defined for this system
            numAtomPerRas=42     # defined for this system
            numProteins=0
            analysisFileName = "sim-{}KRAS-{}-ps.rmsdGdomain".format(numKRAS, i)
            analysisFileName = os.path.join(self.inpath, analysisFileName)
            analysisFileHandle = open(analysisFileName, "w")
            for ppr in pdbPatch.residues:
                if ppr.resname == "THR":
                    proteinStartResidue = ppr.resid
                    break
            while pdbPatch.residues[proteinStartResidue].resname == "THR":
                numProteins += 1
                thisProtein = pdbPatch.select_atoms("resid %d:%d" % (proteinStartResidue, proteinStartResidue + numResPerRas - 1))
                print ("sample has %d residues" % len(thisProtein.residues))
                print ("reference has %d residues" % len(pdbRef.residues))
                R = RMSD(thisProtein, pdbRef, select="name BB", filename="meow.cat", weights=None, tol_mass=999.99)
                R.run()
                analysisFileHandle.write("%s" % R.rmsd[0][2])
                proteinStartResidue += numResPerRas
            analysisFileHandle.close()
</pre>

and the complete output from this code:
<pre>
sample has 186 residues
reference has 186 residues
ERROR:MDAnalysis.analysis.rmsd:Reference and trajectory atom selections do not contain the same number of atoms: N_ref=184, N_traj=920
None
Traceback (most recent call last):
  File "cganalysis.py", line 256, in <module>
    main()
  File "cganalysis.py", line 251, in main
    cg.run()
  File "cganalysis.py", line 187, in run
    R = RMSD(thisProtein, pdbRef, select="name BB", filename="meow.cat", weights=None, tol_mass=999.99)
  File "/home/cneale/.local/lib/python2.7/site-packages/MDAnalysis/analysis/rms.py", line 415, in __init__
    raise SelectionError(err)
MDAnalysis.exceptions.SelectionError: Reference and trajectory atom selections do not contain the same number of atoms: N_ref=184, N_traj=920
</pre>

Thank you for your suggestions,
Chris.

[Oliver Beckstein]

Hi Chris,

There has been at least one fix in RMSD in 0.17.0 that might help (https://github.com/MDAnalysis/mdanalysis/issues/1487 ). Would you be able to try out 0.17.0?

Because of issue #1487 in 0.16.2 I am not sure if it makes sense for me to suggest other ways to make selections because it might not work with the 0.16.2 code anyway… but I’ll try:  the select keyword can take a dict {‘mobile’: selection_string_1, ‘reference’: selection_string_2} that is applied separately and this might work when you provide universes to the RMSD class instead of AtomGroups … maybe.

Oliver

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To post to this group, send email to mdnalysis-...@googlegroups.com.
Visit this group at https://groups.google.com/group/mdnalysis-discussion.
For more options, visit https://groups.google.com/d/optout.

--

Oliver Beckstein

orbe...@gmail.com * orbe...@gmx.net

[Chris Neale]

Hey, yeah, that completely fixes the problem. No errors and RMSD results agree with gromacs gmx rms results to within 10 fm precision.

Thank you Oliver!

To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discussion+unsub...@googlegroups.com.

To post to this group, send email to mdnalysis-...@googlegroups.com.
Visit this group at https://groups.google.com/group/mdnalysis-discussion.
For more options, visit https://groups.google.com/d/optout.