owl/trunk/computeEnergies.java

import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;

import proteinstructure.MaxClusterRunner;
import proteinstructure.MaxClusterRunner.ScoreType;

import tinker.TinkerError;
import tinker.TinkerRunner;

/**
 * Script to calculate energies of a set of decoys by using tinker minimize program.
 * At the moment the force field used for the energy calculation is amber99. But 
 * that could be change to any of the ones supported by tinker. See:
 * ftp://dasher.wustl.edu/pub/tinker/params
 * 
 * @author duarte
 *
 */
public class computeEnergies {

        private static final String TINKER_BIN_DIR = "/project/StruPPi/Software/tinker/bin";
        private static final String FORCEFIELD_PRM_FILE = "/project/StruPPi/Software/tinker/amber/amber99.prm";
        private static final String maxClusterExecutable = "/project/StruPPi/bin/maxcluster";
        private static final double DEFAULT_RMSGRADIENT = 1.0;
        
        private class Row {
                
                String id;
                double rmsdOriginal;
                int rmsdOrigRank;
                double rmsdMinimised;
                int rmsdMinRank;
                double gdtMinimised;
                int gdtMinRank;
                double energy;
                int energyRank;
                
                public Row() {
                        this.id = null;
                        this.rmsdOriginal = 0;
                        this.rmsdMinimised = 0;
                        this.gdtMinimised = 0;
                        this.energy = 0;
                }
                
                public Row(String id, double rmsdOriginal, double rmsdMinimised, double gdtMinimised, double energy) {
                        this.id = id;
                        this.rmsdOriginal = rmsdOriginal;
                        this.rmsdMinimised = rmsdMinimised;
                        this.gdtMinimised = gdtMinimised;
                        this.energy = energy;
                }
                
        }
        
        
        public static void main(String[] args) throws IOException {

                File listFile = new File(args[0]); // 1st argument listFile
                File nativeStruct = new File(args[1]); // 2nd argument: native structure pdb file
                double rmsGradient = DEFAULT_RMSGRADIENT;
                if (args.length>=3) rmsGradient = Double.parseDouble(args[2]); // 3rd (optional) argument: rms gradient for minimization

                // getting file names from list file
                ArrayList<File> pdbFiles = new ArrayList<File>();
                
                BufferedReader in = new BufferedReader(new FileReader(listFile));
                String line;
                int countPdbFiles=0;
                while ((line=in.readLine())!=null) {
                        // ignore empty lines or comments
                        if (line.equals("") || line.startsWith("#")) {
                                continue;
                        }
                        pdbFiles.add(new File(line.trim()));
                        countPdbFiles++;
                }
                
                TinkerRunner tr = new TinkerRunner(TINKER_BIN_DIR,FORCEFIELD_PRM_FILE);
                MaxClusterRunner maxCluster = new MaxClusterRunner(maxClusterExecutable);
                PrintWriter log = new PrintWriter(new File(listFile.getParent(),"tinker.log"));
                                
                // getting energy of native
                Row nativeRow = new computeEnergies().new Row();
                try {
                        double energy = tr.minimize(nativeStruct, rmsGradient, log);
                        nativeRow = new computeEnergies().new Row(nativeStruct.getName(), 0.0, 0.0, 0.0, energy);                       
                } catch(TinkerError e) {
                        System.err.println("Tinker error while computing energy for "+nativeStruct+". Skipping.");
                }       
                
                // getting energy of all decoys in list
                ArrayList<Row> rows = new ArrayList<Row>();
                for (File pdbFile:pdbFiles){
                        try {
                                double energy = tr.minimize(pdbFile, rmsGradient, log);
                                double rmsdOriginal = maxCluster.calculatePairwiseScore(pdbFile.getAbsolutePath(), nativeStruct.getAbsolutePath(),  ScoreType.RMSD);
                                String basename = pdbFile.getAbsolutePath();
                                if (basename.contains(".")) {
                                        basename = basename.substring(0,basename.lastIndexOf("."));
                                }
                                double rmsdMinimised = maxCluster.calculatePairwiseScore(basename+".min.pdb", nativeStruct.getAbsolutePath(), ScoreType.RMSD);
                                double gdtMinimised = maxCluster.calculatePairwiseScore(basename+".min.pdb", nativeStruct.getAbsolutePath(), ScoreType.GDT);
                                Row row = new computeEnergies().new Row(pdbFile.getName(), rmsdOriginal, rmsdMinimised, gdtMinimised, energy);
                                rows.add(row);                  
                                
                        } catch(TinkerError e) {
                                System.err.println("Tinker error while computing energy for "+pdbFile+". Skipping. Error: "+e.getMessage());
                                continue;
                        }       
                
                }               
                
                log.close();
                
                System.out.println("Done "+rows.size()+" pdb files out of "+countPdbFiles);
                
                // sorting on rmsdOrig to get rmsdOrig rank
                sortOnRmsdOrig(rows);
                int rmsdOrigRank=1;
                for (Row row: rows) {
                        row.rmsdOrigRank = rmsdOrigRank;        
                        rmsdOrigRank++;
                }

                // sorting on gdtMin to get gdt rank
                sortOnGdtMin(rows);
                int gdtRank=1;
                for (Row row: rows) {
                        row.gdtMinRank = gdtRank;       
                        gdtRank++;
                }
                
                // sorting on energy to get energy rank
                sortOnEnergy(rows);
                int energyRank=1;
                for (Row row: rows) {
                        row.energyRank = energyRank;    
                        energyRank++;
                }               

                // sorting on rmsdMin for rmsd min rank and printing out
                sortOnRmsdMin(rows);
                int rmsdRank=1;
                for (Row row: rows) {
                        row.rmsdMinRank = rmsdRank;     
                        rmsdRank++;
                }
                
                System.out.println("id\trmsdOrig\trmsdOrigRank\trmsdMin\trmsdMinRank\tgdtMin\tgdtRank\tenergy\tenergyRank");
                for (Row row: rows) {
                        System.out.printf("%s %6.3f %3d %6.3f %3d %6.3f %3d %12.4f %3d\n",row.id,row.rmsdOriginal,row.rmsdOrigRank,row.rmsdMinimised,row.rmsdMinRank,row.gdtMinimised,row.gdtMinRank,row.energy,row.energyRank);
                }
                System.out.println("Native: ");
                System.out.printf("%s %6.3f %3d %6.3f %3d %6.3f %3d %12.4f %3d\n",nativeRow.id,nativeRow.rmsdOriginal,nativeRow.rmsdOrigRank,nativeRow.rmsdMinimised,nativeRow.rmsdMinRank,nativeRow.gdtMinimised,nativeRow.gdtMinRank,nativeRow.energy,nativeRow.energyRank);
                
                
        }

