Skip to content

Commit 651d90b

Browse files
josemduartejosemduarte
committed
Full implementation of ASA calculation including a unit test
1 parent e68c0e8 commit 651d90b

5 files changed

Lines changed: 302 additions & 62 deletions

File tree

‎biojava3-structure/src/main/java/demo/DemoAsa.java‎

Lines changed: 30 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,12 @@
11
package demo;
22

33
import java.io.IOException;
4-
import java.util.ArrayList;
5-
import java.util.List;
64

7-
import org.biojava.bio.structure.AsaCalculator;
8-
import org.biojava.bio.structure.Atom;
9-
import org.biojava.bio.structure.Chain;
10-
import org.biojava.bio.structure.Element;
11-
import org.biojava.bio.structure.Group;
12-
import org.biojava.bio.structure.GroupType;
135
import org.biojava.bio.structure.Structure;
146
import org.biojava.bio.structure.StructureException;
15-
import org.biojava.bio.structure.StructureTools;
167
import org.biojava.bio.structure.align.util.AtomCache;
8+
import org.biojava.bio.structure.asa.AsaCalculator;
9+
import org.biojava.bio.structure.asa.GroupAsa;
1710
import org.biojava3.structure.StructureIO;
1811

1912
public class DemoAsa {
@@ -43,44 +36,24 @@ private static void demoAsa(String pdbCode, int numThreads) throws IOException,
4336
AsaCalculator.DEFAULT_PROBE_SIZE,
4437
1000, numThreads, hetAtoms);
4538

46-
double[] asas = asaCalc.calculateAsa();
39+
GroupAsa[] groupAsas = asaCalc.getGroupAsas();
4740

4841
long end = System.currentTimeMillis();
4942

5043

5144
double tot = 0;
5245

53-
//int i = 0;
5446

55-
Atom[] atoms = StructureTools.getAllNonHAtomArray(structure, hetAtoms);
5647

57-
for (int i=0;i<atoms.length;i++) {
58-
tot+=asas[i];
48+
for (GroupAsa groupAsa: groupAsas) {
49+
System.out.printf("%1s\t%5s\t%3s\t%6.2f\n",
50+
groupAsa.getGroup().getChainId(),
51+
groupAsa.getGroup().getResidueNumber(),
52+
groupAsa.getGroup().getPDBName(),
53+
groupAsa.getAsaU());
54+
tot+=groupAsa.getAsaU();
5955
}
6056

61-
// for (Chain chain:structure.getChains()) {
62-
// for (Group res:chain.getAtomGroups()) {
63-
// //List<Group> alllocs = new ArrayList<Group>();// = alt.getAltLocs();
64-
// //alllocs.add(alt);
65-
// //for(Group res: alllocs) {
66-
// if (!hetAtoms && res.getType().equals(GroupType.HETATM)) continue;
67-
//
68-
// //System.out.printf("%1s\t%3s\t%s\n",chain.getChainID(),res.getResidueNumber().toString(),res.getPDBName());
69-
// double perResidue = 0;
70-
// for (Atom atom:res.getAtoms()) {
71-
// if (atom.getElement()!=Element.H) {
72-
// System.out.println(" "+atom.getName()+" "+atom.getAltLoc()+ " "+atom.getPDBserial());
73-
// perResidue += asas[i];
74-
// i++;
75-
// }
76-
//
77-
// }
78-
// System.out.printf("%1s\t%3d\t%s\t%6.2f\n",chain.getChainID(),res.getResidueNumber().getSeqNum(),res.getPDBName(),perResidue);
79-
// tot+=perResidue;
80-
// //}
81-
// }
82-
// }
83-
8457

8558
System.out.printf("Total area: %9.2f\n",tot);
8659
System.out.printf("Time: %4.1fs\n",((end-start)/1000.0));
@@ -95,7 +68,8 @@ private static void demoAsa(String pdbCode, int numThreads) throws IOException,
9568
AsaCalculator.DEFAULT_PROBE_SIZE,
9669
1000, numThreads, hetAtoms);
9770

