addRingClosingBond in molmodel

[Samuel Flores]

Guys,

Not sure if I should put this here, I don't know if anyone monitors the molmodel forum..

I want to be able to make posttranslational modifications to amino acids in a peptide. For example, let's say I have some sort of non beta branched "stub" amino acid, and I want to add arbitrary atoms to represent some non-canonical residue. I want the user to be able to do this at runtime. I can of course apply the usual bondAtom etc commands to amino acids, but the added atoms don't end up getting residue numbers. Residues are not really compounds of their own that can be retrieved. This lack of residue numbers ends up causing problems everywhere, including writing out PDB files (the added atoms are simply not printed).

It occurs to me that addRingClosing bond would be perfect for this, if it could be applied across Compounds, rather than only within a single Compound. Is there some fundamental reason why this wouldn't work? It's not implemented in molmodel currently. The idea is that I would create the side chain and them join it to the peptide using addRingClosingBond.

Sam

[Michael Sherman]

Hi, Sam. Couldn't each of the two Compounds potentially have their own residue numbers? That would be confusing if residue numbers are supposed to transfer from one to the other. I don't think I understand the connection to addRingClosingBond() -- if I'm remembering it correctly it doesn't add any atoms so I wouldn't think there wouldn't be a residue numbering aspect to it.

Would it make sense to add something more straightforward, like addCompoundToResidue() or whatever? Then it wouldn't be ambiguous that the new compound should inherit the residue number.

Sorry if I'm misunderstanding -- it has been a while since I thought about Molmodel!

Regards,
Sherm

[Michael Sherman]

Another possibility that comes to mind (not sure if it makes sense) would be to add a setResidueNumber() call to atoms that would directly supply the missing residue number.

[Samuel Flores]

Hi Sherm,

Thanks for your help. Following what you said, I figured out that ResidueInfo is loaded when the biopolymer is first polymerized, so I needed to explicitly add the prosthetic atoms afterwards. It's not done automatically by the bondCompound method. So I guess I have to make a loop to do that. Not very convenient, but probably unavoidable.

The other inconvenience, which I do think I should be able to avoid, is setting the BiotypeIndex's for all the prosthetic atoms. You see, I am using bondCompound to add an entire compound to a suitably prepared amino acid (one with a free outboard bond). I would think the added compound would take its biotypes with it. But no. I verify (line 144) that just after adding the compound, the atoms of that compound do not remember their BiotypeIndex's. I still have the original compound lying around, so I query the BiotypeIndex of the appropriate atom and then label the corresponding atom on the Biopolymer. That works. I verify (line 146) that the atom on the biopolymer now has the correct BiotypeIndex (output at the bottom of this post).

So I should just set up a loop to copy the BiotypeIndex's from the original compound to the one that has been attached to the Biopolymer. I can do this, but I see that in e.g. Protein.h, bondCompound is not followed by any such transfer of BiotypeIndex's. Is it because all the BiotypeIndex's are set at the end somewhere? Or is there some way to force bondCompound to remember the BiotypeIndex's ?



131 myBiopolymer.bondCompound(
132 residueIndexString , // This should be the residue index. Will be prepended with / to atom name
133 myCompound ,
134 bondName,
135 bondLength,
136 myDihedral);
139 Compound::AtomIndex myAtomIndex = myBiopolymer.getAtomIndex( Compound::AtomPathName (longAtomName) );
140 Compound::AtomPathName tempAtomName("HG");
141 ResidueInfo myResidueInfo = myBiopolymer.updResidue(ResidueInfo::Index(0));
142 myBiopolymer.updResidue(ResidueInfo::Index(0)) .addAtom(
143 myAtomIndex, specificAtomName);
144 cout<<__FILE__<<":"<<__LINE__<<" myBiopolymer.getAtomBiotypeIndex(Compound::AtomIndex( myBiopolymer.getAtomIndex(longAtomName))) : "<< myBiopolymer.getAtomBiotypeIndex(Compound::AtomIndex( myBiopolymer.getAtomIndex(longAtomName))) <<endl;
145 myBiopolymer.setBiotypeIndex(longAtomName, compoundToAdd.getAtomBiotypeIndex(Compound::AtomIndex( 0)));
146 cout<<__FILE__<<":"<<__LINE__<<" myBiopolymer.getAtomBiotypeIndex(Compound::AtomIndex( myBiopolymer.getAtomIndex(longAtomName))) : "<< myBiopolymer.getAtomBiotypeIndex(Compound::AtomIndex( myBiopolymer.getAtomIndex(longAtomName))) <<endl;


output:
/Users/Sam/svn/RNAToolbox/trunk/src/BiopolymerClass.cpp:144 myBiopolymer.getAtomBiotypeIndex(Compound::AtomIndex( myBiopolymer.getAtomIndex(longAtomName))) : -1111111111
/Users/Sam/svn/RNAToolbox/trunk/src/BiopolymerClass.cpp:146 myBiopolymer.getAtomBiotypeIndex(Compound::AtomIndex( myBiopolymer.getAtomIndex(longAtomName))) : 91


Many thanks

Sam

[Michael Sherman]

Hi, Sam. I took a look through the code and didn't spot anywhere that the biotype was being intentionally forgotten. And I see that Biopolymer::assignResidueBiotypes() explicitly preserves existing Biotypes if they are present (but that assumes they made it through the Compound-building process).

I would step through the bondCompound() call in the debugger while attaching a very small compound (one atom?) with a known preassigned Biotype to check the atom-copying code (or you could sprinkle in some judicious printfs if you prefer that approach). As best I can tell it should have copied the Biotype if it was present initially. Maybe there is a copy constructor or copy assignment operator incorrectly implemented?

I suppose in general it is a good idea to finish building the Compound completely and then assign Biotypes at the end, especially if there is any context-sensitive typing going on. But I think if you prefer to carry through the existing Biotypes it should be just a matter of tracing why they didn't get copied in the first place, and then fix that.

Regards,
Sherm