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