98-
asas = asaCalc.calculateAsa();
71+
// only calculating all atom ASAs without keeping the returned value
72+
asaCalc.calculateAsas();
9973

10074
end = System.currentTimeMillis();
10175
runTimes[nThreads-1] = (end-start)/1000.0;
@@ -107,6 +81,24 @@ private static void demoAsa(String pdbCode, int numThreads) throws IOException,
10781

10882

10983

84+
// for (Chain chain:structure.getChains()) {
85+
// for (Group alt:chain.getAtomGroups()) {
86+
// List<Group> alllocs = alt.getAltLocs();
87+
// //alllocs.add(alt);
88+
// for(Group res: alllocs) {
89+
// if (!hetAtoms && res.getType().equals(GroupType.HETATM)) continue;
90+
//
91+
// System.out.printf("%1s\t%3s\t%s\n",chain.getChainID(),res.getResidueNumber().toString(),res.getPDBName());
92+
// for (Atom atom:res.getAtoms()) {
93+
// if (atom.getElement()!=Element.H) {
94+
// System.out.println(" "+atom.getName()+"\t"+atom.getAltLoc()+ "\t"+atom.getPDBserial());
95+
// }
96+
//
97+
// }
98+
// }
99+
// }
100+
// }
101+
110102
}
111103

112104

‎biojava3-structure/src/main/java/org/biojava/bio/structure/StructureTools.java‎

