trunk/proteinstructure/Pdb.java

package proteinstructure;

import java.io.FileOutputStream;
import java.io.PrintStream;
import java.io.IOException;
import java.util.HashMap;
import java.util.TreeMap;
import java.util.ArrayList;
import java.util.Collections;

public class Pdb {
        
        HashMap<String,Integer> resser_atom2atomserial;
        HashMap<Integer,String> resser2restype;
        HashMap<Integer,Double[]> atomser2coord;
        HashMap<Integer,Integer> atomser2resser;
        HashMap<String,ArrayList<String>> aas2atoms = AA.getaas2atoms();
        String sequence="";
        String accode="";
        String chaincode="";
        String db;
        String chain;
        
        public Pdb (String accode, String chaincode) throws PdbaseInconsistencyError, PdbaseAcCodeNotFoundError {
                this(accode,chaincode,PdbaseInfo.pdbaseDB);
                
        }
        
        public Pdb (String accode, String chaincode, String db) throws PdbaseInconsistencyError, PdbaseAcCodeNotFoundError {
                this.accode=accode;
                this.chaincode=chaincode;
                this.db=db;
                this.chain=chaincode; // we initialise it to chaincode, in read_pdb_data_from_pdbase gets reset to the right internal chain id
                read_pdb_data_from_pdbase(db);
        }
        
        public Pdb (String pdbfile) {
                //TODO implement read_pdb_data_from_file
                //read_pdb_data_from_file(pdbfile);
        }
        
        public void read_pdb_data_from_pdbase(String db) throws PdbaseInconsistencyError, PdbaseAcCodeNotFoundError{
                resser_atom2atomserial = new HashMap<String,Integer>();
                resser2restype = new HashMap<Integer,String>();
                atomser2coord = new HashMap<Integer,Double[]>();
                atomser2resser = new HashMap<Integer,Integer>();

                PdbaseInfo mypdbaseinfo = new PdbaseInfo(accode,chaincode,db);
                ArrayList<ArrayList> resultset = mypdbaseinfo.read_atomData();
                sequence = mypdbaseinfo.read_seq();
                mypdbaseinfo.close();
                
                for (ArrayList result:resultset){
                        int atomserial = (Integer) result.get(0);
                        String atom = (String) result.get(1);
                        String res_type = (String) result.get(2);
                        chain=(String) result.get(3);
                        int res_serial = (Integer) result.get(4);
                        double x = (Double) result.get(5);
                        double y = (Double) result.get(6);
                        double z = (Double) result.get(7);
                        Double[] coords = {x, y, z};
                        ArrayList<String> aalist=AA.aas();
                        if (aalist.contains(res_type)) {
                                atomser2coord.put(atomserial, coords);
                                atomser2resser.put(atomserial, res_serial);
                                resser2restype.put(res_serial, res_type);
                                ArrayList<String> atomlist = aas2atoms.get(res_type);
                                if (atomlist.contains(atom)){
                                        resser_atom2atomserial.put(res_serial+"_"+atom, atomserial);
                                }
                        }
                }
        }
        
        public void dump2pdbfile(String outfile) throws IOException {
                String chainstr=chain;
                if (chain.equals("NULL")){
                        chainstr="A";
                }
                PrintStream Out = new PrintStream(new FileOutputStream(outfile));
                Out.println("HEADER  Dumped from "+db+". pdb accession code="+accode+", pdb chain code="+chaincode);
                for (String resser_atom:resser_atom2atomserial.keySet()){
                        int atomserial = resser_atom2atomserial.get(resser_atom);
                        int res_serial = Integer.parseInt(resser_atom.split("_")[0]);
                        String atom = resser_atom.split("_")[1];
                        String res_type = resser2restype.get(res_serial);
                        Double[] coords = atomser2coord.get(atomserial);
                        Object[] fields = {atomserial, atom, res_type, chainstr, res_serial, coords[0], coords[1], coords[2]};
                        Out.printf("ATOM  %5d  %3s %3s %1s%4d    %8.3f%8.3f%8.3f",fields);
                }
                Out.println("END");
                Out.close();
        }
        
