How to Make a Point Mutation in a Protein Structure and Compare It with the Wild Type (PyMOL Mutagenesis + GROMACS)

The most common MSc project design in structural bioinformatics is also the one this blog has never taught end to end: take a structure, introduce a reported point mutation in silico, simulate the mutant and the wild type, and compare them. Every individual step is already covered here, from downloading the structure through pdb2gmx to RMSD and RMSF analysis. The mutation itself was the missing link. With more than 258,000 released entries in the RCSB Protein Data Bank as of August 2026, the wild-type coordinates you need almost certainly exist. The mutant almost certainly does not, so you have to build it.

Answer: To mutate a residue in PyMOL, load the structure, open Wizard > Mutagenesis, click the target residue, choose the new amino acid from the “No Mutation” dropdown, pick a rotamer with the frame controls, then click Apply and Done. The wizard swaps the side chain onto the existing backbone and does not relax the surroundings, so the mutant must be minimised and simulated before you claim anything.

What does the PyMOL Mutagenesis wizard actually do?

It replaces one side chain with another on the backbone that is already there, and it offers you a library of pre-computed side-chain conformations to choose from. That is the whole operation. The backbone atoms N, CA, C and O stay where the crystallographer or the structure predictor put them, and every neighbouring residue stays exactly where it was.

The rotamers are not random. According to the PyMOL Wiki Mutagenesis page, “The rotamers are ordered according to their frequencies of occurrence in proteins”, and the percentage for the currently displayed rotamer appears as a red label on the preview object. The wizard also lets you switch between backbone-dependent and backbone-independent rotamer sets from the panel, labelled “Backbone Depen. Rotamers” and “Backbone Indep. Rotamers”. Backbone-dependent is the default and is the one you want, because it conditions the side-chain conformations on the local phi and psi angles. If the residue has no phi or psi to read, for example at a chain terminus, PyMOL prints a message in the console saying it fell back to backbone-independent rotamers.

One detail almost no tutorial mentions: the wizard scores every rotamer for steric strain against its neighbours and jumps to the lowest-strain one for you. As you step through rotamers, the PyMOL console prints a line of the form Rotamer 3/9, strain=.... That strain number is a van der Waals overlap score used to rank conformers, not a folding free energy, and it must never be reported as one.

How do you mutate a residue in PyMOL, step by step?

The wizard ships with the open-source build. The mutagenesis wizard source file lives in the schrodinger/pymol-open-source repository, so students on the free or educational build get the same tool as licence holders. Nucleotide mutagenesis was added separately in PyMOL 2.2; if your Wizard menu shows a protein variant and a nucleic acid variant, pick the protein one.

1. Load the structure. Type fetch 1abc in the PyMOL command line, or use File > Open on a local PDB file. Use a structure you have already cleaned. Alternate conformations, missing side chains and unusual protonation all cause problems later, and our guide to preparing a protein for MD covers how to deal with them first.
2. Check the residue number against the file, not against the paper. Author numbering in a PDB entry frequently differs from UniProt numbering, and two chains can carry the same residue number. Run iterate 1abc and chain A and resi 120 and name CA, print(chain, resi, resn) and confirm the printed residue name is the one you expect to mutate away from. Skipping this check is the single most common way a mutation project silently produces the wrong molecule.
3. Open the wizard. Wizard > Mutagenesis. The wizard switches picking to residue level automatically, so a single click selects a whole residue.
4. Click the target residue in the 3D viewer. The panel prompt changes to ask you to select a rotamer or pick a new residue.
5. Choose the new amino acid from the dropdown labelled “No Mutation”. Once you pick one, that button reads “Mutate to ARG” or whichever residue you chose.
6. Cycle the rotamers with the frame arrows in the lower right of the window. Watch the contact discs. The PyMOL Wiki describes small green discs as near contacts and large red discs as significant van der Waals overlap. Prefer a rotamer with no large red discs, and keep an eye on the frequency percentage.
7. Click Apply, then Done. Apply writes the chosen side chain into the object you clicked. The temporary preview object called mutation disappears at that point, which is expected.

If you would rather have a reproducible script than a sequence of clicks, the same wizard is scriptable, which matters when a reviewer asks how the mutant was built:

