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.potential.commands;
39
40 import ffx.potential.bonded.Atom;
41 import ffx.potential.bonded.Residue;
42 import ffx.potential.bonded.Residue.ResidueType;
43 import ffx.potential.cli.PotentialCommand;
44 import ffx.potential.utils.PotentialsUtils;
45 import ffx.utilities.FFXBinding;
46 import picocli.CommandLine.Command;
47 import picocli.CommandLine.Parameters;
48
49 import java.io.File;
50 import java.util.ArrayList;
51 import java.util.List;
52
53 import static ffx.numerics.math.DoubleMath.dist;
54 import static java.lang.String.format;
55 import static org.apache.commons.io.FilenameUtils.getExtension;
56 import static org.apache.commons.io.FilenameUtils.getName;
57 import static org.apache.commons.io.FilenameUtils.removeExtension;
58 import static org.apache.commons.math3.util.FastMath.abs;
59
60
61
62
63
64
65
66 @Command(name = "ChainBreaks", description = " Fix chain breaks in a pdb file.")
67 public class ChainBreaks extends PotentialCommand {
68
69
70 @Parameters(arity = "1", paramLabel = "file",
71 description = "The atomic coordinate file in PDB or XYZ format.")
72 private String filename = null;
73
74 private List<String> chainBreaks = new ArrayList<>();
75 private List<double[]> newCoordinates;
76 private final PotentialsUtils potentialsUtils = new PotentialsUtils();
77
78 public ChainBreaks() {
79 super();
80 }
81
82 public ChainBreaks(FFXBinding binding) {
83 super(binding);
84 }
85
86 public ChainBreaks(String[] args) {
87 super(args);
88 }
89
90 @Override
91 public ChainBreaks run() {
92 if (!init()) {
93 return null;
94 }
95
96
97 activeAssembly = getActiveAssembly(filename);
98 if (activeAssembly == null) {
99 logger.info(helpString());
100 return null;
101 }
102
103
104 String name = getName(filename);
105 String ext = getExtension(name);
106 name = removeExtension(name);
107
108 if (!ext.toLowerCase().contains("pdb")) {
109 logger.info(format(" The file extension does not include 'pdb': %s", filename));
110 return null;
111 }
112
113 List<Residue> residues = activeAssembly.getResidueList();
114 chainBreaks = findChainBreaks(residues, 3.0);
115 logger.info(format(" Fixing Chain Breaks in %s", filename));
116
117
118 String dirString = getBaseDirString(filename);
119 File editedPDBFile = new File(dirString + name + "_edited.pdb");
120
121 logger.info(format(" Saving New Coordinates to:\n %s", editedPDBFile));
122 potentialsUtils.saveAsPDB(activeAssembly, editedPDBFile, false, false);
123
124 return this;
125 }
126
127 private List<String> findChainBreaks(List<Residue> residues, double cutoff) {
128 List<List<Residue>> subChains = new ArrayList<>();
129 List<String> chainBreaks = new ArrayList<>();
130
131
132
133 ResidueType rType = residues.get(0).getResidueType();
134 String startAtName;
135 String endAtName;
136 switch (rType) {
137 case AA:
138 startAtName = "N";
139 endAtName = "C";
140 break;
141 case NA:
142 boolean namedStar = residues.stream()
143 .flatMap((Residue r) -> r.getAtomList().stream())
144 .anyMatch((Atom a) -> a.getName().equals("O5*"));
145 if (namedStar) {
146 startAtName = "O5*";
147 endAtName = "O3*";
148 } else {
149 startAtName = "O5'";
150 endAtName = "O3'";
151 }
152 break;
153 case UNK:
154 default:
155 logger.fine(" Not attempting to find chain breaks for chain with residue " + residues.get(0));
156 return null;
157 }
158
159 List<Residue> subChain = null;
160 Residue previousResidue = null;
161 Atom priorEndAtom = null;
162
163 for (Residue residue : residues) {
164 List<Atom> resAtoms = residue.getAtomList();
165 if (priorEndAtom == null) {
166
167 subChain = new ArrayList<>();
168 subChain.add(residue);
169 subChains.add(subChain);
170 } else {
171
172 Atom startAtom = null;
173 for (Atom a : resAtoms) {
174 if (a.getName().equalsIgnoreCase(startAtName)) {
175 startAtom = a;
176 break;
177 }
178 }
179 if (startAtom == null) {
180 subChain.add(residue);
181 continue;
182 }
183
184 double r = dist(priorEndAtom.getXYZ(null), startAtom.getXYZ(null));
185
186 if (r > cutoff) {
187
188 subChain = new ArrayList<>();
189 subChain.add(residue);
190 subChains.add(subChain);
191 chainBreaks.add("C " + previousResidue + " N " + residue);
192 fixChainBreaks(priorEndAtom.getXYZ(null), startAtom.getXYZ(null));
193 priorEndAtom.setXYZ(newCoordinates.get(0));
194 startAtom.setXYZ(newCoordinates.get(1));
195 } else {
196
197 subChain.add(residue);
198 }
199 }
200
201
202 for (Atom a : resAtoms) {
203 if (a.getName().equalsIgnoreCase(endAtName)) {
204 priorEndAtom = a;
205 break;
206 }
207 }
208 previousResidue = residue;
209 }
210
211 return chainBreaks;
212 }
213
214 private void fixChainBreaks(double[] cCoordinates, double[] nCoordinates) {
215 logger.info(" Generating new coordinates.");
216 newCoordinates = new ArrayList<>();
217 double distance = dist(cCoordinates, nCoordinates);
218 while (distance > 3.0) {
219 double dx = abs((cCoordinates[0] - nCoordinates[0]) / 4.0);
220 double dy = abs((cCoordinates[1] - nCoordinates[1]) / 4.0);
221 double dz = abs((cCoordinates[2] - nCoordinates[2]) / 4.0);
222
223 if (cCoordinates[0] > nCoordinates[0]) {
224 cCoordinates[0] = cCoordinates[0] - dx;
225 nCoordinates[0] = nCoordinates[0] + dx;
226 } else if (cCoordinates[0] < nCoordinates[0]) {
227 cCoordinates[0] = cCoordinates[0] + dx;
228 nCoordinates[0] = nCoordinates[0] - dx;
229 }
230
231 if (cCoordinates[1] > nCoordinates[1]) {
232 cCoordinates[1] = cCoordinates[1] - dy;
233 nCoordinates[1] = nCoordinates[1] + dy;
234 } else if (cCoordinates[1] < nCoordinates[1]) {
235 cCoordinates[1] = cCoordinates[1] + dy;
236 nCoordinates[1] = nCoordinates[1] - dy;
237 }
238
239 if (cCoordinates[2] > nCoordinates[2]) {
240 cCoordinates[2] = cCoordinates[2] - dz;
241 nCoordinates[2] = nCoordinates[2] + dz;
242 } else if (cCoordinates[2] < nCoordinates[0]) {
243 cCoordinates[2] = cCoordinates[2] + dz;
244 nCoordinates[2] = nCoordinates[2] - dz;
245 }
246
247 distance = dist(cCoordinates, nCoordinates);
248 }
249 newCoordinates.add(cCoordinates);
250 newCoordinates.add(nCoordinates);
251 }
252 }