1 |
import gnu.getopt.Getopt; |
2 |
|
3 |
import java.io.File; |
4 |
import java.io.FileNotFoundException; |
5 |
import java.io.FileOutputStream; |
6 |
import java.io.IOException; |
7 |
import java.io.PrintWriter; |
8 |
import java.sql.SQLException; |
9 |
import java.util.Formatter; |
10 |
import java.util.Locale; |
11 |
|
12 |
import proteinstructure.AAinfo; |
13 |
import proteinstructure.ConformationsNotSameSizeError; |
14 |
import proteinstructure.Graph; |
15 |
import proteinstructure.Pdb; |
16 |
import proteinstructure.PdbChainCodeNotFoundError; |
17 |
import proteinstructure.PdbCodeNotFoundError; |
18 |
import proteinstructure.PdbaseInconsistencyError; |
19 |
import proteinstructure.PdbasePdb; |
20 |
|
21 |
import tinker.TinkerError; |
22 |
import tinker.TinkerRunner; |
23 |
|
24 |
|
25 |
public class reconstruct { |
26 |
|
27 |
private static final String TINKERBINDIR = "/project/StruPPi/Software/tinker/bin"; |
28 |
private static final String PRMFILE = "/project/StruPPi/Software/tinker/amber/amber99.prm"; |
29 |
|
30 |
public static void main(String[] args) { |
31 |
|
32 |
String programName = reconstruct.class.getName(); |
33 |
String help = "Usage:\n" + |
34 |
programName+" -p <pdb code> -c <pdb chain code> -t <contact_type> [-r] -d <distance cutoff 1> -D <distance cutoff 2> -i <distance cutoff 3> -b <base name> -o <output dir> [-n <number of models>]\n"; |
35 |
|
36 |
String pdbCode = ""; |
37 |
String pdbChainCode = ""; |
38 |
String ct = ""; |
39 |
double cutoff1 = 0.0; |
40 |
double cutoff2 = 0.0; |
41 |
double cutoff3 = 0.0; |
42 |
String outputDir = ""; |
43 |
String baseName = ""; |
44 |
boolean cross = false; |
45 |
int n = 1; |
46 |
|
47 |
Getopt g = new Getopt(programName, args, "p:c:d:t:rb:o:d:D:i:n:h?"); |
48 |
int c; |
49 |
while ((c = g.getopt()) != -1) { |
50 |
switch(c){ |
51 |
case 'p': |
52 |
pdbCode = g.getOptarg(); |
53 |
break; |
54 |
case 'c': |
55 |
pdbChainCode = g.getOptarg(); |
56 |
break; |
57 |
case 't': |
58 |
ct = g.getOptarg(); |
59 |
break; |
60 |
case 'r': |
61 |
cross = true; |
62 |
break; |
63 |
case 'd': |
64 |
cutoff1 = Double.valueOf(g.getOptarg()); |
65 |
break; |
66 |
case 'D': |
67 |
cutoff2 = Double.valueOf(g.getOptarg()); |
68 |
break; |
69 |
case 'i': |
70 |
cutoff3 = Double.valueOf(g.getOptarg()); |
71 |
break; |
72 |
case 'b': |
73 |
baseName = g.getOptarg(); |
74 |
break; |
75 |
case 'o': |
76 |
outputDir = g.getOptarg(); |
77 |
break; |
78 |
case 'n': |
79 |
n = Integer.parseInt(g.getOptarg()); |
80 |
break; |
81 |
case 'h': |
82 |
case '?': |
83 |
System.out.println(help); |
84 |
System.exit(0); |
85 |
break; // getopt() already printed an error |
86 |
} |
87 |
} |
88 |
|
89 |
if (pdbCode.equals("") || pdbChainCode.equals("") || ct.equals("") || cutoff1==0.0 || outputDir.equals("") || baseName.equals("")){ |
90 |
System.err.println("Must specify at least -p, -c, -t, -d, -o and -b"); |
91 |
System.err.println(help); |
92 |
System.exit(1); |
93 |
} |
94 |
|
95 |
if (baseName.contains(".")) { |
96 |
System.err.println("Basename can't contain a dot (not allowed by tinker). Exiting"); |
97 |
System.exit(1); |
98 |
} |
99 |
|
100 |
if (!AAinfo.isValidContactType(ct) || ct.contains("ALL")) { |
101 |
System.err.println("Invalid contact type specified. Exiting"); |
102 |
System.exit(1); |
103 |
} |
104 |
|
105 |
if (n>999) { |
106 |
System.err.println("Maximum number of models is 999. Specify a lower value. Exiting"); |
107 |
System.exit(1); |
108 |
} |
109 |
|
110 |
boolean doublecm = false; |
111 |
String ct1 = ct; |
112 |
String ct2 = ct; |
113 |
String ct3 = null; |
114 |
if (ct.contains("_")) { |
115 |
ct1 = ct.split("_")[0]; |
116 |
ct2 = ct.split("_")[1]; |
117 |
doublecm = true; |
118 |
} |
119 |
if (cross) { |
120 |
ct3 = ct1+"/"+ct2; |
121 |
} |
122 |
|
123 |
if (cutoff2==0.0) { |
124 |
cutoff2 = cutoff1; |
125 |
} |
126 |
if (cutoff3==0.0) { |
127 |
cutoff3 = cutoff1; |
128 |
} |
129 |
|
130 |
Pdb pdb = null; |
131 |
try { |
132 |
pdb = new PdbasePdb(pdbCode,pdbChainCode); |
133 |
} catch (PdbaseInconsistencyError e) { |
134 |
System.err.println("Pdbase inconsistency for structure "+pdbCode+". Can't continue, exiting"); |
135 |
System.exit(1); |
136 |
} catch (PdbCodeNotFoundError e) { |
137 |
System.err.println("Given pdb code "+pdbCode+" couldn't be found in pdbase. Exiting"); |
138 |
System.exit(1); |
139 |
} catch (SQLException e) { |
140 |
System.err.println("Problems connecting to database for getting pdb data for "+pdbCode+". Exiting"); |
141 |
System.exit(1); |
142 |
} catch (PdbChainCodeNotFoundError e) { |
143 |
System.err.println("Given pdb chain code "+pdbChainCode+" couldn't be found for pdb code "+pdbCode+". Exiting"); |
144 |
System.exit(1); |
145 |
} |
146 |
// we also write the file to the out dir so it can be used later for clustering rmsds etc. |
147 |
File origPdbFile = new File (outputDir,baseName+".native.pdb"); |
148 |
try { |
149 |
pdb.dump2pdbfile(origPdbFile.getAbsolutePath()); |
150 |
} catch (IOException e4) { |
151 |
System.err.println("Couldn't write original pdb file "+origPdbFile.getAbsolutePath()); |
152 |
System.err.println("Continuing without it, this is not needed for the rest of the reconstruction process but only for post processing (e.g. comparing rmsds to original)"); |
153 |
} |
154 |
|
155 |
String sequence = pdb.getSequence(); |
156 |
|
157 |
Graph[] graphs =null; |
158 |
Graph graph1 = pdb.get_graph(ct1, cutoff1); |
159 |
Graph graph2 = null; |
160 |
Graph graph3 = null; |
161 |
if (doublecm) { |
162 |
graph2 = pdb.get_graph(ct2, cutoff2); |
163 |
if (cross) { |
164 |
graph3 = pdb.get_graph(ct3, cutoff3); |
165 |
graphs = new Graph[3]; |
166 |
graphs[0] = graph1; |
167 |
graphs[1] = graph2; |
168 |
graphs[2] = graph3; |
169 |
} else { |
170 |
graphs = new Graph[2]; |
171 |
graphs[0] = graph1; |
172 |
graphs[1] = graph2; |
173 |
} |
174 |
} else { |
175 |
graphs = new Graph[1]; |
176 |
graphs[0] = graph1; |
177 |
} |
178 |
|
179 |
// defining files |
180 |
File reportFile = new File(outputDir,baseName+".report"); |
181 |
|
182 |
// creating TinkerRunner object |
183 |
|
184 |
TinkerRunner tr = null; |
185 |
try { |
186 |
tr = new TinkerRunner(TINKERBINDIR, PRMFILE); |
187 |
} catch (FileNotFoundException e3) { |
188 |
System.err.println("Couldn't find tinker bin dir "+TINKERBINDIR+". Exiting"); |
189 |
System.exit(1); |
190 |
} |
191 |
|
192 |
// call reconstruction |
193 |
|
194 |
try { |
195 |
tr.reconstruct(sequence, graphs, n, outputDir, baseName, false); |
196 |
} catch (IOException e) { |
197 |
System.err.println("Error while running Tinker reconstruction: " + e.getMessage()); |
198 |
} catch (TinkerError e) { |
199 |
System.err.println("Error while running Tinker reconstruction: " + e.getMessage()); |
200 |
} |
201 |
|
202 |
double[] err = tr.getErrorFunctionVal(); |
203 |
double[] mubv = tr.getMaxUpperBoundViol(); |
204 |
double[] mlbv = tr.getMaxLowerBoundViol(); |
205 |
double[] muv = tr.getMaxUpperViol(); |
206 |
double[] mlv = tr.getMaxLowerViol(); |
207 |
int[] nubv = tr.getNumUpperBoundViol(); |
208 |
int[] nlbv = tr.getNumLowerBoundViol(); |
209 |
int[] nuv = tr.getNumUpperViol(); |
210 |
int[] nlv = tr.getNumLowerViol(); |
211 |
double[] rbv = tr.getRmsBoundViol(); |
212 |
double[] rrv = tr.getRmsRestViol(); |
213 |
|
214 |
|
215 |
// calculate rmsds |
216 |
|
217 |
double[] rmsds = new double[n+1]; |
218 |
|
219 |
for (int i = 1; i<=n; i++) { |
220 |
String ext = new Formatter().format(".%03d",i).toString(); |
221 |
File outputPdbFile = new File(outputDir, baseName+ext+".pdb"); |
222 |
try { |
223 |
Pdb outputPdb = tr.getStructure(i); |
224 |
rmsds[i] = pdb.rmsd(outputPdb, "Ca"); |
225 |
} |
226 |
catch (TinkerError e) { |
227 |
System.err.println("Error while trying to retrieve results from Tinker: + e.getMessage()"); |
228 |
} catch (ConformationsNotSameSizeError e) { |
229 |
System.err.println(origPdbFile.getAbsolutePath()+" and "+outputPdbFile.getAbsolutePath()+" don't have the same conformation size, can't calculate rmsd for them."); |
230 |
} |
231 |
} |
232 |
|
233 |
// write report file |
234 |
|
235 |
try { |
236 |
PrintWriter reportOut = new PrintWriter(new FileOutputStream(reportFile)); |
237 |
reportOut.println("run_id\tcutoff\tcutoff2\tcutoff3\tct\tct2\tct3\tnum_res" + |
238 |
"\tresult_id\terror_val" + |
239 |
"\tup_bound_viol\tlow_bound_viol\tmax_bound_up\tmax_bound_low\trms_bound" + |
240 |
"\tup_viol\tlow_viol\tmax_up\tmax_low\trms_viol\trmsd_to_orig"); |
241 |
|
242 |
for (int i=1;i<=n;i++){ |
243 |
String rmsd = String.format(Locale.US,"%6.3f",rmsds[i]); |
244 |
String errStr = String.format(Locale.US, "%6.3f", err[i]); |
245 |
// run_id cutoff cutoff2 cutoff3 ct1 ct2 ct3 num_res |
246 |
reportOut.println(pdbCode+"_"+pdbChainCode+"\t"+cutoff1+"\t"+cutoff2+"\t"+cutoff3+"\t"+ct1+"\t"+ct2+"\t"+ct3+"\t"+sequence.length()+"\t"+ |
247 |
// result_id error_val |
248 |
i + "\t" + errStr + "\t" + |
249 |
// up_bound_viol low_bound_viol max_bound_up max_bound_low rms_bound |
250 |
nubv[i] + "\t" + nlbv[i] + "\t" + mubv[i] + "\t" + mlbv[i] + "\t" + rbv[i] + "\t" + |
251 |
// up_viol low_viol max_up max_low rms_viol |
252 |
nuv[i] + "\t" + nlv[i] + "\t" + muv[i] + "\t" + mlv[i] + "\t"+ rrv[i] + "\t"+ |
253 |
// rmsd_to_orig |
254 |
rmsd); |
255 |
} |
256 |
reportOut.close(); |
257 |
} catch (FileNotFoundException e) { |
258 |
System.err.println("Couldn't write to report file "+reportFile.getAbsolutePath() +". Error: "+e.getMessage()); |
259 |
} |
260 |
|
261 |
} |
262 |
|
263 |
} |