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.quaternary;
022
023import org.biojava.nbio.structure.*;
024
025import java.util.ArrayList;
026import java.util.Arrays;
027import java.util.List;
028
029import javax.vecmath.Point3d;
030
031
032
033
034
035/**
036 * @author Peter Rose
037 *
038 *
039 */
040public class BioAssemblyTools {
041
042
043        /**
044         * Checks if the passed in expression is a unary operator expression
045         * Example: (1,2,3) or (1-60) are unary operator expressions
046         *          (1-60)(61-88) is a binary operator expression, representing
047         *          a cartesian product of the two parenthesised lists
048         *
049         * @param expression
050         * @return true if expression is a unary operator expression
051         */
052        public static boolean isUnaryExpression(String expression) {
053                int first = expression.indexOf("(");
054                int last = expression.lastIndexOf("(");
055                if (first < 0 || last < 0) {
056                        return true;
057                }
058                return ! (first == 0 && last > first);
059        }
060
061        public static List<String> parseUnaryOperatorExpression(String operatorExpression) throws IllegalArgumentException {
062                return parseSubExpression(operatorExpression);
063        }
064
065        private static List<String> parseSubExpression(String expression) throws IllegalArgumentException {
066                // remove parenthesis, if any
067                String tmp = expression.replace("(", "");
068                tmp = tmp.replace(")", "");
069
070                // separate the operators
071                List<String> components = null;
072                try {
073                        components = Arrays.asList(tmp.split(","));
074                } catch (Exception e) {
075                        throw new IllegalArgumentException("Invalid oper_expression: " + expression);
076                }
077
078                // expand ranges if present, i.e. 1-60 -> 1, 2, 3, ..., 60
079                List<String> operators = new ArrayList<String>();
080                for (String component : components) {
081                        if (component.contains("-")) {
082                                operators.addAll(expandRange(component));
083                        } else {
084                                operators.add(component);
085                        }
086                }
087                return operators;
088        }
089
090        /**
091         * Expands a range expression, i.e. (1-6) to a list 1,2,3,4,5,6
092         * @param expression the expression to be expanded
093         * @return list of items in range
094         * @throws IllegalArgumentException
095         */
096        private static List<String> expandRange(String expression) throws IllegalArgumentException {
097                int first = 0;
098                int last = 0;
099                try {
100                        String[] range = expression.split("-");
101                        first = Integer.parseInt(range[0]);
102                        last = Integer.parseInt(range[1]);
103                } catch (Exception e) {
104                        throw new IllegalArgumentException("Invalid range specification in oper_expression: " + expression);
105                }
106
107                List<String> expandedExpression = new ArrayList<String>(last-first+1);
108                for (int i = first; i <= last; i++) {
109                        expandedExpression.add(String.valueOf(i));
110                }
111                return expandedExpression;
112        }
113
114        public static List<OrderedPair<String>> parseBinaryOperatorExpression(String expression)
115                        throws IllegalArgumentException {
116                // split operator expression, i.e. (1,2,3)(4,5) into two subexpressions
117                String[] subExpressions = null;
118                try {
119                        subExpressions = expression.split("\\)\\(");
120                } catch (Exception e) {
121                        throw new IllegalArgumentException("Invalid oper_expression: " + expression);
122                }
123                if (subExpressions.length != 2) {
124                        throw new IllegalArgumentException("Invalid oper_expression: " + expression);
125                }
126                List<String> leftSide = parseSubExpression(subExpressions[0]);
127                List<String> rightSide = parseSubExpression(subExpressions[1]);
128
129                // form the cartesian product of the two lists
130                CartesianProduct<String> product = new CartesianProduct<String>(leftSide, rightSide);
131                return product.getOrderedPairs();
132        }
133
134        public static double[][]  getBiologicalMoleculeBounds(Structure asymStructure,List<BiologicalAssemblyTransformation> transformations) {
135                final double[][] coordinateBounds = new double[2][3];
136                coordinateBounds[0][0] = Double.MAX_VALUE;  // min x
137                coordinateBounds[0][1] = Double.MAX_VALUE;  // min y
138                coordinateBounds[0][2] = Double.MAX_VALUE;  // min z
139                coordinateBounds[1][0] = Double.MIN_VALUE;  // max x
140                coordinateBounds[1][1] = Double.MIN_VALUE;  // max y
141                coordinateBounds[1][2] = Double.MIN_VALUE;  // max z
142
143                double[] transformedCoords = new double[3];
144
145                Atom[] atoms = StructureTools.getAllAtomArray(asymStructure);
146
147                for ( Atom a : atoms) {
148
149                        Chain c = a.getGroup().getChain();
150                        String intChainID = c.getId();
151
152                        for (BiologicalAssemblyTransformation m: transformations) {
153                                if ( ! m.getChainId().equals(intChainID))
154                                        continue;
155                                Point3d coords = a.getCoordsAsPoint3d();
156                                transformedCoords[0] = coords.x;
157                                transformedCoords[1] = coords.y;
158                                transformedCoords[2] = coords.z;
159
160                                if (transformedCoords[0] < coordinateBounds[0][0] ) {
161                                        coordinateBounds[0][0] = transformedCoords[0];  // min x
162                                }
163
164                                if (transformedCoords[1] < coordinateBounds[0][1] ) {
165                                        coordinateBounds[0][1] = transformedCoords[1];  // min y
166                                }
167
168                                if (transformedCoords[2] < coordinateBounds[0][2] ) {
169                                        coordinateBounds[0][2] = transformedCoords[2];  // min z
170                                }
171
172                                if (transformedCoords[0] > coordinateBounds[1][0] ) {
173                                        coordinateBounds[1][0] = transformedCoords[0];  // max x
174                                }
175
176                                if (transformedCoords[1] > coordinateBounds[1][1] ) {
177                                        coordinateBounds[1][1] = transformedCoords[1];  // max y
178                                }
179
180                                if (transformedCoords[2] > coordinateBounds[1][2] ) {
181                                        coordinateBounds[1][2] = transformedCoords[2];  // max z
182                                }
183                        }
184                }
185                return coordinateBounds;
186        }
187        public static double[][] getAtomCoordinateBounds(Structure s){
188
189                Atom[] atoms = StructureTools.getAllAtomArray(s);
190                int atomCount = atoms.length;
191                final double[][] coordinateBounds = new double[2][3];
192                if ( atomCount <= 0 ) {
193                        return coordinateBounds;
194                }
195
196                Atom a = atoms[0];
197
198                coordinateBounds[0][0] = a.getX();  // min x
199                coordinateBounds[0][1] = a.getY();  // min y
200                coordinateBounds[0][2] = a.getZ();  // min z
201                coordinateBounds[1][0] = a.getX();  // max x
202                coordinateBounds[1][1] = a.getY();  // max y
203                coordinateBounds[1][2] = a.getZ();  // max z
204
205                for ( int i=1; i<atomCount; i++ )
206                {
207                        a =atoms[i];
208
209                        if ( a.getX() < coordinateBounds[0][0] ) {
210                                coordinateBounds[0][0] = a.getX();  // min x
211                        }
212
213                        if ( a.getY() < coordinateBounds[0][1] ) {
214                                coordinateBounds[0][1] = a.getY();  // min y
215                        }
216
217                        if ( a.getZ() < coordinateBounds[0][2] ) {
218                                coordinateBounds[0][2] = a.getZ();  // min z
219                        }
220
221                        if ( a.getX() > coordinateBounds[1][0] ) {
222                                coordinateBounds[1][0] = a.getX();  // max x
223                        }
224
225                        if ( a.getY() > coordinateBounds[1][1] ) {
226                                coordinateBounds[1][1] = a.getY();  // max y
227                        }
228
229                        if ( a.getZ() > coordinateBounds[1][2] ) {
230                                coordinateBounds[1][2] = a.getZ();  // max z
231                        }
232                }
233
234                return coordinateBounds;
235        }
236
237        /**
238         * Returns the maximum extend of the structure in the x, y, or z direction.
239         * @param structure
240         * @return maximum extend
241         */
242        public static double getMaximumExtend( final Structure structure ) {
243                double[][] bounds = getAtomCoordinateBounds(structure);
244                double xMax = Math.abs(bounds[0][0] - bounds[1][0]);
245                double yMax = Math.abs(bounds[0][1] - bounds[1][1]);
246                double zMax = Math.abs(bounds[0][2] - bounds[1][2]);
247                return Math.max(xMax, Math.max(yMax, zMax));
248        }
249
250        /**
251         * Returns the maximum extend of the biological molecule in the x, y, or z direction.
252         * @param structure
253         * @return maximum extend
254         */
255        public static double getBiologicalMoleculeMaximumExtend( final Structure structure,List<BiologicalAssemblyTransformation> transformations ) {
256                double[][] bounds = getBiologicalMoleculeBounds(structure, transformations);
257                double xMax = Math.abs(bounds[0][0] - bounds[1][0]);
258                double yMax = Math.abs(bounds[0][1] - bounds[1][1]);
259                double zMax = Math.abs(bounds[0][2] - bounds[1][2]);
260                return Math.max(xMax, Math.max(yMax, zMax));
261        }
262
263        /**
264         * Returns the centroid of the biological molecule.
265         * @param structure
266         * @return centroid
267         * @throws IllegalArgumentException if structure is null
268         */
269
270        public static double[] getBiologicalMoleculeCentroid( final Structure asymUnit,List<BiologicalAssemblyTransformation> transformations ) throws IllegalArgumentException {
271                if ( asymUnit == null ) {
272                        throw new IllegalArgumentException( "null structure" );
273                }
274
275                Atom[] atoms = StructureTools.getAllAtomArray(asymUnit);
276                int atomCount = atoms.length;
277                double[] centroid = new double[3];
278
279                if ( atomCount <= 0 ) {
280                        return centroid;
281                }
282
283                if ( transformations.size() == 0) {
284                        return Calc.getCentroid(atoms).getCoords();
285
286                }
287
288
289
290                int count = 0;
291                double[] transformedCoordinate = new double[3];
292
293                for (int i = 0; i < atomCount; i++)
294                {
295                        Atom atom = atoms[i];
296                        Chain chain = atom.getGroup().getChain();
297                        String intChainID = chain.getId();
298
299
300                        for (BiologicalAssemblyTransformation m: transformations) {
301                                if (!  m.getChainId().equals(intChainID))
302                                        continue;
303
304                                Point3d coords = atom.getCoordsAsPoint3d();
305                                transformedCoordinate[0] = coords.x;
306                                transformedCoordinate[1] = coords.y;
307                                transformedCoordinate[2] = coords.z;
308                                m.transformPoint(transformedCoordinate);
309                                centroid[0] += transformedCoordinate[0];
310                                centroid[1] += transformedCoordinate[1];
311                                centroid[2] += transformedCoordinate[2];
312                                count++;
313                        }
314                }
315
316
317
318                centroid[0] /= count;
319                centroid[1] /= count;
320                centroid[2] /= count;
321
322                return centroid;
323        }
324
325        /** 
326         * Reduce a structure to a single-atom representation (e.g. CA atoms)
327         *
328         * @param orig
329         * @return
330         * @since Biojava 4.1.0
331         */
332        public static Structure getReducedStructure(Structure orig){
333                Structure s = new StructureImpl();
334                s.setPDBHeader(orig.getPDBHeader());
335                for ( Chain c : orig.getChains()){
336
337                        Chain c1 = new ChainImpl();
338                        c1.setId(c.getId());
339                        c1.setName(c.getName());
340                        s.addChain(c1);
341
342                        for (Group g : c.getAtomGroups()){
343
344                                Atom a = null;
345                                switch(g.getType()) {
346                                case AMINOACID:
347                                        a = g.getAtom(StructureTools.CA_ATOM_NAME);
348                                        break;
349                                case NUCLEOTIDE:
350                                        a = g.getAtom(StructureTools.NUCLEOTIDE_REPRESENTATIVE);
351                                        break;
352                                default:
353                                        //omit group
354                                }
355                                if ( a != null){
356
357                                        Group g1 = (Group)g.clone();
358                                        g1.clearAtoms();
359                                        g1.addAtom(a);
360                                        c1.addGroup(g1);
361
362                                }
363
364                        }
365
366                }
367                return s;
368        }
369
370}