How this project works can be separated into three parts:
- Precalculating the data
- Searching the data
- Integrating with ChemAxon
Data is precalculated using the LigandDistanceDataTreeBuilder. This class only generates data for distance interactions between ligands and other sorts of groups, so it does not generate protein-protein data. Right now, it is limited to searching through interactions that are up to 5 angstroms in length. If you want to create your own set of precalculated data, the LigandDistanceDataTreeBuilder is most likely the only file you'll have to modify too much. Add your own rules for what interactions to include and DistanceDataTree will take care of saving the interactions efficiently.
Which brings us to DistanceDataTree. This class is essentially just a tree. It is nested somewhat deeply, but because the nesting is an exact and constant depth, I elected to implement each level of the tree as its own class. So, within the DistanceDataTree class you also have 5 subclasses. They are:
OriginGroupTreeOriginElementTreeOriginAtomNameTreeTargetGroupTreeTargetElementTree
The terms "origin" and "target" are my own, and in my particular case they tend to be replaceable with the words "ligand" and "residue", respectively. However, I did not name them that way to emphasize to future users that any sort of atom will work for the origin and target data. So what does each level of the tree mean? They mean:
DistanceDataTreetakes in origin group names ('HEM', 'PO4', etc.) and returnsOriginGroupTrees.OriginGroupTreetakes in origin elements (Element.O,Element.C, etc.) and returnsOriginElementTrees.OriginElementTreetakes in origin atom names ('OXT', 'CD1', etc.) and returnsOriginAtomNameTrees.OriginAtomNameTreetakes in target group names ('ALA', 'TYR', etc.) and returnsTargetGroupTrees.TargetGroupTreetakes in target elements (Element.O,Element.C, etc.) and returnsTargetElementTrees.TargetElementTreetakes in target atom names ('OXT', 'CD1', etc.) and returns lists ofDistanceResults.
How these classes work is not very complicated but very tedious to explain, so I won't go there. It's basically just a plain old tree.
DistanceResults are the final things stored. They keep track of a distance interaction that took place in some structure. LigandDistanceDataTreeBuilder takes care of generating these guys. Stored in a DistanceResult is the PDB ID where this interaction took place, the serial numbers of the two atoms in the bond, and the distance (rounded down to the tenths place) between the two atoms.
So now that we've covered what data is being stored, here's how we precalculate this stuff.
Our final goal is to save distance results in a very efficient and fast-to-retrieve manner. We will do this by saving them in a directory structure. But first, we must find all distance results in the PDB, and next we must save them to disk. Unfortunately, this first part cannot be done in one fell swoop; the PDB is far too large and you'd run out of memory.
The solution, then, is to calculate all the distance results and save them in serialized distance data trees in groups of a thousand PDB IDs. Then, once you've gone through the whole PDB, you go back and generate the directory structure for the whole PDB, going one thousand-PDB-ID-tree at a time. It takes awhile to do this -- over 24 hours, I'd estimate (I've never done this whole process continuously). Exactly how this is implemented can be seen in the DistanceDataTreeBuilderDriver class.
The directory structure used is not very fancy. It is really the directory counterpart to the DistanceDataTrees we stored them in. The folders are layed out like this:
LigandGroup /
LigandElement /
LigandAtomName /
TargetGroup /
TargetElement /
TargetAtomName.ser
Here, TargetAtomName.ser is a serialized file that compactly saves information about a list of DistanceResults. Each distance result goes on a separate line in that file. Exactly how I save and re-load those DistanceResults is handled by DistanceResult#toSerializedForm() and DistanceResult#parseSerializedResult().
For an example of what the directory looks like, here's a real-life demonstration of what the ligand "ZZZ" has in its directory structure:
ulysse$ tree ZZZ/ | more
ZZZ/
├── C
│ ├── C10
│ │ ├── ASP
│ │ │ ├── N
│ │ │ │ └── N.ser
│ │ │ └── O
│ │ │ └── OD1.ser
│ │ ├── GLY
│ │ │ └── O
│ │ │ └── O.ser
│ │ ├── ILE
│ │ │ └── C
│ │ │ └── CD1.ser
│ │ ├── LEU
│ │ │ └── C
│ │ │ └── CD2.ser
│ │ └── PHE
│ │ └── C
│ │ ├── CD1.ser
│ │ └── CZ.ser
│ ├── C2
│ │ ├── ALA
│ │ │ └── C
│ │ │ └── CB.ser
│ │ ├── ASP
│ │ │ └── O
│ │ │ └── O.ser
│ │ ├── GLY
│ │ │ └── O
│ │ │ └── O.ser
│ │ ├── ILE
│ │ │ ├── C
│ │ │ │ ├── C.ser
│ │ │ │ ├── CA.ser
│ │ │ │ ├── CB.ser
│ │ │ │ ├── CD1.ser
│ │ │ │ └── CG1.ser
│ │ │ ├── N
│ │ │ │ └── N.ser
│ │ │ └── O
The process of saving DistanceDataTrees into directory structures is handled by DistanceDataTreeBuilderDriver. Saving and re-loading DistanceResults is taken care of by DistanceDataTreeSerializer and DistanceResult.
Once all the data has been saved into the directory structure thanks to DistanceDataTreeSerializer, then searching is done with DistanceDataTreeDirSearcher. What this directory searcher does is it takes in a DistanceQuery and looks through the directory structure for folders and files that match the parameters in the distance query. For instance, if you specify that the interaction you want takes place with a ligand called "XYZ", then we will search through the "XYZ/" directory. If you also specify the element of your ligand is Nitrogen, then we'll search through "XYZ/N/". If you don't specify the atom name of your ligand atom, then we'll search through "XYZ/N/*/" (in other words, we'll search through all directories under "XYZ/N/"). If you specify that the interaction happens on an Alanine, then the directories we'll search through must look like this: "XYZ/N/*/ALA/". Nothing terribly magical is going on here.
Once we reach the bottom of the directory structure, we arrive on a bunch "serialized" DistanceResults. Here's an example of what those serialized distance results look like:
ulysse$ cat HEM/C/C1A/ALA/C/CA.ser
19HC~48~4440~61
19HC~48~4835~2268
1CZJ~48~959~786
2WIY~47~3124~2695
4EP6~49~3066~2668
Though that data might seem a little cryptic, it has a very simple meaning. The data is separated by tildes (I'm not sure why I chose to use tildes), and each portion the data means:
PDB ID of the structure where this interaction happened ~ (Distance (in Angstroms) of the interaction) * 10, rounded down ~ Serial number of the origin atom ~ Serial number of the target atom
For example, the first line, which is 19HC~48~4440~61, indicates that there was an interaction between atom #4440 and atom #61 in structure 19HC. The atoms are approximately 4.8 A away from each other.
In case you're wondering how to reproduce these calculated numbers, here's an example of some method that gets the serialized form of a distance result for an interaction, given the two atoms in the interaction and the structure this interaction is occuring in.
public static String getSerializedForm(Structure structure, Atom origin, Atom target) {
String pdbID = structure.getPDBCode();
int distance = (int) (Calc.getDistance(origin, target) * 10);
String originSerial = origin.getPDBserial();
String targetSerial = target.getPDBserial();
return pdbID + "~" + distance + "~" + originSerial + "~" + targetSerial;
}Note that the directory structure does not index by distance; as I've already explained, the information about the distance at which an interaction occured is stored in the DistanceResult. Thus, if the user has specified a particular distance range they want their interactions to fall into, then the directory structure searcher must open up the file containing the distance results, parse each of them, and then make sure their distances are between the user's specified minimum and maximum distances.
This is where things start to get a little hairy, unfortunately. When I say "integrating by ChemAxon", I mean getting this search tool to work with ChemAxon's MarvinView (an applet for drawing molecules), and getting it to work with the substructure search tool they (ChemAxon) provide. The first part of this section, getting it to work with MarvinSketch, is not something I've found an elegant way of doing. I've attached to this gist the code I've been using to try to extract a distance query from a drawing of a molecule, but as you can see it's not at all a very simple thing to do. As complicated as the specification might be, the source code pretty much speaks for itself here, so I won't talk about this too much. Here are a few things to know about ChemAxon's MarvinSketch if you want to try to modify / use that code:
- MarvinSketch exports queries into a type of XML called 'MRV' files. These are human-readable and it helps to look at those files directly when trying to understand what is going on in a query.
- When you "add data" in the MarvinSketch GUI, behind the scenes this creates a
DataSgroupin the exported 'MRV' file. I look at these quite a bit in my code. - It's possible to instantly get the MRV version of your drawing if you want that -- just go to 'File > Save As ...', and save it as an 'MRV' file.
- The final product of this tool would most likely need to have amino acids that have their atom names built-in. To do this, you must first create those molecules in marvin sketch, then save all of your amino acids in an MRV file. Make sure this MRV file is located in a directory where users can simply deliver an XmlHttpRequest to get that data from client-side JavaScript. Once you've put your file in place, then go to chemApplet.jsp, and where there are a bunch of calls to
msketch_paramyou can add a line that goes as
msketch_param("tmpls0", ":My Amino Acids:path/to/your/amino/acids.mrv");(The directory you provide should be relative to the chemApplet.jsp file.) (The "0" in "tmpls0" means "put this in position 0"; you can change this to "1" or "2" or some other number; if you set it to be the same position value as the one the default amino acids are at, then you can replace the amino acids altogether.)
So once we have the query, the next step is to parse the query to extract the information the user provided. This is done in the two files attached to this gist. In particular, the MarvinQueryParser class takes care of extracting information about the target atom, as well as finding the origin atom as well as finding out what the range of the distance interaction is. Right now, these three things are found using these assumptions:
- All atoms without an associated residue are part of the ligand (origin).
- The distance bond will have the field name of 'distance'; the attached data will be text such that after the text is split along spaces, the first and last token represent the minimum and maximum distance, respectively.
- The origin atom will be the atom in the distance bond that is also a ligand.
We find the associated residue of the atoms using the MolAtom#getResidueType() method. In MRV files, this data looks like this:
<!-- This atom has a nonzero residue type. -->
<atom id="a7" elementType="O" residueType="Ser"
x2="5.7750000762939475" y2="1.386846438394310" />
<!-- This atom has a residue type of 0. -->
<atom id="a8" elementType="C"
x2="-3.575000047683716" y2="1.4300000071525574" />To find the distance property, we use ChemAxon's Molecule#getSgroupArray; an Sgroup takes many forms, including the DataSgroup, which is what we're interested in and is what's created when a user attaches data to something. I look for a data Sgroup that has the field name 'Distance' (capitalization is irrelevant) and then from this extract the minimum and maxmimum distance as well as figure out which atom is the origin atom based on which of the two atoms in the interaction has / hasn't an associated residue type.
Once we have extracted all the data that's possible to get from the query, the hard part begins: running the substructure search and figuring out atom names. This is taken care of in the DistanceQuerySearchTool. This program uses a lot of ChemAxon stuff that's a little weird (and by "weird" I mean "I don't entirely understand what's going on"), but it certainly works and is copied from the LigandDAO search tool in the PDB.
In the distance query search tool, we get the data from MarvinQueryParser, then we make a new JChemSearch that will search for the "origin" molecule (this molecule was aquired from MarvinQueryParser). This will return a list of substructure matches. However, we are in particular interested in finding precisely what atoms were matched in the substructure match. To do this, we must use the MolSearch class, which gives us detailed substructure match information for a given pair of molecules. The MolSearch will provide a MolSearch#findAll() method which returns an int[][].
If we let hits be the value returned from MolSearch#findAll(), then hits[i] represents the i-th substructure match between the two passed atoms (it's possible for the query atom to show up within the target atoms multiple times, and each of those gets its own hit). hits[i] provides information about converting indicies in your query to the corresponding indicies in the MolSearch's target. For example, if you have an atom whose index is '5' in your query, then hits[i][5] will give you the index of that atom in your target molecule.
Now that we have a way of finding out what atoms are being matched in a substructure search, the next step is to figure out the names of the atoms the user provided. This is actually a relatively complicated process that requires that we implement Morgan's algorithm.
Morgan's algorithm is a way to identify atoms based on their extended connectivity. If done right, Morgan's algorithm can also allow us to consider most atoms to be different, but in some cases, due to pi-centeredness, consider some atoms to be equivalent (for example, we want to consider the two oxygens in a carboxyl to be equivalent).
Before we can implement Morgan's algorithm, we need to create a graph from a molecule. This needs to be done twice; once for a ChemAxon molecule, and once for a Biojava molecule. Morgan's algorithm will allow us to go back and forth. Implementing the programs to create the graphs from the molecules is not very difficult to do and is done in the files ChemAxonGraphMaker and PDBGraphMaker in the PDBDistances repository.
We then implement Morgan's algorithm. This is done in the Morgan and GraphDriver classes. At a high level, Morgan's algorithm looks like this:
initialize atom values
while the number of unique atom values continues to increase:
for each atom:
atom.value += the sum of the values of its neighbors
Once this algorithm stops (because the number of unique values has stopped increasing), then the atoms have acquired their final values. Two atoms with the same values are equivalent.
Initializing the atom values is not very complicated, and this is where some chemical heuristics have to be made. An atom's initial value is calculated as:
Initial Value = Atomic Number * 10 + Pi-Centeredness
Pi-Centeredness is 1 if the atom is in a pi-center, or 0 if it isn't. We determine an atom to be in a pi center if it meets any of these criteria:
- It forms a bond with another atom, and that bond has an order greater than 1.
- It is electronegative (more on this later), and is single-bonded to an atom that is in turn double-bonded to an electronegative atom.
We determine an atom to be electronegative if its valence electron count is between 5 and 7, inclusive.
That's all there is to it! Now, let's look back at the Distance query pipeline...
You'll recall that at this point we have found our substructure hits and we have a way of getting which atoms matched in a substructure search. So next, we apply morgan's algorithm to both the PDB's version of the match (which will also have atom names) and the ChemAxon's version of the match, and use the atom values to convert MolAtoms to atom names. The code sort of speaks for itself in this case. Once we have the origin atom name(s), we go through each possible name and create a distance query for each. Then, we just run the distance query we've created, and presto magico! the program works.
Attached are the two files that are not part of PDBDistances. These would go in PDB code or something.
As I've mentioned before, the code for setting template libraries goes like this:
msketch_param("tmpls0", ":My Amino Acids:path/to/your/amino/acids.mrv");But what do we put in that MRV file? Well, what you could do is modify the attached allAminoAcidsTemplate.mrv file, which contains an MRV file with all the amino acids (this MRV file is a valid MarvinSketch template library, and you can import it in MarvinSketch by going to "Insert > Template Library" and then clicking on the book with the green "plus" icon). Then, you could open these up in MarvinView and modify the atoms by right-clicking on an atom, going to "Edit Properties", then add a "Name" property with the appropriate name as a value. Then, once you've finished working on an amino acid, save that atom as a MRV file, and replace your template with the new version of the amino acids with the atom names attached.
I haven't actually tried this, so I imagine complications may arise.
Of course, the more proper way to do this would be to run Morgan's algorithm on the amino acid as well, and figure out atom names that way. As a bonus, this method will also identify equivalent atoms in an amino acid too.