owl/trunk/reconstruct.java

import gnu.getopt.Getopt;

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

import proteinstructure.ConformationsNotSameSizeError;
import proteinstructure.Graph;
import proteinstructure.Pdb;
import proteinstructure.PdbChainCodeNotFoundError;
import proteinstructure.PdbCodeNotFoundError;
import proteinstructure.PdbaseInconsistencyError;
import proteinstructure.PdbasePdb;
import proteinstructure.PdbfileFormatError;
import proteinstructure.PdbfilePdb;

import tinker.ConstraintsMaker;
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";
        private static final String PRMTYPE = "amber";
        
        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;

                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 '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);
                }
                
                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);
                }
                String sequence = pdb.getSequence();

                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);
                }

                
                // defining files
                File logFile = new File(outputDir,"tinker.log");
                File prmFile = new File(PRMFILE);
                File xyzFile = new File(outputDir,baseName+".xyz");
                File seqFile = new File(outputDir,baseName+".seq");
                File pdbFile = new File(outputDir,baseName+".pdb");
                File keyFile = new File(outputDir,baseName+".key");
                
                // creating TinkerRunner object
                TinkerRunner tr = null;
                try {
                        tr = new TinkerRunner(TINKERBINDIR,PRMFILE,logFile);
                } catch (FileNotFoundException e3) {
                        System.err.println("Couldn't find tinker bin dir "+TINKERBINDIR+". Exiting");
                        System.exit(1);
                }
                
                // 1. run tinker's protein program      
                try {
                        tr.runProtein(sequence, outputDir, baseName );
                } catch (IOException e2) {
                        System.err.println("Couldn't read file to run 'protein' "+seqFile.getAbsolutePath());
                        System.err.println("Exiting");
                        System.exit(1);                 
                } catch (TinkerError e2) {
                        System.err.println("Tinker error while running 'protein', check log file "+logFile.getAbsolutePath()+". Exiting");
                        System.exit(1);                 
                }

                // 1a. convert xyz file to pdb to be able to map atom serials after
                try {
                        tr.runXyzpdb(xyzFile, seqFile, pdbFile);
                } catch (IOException e1) {
                        System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+" or "+seqFile.getAbsolutePath()+" or write to "+pdbFile.getAbsolutePath()+" for running 'xyzpdb'");
                        System.err.println("Exiting");
                        System.exit(1);
                } catch (TinkerError e1) {
                        System.err.println("Tinker error while running xyzpdb, check log file "+logFile.getAbsolutePath()+". Exiting");
                        System.exit(1);
                }

                // 2. creating constraints into key file
                ConstraintsMaker cm = null;
                try {
                        cm = new ConstraintsMaker(pdbFile,xyzFile,prmFile,PRMTYPE,keyFile);
                } catch (IOException e3) {
                        System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+", "+pdbFile.getAbsolutePath()+" or, "+prmFile.getAbsolutePath()+" write to "+keyFile.getAbsolutePath()+" for creating distance constraints");
                        System.err.println("Exiting");
                        System.exit(1);
                } catch (PdbfileFormatError e3) {
                        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");
                        System.exit(1);
                } 
                try {
                        cm.createConstraints(graph1);
                        if (doublecm) cm.createConstraints(graph2);
                        if (cross) cm.createConstraints(graph3);
                } catch (Exception e2) {
                        System.err.println("Invalid contact type for writing constraints.");
                        System.err.println("Error: "+e2.getMessage());
                        System.err.println("Exiting");
                        System.exit(1);
                }
                cm.closeKeyFile();

                // 3. run tinker's distgeom
                try {
                        System.out.println("Running distgeom...");
                        tr.runDistgeom(xyzFile, outputDir, baseName, n);
                } catch (TinkerError e1) {
                        System.err.println("Tinker error while running 'distgeom', check log file "+logFile.getAbsolutePath()+". Exiting");
                        System.exit(1);                 
                } catch (IOException e1) {
                        System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+" or write 'distgeom' output files to output dir "+outputDir);
                        System.err.println("Exiting");
                        System.exit(1);
                }
                
                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();
                

                // 4. converting xyz output files to pdb files and calculating rmsds

                double[] rmsds = new double[n+1];               
                
                        for (int i = 1; i<=n; i++) {
                                String ext = new Formatter().format(".%03d",i).toString();
                                File outputXyzFile = new File(outputDir, baseName+ext);
                                File outputPdbFile = new File(outputDir, baseName+ext+".pdb");
                                try {
                                        tr.runXyzpdb(outputXyzFile, seqFile, outputPdbFile);
                                        
                                        Pdb outputPdb = new PdbfilePdb(outputPdbFile.getAbsolutePath(),"NULL");
                                        rmsds[i] = pdb.rmsd(outputPdb, "Ca");

                                } catch (IOException e) {
                                        System.err.println("Couldn't read file "+outputXyzFile.getAbsolutePath()+", or "+seqFile.getAbsolutePath()+", or write to "+outputPdbFile.getAbsolutePath()+" while converting with 'xyzpdb'");
                                        System.err.println("Can't calculate rmsd for it");
                                } catch (TinkerError e) {
                                        System.err.println("Tinker error while running 'xyzpdb' to convert"+outputXyzFile.getAbsolutePath()+", check log file "+logFile.getAbsolutePath());
                                        System.err.println("Can't calculate rmsd for it");
                                }
                                catch (PdbfileFormatError e) {
                                        System.err.println("Output pdb file "+outputPdbFile.getAbsolutePath()+" doesn't seem to be in the correcet format. Can't calculate rmsd for it");
                                } catch (PdbChainCodeNotFoundError e) {
                                        // this shouldn't happen, chain code is hard coded, we throw stack trace and continue if it happens
                                        e.printStackTrace();
                                } catch (ConformationsNotSameSizeError e) {
                                        System.err.println(pdbFile.getAbsolutePath()+" and "+outputPdbFile.getAbsolutePath()+" don't have the same conformation size, can't calculate rmsd for them.");
                                }                               
                                
                                tr.closeLog();
                        }                                       

                
                // 6. report
                System.out.println("accession_code\tchain_pdb_code\tcutoff\tcutoff2\tcutoff3\tct\tct2\tct3" +
                                "\tresult_id\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++){
                        System.out.println(pdbCode+"\t"+pdbChainCode+"\t"+cutoff1+"\t"+cutoff2+"\t"+cutoff3+"\t"+ct1+"\t"+ct2+"\t"+ct3+"\t"+
                                        i+nubv[i]+"\t"+nlbv[i]+"\t"+mubv[i]+"\t"+mlbv[i]+"\t"+muv[i]+"\t"+mlv[i]+"\t"+rbv[i]+"\t"+
                                        "\t"+nuv[i]+"\t"+nlv[i]+"\t"+rrv[i]+"\t"+rmsds[i]);
                }
                
        }

}
Revision:	338
Committed:	Fri Oct 5 11:56:28 2007 UTC (16 years, 11 months ago) by duarte
File size:	9727 byte(s)
Log Message:	Initial commit of reconstruct script. Works only sometimes!! Some others doesn't run distgeom without apparent reason! Throwing less exceptions 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			import java.io.IOException;
6			import java.sql.SQLException;
7			import java.util.Formatter;
8
9			import proteinstructure.ConformationsNotSameSizeError;
10			import proteinstructure.Graph;
11			import proteinstructure.Pdb;
12			import proteinstructure.PdbChainCodeNotFoundError;
13			import proteinstructure.PdbCodeNotFoundError;
14			import proteinstructure.PdbaseInconsistencyError;
15			import proteinstructure.PdbasePdb;
16			import proteinstructure.PdbfileFormatError;
17			import proteinstructure.PdbfilePdb;
18
19			import tinker.ConstraintsMaker;
20			import tinker.TinkerError;
21			import tinker.TinkerRunner;
22
23
24			public class reconstruct {
25
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			private static final String PRMTYPE = "amber";
30
31			public static void main(String[] args) {
32
33			String programName = reconstruct.class.getName();
34			String help = "Usage:\n" +
35			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";
36
37			String pdbCode = "";
38			String pdbChainCode = "";
39			String ct = "";
40			double cutoff1 = 0.0;
41			double cutoff2 = 0.0;
42			double cutoff3 = 0.0;
43			String outputDir = "";
44			String baseName = "";
45			boolean cross = false;
46			int n = 1;
47
48			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			break;
82			case 'h':
83			case '?':
84			System.out.println(help);
85			System.exit(0);
86			break; // getopt() already printed an error
87			}
88			}
89
90			if (pdbCode.equals("") \|\| pdbChainCode.equals("") \|\| ct.equals("") \|\| cutoff1==0.0 \|\| outputDir.equals("") \|\| baseName.equals("")){
91			System.err.println("Must specify at least -p, -c, -t, -d, -o and -b");
92			System.err.println(help);
93			System.exit(1);
94			}
95
96			boolean doublecm = false;
97			String ct1 = ct;
98			String ct2 = ct;
99			String ct3 = null;
100			if (ct.contains("_")) {
101			ct1 = ct.split("_")[0];
102			ct2 = ct.split("_")[1];
103			doublecm = true;
104			}
105			if (cross) {
106			ct3 = ct1+"/"+ct2;
107			}
108
109			if (cutoff2==0.0) {
110			cutoff2 = cutoff1;
111			}
112			if (cutoff3==0.0) {
113			cutoff3 = cutoff1;
114			}
115
116
117			Pdb pdb = null;
118			try {
119			pdb = new PdbasePdb(pdbCode,pdbChainCode);
120			} catch (PdbaseInconsistencyError e) {
121			System.err.println("Pdbase inconsistency for structure "+pdbCode+". Can't continue, exiting");
122			System.exit(1);
123			} catch (PdbCodeNotFoundError e) {
124			System.err.println("Given pdb code "+pdbCode+" couldn't be found in pdbase. Exiting");
125			System.exit(1);
126			} catch (SQLException e) {
127			System.err.println("Problems connecting to database for getting pdb data for "+pdbCode+". Exiting");
128			System.exit(1);
129			} catch (PdbChainCodeNotFoundError e) {
130			System.err.println("Given pdb chain code "+pdbChainCode+" couldn't be found for pdb code "+pdbCode+". Exiting");
131			System.exit(1);
132			}
133			String sequence = pdb.getSequence();
134
135			Graph graph1 = pdb.get_graph(ct1, cutoff1);
136			Graph graph2 = null;
137			Graph graph3 = null;
138			if (doublecm) {
139			graph2 = pdb.get_graph(ct2, cutoff2);
140			}
141			if (cross) {
142			graph3 = pdb.get_graph(ct3, cutoff3);
143			}
144
145
146			// defining files
147			File logFile = new File(outputDir,"tinker.log");
148			File prmFile = new File(PRMFILE);
149			File xyzFile = new File(outputDir,baseName+".xyz");
150			File seqFile = new File(outputDir,baseName+".seq");
151			File pdbFile = new File(outputDir,baseName+".pdb");
152			File keyFile = new File(outputDir,baseName+".key");
153
154			// creating TinkerRunner object
155			TinkerRunner tr = null;
156			try {
157			tr = new TinkerRunner(TINKERBINDIR,PRMFILE,logFile);
158			} catch (FileNotFoundException e3) {
159			System.err.println("Couldn't find tinker bin dir "+TINKERBINDIR+". Exiting");
160			System.exit(1);
161			}
162
163			// 1. run tinker's protein program
164			try {
165			tr.runProtein(sequence, outputDir, baseName );
166			} catch (IOException e2) {
167			System.err.println("Couldn't read file to run 'protein' "+seqFile.getAbsolutePath());
168			System.err.println("Exiting");
169			System.exit(1);
170			} catch (TinkerError e2) {
171			System.err.println("Tinker error while running 'protein', check log file "+logFile.getAbsolutePath()+". Exiting");
172			System.exit(1);
173			}
174
175			// 1a. convert xyz file to pdb to be able to map atom serials after
176			try {
177			tr.runXyzpdb(xyzFile, seqFile, pdbFile);
178			} catch (IOException e1) {
179			System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+" or "+seqFile.getAbsolutePath()+" or write to "+pdbFile.getAbsolutePath()+" for running 'xyzpdb'");
180			System.err.println("Exiting");
181			System.exit(1);
182			} catch (TinkerError e1) {
183			System.err.println("Tinker error while running xyzpdb, check log file "+logFile.getAbsolutePath()+". Exiting");
184			System.exit(1);
185			}
186
187			// 2. creating constraints into key file
188			ConstraintsMaker cm = null;
189			try {
190			cm = new ConstraintsMaker(pdbFile,xyzFile,prmFile,PRMTYPE,keyFile);
191			} catch (IOException e3) {
192			System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+", "+pdbFile.getAbsolutePath()+" or, "+prmFile.getAbsolutePath()+" write to "+keyFile.getAbsolutePath()+" for creating distance constraints");
193			System.err.println("Exiting");
194			System.exit(1);
195			} catch (PdbfileFormatError e3) {
196			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");
197			System.exit(1);
198			}
199			try {
200			cm.createConstraints(graph1);
201			if (doublecm) cm.createConstraints(graph2);
202			if (cross) cm.createConstraints(graph3);
203			} catch (Exception e2) {
204			System.err.println("Invalid contact type for writing constraints.");
205			System.err.println("Error: "+e2.getMessage());
206			System.err.println("Exiting");
207			System.exit(1);
208			}
209			cm.closeKeyFile();
210
211			// 3. run tinker's distgeom
212			try {
213			System.out.println("Running distgeom...");
214			tr.runDistgeom(xyzFile, outputDir, baseName, n);
215			} catch (TinkerError e1) {
216			System.err.println("Tinker error while running 'distgeom', check log file "+logFile.getAbsolutePath()+". Exiting");
217			System.exit(1);
218			} catch (IOException e1) {
219			System.err.println("Couldn't read files "+xyzFile.getAbsolutePath()+" or write 'distgeom' output files to output dir "+outputDir);
220			System.err.println("Exiting");
221			System.exit(1);
222			}
223
224			double[] mubv = tr.getMaxUpperBoundViol();
225			double[] mlbv = tr.getMaxLowerBoundViol();
226			double[] muv = tr.getMaxUpperViol();
227			double[] mlv = tr.getMaxLowerViol();
228			int[] nubv = tr.getNumUpperBoundViol();
229			int[] nlbv = tr.getNumLowerBoundViol();
230			int[] nuv = tr.getNumUpperViol();
231			int[] nlv = tr.getNumLowerViol();
232			double[] rbv = tr.getRmsBoundViol();
233			double[] rrv = tr.getRmsRestViol();
234
235
236			// 4. converting xyz output files to pdb files and calculating rmsds
237
238			double[] rmsds = new double[n+1];
239
240			for (int i = 1; i<=n; i++) {
241			String ext = new Formatter().format(".%03d",i).toString();
242			File outputXyzFile = new File(outputDir, baseName+ext);
243			File outputPdbFile = new File(outputDir, baseName+ext+".pdb");
244			try {
245			tr.runXyzpdb(outputXyzFile, seqFile, outputPdbFile);
246
247			Pdb outputPdb = new PdbfilePdb(outputPdbFile.getAbsolutePath(),"NULL");
248			rmsds[i] = pdb.rmsd(outputPdb, "Ca");
249
250			} catch (IOException e) {
251			System.err.println("Couldn't read file "+outputXyzFile.getAbsolutePath()+", or "+seqFile.getAbsolutePath()+", or write to "+outputPdbFile.getAbsolutePath()+" while converting with 'xyzpdb'");
252			System.err.println("Can't calculate rmsd for it");
253			} catch (TinkerError e) {
254			System.err.println("Tinker error while running 'xyzpdb' to convert"+outputXyzFile.getAbsolutePath()+", check log file "+logFile.getAbsolutePath());
255			System.err.println("Can't calculate rmsd for it");
256			}
257			catch (PdbfileFormatError e) {
258			System.err.println("Output pdb file "+outputPdbFile.getAbsolutePath()+" doesn't seem to be in the correcet format. Can't calculate rmsd for it");
259			} catch (PdbChainCodeNotFoundError e) {
260			// this shouldn't happen, chain code is hard coded, we throw stack trace and continue if it happens
261			e.printStackTrace();
262			} catch (ConformationsNotSameSizeError e) {
263			System.err.println(pdbFile.getAbsolutePath()+" and "+outputPdbFile.getAbsolutePath()+" don't have the same conformation size, can't calculate rmsd for them.");
264			}
265
266			tr.closeLog();
267			}
268
269
270			// 6. report
271			System.out.println("accession_code\tchain_pdb_code\tcutoff\tcutoff2\tcutoff3\tct\tct2\tct3" +
272			"\tresult_id\tup_bound_viol\tlow_bound_viol\tmax_bound_up\tmax_bound_low\trms_bound" +
273			"\tup_viol\tlow_viol\tmax_up\tmax_low\trms_viol\trmsd_to_orig");
274
275			for (int i=1;i<=n;i++){
276			System.out.println(pdbCode+"\t"+pdbChainCode+"\t"+cutoff1+"\t"+cutoff2+"\t"+cutoff3+"\t"+ct1+"\t"+ct2+"\t"+ct3+"\t"+
277			i+nubv[i]+"\t"+nlbv[i]+"\t"+mubv[i]+"\t"+mlbv[i]+"\t"+muv[i]+"\t"+mlv[i]+"\t"+rbv[i]+"\t"+
278			"\t"+nuv[i]+"\t"+nlv[i]+"\t"+rrv[i]+"\t"+rmsds[i]);
279			}
280
281			}
282
283			}