Lines changed: 0 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -297,23 +297,6 @@ public static final Atom[] getAllNonHAtomArray(Structure s, boolean hetAtoms) {
297297
List<Atom> atoms = new ArrayList<Atom>();
298298

299299

300-
// for (Chain chain:s.getChains()) {
301-
// for (Group alt:chain.getAtomGroups()) {
302-
// if (!hetAtoms && alt.getType().equals(GroupType.HETATM)) continue;
303-
// List<Group> alllocs = alt.getAltLocs();
304-
// alllocs.add(alt);
305-
// for(Group res: alllocs) {
306-
//
307-
// for (Atom atom:res.getAtoms()) {
308-
// if (atom.getAltLoc()==' ' || atom.getAltLoc()=='A'){
309-
// atoms.add(atom);
310-
// }
311-
// }
312-
// }
313-
// }
314-
// }
315-
316-
317300
AtomIterator iter = new AtomIterator(s);
318301
while (iter.hasNext()){
319302
Atom a = iter.next();

biojava3-structure/src/main/java/org/biojava/bio/structure/AsaCalculator.java renamed to biojava3-structure/src/main/java/org/biojava/bio/structure/asa/AsaCalculator.java

Lines changed: 50 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,24 @@
1-
package org.biojava.bio.structure;
1+
package org.biojava.bio.structure.asa;
22

33
import java.util.ArrayList;
4+
import java.util.TreeMap;
45
import java.util.concurrent.ExecutorService;
56
import java.util.concurrent.Executors;
67

78
import javax.vecmath.Point3d;
89

10+
import org.biojava.bio.structure.AminoAcid;
11+
import org.biojava.bio.structure.Atom;
12+
import org.biojava.bio.structure.Calc;
13+
import org.biojava.bio.structure.Element;
14+
import org.biojava.bio.structure.Group;
15+
import org.biojava.bio.structure.GroupType;
16+
import org.biojava.bio.structure.NucleotideImpl;
17+
import org.biojava.bio.structure.ResidueNumber;
18+
import org.biojava.bio.structure.Structure;
19+
import org.biojava.bio.structure.StructureException;
20+
import org.biojava.bio.structure.StructureTools;
21+
922

1023

1124

@@ -76,14 +89,16 @@ public void run() {
7689
private double cons;
7790

7891
/**
79-
* Constructs a new AsaCalculator. Subsequently call {@link #calculateAsa()}
80-
* to calculate the ASAs.
92+
* Constructs a new AsaCalculator. Subsequently call {@link #calculateAsas()}
93+
* or {@link #getGroupAsas()} to calculate the ASAs
94+
* Only non-Hydrogen atoms are considered in the calculation.
8195
* @param structure
8296
* @param probe
8397
* @param nSpherePoints
8498
* @param nThreads
8599
* @param hetAtoms if true HET residues are considered, if false they aren't, equivalent to
86100
* NACCESS' -h option
101+
* @see StructureTools.getAllNonHAtomArray
87102
*/
88103
public AsaCalculator(Structure structure, double probe, int nSpherePoints, int nThreads, boolean hetAtoms) {
89104
this.atoms = StructureTools.getAllNonHAtomArray(structure, hetAtoms);
@@ -103,8 +118,8 @@ public AsaCalculator(Structure structure, double probe, int nSpherePoints, int n
103118
}
104119

105120
/**
106-
* Constructs a new AsaCalculator. Subsequently call {@link #calculateAsa()}
107-
* to calculate the ASAs.
121+
* Constructs a new AsaCalculator. Subsequently call {@link #calculateAsas()}
122+
* or {@link #getGroupAsas()} to calculate the ASAs.
108123
* @param atoms an array of atoms not containing Hydrogen atoms
109124
* @param probe the probe size
110125
* @param nSpherePoints the number of points to be used in generating the spherical
@@ -134,13 +149,41 @@ public AsaCalculator(Atom[] atoms, double probe, int nSpherePoints, int nThreads
134149
cons = 4.0 * Math.PI / (double)nSpherePoints;
135150
}
136151

152+
/**
153+
* Calculates ASA for all atoms and return them as a GroupAsa
154+
* array (one element per residue in structure) containing ASAs per residue
155+
* and per atom.
156+
* The sorting of Groups in returned array is as specified by {@link org.biojava.bio.structure.ResidueNumber}
157+
* @return
158+
*/
159+
public GroupAsa[] getGroupAsas() {
160+
161+
TreeMap<ResidueNumber, GroupAsa> asas = new TreeMap<ResidueNumber, GroupAsa>();
162+
163+
double[] asasPerAtom = calculateAsas();
164+
165+
for (int i=0;i<atoms.length;i++) {
166+
Group g = atoms[i].getGroup();
167+
if (!asas.containsKey(g.getResidueNumber())) {
168+
GroupAsa groupAsa = new GroupAsa(g);
169+
groupAsa.addAtomAsaU(asasPerAtom[i]);
170+
asas.put(g.getResidueNumber(), groupAsa);
171+
} else {
172+
GroupAsa groupAsa = asas.get(g.getResidueNumber());
173+
groupAsa.addAtomAsaU(asasPerAtom[i]);
174+
}
175+
}
176+
177+
return (GroupAsa[]) asas.values().toArray(new GroupAsa[asas.size()]);
178+
}
179+
137180
/**
138181
* Calculates the Accessible Surface Areas for the atoms given in constructor and with parameters given.
139182
* Beware that the parallel implementation is quite memory hungry. It scales well as long as there is
140183
* enough memory available.
141184
* @return an array with asa values corresponding to each atom of the input array
142185
*/
143-
public double[] calculateAsa() {
186+
public double[] calculateAsas() {
144187

145188
double[] asas = new double[atoms.length];
146189

@@ -333,7 +376,7 @@ else if (atomCode.equals("CA") || atomCode.equals("CB") ||
333376
atomCode.equals("CG1") || atomCode.equals("CG2")) {
334377
return TETRAHEDRAL_CARBON_VDW; // tetrahedral Carbon
335378
}
336-
// left cases depend on amino acid: CD, CD1, CD2, CG
379+
// the rest of the cases (CD, CD1, CD2, CG) depend on amino acid
337380
else {
338381
switch (aa) {
339382
case 'F':

0 commit comments

Comments
 (0)