Residue mutations

So far, perturbations have involved changing an entire molecule. However, we can also change individual residues in a molecule. This can be useful to, e.g. study the effect of residue mutation on a protein, and how this could impact ligand binding or protein folding.

Sub-structure matching

To start, we need to define the sub-structure that we want to change.

Let’s first load the kigaki protein system.

>>> import sire as sr
>>> mols = sr.load_test_files("kigaki.gro", "kigaki.top")
>>> protein = mols[0]
>>> print(protein)
Molecule( Protein:6 num_atoms=302 num_residues=19 )
>>> print(protein.residues().names())
[ResName('LYS'), ResName('ILE'), ResName('GLY'), ResName('ALA'),
 ResName('LYS'), ResName('ILE'), ResName('LYS'), ResName('ILE'),
 ResName('GLY'), ResName('ALA'), ResName('LYS'), ResName('ILE'),
 ResName('LYS'), ResName('ILE'), ResName('GLY'), ResName('ALA'),
 ResName('LYS'), ResName('ILE'), ResName('NH2')]

To start with, let’s mutate the first alanine residue that we find in the protein into a lysine. First, let’s select the first alanine…

>>> ala = protein["resname ALA"][0]
>>> print(ala)
Residue( ALA:4   num_atoms=10 )

To mutate this into a lysine, we need to tell sire what a lysine looks like. Fortunately, our protein already contains lysine residues, so let’s select on of those…

>>> lys = protein["resname LYS"][1]
>>> print(lys)
Residue( LYS:5   num_atoms=22 )

Note

We have selected the second lysine residue in the protein, as the first lysine residue is the N-terminal residue. Since our alanine is a mid-chain residue, it makes more sense to use a mid-chain lysine.

Next, we need to match the alanine residue to the lysine residue, just as we did when perturbing entire molecules. Calling sire.morph.match() on sub-views of molecules will only match the atoms in those views.

>>> mapping = sr.morph.match(ala, lys, match_light_atoms=True)
>>> print(mapping)
AtomMapping( size=9, unmapped0=1, unmapped1=13
0: MolNum(6) Atom( HA:54 ) <=> MolNum(6) Atom( HA:64 )
1: MolNum(6) Atom( CA:53 ) <=> MolNum(6) Atom( CA:63 )
2: MolNum(6) Atom( O:60 ) <=> MolNum(6) Atom( O:82 )
3: MolNum(6) Atom( C:59 ) <=> MolNum(6) Atom( C:81 )
4: MolNum(6) Atom( HB2:57 ) <=> MolNum(6) Atom( HB2:67 )
5: MolNum(6) Atom( H:52 ) <=> MolNum(6) Atom( H:62 )
6: MolNum(6) Atom( N:51 ) <=> MolNum(6) Atom( N:61 )
7: MolNum(6) Atom( HB1:56 ) <=> MolNum(6) Atom( HB1:66 )
8: MolNum(6) Atom( CB:55 ) <=> MolNum(6) Atom( CB:65 )
)

Warning

This used the default MCS matching algorithm built into sire. This is not supported on Windows. You may want to use a different matching algorithm, e.g. using Kartograf

We can see that this mapping has nicely matched the atoms in alanine to the corresponding atoms in lysine.

Next, we will align the lysine in the mapping onto the alanine.

>>> mapping = mapping.align()

And finally, we will create the merged molecule.

>>> merged = mapping.merge(as_new_molecule=False)
>>> print(merged)
Molecule( Protein:6 num_atoms=315 num_residues=19 )

This has created the merged molecule as before. To see what has changed, we can check the perturbation…

>>> p_omm = merged.perturbation().to_openmm()
>>> print(p_omm.changed_atoms())
       atom  charge0  charge1    sigma0    sigma1  epsilon0  epsilon1  alpha0  alpha1  kappa0  kappa1
