View Javadoc
1   //******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  //******************************************************************************
38  package ffx.algorithms.commands;
39  
40  import ffx.algorithms.cli.AlgorithmsCommand;
41  import ffx.numerics.Potential;
42  import ffx.potential.ForceFieldEnergy;
43  import ffx.potential.MolecularAssembly;
44  import ffx.potential.bonded.Polymer;
45  import ffx.potential.bonded.Residue;
46  import ffx.potential.bonded.Rotamer;
47  import ffx.potential.bonded.RotamerLibrary;
48  import ffx.potential.parameters.ForceField;
49  import ffx.potential.parsers.ForceFieldFilter;
50  import ffx.potential.parsers.PDBFilter;
51  import ffx.utilities.FFXBinding;
52  import ffx.utilities.Keyword;
53  import org.apache.commons.configuration2.CompositeConfiguration;
54  import picocli.CommandLine.Command;
55  import picocli.CommandLine.Option;
56  import picocli.CommandLine.Parameters;
57  
58  import java.io.File;
59  import java.util.ArrayList;
60  import java.util.Arrays;
61  import java.util.Collections;
62  import java.util.List;
63  
64  import static ffx.utilities.StringUtils.parseAtomRanges;
65  
66  /**
67   * The MutatePDB script mutates residue(s) of a PDB file.
68   * <br>
69   * Usage:
70   * <br>
71   * ffxc MutatePDB [options] &lt;pdb&gt;
72   */
73  @Command(description = " Mutate PDB residue(s).", name = "MutatePDB")
74  public class MutatePDB extends AlgorithmsCommand {
75  
76    /**
77     * -r or --resid Residue numbers.
78     */
79    @Option(names = {"--resid", "-r"}, paramLabel = "1", defaultValue = "1",
80        description = "Residue number(s).")
81    private String resIDs;
82  
83    /**
84     * -n or --resname New residue names.
85     */
86    @Option(names = {"--resname", "-n"}, paramLabel = "ALA", defaultValue = "ALA",
87        description = "New residue name(s).")
88    private String resNameString;
89  
90    /**
91     * -ch or --chain Single character chain name (default is ' '). If only one chain exists, that chain will be mutated.
92     */
93    @Option(names = {"--chain", "--ch"}, paramLabel = " ", defaultValue = " ",
94        description = "Single character chain name (default is ' ').")
95    private String chainString;
96  
97    /**
98     * -R or --rotamer Rotamer number to apply.
99     */
100   @Option(names = {"--rotamer", "-R"}, paramLabel = "-1", defaultValue = "-1",
101       description = "Rotamer number to apply.")
102   private int rotamer;
103 
104   /**
105    * --allChains  Mutate all copies of a chain in a multimeric protein.
106    */
107   @Option(names = {"--allChains"}, paramLabel = "false", defaultValue = "false",
108       description = "Mutate all copies of a chains in a multimeric protein.")
109   private boolean allChains;
110 
111   /**
112    * A PDB filename.
113    */
114   @Parameters(arity = "1", paramLabel = "file",
115       description = "A PDB input file.")
116   private String filename;
117 
118   private ForceFieldEnergy forceFieldEnergy;
119 
120   /**
121    * MutatePDB Constructor.
122    */
123   public MutatePDB() {
124     super();
125   }
126 
127   /**
128    * MutatePDB Constructor.
129    *
130    * @param binding The Binding to use.
131    */
132   public MutatePDB(FFXBinding binding) {
133     super(binding);
134   }
135 
136   /**
137    * MutatePDB constructor that sets the command line arguments.
138    *
139    * @param args Command line arguments.
140    */
141   public MutatePDB(String[] args) {
142     super(args);
143   }
144 
145   /**
146    * Execute the script.
147    */
148   @Override
149   public MutatePDB run() {
150 
151     if (!init()) {
152       return this;
153     }
154 
155     // The "false" assembly provides access to the chainIDs without compromising the mutated molecular assembly.
156     // Used if --allChains is true.
157 
158     // Load the MolecularAssembly.
159     MolecularAssembly falseAssembly = getActiveAssembly(filename);
160     if (falseAssembly == null) {
161       logger.info(helpString());
162       return this;
163     }
164 
165     List<Integer> resIDint = parseAtomRanges("residueID", resIDs, falseAssembly.getAtomList().size());
166 
167     String[] resNameArr = Arrays.stream(resNameString.split("\\.|,|;")).map(String::trim).toArray(String[]::new);
168 
169     String[] chainStringArr = Arrays.stream(chainString.split("\\.|,|;")).map(String::trim).toArray(String[]::new);
170     String s = "";
171     for (String sub : chainStringArr) {
172       if (sub.length() >= 2) {
173         logger.warning("Chain ID's have to be single characters separated by commas!");
174         return this;
175       } else {
176         s += sub;
177       }
178     }
179     char[] chainArr = s.toCharArray();
180 
181     // For every chain, mutate the residue.
182     Polymer[] chains = falseAssembly.getChains();
183 
184     if (chains.length == 1 && chainArr.length == 0) {
185       chainArr = new char[]{chains[0].getChainID()};
186     }
187 
188     if (resIDint.size() != resNameArr.length) { // || resIDint.size() != chainArr.length) {
189       logger.warning("The number of chains, residue names, and residue ids must be the same!");
190       return this;
191     }
192 
193     int destRotamer = 0;
194     if (rotamer > -1) {
195       destRotamer = rotamer;
196     }
197 
198     // Read in command line.
199     File structureFile = new File(filename);
200     int index = filename.lastIndexOf(".");
201     String name = filename.substring(0, index);
202     MolecularAssembly molecularAssembly = new MolecularAssembly(name);
203     molecularAssembly.setFile(structureFile);
204 
205     CompositeConfiguration properties = Keyword.loadProperties(structureFile);
206     ForceFieldFilter forceFieldFilter = new ForceFieldFilter(properties);
207     ForceField forceField = forceFieldFilter.parse();
208     molecularAssembly.setForceField(forceField);
209 
210     PDBFilter pdbFilter = new PDBFilter(structureFile, molecularAssembly, forceField, properties);
211 
212     List<PDBFilter.Mutation> mutations = new ArrayList<>();
213     for (int i = 0; i <= resIDint.size() - 1; i++) {
214       // need to add one to get the correct resID because parseAtomRanges removes one
215       if (allChains) {
216         for (Polymer currentChain : chains) {
217           logger.info("\n Mutating residue number " + (resIDint.get(i) + 1) + " of chain " + currentChain.getChainID() + " to " + resNameArr[i]);
218           mutations.add(new PDBFilter.Mutation(resIDint.get(i) + 1, currentChain.getChainID(), resNameArr[i]));
219         }
220       } else {
221         logger.info("\n Mutating residue number " + (resIDint.get(i) + 1) + " of chain " + chainArr[i] + " to " + resNameArr[i]);
222         mutations.add(new PDBFilter.Mutation(resIDint.get(i) + 1, chainArr[i], resNameArr[i]));
223       }
224     }
225     pdbFilter.mutate(mutations);
226 
227     pdbFilter.readFile();
228     pdbFilter.applyAtomProperties();
229     molecularAssembly.finalize(true, forceField);
230 
231     if (destRotamer > -1) {
232       if (allChains) {
233         chains = molecularAssembly.getChains();
234         for (Polymer currentChain : chains) {
235           for (int i = 0; i <= resIDint.size() - 1; i++) {
236             Residue residue = currentChain.getResidue(resIDint.get(i) + 1);
237             Rotamer[] rotamers = residue.getRotamers();
238             if (rotamers != null && rotamers.length > 0) {
239               RotamerLibrary.applyRotamer(residue, rotamers[destRotamer]);
240             } else {
241               logger.info(" No rotamer to apply.");
242             }
243           }
244         }
245       } else {
246         for (int i = 0; i <= resIDint.size() - 1; i++) {
247           Polymer polymer = molecularAssembly.getChain(String.valueOf(chainArr[i]));
248           Residue residue = polymer.getResidue(resIDint.get(i) + 1);
249           Rotamer[] rotamers = residue.getRotamers();
250           if (rotamers != null && rotamers.length > destRotamer) {
251             RotamerLibrary.applyRotamer(residue, rotamers[destRotamer]);
252           } else {
253             logger.info(" No rotamer to apply.");
254           }
255         }
256       }
257     }
258     pdbFilter.writeFile(structureFile, false);
259 
260     for (PDBFilter.Mutation mut : mutations) {
261       mut.calculateTorsion();
262     }
263 
264     forceFieldEnergy = molecularAssembly.getPotentialEnergy();
265 
266     return this;
267   }
268 
269   /**
270    * {@inheritDoc}
271    */
272   @Override
273   public List<Potential> getPotentials() {
274     List<Potential> potentials;
275     if (forceFieldEnergy == null) {
276       potentials = Collections.emptyList();
277     } else {
278       potentials = Collections.singletonList(forceFieldEnergy);
279     }
280     return potentials;
281   }
282 
283 }