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.crystal.Crystal;
41 import ffx.numerics.math.DoubleMath;
42 import ffx.potential.ForceFieldEnergy;
43 import ffx.potential.bonded.Atom;
44 import ffx.potential.bonded.Residue;
45 import ffx.potential.cli.PotentialCommand;
46 import ffx.potential.parsers.PDBFilter;
47 import ffx.potential.parsers.SystemFilter;
48 import ffx.potential.parsers.XYZFilter;
49 import ffx.utilities.FFXBinding;
50 import picocli.CommandLine.Command;
51 import picocli.CommandLine.Option;
52 import picocli.CommandLine.Parameters;
53
54 import java.util.Arrays;
55 import java.util.List;
56 import java.util.stream.Collectors;
57
58 import static java.lang.String.format;
59 import static org.apache.commons.math3.util.FastMath.toDegrees;
60
61
62
63
64
65
66
67 @Command(description = "Calculates nucleic acid torsions as well as information regarding the sugar pucker.", name = "NucleicAcidAnalysis")
68 public class NucleicAcidAnalysis extends PotentialCommand {
69
70 @Option(names = {"--at", "--allTorsions"}, paramLabel = "false", defaultValue = "false",
71 description = "Print all torsions and information.")
72 private boolean allTorsions = false;
73
74 @Parameters(arity = "1", paramLabel = "file",
75 description = "The atomic coordinate file in PDB or XYZ or ARC format.")
76 private String filename = null;
77
78 private List<Residue> residues;
79
80 public NucleicAcidAnalysis() {
81 super();
82 }
83
84 public NucleicAcidAnalysis(FFXBinding binding) {
85 super(binding);
86 }
87
88 public NucleicAcidAnalysis(String[] args) {
89 super(args);
90 }
91
92 @Override
93 public NucleicAcidAnalysis run() {
94
95 if (!init()) {
96 return this;
97 }
98
99 activeAssembly = getActiveAssembly(filename);
100 if (activeAssembly == null) {
101 logger.info(helpString());
102 return this;
103 }
104
105 ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy();
106 double[] x = new double[forceFieldEnergy.getNumberOfVariables()];
107 forceFieldEnergy.getCoordinates(x);
108 forceFieldEnergy.energy(x);
109
110
111 List<String> nucleicAcidNames = Arrays.asList(
112 "DA", "DC", "DG", "DT", "DU", "A", "C", "G", "T", "U",
113 "DAD", "DCY", "DGU", "DTY", "URA");
114
115 residues = activeAssembly.getResidueList().stream()
116 .filter(residue -> nucleicAcidNames.contains(residue.getName().toUpperCase()))
117 .collect(Collectors.toList());
118
119 SystemFilter systemFilter = potentialFunctions.getFilter();
120 if (systemFilter instanceof XYZFilter || systemFilter instanceof PDBFilter) {
121 int numModels = systemFilter.countNumModels();
122
123 if (numModels == 1 && allTorsions) {
124
125 printHeader();
126 analyzeAndPrintResidues();
127 } else {
128
129 int frameNumber = 1;
130 logger.info(format("\nFrame %d:", frameNumber));
131 logger.info("Residue Name χ P Sugar pucker");
132 logger.info("-----------------------------------------------------");
133
134 for (Residue residue : residues) {
135 Double chiTemp = getDihedral(residue, "O4'", "C1'", "N9", "C4");
136 double chi = (chiTemp != null) ? chiTemp : getDihedral(residue, "O4'", "C1'", "N1", "C2");
137 double P = calculateP(residue);
138 String sugarPucker = determineSugarPucker(P);
139
140 logger.info(format("%-10s %-8s %-10s %-10s %-14s",
141 residue.getResidueNumber(),
142 residue.getName(),
143 formatValue(chi),
144 formatValue(P),
145 sugarPucker
146 ));
147 }
148 frameNumber++;
149 while (systemFilter.readNext()) {
150 Crystal crystal = activeAssembly.getCrystal();
151 forceFieldEnergy.setCrystal(crystal);
152 forceFieldEnergy.getCoordinates(x);
153 forceFieldEnergy.energy(x);
154
155
156 residues = activeAssembly.getResidueList().stream()
157 .filter(residue -> nucleicAcidNames.contains(residue.getName().toUpperCase()))
158 .collect(Collectors.toList());
159
160 logger.info(format("\nFrame %d:", frameNumber));
161 logger.info("Residue Name χ P Sugar pucker");
162 logger.info("-----------------------------------------------------");
163
164 for (Residue residue : residues) {
165 Double chiTemp = getDihedral(residue, "O4'", "C1'", "N9", "C4");
166 double chi = (chiTemp != null) ? chiTemp : getDihedral(residue, "O4'", "C1'", "N1", "C2");
167 double P = calculateP(residue);
168 String sugarPucker = determineSugarPucker(P);
169
170 logger.info(format("%-10s %-8s %-10s %-10s %-14s",
171 residue.getResidueNumber(),
172 residue.getName(),
173 formatValue(chi),
174 formatValue(P),
175 sugarPucker
176 ));
177 }
178 frameNumber++;
179 }
180 }
181 }
182
183 return this;
184 }
185
186 private void printHeader() {
187 logger.info("Residue Name V0 V1 V2 V3 V4 P νmax χ γ α β δ ε ζ Sugar pucker Stage");
188 logger.info("------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------");
189 }
190
191 private void analyzeAndPrintResidues() {
192 for (Residue residue : residues) {
193 Double v0 = getDihedral(residue, "C2'", "C1'", "O4'", "C4'");
194 Double v1 = getDihedral(residue, "O4'", "C1'", "C2'", "C3'");
195 Double v2 = getDihedral(residue, "C1'", "C2'", "C3'", "C4'");
196 Double v3 = getDihedral(residue, "O4'", "C4'", "C3'", "C2'");
197 Double v4 = getDihedral(residue, "C1'", "O4'", "C4'", "C3'");
198 Double chiTemp = getDihedral(residue, "O4'", "C1'", "N9", "C4");
199 Double chi = (chiTemp != null) ? chiTemp : getDihedral(residue, "O4'", "C1'", "N1", "C2");
200 Double gamma = getDihedral(residue, "O5'", "C5'", "C4'", "C3'");
201 Double alpha = getDihedral(residue, "O3'", "P", "O5'", "C5'");
202 Double beta = getDihedral(residue, "P", "O5'", "C5'", "C4'");
203 Double delta = getDihedral(residue, "C5'", "C4'", "C3'", "O3'");
204 Double epsilon = getDihedral(residue, "C4'", "C3'", "O3'", "P");
205 Double zeta = getDihedral(residue, "C3'", "O3'", "P", "O5'");
206
207 Double P = null;
208 if (v0 != null && v1 != null && v3 != null && v4 != null && v2 != null) {
209 P = calculateP(v0, v1, v2, v3, v4);
210 }
211
212 Double nuMax = null;
213 if (v2 != null && P != null) {
214 nuMax = Math.abs(v2 / Math.cos(Math.toRadians(P)));
215 }
216
217 String type = determineType(residue, P);
218 String stage = determineStage(delta, chi, P);
219
220 logger.info(format("%-10s %-8s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-10s %-14s %-18s",
221 residue.getResidueNumber(),
222 residue.getName(),
223 formatValue(v0), formatValue(v1), formatValue(v2),
224 formatValue(v3), formatValue(v4),
225 formatValue(P), formatValue(nuMax),
226 formatValue(chi), formatValue(gamma),
227 formatValue(alpha), formatValue(beta),
228 formatValue(delta), formatValue(epsilon),
229 formatValue(zeta),
230 type,
231 stage
232 ));
233 }
234 }
235
236 private static Double getDihedral(Residue residue, String atom1, String atom2, String atom3, String atom4) {
237 try {
238 Atom a1 = residue.getAtomByName(atom1, true);
239 Atom a2 = residue.getAtomByName(atom2, true);
240 Atom a3 = residue.getAtomByName(atom3, true);
241 Atom a4 = residue.getAtomByName(atom4, true);
242
243 if (a1 != null && a2 != null && a3 != null && a4 != null) {
244 return toDegrees(DoubleMath.dihedralAngle(a1.getXYZ(null), a2.getXYZ(null), a3.getXYZ(null), a4.getXYZ(null)));
245 }
246 } catch (Exception e) {
247 logger.warning(format("Could not calculate dihedral for %s atoms %s-%s-%s-%s", residue, atom1, atom2, atom3, atom4));
248 }
249 return null;
250 }
251
252 private static String formatValue(Double value) {
253 return value != null ? format("%.2f", value) : "N/A";
254 }
255
256 private static double calculateP(Residue residue) {
257 Double v0 = getDihedral(residue, "C2'", "C1'", "O4'", "C4'");
258 Double v1 = getDihedral(residue, "O4'", "C1'", "C2'", "C3'");
259 Double v2 = getDihedral(residue, "C1'", "C2'", "C3'", "C4'");
260 Double v3 = getDihedral(residue, "O4'", "C4'", "C3'", "C2'");
261 Double v4 = getDihedral(residue, "C1'", "O4'", "C4'", "C3'");
262
263 if (v0 != null && v1 != null && v3 != null && v4 != null && v2 != null) {
264 return calculateP(v0, v1, v2, v3, v4);
265 }
266 return Double.NaN;
267 }
268
269 private static double calculateP(double v0, double v1, double v2, double v3, double v4) {
270 try {
271 double sin36 = Math.sin(Math.toRadians(36));
272 double sin72 = Math.sin(Math.toRadians(72));
273 double denominator = 2 * v2 * (sin36 + sin72);
274 if (Math.abs(denominator) < 1e-10) {
275 return Double.NaN;
276 }
277 double p = ((v4 - v0) - (v3 - v1)) / denominator;
278
279 double P;
280 if (v2 < 0) {
281 P = Math.toDegrees(Math.atan(p)) + 180.0;
282 } else if (p < 0) {
283 P = Math.toDegrees(Math.atan(p)) + 360.0;
284 } else {
285 P = Math.toDegrees(Math.atan(p));
286 }
287 return P;
288 } catch (Exception e) {
289 return Double.NaN;
290 }
291 }
292
293 private static String determineType(Residue residue, Double P) {
294 if (P == null) return "Unknown";
295
296 String base;
297 switch (residue.getName()) {
298 case "DAD":
299 base = "Ade";
300 break;
301 case "DGU":
302 base = "Gua";
303 break;
304 case "DCY":
305 base = "Cyt";
306 break;
307 case "DTY":
308 base = "Thy";
309 break;
310 case "URA":
311 base = "Ura";
312 break;
313 default:
314 base = "Unknown";
315 break;
316 }
317
318 String sugarPucker = determineSugarPucker(P);
319 return base + ", " + sugarPucker;
320 }
321
322 private static String determineSugarPucker(Double P) {
323 if (P == null) return "Unknown";
324
325 if (P >= 0 && P < 36) {
326 return "C3'-endo";
327 } else if (P >= 36 && P < 72) {
328 return "C4'-endo";
329 } else if (P >= 72 && P < 108) {
330 return "O4'-endo";
331 } else if (P >= 108 && P < 144) {
332 return "C1'-exo";
333 } else if (P >= 144 && P < 180) {
334 return "C2'-endo";
335 } else if (P >= 180 && P < 216) {
336 return "C3'-exo";
337 } else if (P >= 216 && P < 252) {
338 return "C4'-exo";
339 } else if (P >= 252 && P < 288) {
340 return "O4'-exo";
341 } else if (P >= 288 && P < 324) {
342 return "C1'-endo";
343 } else if (P >= 324 || P < 0) {
344 return "C2'-endo";
345 }
346 return "Unknown";
347 }
348
349 private static String determineStage(Double delta, Double chi, Double P) {
350 if (delta == null || chi == null || P == null) return "Unknown";
351
352 if (delta > 120.0 && chi > -160) {
353 return "Stage 1: Slide and Roll Change";
354 } else if (delta < 120.0 && P >= 0 && P < 180.0) {
355 return "Stage 2: Backbone and Sugar Pucker Change";
356 } else if (delta < 120.0 && P >= 180.0 && P < 360.0) {
357 return "Stage 3: Roll and Inclination to Finish Transition";
358 }
359 return "Unknown";
360 }
361 }