        private static void sortOnRmsdOrig(ArrayList<Row> rows) {
                Collections.sort(rows, new Comparator<Row>() {
                        public int compare(Row o1, Row o2) {
                                return new Double(o1.rmsdOriginal).compareTo(o2.rmsdOriginal);
                        }
                });
        }
        
        private static void sortOnGdtMin(ArrayList<Row> rows) {
                Collections.sort(rows, new Comparator<Row>() {
                        public int compare(Row o1, Row o2) {
                                // ordered is inversed here (in gdt maximum score is best)
                                return new Double(o2.gdtMinimised).compareTo(o1.gdtMinimised);
                        }
                });
        }
        
        private static void sortOnEnergy(ArrayList<Row> rows) {
                Collections.sort(rows, new Comparator<Row>() {
                        public int compare(Row o1, Row o2) {
                                return new Double(o1.energy).compareTo(o2.energy);
                        }
                });
        }
        
        private static void sortOnRmsdMin(ArrayList<Row> rows) {
                Collections.sort(rows, new Comparator<Row>() {
                        public int compare(Row o1, Row o2) {
                                return new Double(o1.rmsdMinimised).compareTo(o2.rmsdMinimised);
                        }
                });
        }
        
}
Revision:	573
Committed:	Fri Mar 28 15:56:24 2008 UTC (16 years, 6 months ago) by duarte
File size:	6520 byte(s)
Log Message:	Fixed script: was exiting after computing native energy!
Line	User	Rev	File contents
1	duarte	572	import java.io.BufferedReader;
2			import java.io.File;
3			import java.io.FileReader;
4			import java.io.IOException;
5			import java.io.PrintWriter;
6			import java.util.ArrayList;
7			import java.util.Collections;
8			import java.util.Comparator;
9
10			import proteinstructure.MaxClusterRunner;
11			import proteinstructure.MaxClusterRunner.ScoreType;
12
13			import tinker.TinkerError;
14			import tinker.TinkerRunner;
15
16			/**
17			* Script to calculate energies of a set of decoys by using tinker minimize program.
18			* At the moment the force field used for the energy calculation is amber99. But
19			* that could be change to any of the ones supported by tinker. See:
20			* ftp://dasher.wustl.edu/pub/tinker/params
21			*
22			* @author duarte
23			*
24			*/
25			public class computeEnergies {
26
27			private static final String TINKER_BIN_DIR = "/project/StruPPi/Software/tinker/bin";
28			private static final String FORCEFIELD_PRM_FILE = "/project/StruPPi/Software/tinker/amber/amber99.prm";
29			private static final String maxClusterExecutable = "/project/StruPPi/bin/maxcluster";
30			private static final double DEFAULT_RMSGRADIENT = 1.0;
31
32			private class Row {
33
34			String id;
35			double rmsdOriginal;
36			int rmsdOrigRank;
37			double rmsdMinimised;
38			int rmsdMinRank;
39			double gdtMinimised;
40			int gdtMinRank;
41			double energy;
42			int energyRank;
43
44			public Row() {
45			this.id = null;
46			this.rmsdOriginal = 0;
47			this.rmsdMinimised = 0;
48			this.gdtMinimised = 0;
49			this.energy = 0;
50			}
51
52			public Row(String id, double rmsdOriginal, double rmsdMinimised, double gdtMinimised, double energy) {
53			this.id = id;
54			this.rmsdOriginal = rmsdOriginal;
55			this.rmsdMinimised = rmsdMinimised;
56			this.gdtMinimised = gdtMinimised;
57			this.energy = energy;
58			}
59
60			}
61
62
63			public static void main(String[] args) throws IOException {
64
65			File listFile = new File(args[0]); // 1st argument listFile
66			File nativeStruct = new File(args[1]); // 2nd argument: native structure pdb file
67			double rmsGradient = DEFAULT_RMSGRADIENT;
68			if (args.length>=3) rmsGradient = Double.parseDouble(args[2]); // 3rd (optional) argument: rms gradient for minimization
69
70			// getting file names from list file
71			ArrayList<File> pdbFiles = new ArrayList<File>();
72
73			BufferedReader in = new BufferedReader(new FileReader(listFile));
74			String line;
75			int countPdbFiles=0;
76			while ((line=in.readLine())!=null) {
77			// ignore empty lines or comments
78			if (line.equals("") \|\| line.startsWith("#")) {
79			continue;
80			}
81			pdbFiles.add(new File(line.trim()));
82			countPdbFiles++;
83			}
84
85			TinkerRunner tr = new TinkerRunner(TINKER_BIN_DIR,FORCEFIELD_PRM_FILE);
86			MaxClusterRunner maxCluster = new MaxClusterRunner(maxClusterExecutable);
87			PrintWriter log = new PrintWriter(new File(listFile.getParent(),"tinker.log"));
88
89			// getting energy of native
90			Row nativeRow = new computeEnergies().new Row();
91			try {
92			double energy = tr.minimize(nativeStruct, rmsGradient, log);
93			nativeRow = new computeEnergies().new Row(nativeStruct.getName(), 0.0, 0.0, 0.0, energy);
94			} catch(TinkerError e) {
95			System.err.println("Tinker error while computing energy for "+nativeStruct+". Skipping.");
96			}
97	duarte	573
98	duarte	572	// getting energy of all decoys in list
99			ArrayList<Row> rows = new ArrayList<Row>();
100			for (File pdbFile:pdbFiles){
101			try {
102			double energy = tr.minimize(pdbFile, rmsGradient, log);
103			double rmsdOriginal = maxCluster.calculatePairwiseScore(pdbFile.getAbsolutePath(), nativeStruct.getAbsolutePath(), ScoreType.RMSD);
104			String basename = pdbFile.getAbsolutePath();
105			if (basename.contains(".")) {
106			basename = basename.substring(0,basename.lastIndexOf("."));
107			}
108			double rmsdMinimised = maxCluster.calculatePairwiseScore(basename+".min.pdb", nativeStruct.getAbsolutePath(), ScoreType.RMSD);
109			double gdtMinimised = maxCluster.calculatePairwiseScore(basename+".min.pdb", nativeStruct.getAbsolutePath(), ScoreType.GDT);
110			Row row = new computeEnergies().new Row(pdbFile.getName(), rmsdOriginal, rmsdMinimised, gdtMinimised, energy);
111			rows.add(row);
112
113			} catch(TinkerError e) {
114			System.err.println("Tinker error while computing energy for "+pdbFile+". Skipping. Error: "+e.getMessage());
115			continue;
116			}
117
118			}
119
120			log.close();
121
122			System.out.println("Done "+rows.size()+" pdb files out of "+countPdbFiles);
123
124			// sorting on rmsdOrig to get rmsdOrig rank
125			sortOnRmsdOrig(rows);
126			int rmsdOrigRank=1;
127			for (Row row: rows) {
128			row.rmsdOrigRank = rmsdOrigRank;
129			rmsdOrigRank++;
130			}
131
132			// sorting on gdtMin to get gdt rank
133			sortOnGdtMin(rows);
134			int gdtRank=1;
135			for (Row row: rows) {
136			row.gdtMinRank = gdtRank;
137			gdtRank++;
138			}
139
140			// sorting on energy to get energy rank
141			sortOnEnergy(rows);
142			int energyRank=1;
143			for (Row row: rows) {
144			row.energyRank = energyRank;
145			energyRank++;
146			}
147
148			// sorting on rmsdMin for rmsd min rank and printing out
149			sortOnRmsdMin(rows);
150			int rmsdRank=1;
151			for (Row row: rows) {
152			row.rmsdMinRank = rmsdRank;
153			rmsdRank++;
154			}
155
156			System.out.println("id\trmsdOrig\trmsdOrigRank\trmsdMin\trmsdMinRank\tgdtMin\tgdtRank\tenergy\tenergyRank");
157			for (Row row: rows) {
158			System.out.printf("%s %6.3f %3d %6.3f %3d %6.3f %3d %12.4f %3d\n",row.id,row.rmsdOriginal,row.rmsdOrigRank,row.rmsdMinimised,row.rmsdMinRank,row.gdtMinimised,row.gdtMinRank,row.energy,row.energyRank);
159			}
160			System.out.println("Native: ");
161			System.out.printf("%s %6.3f %3d %6.3f %3d %6.3f %3d %12.4f %3d\n",nativeRow.id,nativeRow.rmsdOriginal,nativeRow.rmsdOrigRank,nativeRow.rmsdMinimised,nativeRow.rmsdMinRank,nativeRow.gdtMinimised,nativeRow.gdtMinRank,nativeRow.energy,nativeRow.energyRank);
162
163
164			}
165
166			private static void sortOnRmsdOrig(ArrayList<Row> rows) {
167			Collections.sort(rows, new Comparator<Row>() {
168			public int compare(Row o1, Row o2) {
169			return new Double(o1.rmsdOriginal).compareTo(o2.rmsdOriginal);
170			}
171			});
172			}
173
174			private static void sortOnGdtMin(ArrayList<Row> rows) {
175			Collections.sort(rows, new Comparator<Row>() {
176			public int compare(Row o1, Row o2) {
177			// ordered is inversed here (in gdt maximum score is best)
178			return new Double(o2.gdtMinimised).compareTo(o1.gdtMinimised);
179			}
180			});
181			}
182
183			private static void sortOnEnergy(ArrayList<Row> rows) {
184			Collections.sort(rows, new Comparator<Row>() {
185			public int compare(Row o1, Row o2) {
186			return new Double(o1.energy).compareTo(o2.energy);
187			}
188			});
189			}
190
191			private static void sortOnRmsdMin(ArrayList<Row> rows) {
192			Collections.sort(rows, new Comparator<Row>() {
193			public int compare(Row o1, Row o2) {
194			return new Double(o1.rmsdMinimised).compareTo(o2.rmsdMinimised);
195			}
196			});
197			}
198
199			}