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 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
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
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 }