0      N:51  -0.4157  -0.3479  0.325000  0.325000  0.711280  0.711280     0.0     0.0     0.0     0.0
1      H:52   0.2719   0.2747  0.106908  0.106908  0.065689  0.065689     0.0     0.0     0.0     0.0
2     CA:53   0.0337  -0.2400  0.339967  0.339967  0.457730  0.457730     0.0     0.0     0.0     0.0
3     HA:54   0.0823   0.1426  0.247135  0.247135  0.065689  0.065689     0.0     0.0     0.0     0.0
4     CB:55  -0.1825  -0.0094  0.339967  0.339967  0.457730  0.457730     0.0     0.0     0.0     0.0
5    HB1:56   0.0603   0.0362  0.264953  0.264953  0.065689  0.065689     0.0     0.0     0.0     0.0
6    HB2:57   0.0603   0.0362  0.264953  0.264953  0.065689  0.065689     0.0     0.0     0.0     0.0
7    HB3:58   0.0603   0.0000  0.264953  0.264953  0.065689  0.000000     0.0     1.0     1.0     1.0
8      C:59   0.5973   0.7341  0.339967  0.339967  0.359824  0.359824     0.0     0.0     0.0     0.0
9      O:60  -0.5679  -0.5894  0.295992  0.295992  0.878640  0.878640     0.0     0.0     0.0     0.0
10  Xxx:303   0.0000   0.0187  0.339967  0.339967  0.000000  0.457730     1.0     0.0     1.0     1.0
11  Xxx:304   0.0000   0.0103  0.264953  0.264953  0.000000  0.065689     1.0     0.0     1.0     1.0
12  Xxx:305   0.0000   0.0103  0.264953  0.264953  0.000000  0.065689     1.0     0.0     1.0     1.0
13  Xxx:306   0.0000  -0.0479  0.339967  0.339967  0.000000  0.457730     1.0     0.0     1.0     1.0
14  Xxx:307   0.0000   0.0621  0.264953  0.264953  0.000000  0.065689     1.0     0.0     1.0     1.0
15  Xxx:308   0.0000   0.0621  0.264953  0.264953  0.000000  0.065689     1.0     0.0     1.0     1.0
16  Xxx:309   0.0000  -0.0143  0.339967  0.339967  0.000000  0.457730     1.0     0.0     1.0     1.0
17  Xxx:310   0.0000   0.1135  0.195998  0.195998  0.000000  0.065689     1.0     0.0     1.0     1.0
18  Xxx:311   0.0000   0.1135  0.195998  0.195998  0.000000  0.065689     1.0     0.0     1.0     1.0
19  Xxx:312   0.0000  -0.3854  0.325000  0.325000  0.000000  0.711280     1.0     0.0     1.0     1.0
20  Xxx:313   0.0000   0.3400  0.106908  0.106908  0.000000  0.065689     1.0     0.0     1.0     1.0
21  Xxx:314   0.0000   0.3400  0.106908  0.106908  0.000000  0.065689     1.0     0.0     1.0     1.0
22  Xxx:315   0.0000   0.3400  0.106908  0.106908  0.000000  0.065689     1.0     0.0     1.0     1.0

This shows how the atoms in alanine have been perturbed into the equivalent atoms in lysine, plus how a number of ghost atoms have been added that correspond to the additional atoms in lysine. In addition, the HB3 atom of alanine is turned into a ghost.

One-stop function

Just as before, the entire set of steps above can be performed in one step via the sire.morph.merge() function.

>>> merged = sr.morph.merge(ala, lys)
>>> print(merged)
Molecule( Protein:3627 num_atoms=315 num_residues=19 )
>>> p_omm = merged.perturbation().to_openmm()
>>> print(p_omm.changed_atoms())
       atom  charge0  charge1    sigma0    sigma1  epsilon0  epsilon1  alpha0  alpha1  kappa0  kappa1
