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>]\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);
                        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);
                RIGraph graph2 = null;
                RIGraph graph3 = null;
                if (doublecm) {
                        graph2 = pdb.get_graph(ct2, cutoff2);           
                        if (cross) {
                                graph3 = pdb.get_graph(ct3, cutoff3);
                                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 {
                        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:	441
Committed:	Thu Nov 29 15:19:37 2007 UTC (16 years, 10 months ago) by duarte
Original Path:	*branches/aglappe-jung/reconstruct.java*
File size:	8199 byte(s)
Log Message:	Now Pdb constructors don't load data but rather intialise the pdbCode, loading occurs upon call of load(pdbChainCode,modelSerial) New methods getChains() and getModels() in all Pdb classes New Exception PdbLoadError New tester for getChains and getModels: testGetChains Changed all calls to Pdb construction accordingly (includes changing excpetions)
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	duarte	441	import proteinstructure.PdbLoadError;
15	duarte	420	import proteinstructure.RIGraph;
16	duarte	338	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>]\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	duarte	393	//double forceConstant = 100.0;
46
47	duarte	338	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	duarte	393	break;
81			//case 'f':
82			// forceConstant = Double.valueOf(g.getOptarg());
83			// break;
84	duarte	338	case 'h':
85			case '?':
86			System.out.println(help);
87			System.exit(0);
88			break; // getopt() already printed an error
89			}
90			}
91
92			if (pdbCode.equals("") \|\| pdbChainCode.equals("") \|\| ct.equals("") \|\| cutoff1==0.0 \|\| outputDir.equals("") \|\| baseName.equals("")){
93			System.err.println("Must specify at least -p, -c, -t, -d, -o and -b");
94			System.err.println(help);
95			System.exit(1);
96			}
97
98	duarte	346	if (baseName.contains(".")) {
99			System.err.println("Basename can't contain a dot (not allowed by tinker). Exiting");
100			System.exit(1);
101			}
102
103			if (!AAinfo.isValidContactType(ct) \|\| ct.contains("ALL")) {
104			System.err.println("Invalid contact type specified. Exiting");
105			System.exit(1);
106			}
107
108			if (n>999) {
109			System.err.println("Maximum number of models is 999. Specify a lower value. Exiting");
110			System.exit(1);
111			}
112
113	duarte	338	boolean doublecm = false;
114			String ct1 = ct;
115			String ct2 = ct;
116			String ct3 = null;
117			if (ct.contains("_")) {
118			ct1 = ct.split("_")[0];
119			ct2 = ct.split("_")[1];
120			doublecm = true;
121			}
122			if (cross) {
123			ct3 = ct1+"/"+ct2;
124			}
125
126			if (cutoff2==0.0) {
127			cutoff2 = cutoff1;
128			}
129			if (cutoff3==0.0) {
130			cutoff3 = cutoff1;
131			}
132
133			Pdb pdb = null;
134			try {
135	duarte	441	pdb = new PdbasePdb(pdbCode);
136			pdb.load(pdbChainCode);
137			} catch (PdbLoadError e) {
138			System.err.println("Error while loading pdb data. Specific error "+e.getMessage());
139	duarte	338	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			}
147	duarte	349	// we also write the file to the out dir so it can be used later for clustering rmsds etc.
148	duarte	367	File origPdbFile = new File (outputDir,baseName+".native.pdb");
149	duarte	349	try {
150			pdb.dump2pdbfile(origPdbFile.getAbsolutePath());
151			} catch (IOException e4) {
152			System.err.println("Couldn't write original pdb file "+origPdbFile.getAbsolutePath());
153			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)");
154			}
155
156	duarte	338	String sequence = pdb.getSequence();
157
158	duarte	420	RIGraph[] graphs =null;
159			RIGraph graph1 = pdb.get_graph(ct1, cutoff1);
160			RIGraph graph2 = null;
161			RIGraph graph3 = null;
162	duarte	338	if (doublecm) {
163	stehr	379	graph2 = pdb.get_graph(ct2, cutoff2);
164			if (cross) {
165			graph3 = pdb.get_graph(ct3, cutoff3);
166	duarte	420	graphs = new RIGraph[3];
167	stehr	379	graphs[0] = graph1;
168			graphs[1] = graph2;
169			graphs[2] = graph3;
170			} else {
171	duarte	420	graphs = new RIGraph[2];
172	stehr	379	graphs[0] = graph1;
173			graphs[1] = graph2;
174			}
175			} else {
176	duarte	420	graphs = new RIGraph[1];
177	stehr	379	graphs[0] = graph1;
178	duarte	338	}
179
180			// defining files
181	duarte	347	File reportFile = new File(outputDir,baseName+".report");
182	duarte	338
183			// creating TinkerRunner object
184	stehr	379
185	duarte	338	TinkerRunner tr = null;
186			try {
187	stehr	379	tr = new TinkerRunner(TINKERBINDIR, PRMFILE);
188	duarte	338	} catch (FileNotFoundException e3) {
189			System.err.println("Couldn't find tinker bin dir "+TINKERBINDIR+". Exiting");
190			System.exit(1);
191			}
192
193	stehr	379	// call reconstruction
194
195	duarte	338	try {
196	stehr	379	tr.reconstruct(sequence, graphs, n, outputDir, baseName, false);
197			} catch (IOException e) {
198			System.err.println("Error while running Tinker reconstruction: " + e.getMessage());
199			} catch (TinkerError e) {
200			System.err.println("Error while running Tinker reconstruction: " + e.getMessage());
201	duarte	338	}
202	stehr	379
203	duarte	368	double[] err = tr.getErrorFunctionVal();
204	duarte	338	double[] mubv = tr.getMaxUpperBoundViol();
205			double[] mlbv = tr.getMaxLowerBoundViol();
206			double[] muv = tr.getMaxUpperViol();
207			double[] mlv = tr.getMaxLowerViol();
208			int[] nubv = tr.getNumUpperBoundViol();
209			int[] nlbv = tr.getNumLowerBoundViol();
210			int[] nuv = tr.getNumUpperViol();
211			int[] nlv = tr.getNumLowerViol();
212			double[] rbv = tr.getRmsBoundViol();
213			double[] rrv = tr.getRmsRestViol();
214
215
216	stehr	379	// calculate rmsds
217	duarte	338
218			double[] rmsds = new double[n+1];
219
220	duarte	347	for (int i = 1; i<=n; i++) {
221			String ext = new Formatter().format(".%03d",i).toString();
222			File outputPdbFile = new File(outputDir, baseName+ext+".pdb");
223			try {
224	stehr	379	Pdb outputPdb = tr.getStructure(i);
225	duarte	347	rmsds[i] = pdb.rmsd(outputPdb, "Ca");
226			}
227	stehr	379	catch (TinkerError e) {
228			System.err.println("Error while trying to retrieve results from Tinker: + e.getMessage()");
229	duarte	347	} catch (ConformationsNotSameSizeError e) {
230	stehr	379	System.err.println(origPdbFile.getAbsolutePath()+" and "+outputPdbFile.getAbsolutePath()+" don't have the same conformation size, can't calculate rmsd for them.");
231	duarte	347	}
232			}
233
234	stehr	379	// write report file
235	duarte	338
236	duarte	347	try {
237			PrintWriter reportOut = new PrintWriter(new FileOutputStream(reportFile));
238	duarte	373	reportOut.println("run_id\tcutoff\tcutoff2\tcutoff3\tct\tct2\tct3\tnum_res" +
239			"\tresult_id\terror_val" +
240			"\tup_bound_viol\tlow_bound_viol\tmax_bound_up\tmax_bound_low\trms_bound" +
241	duarte	347	"\tup_viol\tlow_viol\tmax_up\tmax_low\trms_viol\trmsd_to_orig");
242
243			for (int i=1;i<=n;i++){
244	duarte	368	String rmsd = String.format(Locale.US,"%6.3f",rmsds[i]);
245			String errStr = String.format(Locale.US, "%6.3f", err[i]);
246	duarte	373	// run_id cutoff cutoff2 cutoff3 ct1 ct2 ct3 num_res
247			reportOut.println(pdbCode+"_"+pdbChainCode+"\t"+cutoff1+"\t"+cutoff2+"\t"+cutoff3+"\t"+ct1+"\t"+ct2+"\t"+ct3+"\t"+sequence.length()+"\t"+
248			// result_id error_val
249			i + "\t" + errStr + "\t" +
250			// up_bound_viol low_bound_viol max_bound_up max_bound_low rms_bound
251			nubv[i] + "\t" + nlbv[i] + "\t" + mubv[i] + "\t" + mlbv[i] + "\t" + rbv[i] + "\t" +
252			// up_viol low_viol max_up max_low rms_viol
253			nuv[i] + "\t" + nlv[i] + "\t" + muv[i] + "\t" + mlv[i] + "\t"+ rrv[i] + "\t"+
254			// rmsd_to_orig
255			rmsd);
256	duarte	347	}
257			reportOut.close();
258			} catch (FileNotFoundException e) {
259			System.err.println("Couldn't write to report file "+reportFile.getAbsolutePath() +". Error: "+e.getMessage());
260	duarte	338	}
261
262			}
263
264			}