fetch 1abc, async=0
wizard mutagenesis
cmd.get_wizard().set_mode("GLY")
cmd.get_wizard().do_select("/1abc//A/120")
cmd.get_wizard().apply()
set_wizard
save 1abc_R120G.pdb, 1abc

Save that as a .pml file next to your data. It is the honest answer to “how did you generate the mutant”, and it takes the ambiguity out of your methods section.

How do you save only the mutated structure?

This is where a surprising number of projects go wrong. Apply merges the new side chain into the object that contained the residue you clicked, so after Apply the mutation lives inside your original object. There is no separate mutant object waiting to be saved.

Save the object by name: save 1abc_R120G.pdb, 1abc. Three failure modes to avoid. Saving a .pse file gives you a PyMOL session, not a structure file, and pdb2gmx cannot read it. Running save mutant.pdb with no selection writes every enabled object, so if the wild type is still loaded you get both molecules in one file. And saving before clicking Apply writes the unmutated structure, because the preview object is not part of your molecule yet.

After saving, reopen the file and confirm the mutation is present with iterate 1abc_R120G and chain A and resi 120 and name CA, print(resn). Then run the mutant through structure validation. A Ramachandran and clash check on the mutant costs two minutes and catches a bad rotamer before it costs you a week of simulation time.

What can the Mutagenesis wizard not tell you?

Be clear about this in your thesis, because an examiner will ask. The wizard grafts a side chain onto a fixed backbone and selects a rotamer. It does not move the backbone. It does not repack the neighbouring side chains to accommodate the new one. It does not run any energy minimisation. It does not produce a folding free energy, and it cannot tell you whether a mutation is stabilising or destabilising.

That limitation is not a flaw, it is the reason the workflow continues. A mutant built this way is a starting geometry, nothing more. Energy minimisation relieves the local strain that the graft introduced, and the equilibration and production run let the backbone and the neighbours respond to the substitution. Any structural claim you make has to come from the simulation, not from the wizard.

Should you use PyMOL, a stability predictor, or rebuild the model?

The right tool depends on the question. If you want a starting structure for MD, the wizard is enough. If your research question is the stability change itself, you need a method designed to estimate it.

Rosetta sits alongside FoldX for this purpose and is free for academic users, at the cost of a heavier installation. Whichever predictor you use, report the tool version and treat the number as a prediction, not a measurement.

How do you set up a fair wild-type vs mutant comparison in GROMACS?

A comparison is only as good as its controls, and this is where most student projects lose marks. The rule is simple: everything except the mutation must be identical. Build both systems in the same session, with the same commands, changing only the input coordinate file.

Treat these as a fixed protocol applied twice. Use the same force field and water model in pdb2gmx. Use the same box type and the same solute-to-box distance in editconf. Solvate identically and add ions to the same concentration, accepting that the neutralising counter-ion count can differ by one if the mutation changes the net charge, which is expected and should be stated. Run the same energy minimisation and the same NVT and NPT equilibration. Use the same production length and the same .mdp settings, with a different random seed for each replicate.

The GROMACS manual describes gmx pdb2gmx as a tool that reads a coordinate file and database files, adds hydrogens, and generates GROMACS coordinates and a topology. That hydrogen-adding step is why you should hand pdb2gmx the mutant with the wizard’s hydrogens removed and let the force field name them, rather than importing PyMOL’s naming into your topology.

Then compare. RMSD of the backbone for global drift, per-residue RMSF for local flexibility with attention to the window around the mutation site, and radius of gyration, SASA and hydrogen bonds for compactness and exposure. Both analyses are already covered in detail in our RMSD and RMSF guide and our radius of gyration, SASA and hydrogen bond guide, so run the same commands on both trajectories rather than reinventing them. For the figure, superimpose the mutant on the wild type in PyMOL and show the two side chains side by side.

Three controls a reviewer will look for. Run the wild type yourself under your own protocol instead of reusing a published trajectory. Run at least two independent replicates per system, because a single trajectory cannot separate a real difference from normal sampling noise. And discard the equilibration window before comparing anything, using the criteria in our guide to how long to run an MD simulation.

What can you conclude from a wild-type vs mutant comparison?

You can say that under your force field, your protocol and your sampling, the mutant showed higher or lower backbone RMSD, that a specific loop region showed increased RMSF, that a particular hydrogen bond present in the wild type was absent or short-lived in the mutant, or that the mutant sampled a more expanded conformation by radius of gyration. Those are observations about your simulations, and they are defensible if the protocol was identical and the replicates agree.