        public void dumpseq(String seqfile) throws IOException {
                PrintStream Out = new PrintStream(new FileOutputStream(seqfile));
                Out.println(">"+accode+"_"+chaincode);
                Out.println(sequence);
                Out.close();
        }
        
        public int get_length(){
                return resser2restype.size();
        }
        
        public HashMap<Integer,Double[]> get_coords_for_ct(String ct) {
                HashMap<Integer,Double[]> coords = new HashMap<Integer,Double[]>(); 
                HashMap<String,String[]> restype2atoms = AA.ct2atoms(ct);
                for (int resser:resser2restype.keySet()){
                        String[] atoms = restype2atoms.get(resser2restype.get(resser));
                        for (String atom:atoms){
                                if (resser_atom2atomserial.containsKey(resser+"_"+atom)){
                                        int atomser = resser_atom2atomserial.get(resser+"_"+atom);
                                        Double[] coord = atomser2coord.get(atomser);
                                        coords.put(atomser, coord);
                                }
                                else if (atom.equals("O") && resser_atom2atomserial.containsKey(resser+"_"+"OXT")){
                                        int atomser = resser_atom2atomserial.get(resser+"_"+"OXT");
                                        Double[] coord = atomser2coord.get(atomser);
                                        coords.put(atomser, coord);
                                }
                                else {
                                        System.err.println("Couldn't find "+atom+" atom for resser="+resser+". Continuing without that atom for this resser.");
                                }
                        }
                }
                return coords;
        }
        
        public TreeMap<Contact, Double> calculate_dist_matrix(String ct){
                TreeMap<Contact,Double> dist_matrix = new TreeMap<Contact,Double>();
                if (!ct.contains("/")){
                        HashMap<Integer,Double[]> coords = get_coords_for_ct(ct);
                        for (int i_atomser:coords.keySet()){
                                for (int j_atomser:coords.keySet()){
                                        if (j_atomser>i_atomser) {
                                                Contact pair = new Contact(i_atomser,j_atomser);
                                                dist_matrix.put(pair, distance(coords.get(i_atomser),coords.get(j_atomser)));
                                        }
                                }
                        }
                } else {
                        String i_ct = ct.split("/")[0];
                        String j_ct = ct.split("/")[1];
                        HashMap<Integer,Double[]> i_coords = get_coords_for_ct(i_ct);
                        HashMap<Integer,Double[]> j_coords = get_coords_for_ct(j_ct);
                        for (int i_atomser:i_coords.keySet()){
                                for (int j_atomser:j_coords.keySet()){
                                        if (j_atomser!=i_atomser){
                                                Contact pair = new Contact(i_atomser,j_atomser);
                                                dist_matrix.put(pair, distance(i_coords.get(i_atomser),j_coords.get(j_atomser)));
                                        }
                                }
                        }
                }
                return dist_matrix;
        }
        
        public Graph get_graph(String ct, double cutoff){
                TreeMap<Contact,Double> dist_matrix = calculate_dist_matrix(ct);
                ArrayList<Contact> contacts = new ArrayList<Contact>();
                for (Contact pair:dist_matrix.keySet()){
                        int i_atomser=pair.i;
                        int j_atomser=pair.j;
                        if (dist_matrix.get(pair)<=cutoff){
                                int i_resser = atomser2resser.get(i_atomser);
                                int j_resser = atomser2resser.get(j_atomser);
                                Contact resser_pair = new Contact(i_resser,j_resser);
                                if (i_resser!=j_resser && (! contacts.contains(resser_pair))){
                                        contacts.add(resser_pair);
                                }
                        }
                }
                Collections.sort(contacts);
                Graph graph = new Graph (contacts,cutoff,ct);
                return graph;
        }
        
