ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/owl/trunk/reconstruct.java
Revision: 346
Committed: Thu Oct 11 08:38:36 2007 UTC (17 years ago) by duarte
File size: 10094 byte(s)
Log Message:
Doing more input error checks
In ConstraintsMaker createConstraints won't throw Exception anymore, we only assery that ALL can't be used
Line User Rev File contents
1 duarte 338 import gnu.getopt.Getopt;
2    
3     import java.io.File;
4     import java.io.FileNotFoundException;
5     import java.io.IOException;
6     import java.sql.SQLException;
7     import java.util.Formatter;
8    
9 duarte 346 import proteinstructure.AAinfo;
10 duarte 338 import proteinstructure.ConformationsNotSameSizeError;
11     import proteinstructure.Graph;
12     import proteinstructure.Pdb;
13     import proteinstructure.PdbChainCodeNotFoundError;
14     import proteinstructure.PdbCodeNotFoundError;
15     import proteinstructure.PdbaseInconsistencyError;
16     import proteinstructure.PdbasePdb;
17     import proteinstructure.PdbfileFormatError;
18     import proteinstructure.PdbfilePdb;
19    
20     import tinker.ConstraintsMaker;
21     import tinker.TinkerError;
22     import tinker.TinkerRunner;
23    
24    
25     public class reconstruct {
26    
27    
28     private static final String TINKERBINDIR = "/project/StruPPi/Software/tinker/bin";
29     private static final String PRMFILE = "/project/StruPPi/Software/tinker/amber/amber99.prm";
30     private static final String PRMTYPE = "amber";
31    
32     public static void main(String[] args) {
33    
34     String programName = reconstruct.class.getName();
35     String help = "Usage:\n" +
36     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";
37    
38     String pdbCode = "";
39     String pdbChainCode = "";
40     String ct = "";
41     double cutoff1 = 0.0;
42     double cutoff2 = 0.0;
43     double cutoff3 = 0.0;
44     String outputDir = "";
45     String baseName = "";
46     boolean cross = false;
47     int n = 1;
48    
49     Getopt g = new Getopt(programName, args, "p:c:d:t:rb:o:d:D:i:n:h?");
50     int c;
51     while ((c = g.getopt()) != -1) {
52     switch(c){
53     case 'p':
54     pdbCode = g.getOptarg();
55     break;
56     case 'c':
57     pdbChainCode = g.getOptarg();
58     break;
59     case 't':
60     ct = g.getOptarg();
61     break;
62     case 'r':
63     cross = true;
64     break;
65     case 'd':
66     cutoff1 = Double.valueOf(g.getOptarg());
67     break;
68     case 'D':
69     cutoff2 = Double.valueOf(g.getOptarg());
70     break;
71     case 'i':
72     cutoff3 = Double.valueOf(g.getOptarg());
73     break;
74     case 'b':
75     baseName = g.getOptarg();
76     break;
77     case 'o':
78     outputDir = g.getOptarg();
79     break;
80     case 'n':
81     n = Integer.parseInt(g.getOptarg());
82     break;
83     case 'h':
84     case '?':
85     System.out.println(help);
86     System.exit(0);
87     break; // getopt() already printed an error
88     }
89     }
90    
91     if (pdbCode.equals("") || pdbChainCode.equals("") || ct.equals("") || cutoff1==0.0 || outputDir.equals("") || baseName.equals("")){
92     System.err.println("Must specify at least -p, -c, -t, -d, -o and -b");
93     System.err.println(help);
94     System.exit(1);
95     }
96    
97 duarte 346 if (baseName.contains(".")) {
98     System.err.println("Basename can't contain a dot (not allowed by tinker). Exiting");
99     System.exit(1);
100     }
101    
102     if (!AAinfo.isValidContactType(ct) || ct.contains("ALL")) {
103     System.err.println("Invalid contact type specified. Exiting");
104     System.exit(1);
105     }
106    
107     if (n>999) {
108     System.err.println("Maximum number of models is 999. Specify a lower value. Exiting");
109     System.exit(1);
110     }
111    
112 duarte 338 boolean doublecm = false;
113     String ct1 = ct;
114     String ct2 = ct;
115     String ct3 = null;
116     if (ct.contains("_")) {
117     ct1 = ct.split("_")[0];
118     ct2 = ct.split("_")[1];
119     doublecm = true;
120     }
121     if (cross) {
122     ct3 = ct1+"/"+ct2;
123     }
124    
125     if (cutoff2==0.0) {
126     cutoff2 = cutoff1;
127     }
128     if (cutoff3==0.0) {
129     cutoff3 = cutoff1;
130     }
131    
132    
133     Pdb pdb = null;
134     try {
135     pdb = new PdbasePdb(pdbCode,pdbChainCode);
136     } catch (PdbaseInconsistencyError e) {
137     System.err.println("Pdbase inconsistency for structure "+pdbCode+". Can't continue, exiting");
138     System.exit(1);
139     } catch (PdbCodeNotFoundError e) {
140     System.err.println("Given pdb code "+pdbCode+" couldn't be found in pdbase. Exiting");
141     System.exit(1);
142     } catch (SQLException e) {
143     System.err.println("Problems connecting to database for getting pdb data for "+pdbCode+". Exiting");
144     System.exit(1);
145     } catch (PdbChainCodeNotFoundError e) {
146     System.err.println("Given pdb chain code "+pdbChainCode+" couldn't be found for pdb code "+pdbCode+". Exiting");
147     System.exit(1);
148     }
149     String sequence = pdb.getSequence();
150    
151     Graph graph1 = pdb.get_graph(ct1, cutoff1);
152     Graph graph2 = null;
153     Graph graph3 = null;
154     if (doublecm) {
155     graph2 = pdb.get_graph(ct2, cutoff2);
156     }
157     if (cross) {
158     graph3 = pdb.get_graph(ct3, cutoff3);
159     }
160    
161    
162     // defining files
163     File logFile = new File(outputDir,"tinker.log");
164     File prmFile = new File(PRMFILE);
165     File xyzFile = new File(outputDir,baseName+".xyz");
166     File seqFile = new File(outputDir,baseName+".seq");
167     File pdbFile = new File(outputDir,baseName+".pdb");
168     File keyFile = new File(outputDir,baseName+".key");
169    
170     // creating TinkerRunner object
171     TinkerRunner tr = null;
172     try {
173     tr = new TinkerRunner(TINKERBINDIR,PRMFILE,logFile);
174     } catch (FileNotFoundException e3) {
175     System.err.println("Couldn't find tinker bin dir "+TINKERBINDIR+". Exiting");
176     System.exit(1);
177     }
178    
179     // 1. run tinker's protein program
180     try {
181     tr.runProtein(sequence, outputDir, baseName );
182     } catch (IOException e2) {
183     System.err.println("Couldn't read file to run 'protein' "+seqFile.getAbsolutePath());
184     System.err.println("Exiting");
185     System.exit(1);
186     } catch (TinkerError e2) {
187     System.err.println("Tinker error while running 'protein', check log file "+logFile.getAbsolutePath()+". Exiting");
188     System.exit(1);
189     }
190    
191     // 1a. convert xyz file to pdb to be able to map atom serials after
192     try {
193     tr.runXyzpdb(xyzFile, seqFile, pdbFile);
194     } catch (IOException e1) {
195     System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+" or "+seqFile.getAbsolutePath()+" or write to "+pdbFile.getAbsolutePath()+" for running 'xyzpdb'");
196     System.err.println("Exiting");
197     System.exit(1);
198     } catch (TinkerError e1) {
199     System.err.println("Tinker error while running xyzpdb, check log file "+logFile.getAbsolutePath()+". Exiting");
200     System.exit(1);
201     }
202    
203     // 2. creating constraints into key file
204     ConstraintsMaker cm = null;
205     try {
206     cm = new ConstraintsMaker(pdbFile,xyzFile,prmFile,PRMTYPE,keyFile);
207     } catch (IOException e3) {
208     System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+", "+pdbFile.getAbsolutePath()+" or, "+prmFile.getAbsolutePath()+" write to "+keyFile.getAbsolutePath()+" for creating distance constraints");
209     System.err.println("Exiting");
210     System.exit(1);
211     } catch (PdbfileFormatError e3) {
212     System.err.println("pdb file "+pdbFile.getAbsolutePath()+" converted from "+xyzFile.getAbsolutePath()+" doesn't seem to be in the right format. Check log? ("+logFile.getAbsolutePath()+"). Exiting");
213     System.exit(1);
214     }
215 duarte 346
216     cm.createConstraints(graph1);
217     if (doublecm) cm.createConstraints(graph2);
218     if (cross) cm.createConstraints(graph3);
219    
220 duarte 338 cm.closeKeyFile();
221    
222     // 3. run tinker's distgeom
223     try {
224     System.out.println("Running distgeom...");
225     tr.runDistgeom(xyzFile, outputDir, baseName, n);
226     } catch (TinkerError e1) {
227 duarte 340 System.err.println(e1.getMessage());
228     System.err.println("Exiting");
229 duarte 338 System.exit(1);
230     } catch (IOException e1) {
231     System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+" or write 'distgeom' output files to output dir "+outputDir);
232     System.err.println("Exiting");
233     System.exit(1);
234 stehr 344 } catch(InterruptedException e) {
235     System.err.println("Distgeom was interrupted:" + e.getMessage());
236     System.err.println("Exiting.");
237     System.exit(1);
238 duarte 338 }
239    
240     double[] mubv = tr.getMaxUpperBoundViol();
241     double[] mlbv = tr.getMaxLowerBoundViol();
242     double[] muv = tr.getMaxUpperViol();
243     double[] mlv = tr.getMaxLowerViol();
244     int[] nubv = tr.getNumUpperBoundViol();
245     int[] nlbv = tr.getNumLowerBoundViol();
246     int[] nuv = tr.getNumUpperViol();
247     int[] nlv = tr.getNumLowerViol();
248     double[] rbv = tr.getRmsBoundViol();
249     double[] rrv = tr.getRmsRestViol();
250    
251    
252     // 4. converting xyz output files to pdb files and calculating rmsds
253    
254     double[] rmsds = new double[n+1];
255    
256     for (int i = 1; i<=n; i++) {
257     String ext = new Formatter().format(".%03d",i).toString();
258     File outputXyzFile = new File(outputDir, baseName+ext);
259     File outputPdbFile = new File(outputDir, baseName+ext+".pdb");
260     try {
261     tr.runXyzpdb(outputXyzFile, seqFile, outputPdbFile);
262    
263     Pdb outputPdb = new PdbfilePdb(outputPdbFile.getAbsolutePath(),"NULL");
264     rmsds[i] = pdb.rmsd(outputPdb, "Ca");
265    
266     } catch (IOException e) {
267     System.err.println("Couldn't read file "+outputXyzFile.getAbsolutePath()+", or "+seqFile.getAbsolutePath()+", or write to "+outputPdbFile.getAbsolutePath()+" while converting with 'xyzpdb'");
268     System.err.println("Can't calculate rmsd for it");
269     } catch (TinkerError e) {
270     System.err.println("Tinker error while running 'xyzpdb' to convert"+outputXyzFile.getAbsolutePath()+", check log file "+logFile.getAbsolutePath());
271     System.err.println("Can't calculate rmsd for it");
272     }
273     catch (PdbfileFormatError e) {
274     System.err.println("Output pdb file "+outputPdbFile.getAbsolutePath()+" doesn't seem to be in the correcet format. Can't calculate rmsd for it");
275     } catch (PdbChainCodeNotFoundError e) {
276     // this shouldn't happen, chain code is hard coded, we throw stack trace and continue if it happens
277     e.printStackTrace();
278     } catch (ConformationsNotSameSizeError e) {
279     System.err.println(pdbFile.getAbsolutePath()+" and "+outputPdbFile.getAbsolutePath()+" don't have the same conformation size, can't calculate rmsd for them.");
280     }
281    
282     tr.closeLog();
283     }
284    
285    
286     // 6. report
287     System.out.println("accession_code\tchain_pdb_code\tcutoff\tcutoff2\tcutoff3\tct\tct2\tct3" +
288     "\tresult_id\tup_bound_viol\tlow_bound_viol\tmax_bound_up\tmax_bound_low\trms_bound" +
289     "\tup_viol\tlow_viol\tmax_up\tmax_low\trms_viol\trmsd_to_orig");
290    
291     for (int i=1;i<=n;i++){
292     System.out.println(pdbCode+"\t"+pdbChainCode+"\t"+cutoff1+"\t"+cutoff2+"\t"+cutoff3+"\t"+ct1+"\t"+ct2+"\t"+ct3+"\t"+
293 duarte 340 i+"\t"+nubv[i]+"\t"+nlbv[i]+"\t"+mubv[i]+"\t"+mlbv[i]+"\t"+muv[i]+"\t"+mlv[i]+"\t"+rbv[i]+"\t"+
294 duarte 338 "\t"+nuv[i]+"\t"+nlv[i]+"\t"+rrv[i]+"\t"+rmsds[i]);
295     }
296    
297     }
298    
299     }