1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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
68
69
70
71
72
73 @Command(description = " Mutate PDB residue(s).", name = "MutatePDB")
74 public class MutatePDB extends AlgorithmsCommand {
75
76
77
78
79 @Option(names = {"--resid", "-r"}, paramLabel = "1", defaultValue = "1",
80 description = "Residue number(s).")
81 private String resIDs;
82
83
84
85
86 @Option(names = {"--resname", "-n"}, paramLabel = "ALA", defaultValue = "ALA",
87 description = "New residue name(s).")
88 private String resNameString;
89
90
91
92
93 @Option(names = {"--chain", "--ch"}, paramLabel = " ", defaultValue = " ",
94 description = "Single character chain name (default is ' ').")
95 private String chainString;
96
97
98
99
100 @Option(names = {"--rotamer", "-R"}, paramLabel = "-1", defaultValue = "-1",
101 description = "Rotamer number to apply.")
102 private int rotamer;
103
104
105
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
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
122
123 public MutatePDB() {
124 super();
125 }
126
127
128
129
130
131
132 public MutatePDB(FFXBinding binding) {
133 super(binding);
134 }
135
136
137
138
139
140
141 public MutatePDB(String[] args) {
142 super(args);
143 }
144
145
146
147
148 @Override
149 public MutatePDB run() {
150
151 if (!init()) {
152 return this;
153 }
154
155
156
157
158
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
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) {
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
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 String[] patches = properties.getStringArray("patch");
209 if (patches != null) {
210 for (String patch : patches) {
211 logger.info(" Attempting to read force field patch from " + patch + ".");
212 CompositeConfiguration patchConfiguration = new CompositeConfiguration();
213 patchConfiguration.addProperty("parameters", patch);
214 forceFieldFilter = new ForceFieldFilter(patchConfiguration);
215 ForceField patchForceField = forceFieldFilter.parse();
216 forceField.append(patchForceField);
217 }
218 }
219 molecularAssembly.setForceField(forceField);
220
221 PDBFilter pdbFilter = new PDBFilter(structureFile, molecularAssembly, forceField, properties);
222
223 List<PDBFilter.Mutation> mutations = new ArrayList<>();
224 for (int i = 0; i <= resIDint.size() - 1; i++) {
225
226 if (allChains) {
227 for (Polymer currentChain : chains) {
228 logger.info("\n Mutating residue number " + (resIDint.get(i) + 1) + " of chain " + currentChain.getChainID() + " to " + resNameArr[i]);
229 mutations.add(new PDBFilter.Mutation(resIDint.get(i) + 1, currentChain.getChainID(), resNameArr[i]));
230 }
231 } else {
232 logger.info("\n Mutating residue number " + (resIDint.get(i) + 1) + " of chain " + chainArr[i] + " to " + resNameArr[i]);
233 mutations.add(new PDBFilter.Mutation(resIDint.get(i) + 1, chainArr[i], resNameArr[i]));
234 }
235 }
236 pdbFilter.mutate(mutations);
237
238 pdbFilter.readFile();
239 pdbFilter.applyAtomProperties();
240 molecularAssembly.finalize(true, forceField);
241
242 if (destRotamer > -1) {
243 if (allChains) {
244 chains = molecularAssembly.getChains();
245 for (Polymer currentChain : chains) {
246 for (int i = 0; i <= resIDint.size() - 1; i++) {
247 Residue residue = currentChain.getResidue(resIDint.get(i) + 1);
248 Rotamer[] rotamers = residue.getRotamers();
249 if (rotamers != null && rotamers.length > 0) {
250 RotamerLibrary.applyRotamer(residue, rotamers[destRotamer]);
251 } else {
252 logger.info(" No rotamer to apply.");
253 }
254 }
255 }
256 } else {
257 for (int i = 0; i <= resIDint.size() - 1; i++) {
258 Polymer polymer = molecularAssembly.getChain(String.valueOf(chainArr[i]));
259 Residue residue = polymer.getResidue(resIDint.get(i) + 1);
260 Rotamer[] rotamers = residue.getRotamers();
261 if (rotamers != null && rotamers.length > destRotamer) {
262 RotamerLibrary.applyRotamer(residue, rotamers[destRotamer]);
263 } else {
264 logger.info(" No rotamer to apply.");
265 }
266 }
267 }
268 }
269 pdbFilter.writeFile(structureFile, false);
270
271 for (PDBFilter.Mutation mut : mutations) {
272 mut.calculateTorsion();
273 }
274
275 forceFieldEnergy = molecularAssembly.getPotentialEnergy();
276
277 return this;
278 }
279
280
281
282
283 @Override
284 public List<Potential> getPotentials() {
285 List<Potential> potentials;
286 if (forceFieldEnergy == null) {
287 potentials = Collections.emptyList();
288 } else {
289 potentials = Collections.singletonList(forceFieldEnergy);
290 }
291 return potentials;
292 }
293
294 }