Source code

001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.structure.symmetry.core;
022
023import org.biojava.nbio.structure.geometry.CalcPoint;
024import javax.vecmath.Point3d;
025
026import java.util.*;
027import java.util.Map.Entry;
028
029public class HelicalRepeatUnit {
030        private QuatSymmetrySubunits subunits = null;
031        private List<Point3d> repeatUnitCenters = new ArrayList<Point3d>();
032        private List<Point3d[]> repeatUnits = new ArrayList<Point3d[]>();
033        private List<List<Integer>> repeatUnitIndices = new ArrayList<List<Integer>>();
034        private Map<Integer[], Integer> interactingNeighbors = Collections.emptyMap();
035
036public HelicalRepeatUnit(QuatSymmetrySubunits subunits) {
037        this.subunits = subunits;
038}
039
040public List<Point3d> getRepeatUnitCenters() {
041        if (repeatUnitCenters.isEmpty()) {
042                run();
043        }
044        return repeatUnitCenters;
045}
046
047public List<Point3d[]> getRepeatUnits() {
048        if (repeatUnits.isEmpty()) {
049                run();
050        }
051        return repeatUnits;
052}
053
054public List<List<Integer>> getRepeatUnitIndices() {
055        return repeatUnitIndices;
056}
057
058public Map<Integer[], Integer> getInteractingRepeatUnits() {
059        if (interactingNeighbors.isEmpty()) {
060                run();
061        }
062        return interactingNeighbors;
063}
064
065private void run() {
066        this.repeatUnitCenters = calcRepeatUnitCenters();
067        if (this.repeatUnitCenters.size() == 0) {
068                return;
069        }
070        this.repeatUnits = calcRepeatUnits();
071        this.interactingNeighbors = findInteractingNeigbors();
072}
073
074private List<Point3d> calcRepeatUnitCenters() {
075
076        // TODO why do we use models here? it should not matter. Setting to 0 all
077        List<Integer> models = new ArrayList<Integer>(subunits.getSubunitCount());
078        for (int s = 0; s <subunits.getSubunitCount(); s++)
079                models.add(0);
080        Set<Integer> uniqueModels = new HashSet<Integer>(Arrays.asList(1));
081
082        int modelCount = uniqueModels.size();
083        List<Integer> folds = this.subunits.getFolds();
084        int maxFold = folds.get(folds.size()-1);
085
086        List<Point3d> repeatCenters = new ArrayList<Point3d>();
087        List<Point3d> centers = subunits.getCenters();
088
089//      if (modelCount == maxFold && subunits.getSubunitCount() > 3) {
090        if (maxFold%modelCount == 0 && modelCount > 1 && subunits.getSubunitCount() > 3) {
091//              System.out.println("calcRepeatUnitCenters case 1");
092                for (int i = 0; i < modelCount; i++) {
093                        List<Integer> subunitIndices = new ArrayList<Integer>();
094                        Point3d p = new Point3d();
095                        int count = 0;
096//                      System.out.println("Models: " + models.size());
097                        for (int j = 0; j < models.size(); j++) {
098                                if (models.get(j) == i) {
099                                        p.add(centers.get(j));
100                                        subunitIndices.add(j);
101                                        count++;
102                                }
103                        }
104//                      System.out.println("count: " + count);
105                        p.scale(1.0/count);
106//                      System.out.println("Orig Repeat unit: " + p);
107                        repeatCenters.add(p);
108                        repeatUnitIndices.add(subunitIndices);
109                }
110        } else {
111//              System.out.println("calcRepeatUnitCenters case21");
112                // TODO need to group subunits into repeating groups
113                // Case of 3B5U: A14, but seems to form (A2)*7 and symmetry related subunits don't have direct contact
114                List<Integer> sequenceClusterIds = subunits.getClusterIds();
115                for (int i = 0; i < subunits.getSubunitCount(); i++) {
116                        List<Integer> subunitIndices = new ArrayList<Integer>(1);
117                        if (sequenceClusterIds.get(i) == 0) {
118                                repeatCenters.add(new Point3d(centers.get(i)));
119//                              System.out.println("Orig Repeat unit: " + centers.get(i));
120                                subunitIndices.add(i);
121                                repeatUnitIndices.add(subunitIndices);
122                        }
123                }
124        }
125
126        // helixes should have at least 3 repeat centers
127//      System.out.println("Number of repeat centers: " + repeatCenters.size());
128        if (repeatCenters.size() < 3) {
129                repeatCenters.clear();
130        }
131
132        return repeatCenters;
133}
134
135private List<Point3d[]> calcRepeatUnits() {
136
137        // TODO why do we use models here? it should not matter. Setting to 0 all
138        List<Integer> models = new ArrayList<Integer>(
139                        subunits.getSubunitCount());
140        for (int s = 0; s < subunits.getSubunitCount(); s++)
141                models.add(0);
142        Set<Integer> uniqueModels = new HashSet<Integer>(Arrays.asList(1));
143
144        int modelCount = uniqueModels.size();
145        List<Integer> folds = this.subunits.getFolds();
146        int maxFold = folds.get(folds.size()-1);
147
148        List<Point3d[]> repeatTraces = new ArrayList<Point3d[]>();
149        List<Point3d[]> traces = subunits.getTraces();
150
151//      if (modelCount == maxFold && subunitCount > 3) {
152        if (maxFold%modelCount == 0 && modelCount > 1 && subunits.getSubunitCount() > 3) {
153                for (int i = 0; i < modelCount; i++) {
154                        List<Point3d> coords = new ArrayList<Point3d>();
155                        for (int j = 0; j < models.size(); j++) {
156                                if (models.get(j) == i) {
157                                        coords.addAll(Arrays.asList(traces.get(j)));
158                                }
159                        }
160                        Point3d[] x = new Point3d[coords.size()];
161                        coords.toArray(x);
162//                      repeatTraces.add(x); // make sure we don't accidently change the original coordinates
163                        repeatTraces.add(CalcPoint.clonePoint3dArray(x));
164                }
165        } else {
166                List<Integer> sequenceClusterIds = subunits.getClusterIds();
167                for (int i = 0; i < subunits.getSubunitCount(); i++) {
168                        if (sequenceClusterIds.get(i) == 0) {
169                                Point3d[] x = subunits.getTraces().get(i);
170                                repeatTraces.add(CalcPoint.clonePoint3dArray(x));
171                        }
172                }
173        }
174
175//      for (int i = 0; i < repeatTraces.size(); i++) {
176//              System.out.println("Repeat " + i);
177//              System.out.println(Arrays.toString(repeatTraces.get(i)));
178//      }
179        return repeatTraces;
180}
181
182private Map<Integer[], Integer> findInteractingNeigbors() {
183        Map<Integer[], Integer>  contactMap = new HashMap<Integer[], Integer>();
184
185        Map<Integer, List<Integer[]>> distanceMap = findClosestPairs(8);
186        for (List<Integer[]> pairs: distanceMap.values())
187        for (Integer[] pair: pairs) {
188                Integer contacts = calcContactNumber(repeatUnits.get(pair[0]), repeatUnits.get(pair[1]));
189//              System.out.println("contacts: " + pair[0] + "-" + pair[1] + ": " + contacts);
190                if (contacts > 0) {
191                        contactMap.put(pair, contacts);
192                }
193        }
194
195        return contactMap;
196}
197
198private Map<Integer, List<Integer[]>> findClosestPairs(int maxNeighbors) {
199        Map<Integer, List<Integer[]>>  reducedMap = new TreeMap<Integer, List<Integer[]>>();
200
201        Map<Integer, List<Integer[]>>  distanceMap = new TreeMap<Integer, List<Integer[]>>();
202        int nCenters = repeatUnitCenters.size();
203//      System.out.println("repeatUnitCenters: " + repeatUnitCenters);
204
205        for (int i = 0; i < nCenters-1; i++) {
206                for (int j = i+1; j < nCenters; j++) {
207                        float dist = (float)repeatUnitCenters.get(i).distance(repeatUnitCenters.get(j));
208                        // round to 2 digits precision
209//                      System.out.println("dist pair: " + i + "-" + j + ": " + dist);
210                        dist *= 100;
211                        int intDist = Math.round(dist);
212                        List<Integer[]> pairs = distanceMap.get(intDist);
213                        // save only one representative pair for each distance
214                        if (pairs == null) {
215                                pairs = new ArrayList<Integer[]>();
216                        }
217                        Integer[] pair = new Integer[2];
218                        pair[0] = i;
219                        pair[1] = j;
220                        pairs.add(pair);
221                        distanceMap.put(intDist, pairs);
222                }
223        }
224
225        int count = 0;
226        for (Entry<Integer, List<Integer[]>> entry: distanceMap.entrySet()) {
227                if (! (reducedMap.containsKey(entry.getKey())) ) {
228                        reducedMap.put(entry.getKey(), entry.getValue());
229//                      System.out.println("dist pair: " + entry.getValue() + ": " + entry.getKey());
230                        count++;
231                        if (count == maxNeighbors) {
232                                break;
233                        }
234                }
235        }
236        distanceMap.clear();
237
238        return reducedMap;
239}
240
241private static int calcContactNumber(Point3d[] a, Point3d[] b) {
242        int contacts = 0;
243        for (Point3d pa : a) {
244                for (Point3d pb : b) {
245                        if (pa.distance(pb) < 10) {
246                                contacts++;
247                        }
248                }
249        }
250        return contacts;
251}
252}