        public Double distance(Double[] coords1, Double[] coords2){
                return Math.sqrt(Math.pow((coords1[0]-coords2[0]),2)+Math.pow((coords1[1]-coords2[1]),2)+Math.pow((coords1[2]-coords2[2]),2));
        }
}
Revision:	124
Committed:	Wed May 9 11:57:24 2007 UTC (17 years, 4 months ago) by duarte
File size:	6497 byte(s)
Log Message:	Added error handling for PdbaseInfo: 2 new exception classes
Line	User	Rev	File contents
1	duarte	123	package proteinstructure;
2
3			import java.io.FileOutputStream;
4			import java.io.PrintStream;
5			import java.io.IOException;
6			import java.util.HashMap;
7			import java.util.TreeMap;
8			import java.util.ArrayList;
9			import java.util.Collections;
10
11			public class Pdb {
12
13			HashMap<String,Integer> resser_atom2atomserial;
14			HashMap<Integer,String> resser2restype;
15			HashMap<Integer,Double[]> atomser2coord;
16			HashMap<Integer,Integer> atomser2resser;
17			HashMap<String,ArrayList<String>> aas2atoms = AA.getaas2atoms();
18			String sequence="";
19			String accode="";
20			String chaincode="";
21			String db;
22			String chain;
23
24	duarte	124	public Pdb (String accode, String chaincode) throws PdbaseInconsistencyError, PdbaseAcCodeNotFoundError {
25	duarte	123	this(accode,chaincode,PdbaseInfo.pdbaseDB);
26
27			}
28
29	duarte	124	public Pdb (String accode, String chaincode, String db) throws PdbaseInconsistencyError, PdbaseAcCodeNotFoundError {
30	duarte	123	this.accode=accode;
31			this.chaincode=chaincode;
32			this.db=db;
33			this.chain=chaincode; // we initialise it to chaincode, in read_pdb_data_from_pdbase gets reset to the right internal chain id
34			read_pdb_data_from_pdbase(db);
35			}
36
37			public Pdb (String pdbfile) {
38			//TODO implement read_pdb_data_from_file
39			//read_pdb_data_from_file(pdbfile);
40			}
41
42	duarte	124	public void read_pdb_data_from_pdbase(String db) throws PdbaseInconsistencyError, PdbaseAcCodeNotFoundError{
43	duarte	123	resser_atom2atomserial = new HashMap<String,Integer>();
44			resser2restype = new HashMap<Integer,String>();
45			atomser2coord = new HashMap<Integer,Double[]>();
46			atomser2resser = new HashMap<Integer,Integer>();
47
48			PdbaseInfo mypdbaseinfo = new PdbaseInfo(accode,chaincode,db);
49			ArrayList<ArrayList> resultset = mypdbaseinfo.read_atomData();
50			sequence = mypdbaseinfo.read_seq();
51			mypdbaseinfo.close();
52
53			for (ArrayList result:resultset){
54			int atomserial = (Integer) result.get(0);
55			String atom = (String) result.get(1);
56			String res_type = (String) result.get(2);
57			chain=(String) result.get(3);
58			int res_serial = (Integer) result.get(4);
59			double x = (Double) result.get(5);
60			double y = (Double) result.get(6);
61			double z = (Double) result.get(7);
62			Double[] coords = {x, y, z};
63			ArrayList<String> aalist=AA.aas();
64			if (aalist.contains(res_type)) {
65			atomser2coord.put(atomserial, coords);
66			atomser2resser.put(atomserial, res_serial);
67			resser2restype.put(res_serial, res_type);
68			ArrayList<String> atomlist = aas2atoms.get(res_type);
69			if (atomlist.contains(atom)){
70			resser_atom2atomserial.put(res_serial+"_"+atom, atomserial);
71			}
72			}
73			}
74			}
75
76			public void dump2pdbfile(String outfile) throws IOException {
77			String chainstr=chain;
78			if (chain.equals("NULL")){
79			chainstr="A";
80			}
81			PrintStream Out = new PrintStream(new FileOutputStream(outfile));
82			Out.println("HEADER Dumped from "+db+". pdb accession code="+accode+", pdb chain code="+chaincode);
83			for (String resser_atom:resser_atom2atomserial.keySet()){
84			int atomserial = resser_atom2atomserial.get(resser_atom);
85			int res_serial = Integer.parseInt(resser_atom.split("_")[0]);
86			String atom = resser_atom.split("_")[1];
87			String res_type = resser2restype.get(res_serial);
88			Double[] coords = atomser2coord.get(atomserial);
89			Object[] fields = {atomserial, atom, res_type, chainstr, res_serial, coords[0], coords[1], coords[2]};
90			Out.printf("ATOM %5d %3s %3s %1s%4d %8.3f%8.3f%8.3f",fields);
91			}
92			Out.println("END");
93			Out.close();
94			}
95
96			public void dumpseq(String seqfile) throws IOException {
97			PrintStream Out = new PrintStream(new FileOutputStream(seqfile));
98			Out.println(">"+accode+"_"+chaincode);
99			Out.println(sequence);
100			Out.close();
101			}
102
103			public int get_length(){
104			return resser2restype.size();
105			}
106
107			public HashMap<Integer,Double[]> get_coords_for_ct(String ct) {
108			HashMap<Integer,Double[]> coords = new HashMap<Integer,Double[]>();
109			HashMap<String,String[]> restype2atoms = AA.ct2atoms(ct);
110			for (int resser:resser2restype.keySet()){
111			String[] atoms = restype2atoms.get(resser2restype.get(resser));
112			for (String atom:atoms){
113			if (resser_atom2atomserial.containsKey(resser+"_"+atom)){
114			int atomser = resser_atom2atomserial.get(resser+"_"+atom);
115			Double[] coord = atomser2coord.get(atomser);
116			coords.put(atomser, coord);
117			}
118			else if (atom.equals("O") && resser_atom2atomserial.containsKey(resser+"_"+"OXT")){
119			int atomser = resser_atom2atomserial.get(resser+"_"+"OXT");
120			Double[] coord = atomser2coord.get(atomser);
121			coords.put(atomser, coord);
122			}
123			else {
124			System.err.println("Couldn't find "+atom+" atom for resser="+resser+". Continuing without that atom for this resser.");
125			}
126			}
127			}
128			return coords;
129			}
130
131			public TreeMap<Contact, Double> calculate_dist_matrix(String ct){
132			TreeMap<Contact,Double> dist_matrix = new TreeMap<Contact,Double>();
133			if (!ct.contains("/")){
134			HashMap<Integer,Double[]> coords = get_coords_for_ct(ct);
135			for (int i_atomser:coords.keySet()){
136			for (int j_atomser:coords.keySet()){
137			if (j_atomser>i_atomser) {
138			Contact pair = new Contact(i_atomser,j_atomser);
139			dist_matrix.put(pair, distance(coords.get(i_atomser),coords.get(j_atomser)));
140			}
141			}
142			}
143			} else {
144			String i_ct = ct.split("/")[0];
145			String j_ct = ct.split("/")[1];
146			HashMap<Integer,Double[]> i_coords = get_coords_for_ct(i_ct);
147			HashMap<Integer,Double[]> j_coords = get_coords_for_ct(j_ct);
148			for (int i_atomser:i_coords.keySet()){
149			for (int j_atomser:j_coords.keySet()){
150			if (j_atomser!=i_atomser){
151			Contact pair = new Contact(i_atomser,j_atomser);
152			dist_matrix.put(pair, distance(i_coords.get(i_atomser),j_coords.get(j_atomser)));
153			}
154			}
155			}
156			}
157			return dist_matrix;
158			}
159
160			public Graph get_graph(String ct, double cutoff){
161			TreeMap<Contact,Double> dist_matrix = calculate_dist_matrix(ct);
162			ArrayList<Contact> contacts = new ArrayList<Contact>();
163			for (Contact pair:dist_matrix.keySet()){
164			int i_atomser=pair.i;
165			int j_atomser=pair.j;
166			if (dist_matrix.get(pair)<=cutoff){
167			int i_resser = atomser2resser.get(i_atomser);
168			int j_resser = atomser2resser.get(j_atomser);
169			Contact resser_pair = new Contact(i_resser,j_resser);
170			if (i_resser!=j_resser && (! contacts.contains(resser_pair))){
171			contacts.add(resser_pair);
172			}
173			}
174			}
175			Collections.sort(contacts);
176			Graph graph = new Graph (contacts,cutoff,ct);
177			return graph;
178			}
179
180			public Double distance(Double[] coords1, Double[] coords2){
181			return Math.sqrt(Math.pow((coords1[0]-coords2[0]),2)+Math.pow((coords1[1]-coords2[1]),2)+Math.pow((coords1[2]-coords2[2]),2));
182			}
183			}