owl/trunk/reconstruct.java

import gnu.getopt.Getopt;

import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.PrintWriter;
import java.sql.SQLException;
import java.util.Formatter;
import java.util.Locale;

import proteinstructure.AAinfo;
import proteinstructure.ConformationsNotSameSizeError;
import proteinstructure.PdbLoadError;
import proteinstructure.RIGraph;
import proteinstructure.Pdb;
import proteinstructure.PdbCodeNotFoundError;
import proteinstructure.PdbasePdb;

import tinker.TinkerError;
import tinker.TinkerRunner;


public class reconstruct {

        private static final String TINKERBINDIR = "/project/StruPPi/Software/tinker/bin";
        private static final String PRMFILE = "/project/StruPPi/Software/tinker/amber/amber99.prm";
        
        public static void main(String[] args) {
                
                String programName = reconstruct.class.getName();
                String help = "Usage:\n" +
                        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>] [-m <min range>] [-M <max range>] [-f <force constant>]\n"; 

                String pdbCode = "";
                String pdbChainCode = "";
                String ct = "";
                double cutoff1 = 0.0;
                double cutoff2 = 0.0;
                double cutoff3 = 0.0;
                String outputDir = "";
                String baseName = "";
                boolean cross = false;
                int n = 1;
                double forceConstant = 0;
                int minRange = 0;
                int maxRange = 0;
                
                Getopt g = new Getopt(programName, args, "p:c:d:t:rb:o:d:D:i:n:m:M:f:h?");
                int c;
                while ((c = g.getopt()) != -1) {
                        switch(c){
                        case 'p':
                                pdbCode = g.getOptarg();
                                break;
                        case 'c':
                                pdbChainCode = g.getOptarg();
                                break;
                        case 't':
                                ct = g.getOptarg();
                                break;
                        case 'r':
                                cross = true;
                                break;
                        case 'd':
                                cutoff1 = Double.valueOf(g.getOptarg());
                                break;
                        case 'D':
                                cutoff2 = Double.valueOf(g.getOptarg());
                                break;
                        case 'i':
                                cutoff3 = Double.valueOf(g.getOptarg());
                                break;
                        case 'b':
                                baseName = g.getOptarg();
                                break;
                        case 'o':
                                outputDir = g.getOptarg();
                                break;
                        case 'n':
                                n = Integer.parseInt(g.getOptarg());
                                break;
                        case 'm':
                                minRange = Integer.parseInt(g.getOptarg());
                                break;
                        case 'M':
                                maxRange = Integer.parseInt(g.getOptarg());
                                break;                          
                        case 'f':
                                forceConstant = Double.valueOf(g.getOptarg());
                                break;                          
                        case 'h':
                        case '?':
                                System.out.println(help);
                                System.exit(0);
                                break; // getopt() already printed an error
                        }
                }

                if (pdbCode.equals("") || pdbChainCode.equals("") || ct.equals("") || cutoff1==0.0 || outputDir.equals("") || baseName.equals("")){
                        System.err.println("Must specify at least -p, -c, -t, -d, -o and -b");
                        System.err.println(help);
                        System.exit(1);
                }
                
                if (baseName.contains(".")) {
                        System.err.println("Basename can't contain a dot (not allowed by tinker). Exiting");
                        System.exit(1);
                }
                
                if (!AAinfo.isValidContactType(ct) || ct.contains("ALL")) {
                        System.err.println("Invalid contact type specified. Exiting");
                        System.exit(1);
                }
                
                if (n>999) {
                        System.err.println("Maximum number of models is 999. Specify a lower value. Exiting");
                        System.exit(1);
                }
                
                boolean doublecm = false;
                String ct1 = ct;
                String ct2 = ct;
                String ct3 = null;
                if (ct.contains("_")) {
                        ct1 = ct.split("_")[0];
                        ct2 = ct.split("_")[1];
                        doublecm = true;
                }
                if (cross) {
                        ct3 = ct1+"/"+ct2;
                }
                
                if (cutoff2==0.0) {
                        cutoff2 = cutoff1;
                }
                if (cutoff3==0.0) {
                        cutoff3 = cutoff1;
                }
                
                Pdb pdb = null;
                try {
                        pdb = new PdbasePdb(pdbCode);
                        pdb.load(pdbChainCode);
                } catch (PdbLoadError e) {
                        System.err.println("Error while loading pdb data. Specific error "+e.getMessage());
                        System.exit(1);
                } catch (PdbCodeNotFoundError e) {
                        System.err.println("Given pdb code "+pdbCode+" couldn't be found in pdbase. Exiting");
                        System.exit(1);
                } catch (SQLException e) {
                        System.err.println("Problems connecting to database for getting pdb data for "+pdbCode+". Exiting");
                        System.exit(1);
                }
                // we also write the file to the out dir so it can be used later for clustering rmsds etc.
                File origPdbFile = new File (outputDir,baseName+".native.pdb");
                try {
                        pdb.dump2pdbfile(origPdbFile.getAbsolutePath());
                } catch (IOException e4) {
                        System.err.println("Couldn't write original pdb file "+origPdbFile.getAbsolutePath());
                        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)");
                }

                String sequence = pdb.getSequence();

                RIGraph[] graphs =null;
                RIGraph graph1 = pdb.get_graph(ct1, cutoff1);
                if (maxRange>0) graph1.restrictContactsToMaxRange(maxRange);
                if (minRange>0) graph1.restrictContactsToMinRange(minRange);
                RIGraph graph2 = null;
                RIGraph graph3 = null;
                if (doublecm) {
                        graph2 = pdb.get_graph(ct2, cutoff2);
                        if (maxRange>0) graph2.restrictContactsToMaxRange(maxRange);
                        if (minRange>0) graph2.restrictContactsToMinRange(minRange);
                        if (cross) {
                                graph3 = pdb.get_graph(ct3, cutoff3);
                                if (maxRange>0) graph3.restrictContactsToMaxRange(maxRange);
                                if (minRange>0) graph3.restrictContactsToMinRange(minRange);
                                graphs = new RIGraph[3];
                                graphs[0] = graph1;
                                graphs[1] = graph2;
                                graphs[2] = graph3;
                        } else {
                                graphs = new RIGraph[2];
                                graphs[0] = graph1;
                                graphs[1] = graph2;
                        }
                } else {
                        graphs = new RIGraph[1];
                        graphs[0] = graph1;
                }
                
                // defining files
                File reportFile = new File(outputDir,baseName+".report");
                
                // creating TinkerRunner object
                
                TinkerRunner tr = null;
                try {
                        if (forceConstant!=0) {
                                tr = new TinkerRunner(TINKERBINDIR, PRMFILE, forceConstant);
                        } else {
                                // force constant not given in command line: take the default force constant
                                tr = new TinkerRunner(TINKERBINDIR, PRMFILE); 
                        }
                } catch (FileNotFoundException e3) {
                        System.err.println("Couldn't find tinker bin dir "+TINKERBINDIR+". Exiting");
                        System.exit(1);
                }
                
                // call reconstruction
                
                try {
                        tr.reconstruct(sequence, graphs, n, outputDir, baseName, false);
                } catch (IOException e) {
                        System.err.println("Error while running Tinker reconstruction: " + e.getMessage());
                        System.exit(1);
                } catch (TinkerError e) {
                        System.err.println("Error while running Tinker reconstruction: " + e.getMessage());
                        System.exit(1);
                }
                                
                double[] err = tr.getErrorFunctionVal();
                double[] mubv = tr.getMaxUpperBoundViol();
                double[] mlbv = tr.getMaxLowerBoundViol();
                double[] muv = tr.getMaxUpperViol();
                double[] mlv = tr.getMaxLowerViol();
                int[] nubv = tr.getNumUpperBoundViol();
                int[] nlbv = tr.getNumLowerBoundViol();
                int[] nuv = tr.getNumUpperViol();
                int[] nlv = tr.getNumLowerViol();
                double[] rbv = tr.getRmsBoundViol();
                double[] rrv = tr.getRmsRestViol();
                

                // calculate rmsds

                double[] rmsds = new double[n+1];               
                
                for (int i = 1; i<=n; i++) {
                        String ext = new Formatter().format(".%03d",i).toString();
                        File outputPdbFile = new File(outputDir, baseName+ext+".pdb");
                        try {
                                Pdb outputPdb = tr.getStructure(i);
                                rmsds[i] = pdb.rmsd(outputPdb, "Ca");
                        }
                        catch (TinkerError e) {
                                System.err.println("Error while trying to retrieve results from Tinker: "+ e.getMessage());
                        } catch (ConformationsNotSameSizeError e) {
                                System.err.println(origPdbFile.getAbsolutePath()+" and "+outputPdbFile.getAbsolutePath()+" don't have the same conformation size, can't calculate rmsd for them.");
                        }                               
                }                                       

                // write report file
                
                try {
                        PrintWriter reportOut = new PrintWriter(new FileOutputStream(reportFile));
                        reportOut.println("run_id\tcutoff\tcutoff2\tcutoff3\tct\tct2\tct3\tnum_res" +
                                                "\tresult_id\terror_val" +
                                                "\tup_bound_viol\tlow_bound_viol\tmax_bound_up\tmax_bound_low\trms_bound" +
                                                "\tup_viol\tlow_viol\tmax_up\tmax_low\trms_viol\trmsd_to_orig");

                        for (int i=1;i<=n;i++){
                                String rmsd = String.format(Locale.US,"%6.3f",rmsds[i]);
                                String errStr = String.format(Locale.US, "%6.3f", err[i]);
                                //                         run_id                cutoff      cutoff2      cutoff3      ct1      ct2      ct3         num_res
                                reportOut.println(pdbCode+"_"+pdbChainCode+"\t"+cutoff1+"\t"+cutoff2+"\t"+cutoff3+"\t"+ct1+"\t"+ct2+"\t"+ct3+"\t"+sequence.length()+"\t"+
                                //  result_id     error_val         
                                                i + "\t" + errStr + "\t" +
                                //   up_bound_viol    low_bound_viol   max_bound_up    max_bound_low      rms_bound
                                                nubv[i] + "\t" + nlbv[i] + "\t" + mubv[i] + "\t" + mlbv[i] + "\t" + rbv[i] + "\t" +
                                //      up_viol         low_viol       max_up         max_low         rms_viol
                                                nuv[i] + "\t" + nlv[i] + "\t" + muv[i] + "\t" + mlv[i] + "\t"+ rrv[i] + "\t"+
                                //      rmsd_to_orig
                                                rmsd);
                        }
                        reportOut.close();
                } catch (FileNotFoundException e) {
                        System.err.println("Couldn't write to report file "+reportFile.getAbsolutePath() +". Error: "+e.getMessage());
                }
                
        }

}
Revision:	538
Committed:	Thu Feb 14 11:56:26 2008 UTC (17 years ago) by duarte
File size:	9049 byte(s)
Log Message:	New option in reconstruct: -f force constant. New constructor in TinkerRunner that takes the force constant, previous constructor will work as before with the default constant.
Line	File contents
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.PdbLoadError;
15	import proteinstructure.RIGraph;
16	import proteinstructure.Pdb;
17	import proteinstructure.PdbCodeNotFoundError;
18	import proteinstructure.PdbasePdb;
19
20	import tinker.TinkerError;
21	import tinker.TinkerRunner;
22
23
24	public class reconstruct {
25
26	private static final String TINKERBINDIR = "/project/StruPPi/Software/tinker/bin";
27	private static final String PRMFILE = "/project/StruPPi/Software/tinker/amber/amber99.prm";
28
29	public static void main(String[] args) {
30
31	String programName = reconstruct.class.getName();
32	String help = "Usage:\n" +
33	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>] [-m <min range>] [-M <max range>] [-f <force constant>]\n";
34
35	String pdbCode = "";
36	String pdbChainCode = "";
37	String ct = "";
38	double cutoff1 = 0.0;
39	double cutoff2 = 0.0;
40	double cutoff3 = 0.0;
41	String outputDir = "";
42	String baseName = "";
43	boolean cross = false;
44	int n = 1;
45	double forceConstant = 0;
46	int minRange = 0;
47	int maxRange = 0;
48
49	Getopt g = new Getopt(programName, args, "p:c:d:t:rb:o:d:D:i:n:m:M:f: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 'm':
84	minRange = Integer.parseInt(g.getOptarg());
85	break;
86	case 'M':
87	maxRange = Integer.parseInt(g.getOptarg());
88	break;
89	case 'f':
90	forceConstant = Double.valueOf(g.getOptarg());
91	break;
92	case 'h':
93	case '?':
94	System.out.println(help);
95	System.exit(0);
96	break; // getopt() already printed an error
97	}
98	}
99
100	if (pdbCode.equals("") \|\| pdbChainCode.equals("") \|\| ct.equals("") \|\| cutoff1==0.0 \|\| outputDir.equals("") \|\| baseName.equals("")){
101	System.err.println("Must specify at least -p, -c, -t, -d, -o and -b");
102	System.err.println(help);
103	System.exit(1);
104	}
105
106	if (baseName.contains(".")) {
107	System.err.println("Basename can't contain a dot (not allowed by tinker). Exiting");
108	System.exit(1);
109	}
110
111	if (!AAinfo.isValidContactType(ct) \|\| ct.contains("ALL")) {
112	System.err.println("Invalid contact type specified. Exiting");
113	System.exit(1);
114	}
115
116	if (n>999) {
117	System.err.println("Maximum number of models is 999. Specify a lower value. Exiting");
118	System.exit(1);
119	}
120
121	boolean doublecm = false;
122	String ct1 = ct;
123	String ct2 = ct;
124	String ct3 = null;
125	if (ct.contains("_")) {
126	ct1 = ct.split("_")[0];
127	ct2 = ct.split("_")[1];
128	doublecm = true;
129	}
130	if (cross) {
131	ct3 = ct1+"/"+ct2;
132	}
133
134	if (cutoff2==0.0) {
135	cutoff2 = cutoff1;
136	}
137	if (cutoff3==0.0) {
138	cutoff3 = cutoff1;
139	}
140
141	Pdb pdb = null;
142	try {
143	pdb = new PdbasePdb(pdbCode);
144	pdb.load(pdbChainCode);
145	} catch (PdbLoadError e) {
146	System.err.println("Error while loading pdb data. Specific error "+e.getMessage());
147	System.exit(1);
148	} catch (PdbCodeNotFoundError e) {
149	System.err.println("Given pdb code "+pdbCode+" couldn't be found in pdbase. Exiting");
150	System.exit(1);
151	} catch (SQLException e) {
152	System.err.println("Problems connecting to database for getting pdb data for "+pdbCode+". Exiting");
153	System.exit(1);
154	}
155	// we also write the file to the out dir so it can be used later for clustering rmsds etc.
156	File origPdbFile = new File (outputDir,baseName+".native.pdb");
157	try {
158	pdb.dump2pdbfile(origPdbFile.getAbsolutePath());
159	} catch (IOException e4) {
160	System.err.println("Couldn't write original pdb file "+origPdbFile.getAbsolutePath());
161	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)");
162	}
163
164	String sequence = pdb.getSequence();
165
166	RIGraph[] graphs =null;
167	RIGraph graph1 = pdb.get_graph(ct1, cutoff1);
168	if (maxRange>0) graph1.restrictContactsToMaxRange(maxRange);
169	if (minRange>0) graph1.restrictContactsToMinRange(minRange);
170	RIGraph graph2 = null;
171	RIGraph graph3 = null;
172	if (doublecm) {
173	graph2 = pdb.get_graph(ct2, cutoff2);
174	if (maxRange>0) graph2.restrictContactsToMaxRange(maxRange);
175	if (minRange>0) graph2.restrictContactsToMinRange(minRange);
176	if (cross) {
177	graph3 = pdb.get_graph(ct3, cutoff3);
178	if (maxRange>0) graph3.restrictContactsToMaxRange(maxRange);
179	if (minRange>0) graph3.restrictContactsToMinRange(minRange);
180	graphs = new RIGraph[3];
181	graphs[0] = graph1;
182	graphs[1] = graph2;
183	graphs[2] = graph3;
184	} else {
185	graphs = new RIGraph[2];
186	graphs[0] = graph1;
187	graphs[1] = graph2;
188	}
189	} else {
190	graphs = new RIGraph[1];
191	graphs[0] = graph1;
192	}
193
194	// defining files
195	File reportFile = new File(outputDir,baseName+".report");
196
197	// creating TinkerRunner object
198
199	TinkerRunner tr = null;
200	try {
201	if (forceConstant!=0) {
202	tr = new TinkerRunner(TINKERBINDIR, PRMFILE, forceConstant);
203	} else {
204	// force constant not given in command line: take the default force constant
205	tr = new TinkerRunner(TINKERBINDIR, PRMFILE);
206	}
207	} catch (FileNotFoundException e3) {
208	System.err.println("Couldn't find tinker bin dir "+TINKERBINDIR+". Exiting");
209	System.exit(1);
210	}
211
212	// call reconstruction
213
214	try {
215	tr.reconstruct(sequence, graphs, n, outputDir, baseName, false);
216	} catch (IOException e) {
217	System.err.println("Error while running Tinker reconstruction: " + e.getMessage());
218	System.exit(1);
219	} catch (TinkerError e) {
220	System.err.println("Error while running Tinker reconstruction: " + e.getMessage());
221	System.exit(1);
222	}
223
224	double[] err = tr.getErrorFunctionVal();
225	double[] mubv = tr.getMaxUpperBoundViol();
226	double[] mlbv = tr.getMaxLowerBoundViol();
227	double[] muv = tr.getMaxUpperViol();
228	double[] mlv = tr.getMaxLowerViol();
229	int[] nubv = tr.getNumUpperBoundViol();
230	int[] nlbv = tr.getNumLowerBoundViol();
231	int[] nuv = tr.getNumUpperViol();
232	int[] nlv = tr.getNumLowerViol();
233	double[] rbv = tr.getRmsBoundViol();
234	double[] rrv = tr.getRmsRestViol();
235
236
237	// calculate rmsds
238
239	double[] rmsds = new double[n+1];
240
241	for (int i = 1; i<=n; i++) {
242	String ext = new Formatter().format(".%03d",i).toString();
243	File outputPdbFile = new File(outputDir, baseName+ext+".pdb");
244	try {
245	Pdb outputPdb = tr.getStructure(i);
246	rmsds[i] = pdb.rmsd(outputPdb, "Ca");
247	}
248	catch (TinkerError e) {
249	System.err.println("Error while trying to retrieve results from Tinker: "+ e.getMessage());
250	} catch (ConformationsNotSameSizeError e) {
251	System.err.println(origPdbFile.getAbsolutePath()+" and "+outputPdbFile.getAbsolutePath()+" don't have the same conformation size, can't calculate rmsd for them.");
252	}
253	}
254
255	// write report file
256
257	try {
258	PrintWriter reportOut = new PrintWriter(new FileOutputStream(reportFile));
259	reportOut.println("run_id\tcutoff\tcutoff2\tcutoff3\tct\tct2\tct3\tnum_res" +
260	"\tresult_id\terror_val" +
261	"\tup_bound_viol\tlow_bound_viol\tmax_bound_up\tmax_bound_low\trms_bound" +
262	"\tup_viol\tlow_viol\tmax_up\tmax_low\trms_viol\trmsd_to_orig");
263
264	for (int i=1;i<=n;i++){
265	String rmsd = String.format(Locale.US,"%6.3f",rmsds[i]);
266	String errStr = String.format(Locale.US, "%6.3f", err[i]);
267	// run_id cutoff cutoff2 cutoff3 ct1 ct2 ct3 num_res
268	reportOut.println(pdbCode+"_"+pdbChainCode+"\t"+cutoff1+"\t"+cutoff2+"\t"+cutoff3+"\t"+ct1+"\t"+ct2+"\t"+ct3+"\t"+sequence.length()+"\t"+
269	// result_id error_val
270	i + "\t" + errStr + "\t" +
271	// up_bound_viol low_bound_viol max_bound_up max_bound_low rms_bound
272	nubv[i] + "\t" + nlbv[i] + "\t" + mubv[i] + "\t" + mlbv[i] + "\t" + rbv[i] + "\t" +
273	// up_viol low_viol max_up max_low rms_viol
274	nuv[i] + "\t" + nlv[i] + "\t" + muv[i] + "\t" + mlv[i] + "\t"+ rrv[i] + "\t"+
275	// rmsd_to_orig
276	rmsd);
277	}
278	reportOut.close();
279	} catch (FileNotFoundException e) {
280	System.err.println("Couldn't write to report file "+reportFile.getAbsolutePath() +". Error: "+e.getMessage());
281	}
282
283	}
284
285	}