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.Graph;
import proteinstructure.Pdb;
import proteinstructure.PdbChainCodeNotFoundError;
import proteinstructure.PdbCodeNotFoundError;
import proteinstructure.PdbaseInconsistencyError;
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>]\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 = 100.0;
                
                Getopt g = new Getopt(programName, args, "p:c:d:t:rb:o:d:D:i:n: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 '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,pdbChainCode);
                } catch (PdbaseInconsistencyError e) {
                        System.err.println("Pdbase inconsistency for structure "+pdbCode+". Can't continue, exiting");
                        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);
                } catch (PdbChainCodeNotFoundError e) {
                        System.err.println("Given pdb chain code "+pdbChainCode+" couldn't be found for pdb code "+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();

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