trunk/src/writeSASAtoBFactor.java

import java.io.*;
import java.sql.SQLException;
import java.util.HashMap;

import owl.core.runners.NaccessRunner;
import owl.core.structure.Pdb;
import owl.core.structure.PdbCodeNotFoundError;
import owl.core.structure.PdbLoadError;
import owl.core.structure.PdbasePdb;
import owl.core.structure.PdbfilePdb;


/**
 * Loads a structure from a file or from PDBase and saves it as a pdb file where the b-factor column is the
 * relative per-residue solvent accessible surface area (SASA) as calculated by NACCESS. These values can be
 * displayed in e.g. Pymol with the command 'spectrum b'.
 * If a cutoff is given, the b-factor is set to 1 (exposed) or 0 (buried) depending on whether the SASA exceeds
 * the given cutoff or not. Otherwise the SASA value itself (in percent but can sometimes exceed 100) is written.
 * @author stehr
 */
public class writeSASAtoBFactor {

        private static final String NACCESS_EXECUTABLE = "/project/StruPPi/bin/naccess";
        
        /**
         * @param args
         */
        public static void main(String[] args) {
                if(args.length < 2) {
                        System.out.println(writeSASAtoBFactor.class.getName() + " <pdbcode or filename> <outfile> [<cutoff for buried/exposed>]");
                        System.exit(1);
                }
                String arg = args[0];
                String outFile = args[1];
                double cutoff = -1;
                if(args.length > 2) {
                        cutoff = Double.parseDouble(args[2]);
                }

                Pdb pdb = readStructureOrExit(arg);
                                
                // now we have a pdb structure
                try {
                        NaccessRunner naccRunner = new NaccessRunner(new File(NACCESS_EXECUTABLE), "");
                        naccRunner.runNaccess(pdb);
                        HashMap<Integer, Double> sasas = pdb.getSurfaceAccessibilities();
                        if(sasas == null) {
                                System.err.println("Error: no surface accessibility values found.");
                                System.exit(1);
                        }
                        if(cutoff >= 0) {
                                System.out.println("Discretizing values");
                                sasas = discretizeValues(sasas, cutoff);
                        }
                        pdb.setBFactorsPerResidue(sasas);
                        try {
                                System.out.println("Writing " + outFile);
                                pdb.writeToPDBFile(outFile);
                        } catch (IOException e) {
                                System.err.println("Error writing to file " + outFile + ": " + e.getMessage());
                                System.exit(1);
                        }
                } catch (IOException e) {
                        System.err.println("Error running NACCESS: " + e.getMessage());
                        System.exit(1);
                }
                
                System.out.println("done.");
        }

        /**
         * Loads a pdb structure where arg can be a pdbcode+chaincode or a pdb file name.
         * If something goes wrong, prints an error message and exits.
         * @param arg a pdbcode+chaincode (e.g. 1tdrB) or a pdb file name
         * @return the structure object
         */
        public static Pdb readStructureOrExit(String arg) {
                Pdb pdb = null;
                
                // check if argument is a filename
                File inFile = new File(arg);
                if(inFile.canRead()) {
                        System.out.println("Reading file " + arg);
                        pdb = new PdbfilePdb(arg);
                        try {
                                String[] chains = pdb.getChains();
                                System.out.println("Loading chain " + chains[0]);
                                pdb.load(chains[0]);
                        } catch (PdbLoadError e) {
                                System.err.println("Error loading file " + arg + ":" + e.getMessage());
                        }
                } else {
                        // check if argument is a pdb code
                        if(arg.length() < 4 || arg.length() > 5) {
                                System.err.println(arg + "is neither a valid file name nor a valid pdb code");
                                System.exit(1);
                        } else {
                                String pdbCode = arg.substring(0,4);
                                String chainCode = arg.substring(4,5);
                                try {
                                        System.out.println("Loading pdb code " + pdbCode);
                                        pdb = new PdbasePdb(pdbCode);
                                        if(chainCode.length() == 0) {
                                                try {
                                                        chainCode = pdb.getChains()[0];
                                                } catch (PdbLoadError e) {
                                                        System.err.println("Error loading pdb structure:" + e.getMessage());
                                                        System.exit(1);
                                                }
                                        }
                                        try {
                                                System.out.println("Loading chain " + chainCode);
                                                pdb.load(pdb.getChains()[0]);
                                        } catch (PdbLoadError e) {
                                                System.err.println("Error loading pdb structure:" + e.getMessage());
                                                System.exit(1);
                                        }
                                } catch (SQLException e) {
                                        System.err.println("Database error: " + e.getMessage());
                                        System.exit(1);
                                } catch (PdbCodeNotFoundError e) {
                                        System.err.println("Pdb code " + pdbCode + " not found in database.");
                                        System.exit(1);
                                }

                        }
                }
                return pdb;
        }
        
        /**
         * Transform all values above the threshold to 1 and all below or equal to 0. 
         * @param map
         * @param cutoff
         * @return the map with the new values
         */
        public static HashMap<Integer, Double> discretizeValues(HashMap<Integer, Double> map, double cutoff) {
                HashMap<Integer, Double> map2 = new HashMap<Integer, Double>();
                for(int key:map.keySet()) {
                        if(map.get(key) > cutoff) {
                                map2.put(key, 1.0);
                        } else {
                                map2.put(key, 0.0);
                        }
                }
                return map2; 
        }
        
}
Revision:	1005
Committed:	Wed Mar 31 12:29:26 2010 UTC (14 years, 5 months ago) by hstehr
File size:	4641 byte(s)
Log Message:	refactoring: renaming proteinstructure to structure and tools to util; moving connections,features,runners,sequence,structure,util to owl.core
Line	File contents
1	import java.io.*;
2	import java.sql.SQLException;
3	import java.util.HashMap;
4
5	import owl.core.runners.NaccessRunner;
6	import owl.core.structure.Pdb;
7	import owl.core.structure.PdbCodeNotFoundError;
8	import owl.core.structure.PdbLoadError;
9	import owl.core.structure.PdbasePdb;
10	import owl.core.structure.PdbfilePdb;
11
12
13	/**
14	* Loads a structure from a file or from PDBase and saves it as a pdb file where the b-factor column is the
15	* relative per-residue solvent accessible surface area (SASA) as calculated by NACCESS. These values can be
16	* displayed in e.g. Pymol with the command 'spectrum b'.
17	* If a cutoff is given, the b-factor is set to 1 (exposed) or 0 (buried) depending on whether the SASA exceeds
18	* the given cutoff or not. Otherwise the SASA value itself (in percent but can sometimes exceed 100) is written.
19	* @author stehr
20	*/
21	public class writeSASAtoBFactor {
22
23	private static final String NACCESS_EXECUTABLE = "/project/StruPPi/bin/naccess";
24
25	/**
26	* @param args
27	*/
28	public static void main(String[] args) {
29	if(args.length < 2) {
30	System.out.println(writeSASAtoBFactor.class.getName() + " <pdbcode or filename> <outfile> [<cutoff for buried/exposed>]");
31	System.exit(1);
32	}
33	String arg = args[0];
34	String outFile = args[1];
35	double cutoff = -1;
36	if(args.length > 2) {
37	cutoff = Double.parseDouble(args[2]);
38	}
39
40	Pdb pdb = readStructureOrExit(arg);
41
42	// now we have a pdb structure
43	try {
44	NaccessRunner naccRunner = new NaccessRunner(new File(NACCESS_EXECUTABLE), "");
45	naccRunner.runNaccess(pdb);
46	HashMap<Integer, Double> sasas = pdb.getSurfaceAccessibilities();
47	if(sasas == null) {
48	System.err.println("Error: no surface accessibility values found.");
49	System.exit(1);
50	}
51	if(cutoff >= 0) {
52	System.out.println("Discretizing values");
53	sasas = discretizeValues(sasas, cutoff);
54	}
55	pdb.setBFactorsPerResidue(sasas);
56	try {
57	System.out.println("Writing " + outFile);
58	pdb.writeToPDBFile(outFile);
59	} catch (IOException e) {
60	System.err.println("Error writing to file " + outFile + ": " + e.getMessage());
61	System.exit(1);
62	}
63	} catch (IOException e) {
64	System.err.println("Error running NACCESS: " + e.getMessage());
65	System.exit(1);
66	}
67
68	System.out.println("done.");
69	}
70
71	/**
72	* Loads a pdb structure where arg can be a pdbcode+chaincode or a pdb file name.
73	* If something goes wrong, prints an error message and exits.
74	* @param arg a pdbcode+chaincode (e.g. 1tdrB) or a pdb file name
75	* @return the structure object
76	*/
77	public static Pdb readStructureOrExit(String arg) {
78	Pdb pdb = null;
79
80	// check if argument is a filename
81	File inFile = new File(arg);
82	if(inFile.canRead()) {
83	System.out.println("Reading file " + arg);
84	pdb = new PdbfilePdb(arg);
85	try {
86	String[] chains = pdb.getChains();
87	System.out.println("Loading chain " + chains[0]);
88	pdb.load(chains[0]);
89	} catch (PdbLoadError e) {
90	System.err.println("Error loading file " + arg + ":" + e.getMessage());
91	}
92	} else {
93	// check if argument is a pdb code
94	if(arg.length() < 4 \|\| arg.length() > 5) {
95	System.err.println(arg + "is neither a valid file name nor a valid pdb code");
96	System.exit(1);
97	} else {
98	String pdbCode = arg.substring(0,4);
99	String chainCode = arg.substring(4,5);
100	try {
101	System.out.println("Loading pdb code " + pdbCode);
102	pdb = new PdbasePdb(pdbCode);
103	if(chainCode.length() == 0) {
104	try {
105	chainCode = pdb.getChains()[0];
106	} catch (PdbLoadError e) {
107	System.err.println("Error loading pdb structure:" + e.getMessage());
108	System.exit(1);
109	}
110	}
111	try {
112	System.out.println("Loading chain " + chainCode);
113	pdb.load(pdb.getChains()[0]);
114	} catch (PdbLoadError e) {
115	System.err.println("Error loading pdb structure:" + e.getMessage());
116	System.exit(1);
117	}
118	} catch (SQLException e) {
119	System.err.println("Database error: " + e.getMessage());
120	System.exit(1);
121	} catch (PdbCodeNotFoundError e) {
122	System.err.println("Pdb code " + pdbCode + " not found in database.");
123	System.exit(1);
124	}
125
126	}
127	}
128	return pdb;
129	}
130
131	/**
132	* Transform all values above the threshold to 1 and all below or equal to 0.
133	* @param map
134	* @param cutoff
135	* @return the map with the new values
136	*/
137	public static HashMap<Integer, Double> discretizeValues(HashMap<Integer, Double> map, double cutoff) {
138	HashMap<Integer, Double> map2 = new HashMap<Integer, Double>();
139	for(int key:map.keySet()) {
140	if(map.get(key) > cutoff) {
141	map2.put(key, 1.0);
142	} else {
143	map2.put(key, 0.0);
144	}
145	}
146	return map2;
147	}
148
149	}