You cannot say the mutation is destabilising, that it reduces binding affinity, or that it causes a disease phenotype. Unbiased MD does not measure folding free energy, and a few hundred nanoseconds of a folded protein does not sample unfolding. Free energy claims need a method built for them. If replicates disagree, report that honestly and treat it as insufficient sampling rather than picking the run that supports your hypothesis.

Troubleshooting: real errors and their fixes

Every rotamer shows large red discs and none look clean. The site is too tight for that residue on a fixed backbone, which is informative in itself. Take the lowest-strain rotamer, accept it as a starting geometry, and let energy minimisation and equilibration resolve the clash. If minimisation cannot converge, rebuild the mutant with a tool that repacks neighbours, such as FoldX, or rebuild the model from the mutant sequence.

You mutated the wrong residue because two chains share the number. Author numbering is not unique across chains. Before mutating, run the iterate check above with an explicit chain selector, and in the scripted route always use a full selection such as /1abc//A/120 rather than resi 120 alone.

The residue has alternate conformations in the crystal file. Mutating a residue that carries altloc A and altloc B leaves duplicated atoms that pdb2gmx will reject. Strip the alternates first with remove not (alt ''+A) followed by alter all, alt='', then mutate.

pdb2gmx fails on the mutated residue with an atom-not-found or hydrogen-naming error. PyMOL’s hydrogen names do not always match the force field’s residue topology entries. Add -ignh, documented in the GROMACS manual as “Ignore hydrogen atoms that are in the coordinate file”, so pdb2gmx strips the incoming hydrogens and adds its own with the correct names. If a heavy atom is genuinely missing, that is a preparation problem, not a naming problem, and -missing is not the fix. Our list of common GROMACS errors covers the rest.

The two systems cannot be compared because the boxes differ. If you built the wild type last week and the mutant today, the box dimensions, water count and ion count will not match, and SASA and radius of gyration comparisons become hard to defend. Rebuild both from scratch in one session using a single script with the input filename as the only variable.

The saved file does not contain the mutation. Confirm you clicked Apply before saving, that you saved a .pdb and not a .pse, and that you named the object in the save command. Reopen the saved file in a fresh PyMOL session and check the residue name directly rather than trusting the viewer you just worked in.

Frequently asked questions

Can you mutate a residue in the open-source PyMOL?

Yes. The mutagenesis wizard is part of the open-source PyMOL codebase and behaves the same way as in the licensed build. Students using the free or educational version get the same rotamer library, the same strain-based rotamer ranking, and the same Apply behaviour.

Does the PyMOL Mutagenesis wizard minimise the structure?

No. It places a side-chain rotamer on the existing backbone and stops there. No minimisation, no backbone movement, no repacking of neighbours. The strain score it prints ranks rotamers against each other and is not an energy you can report.

How many replicates do you need for a mutant vs wild-type comparison?

At least two independent runs per system, and three is a safer default for a thesis. Independent means a different random seed, ideally from a separately equilibrated starting point. A difference that appears in one replicate and not the other is sampling noise, not a finding.

Should you mutate before or after preparing the structure for MD?

Prepare first, mutate second. Fix missing residues, remove alternate conformations, decide protonation states and strip unwanted heteroatoms on the wild type, then mutate the cleaned file. That way the wild type and mutant differ only in the substituted side chain.

Can MD tell you whether a mutation is destabilising?

Not on its own. Standard unbiased MD does not compute folding free energy and rarely samples unfolding on accessible timescales. Use a dedicated predictor such as FoldX, Rosetta or DynaMut2 for a ΔΔG estimate, and use MD to describe the structural and dynamic consequences.

Putting it together

The mutation itself takes about a minute in PyMOL. Everything that makes the result publishable happens around it: checking the residue number against the file, saving the object rather than the session, validating the mutant, and then running the mutant and the wild type through one identical GROMACS protocol with replicates. State plainly in your methods that the mutant was built by side-chain replacement in PyMOL with no backbone relaxation, and let the simulation carry the argument.

If you are planning this as a thesis project, the computational biology skills roadmap shows where structure preparation, mutation and MD analysis sit relative to everything else you will need.

Written by the StemSkills Lab team, who bring more than a decade of combined work in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling.