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}