0      N:51  -0.4157  -0.3479  0.325000  0.325000  0.711280  0.711280     0.0     0.0     0.0     0.0
1      H:52   0.2719   0.2747  0.106908  0.106908  0.065689  0.065689     0.0     0.0     0.0     0.0
2     CA:53   0.0337  -0.2400  0.339967  0.339967  0.457730  0.457730     0.0     0.0     0.0     0.0
3     HA:54   0.0823   0.1426  0.247135  0.247135  0.065689  0.065689     0.0     0.0     0.0     0.0
4     CB:55  -0.1825  -0.0094  0.339967  0.339967  0.457730  0.457730     0.0     0.0     0.0     0.0
5    HB1:56   0.0603   0.0362  0.264953  0.264953  0.065689  0.065689     0.0     0.0     0.0     0.0
6    HB2:57   0.0603   0.0362  0.264953  0.264953  0.065689  0.065689     0.0     0.0     0.0     0.0
7    HB3:58   0.0603   0.0000  0.264953  0.264953  0.065689  0.000000     0.0     1.0     1.0     1.0
8      C:59   0.5973   0.7341  0.339967  0.339967  0.359824  0.359824     0.0     0.0     0.0     0.0
9      O:60  -0.5679  -0.5894  0.295992  0.295992  0.878640  0.878640     0.0     0.0     0.0     0.0
10  Xxx:303   0.0000   0.0187  0.339967  0.339967  0.000000  0.457730     1.0     0.0     1.0     1.0
11  Xxx:304   0.0000   0.0103  0.264953  0.264953  0.000000  0.065689     1.0     0.0     1.0     1.0
12  Xxx:305   0.0000   0.0103  0.264953  0.264953  0.000000  0.065689     1.0     0.0     1.0     1.0
13  Xxx:306   0.0000  -0.0479  0.339967  0.339967  0.000000  0.457730     1.0     0.0     1.0     1.0
14  Xxx:307   0.0000   0.0621  0.264953  0.264953  0.000000  0.065689     1.0     0.0     1.0     1.0
15  Xxx:308   0.0000   0.0621  0.264953  0.264953  0.000000  0.065689     1.0     0.0     1.0     1.0
16  Xxx:309   0.0000  -0.0143  0.339967  0.339967  0.000000  0.457730     1.0     0.0     1.0     1.0
17  Xxx:310   0.0000   0.1135  0.195998  0.195998  0.000000  0.065689     1.0     0.0     1.0     1.0
18  Xxx:311   0.0000   0.1135  0.195998  0.195998  0.000000  0.065689     1.0     0.0     1.0     1.0
19  Xxx:312   0.0000  -0.3854  0.325000  0.325000  0.000000  0.711280     1.0     0.0     1.0     1.0
20  Xxx:313   0.0000   0.3400  0.106908  0.106908  0.000000  0.065689     1.0     0.0     1.0     1.0
21  Xxx:314   0.0000   0.3400  0.106908  0.106908  0.000000  0.065689     1.0     0.0     1.0     1.0
22  Xxx:315   0.0000   0.3400  0.106908  0.106908  0.000000  0.065689     1.0     0.0     1.0     1.0

This merged molecule can be used in a free energy simulation in the same way as any other merged molecule.

Protein Mutation

Just as for other merged molecules, we can extract the end states using the sire.morph.extract_reference() and sire.morph.extract_perturbed() functions.

>>> ref_prot = sr.morph.extract_reference(merged)
>>> pert_prot = sr.morph.extract_perturbed(merged)
>>> print(ref_prot.residues().names())
[ResName('LYS'), ResName('ILE'), ResName('GLY'), ResName('ALA'),
 ResName('LYS'), ResName('ILE'), ResName('LYS'), ResName('ILE'),
 ResName('GLY'), ResName('ALA'), ResName('LYS'), ResName('ILE'),
 ResName('LYS'), ResName('ILE'), ResName('GLY'), ResName('ALA'),
 ResName('LYS'), ResName('ILE'), ResName('NH2')]
>>> print(pert_prot.residues().names())
[ResName('LYS'), ResName('ILE'), ResName('GLY'), ResName('LYS'),
 ResName('LYS'), ResName('ILE'), ResName('LYS'), ResName('ILE'),
 ResName('GLY'), ResName('ALA'), ResName('LYS'), ResName('ILE'),
 ResName('LYS'), ResName('ILE'), ResName('GLY'), ResName('ALA'),
 ResName('LYS'), ResName('ILE'), ResName('NH2')]

This shows that the first alanine residue has been mutated into a lysine residue.

Interestingly, if you are only interested in mutation, and are not interested in the merged molecule, then mutation is just the process of performing a merge, and then extracting the perturbed end state.

To make this easier, there is a one-stop function for this, sire.morph.mutate().

>>> mutated_protein = sr.morph.mutate(ala, lys)
>>> print(mutated_protein.residues().names())
[ResName('LYS'), ResName('ILE'), ResName('GLY'), ResName('LYS'),
 ResName('LYS'), ResName('ILE'), ResName('LYS'), ResName('ILE'),
 ResName('GLY'), ResName('ALA'), ResName('LYS'), ResName('ILE'),
 ResName('LYS'), ResName('ILE'), ResName('GLY'), ResName('ALA'),
 ResName('LYS'), ResName('ILE'), ResName('NH2')]

Warning

By default this uses the (potentially slow) internal MCS matching algorithm (which also is not supported on Windows). You may want to use a different matching algorithm, either by passing in the dictionary of the mapping you want to use via the match argument, or by using Kartograf. The match argument works identically in this function as it does in the merge and match functions.

In this case, we copied and pasted the lysine residue from one part of the protein over the alanine. But, you can copy and paste in this way between different molecules. For example, you could have a template library of different residues that you could use for mutation.

Assuming template["lys"] contained your template for a lysine residue, then you could have run the following;

>>> mutated_protein = sr.morph.mutate(ala, template["lys"])

This works for any kind of molecules - not just proteins. In future versions of sire we will add functionality to make it easier to manage libraries of templates, and to perform molecular editing by copying and pasting between molecules, and between fragments constructed via, e.g. smiles strings.