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.contact; 022 023import org.biojava.nbio.structure.Atom; 024import org.biojava.nbio.structure.Chain; 025import org.biojava.nbio.structure.Element; 026import org.biojava.nbio.structure.EntityInfo; 027import org.biojava.nbio.structure.Group; 028import org.biojava.nbio.structure.GroupType; 029import org.biojava.nbio.structure.ResidueNumber; 030import org.biojava.nbio.structure.Structure; 031import org.biojava.nbio.structure.asa.AsaCalculator; 032import org.biojava.nbio.structure.asa.GroupAsa; 033import org.biojava.nbio.structure.chem.ChemComp; 034import org.biojava.nbio.structure.chem.PolymerType; 035import org.biojava.nbio.structure.io.FileConvert; 036import org.biojava.nbio.structure.io.FileParsingParameters; 037import org.biojava.nbio.structure.io.cif.AbstractCifFileSupplier; 038import org.biojava.nbio.structure.xtal.CrystalTransform; 039import org.rcsb.cif.CifBuilder; 040import org.rcsb.cif.CifIO; 041import org.rcsb.cif.model.Category; 042import org.rcsb.cif.schema.StandardSchemata; 043import org.rcsb.cif.schema.mm.MmCifBlockBuilder; 044import org.slf4j.Logger; 045import org.slf4j.LoggerFactory; 046 047import java.io.IOException; 048import java.io.Serializable; 049import java.io.UncheckedIOException; 050import java.util.ArrayList; 051import java.util.List; 052import java.util.Map; 053import java.util.TreeMap; 054 055 056/** 057 * An interface between 2 molecules (2 sets of atoms). 058 * 059 * @author duarte_j 060 * 061 */ 062public class StructureInterface implements Serializable, Comparable<StructureInterface> { 063 064 private static final long serialVersionUID = 1L; 065 066 private static final Logger logger = LoggerFactory.getLogger(StructureInterface.class); 067 068 /** 069 * Interfaces with larger inverse self contact overlap score will be considered isologous 070 */ 071 private static final double SELF_SCORE_FOR_ISOLOGOUS = 0.3; 072 073 private int id; 074 private double totalArea; 075 private AtomContactSet contacts; 076 private GroupContactSet groupContacts; 077 078 private Pair<Atom[]> molecules; 079 080 /** 081 * The identifier for each of the atom arrays (usually a chain identifier, i.e. a single capital letter) 082 * Serves to identify the molecules within the Asymmetric Unit of the crystal 083 */ 084 private Pair<String> moleculeIds; 085 086 /** 087 * The transformations (crystal operators) applied to each molecule (if applicable) 088 */ 089 private Pair<CrystalTransform> transforms; 090 091 private Map<ResidueNumber, GroupAsa> groupAsas1; 092 private Map<ResidueNumber, GroupAsa> groupAsas2; 093 094 private StructureInterfaceCluster cluster; 095 096 /** 097 * Constructs a StructureInterface 098 * @param firstMolecule the atoms of the first molecule 099 * @param secondMolecule the atoms of the second molecule 100 * @param firstMoleculeId an identifier that identifies the first molecule within the Asymmetric Unit 101 * @param secondMoleculeId an identifier that identifies the second molecule within the Asymmetric Unit 102 * @param contacts the contacts between the 2 molecules 103 * @param firstTransf the transformation (crystal operator) applied to first molecule 104 * @param secondTransf the transformation (crystal operator) applied to second molecule 105 */ 106 public StructureInterface( 107 Atom[] firstMolecule, Atom[] secondMolecule, 108 String firstMoleculeId, String secondMoleculeId, 109 AtomContactSet contacts, 110 CrystalTransform firstTransf, CrystalTransform secondTransf) { 111 112 this.molecules = new Pair<>(firstMolecule, secondMolecule); 113 this.moleculeIds = new Pair<>(firstMoleculeId,secondMoleculeId); 114 this.contacts = contacts; 115 this.transforms = new Pair<>(firstTransf, secondTransf); 116 117 this.groupAsas1 = new TreeMap<>(); 118 this.groupAsas2 = new TreeMap<>(); 119 } 120 121 /** 122 * Constructs an empty StructureInterface 123 */ 124 public StructureInterface() { 125 this.groupAsas1 = new TreeMap<>(); 126 this.groupAsas2 = new TreeMap<>(); 127 } 128 129 public int getId() { 130 return id; 131 } 132 133 public void setId(int id) { 134 this.id = id; 135 } 136 137 /** 138 * Returns a pair of identifiers for each of the 2 member molecules that 139 * identify them uniquely in the crystal: 140 * <molecule id (asym unit id)>+<operator id>+<crystal translation> 141 * @return 142 */ 143 public Pair<String> getCrystalIds() { 144 return new Pair<>( 145 moleculeIds.getFirst()+transforms.getFirst().getTransformId()+transforms.getFirst().getCrystalTranslation(), 146 moleculeIds.getSecond()+transforms.getSecond().getTransformId()+transforms.getSecond().getCrystalTranslation()); 147 } 148 149 /** 150 * Returns the total area buried upon formation of this interface, 151 * defined as: 1/2[ (ASA1u-ASA1c) + (ASA2u-ASA2u) ] , with: 152 * <p>ASAxu = ASA of first/second unbound chain</p> 153 * <p>ASAxc = ASA of first/second complexed chain</p> 154 * In the area calculation HETATOM groups not part of the main protein/nucleotide chain 155 * are not included. 156 * @return 157 */ 158 public double getTotalArea() { 159 return totalArea; 160 } 161 162 public void setTotalArea(double totalArea) { 163 this.totalArea = totalArea; 164 } 165 166 public AtomContactSet getContacts() { 167 return contacts; 168 } 169 170 public void setContacts(AtomContactSet contacts) { 171 this.contacts = contacts; 172 } 173 174 public Pair<Atom[]> getMolecules() { 175 return molecules; 176 } 177 178 public void setMolecules(Pair<Atom[]> molecules) { 179 this.molecules = molecules; 180 } 181 182 /** 183 * Return the pair of identifiers identifying each of the 2 molecules of this interface 184 * in the asymmetry unit (usually the chain identifier if this interface is between 2 chains) 185 * @return 186 */ 187 public Pair<String> getMoleculeIds() { 188 return moleculeIds; 189 } 190 191 public void setMoleculeIds(Pair<String> moleculeIds) { 192 this.moleculeIds = moleculeIds; 193 } 194 195 /** 196 * Return the 2 crystal transform operations performed on each of the 197 * molecules of this interface. 198 * @return 199 */ 200 public Pair<CrystalTransform> getTransforms() { 201 return transforms; 202 } 203 204 public void setTransforms(Pair<CrystalTransform> transforms) { 205 this.transforms = transforms; 206 } 207 208 /** 209 * Set ASA annotations by passing the uncomplexed ASA values of the 2 partners. 210 * This will calculate complexed ASA and set the ASA values in the member variables. 211 * @param asas1 ASA values for atoms of partner 1 212 * @param asas2 ASA values for atoms of partner 2 213 * @param nSpherePoints the number of sphere points to be used for complexed ASA calculation 214 * @param nThreads the number of threads to be used for complexed ASA calculation 215 * @param cofactorSizeToUse the minimum size of cofactor molecule (non-chain HET atoms) that will be used in ASA calculation 216 */ 217 void setAsas(double[] asas1, double[] asas2, int nSpherePoints, int nThreads, int cofactorSizeToUse) { 218 219 Atom[] atoms = getAtomsForAsa(cofactorSizeToUse); 220 AsaCalculator asaCalc = new AsaCalculator(atoms, 221 AsaCalculator.DEFAULT_PROBE_SIZE, nSpherePoints, nThreads); 222 223 double[] complexAsas = asaCalc.calculateAsas(); 224 225 if (complexAsas.length!=asas1.length+asas2.length) 226 throw new IllegalArgumentException("The size of ASAs of complex doesn't match that of ASAs 1 + ASAs 2"); 227 228 229 groupAsas1 = new TreeMap<>(); 230 groupAsas2 = new TreeMap<>(); 231 232 this.totalArea = 0; 233 234 for (int i=0;i<asas1.length;i++) { 235 Group g = atoms[i].getGroup(); 236 237 if (!g.getType().equals(GroupType.HETATM) || 238 isInChain(g)) { 239 // interface area should be only for protein/nucleotide but not hetatoms that are not part of the chain 240 this.totalArea += (asas1[i] - complexAsas[i]); 241 } 242 243 if (!groupAsas1.containsKey(g.getResidueNumber())) { 244 GroupAsa groupAsa = new GroupAsa(g); 245 groupAsa.addAtomAsaU(asas1[i]); 246 groupAsa.addAtomAsaC(complexAsas[i]); 247 groupAsas1.put(g.getResidueNumber(), groupAsa); 248 } else { 249 GroupAsa groupAsa = groupAsas1.get(g.getResidueNumber()); 250 groupAsa.addAtomAsaU(asas1[i]); 251 groupAsa.addAtomAsaC(complexAsas[i]); 252 } 253 } 254 255 for (int i=0;i<asas2.length;i++) { 256 Group g = atoms[i+asas1.length].getGroup(); 257 258 if (!g.getType().equals(GroupType.HETATM) || 259 isInChain(g)) { 260 // interface area should be only for protein/nucleotide but not hetatoms that are not part of the chain 261 this.totalArea += (asas2[i] - complexAsas[i+asas1.length]); 262 } 263 264 if (!groupAsas2.containsKey(g.getResidueNumber())) { 265 GroupAsa groupAsa = new GroupAsa(g); 266 groupAsa.addAtomAsaU(asas2[i]); 267 groupAsa.addAtomAsaC(complexAsas[i+asas1.length]); 268 groupAsas2.put(g.getResidueNumber(), groupAsa); 269 } else { 270 GroupAsa groupAsa = groupAsas2.get(g.getResidueNumber()); 271 groupAsa.addAtomAsaU(asas2[i]); 272 groupAsa.addAtomAsaC(complexAsas[i+asas1.length]); 273 } 274 } 275 276 // our interface area definition: average of bsa of both molecules 277 this.totalArea = this.totalArea/2.0; 278 279 } 280 281 protected Atom[] getFirstAtomsForAsa(int cofactorSizeToUse) { 282 283 return getAllNonHAtomArray(molecules.getFirst(), cofactorSizeToUse); 284 } 285 286 protected Atom[] getSecondAtomsForAsa(int cofactorSizeToUse) { 287 288 return getAllNonHAtomArray(molecules.getSecond(), cofactorSizeToUse); 289 } 290 291 protected Atom[] getAtomsForAsa(int cofactorSizeToUse) { 292 Atom[] atoms1 = getFirstAtomsForAsa(cofactorSizeToUse); 293 Atom[] atoms2 = getSecondAtomsForAsa(cofactorSizeToUse); 294 295 Atom[] atoms = new Atom[atoms1.length+atoms2.length]; 296 System.arraycopy(atoms1, 0, atoms, 0, atoms1.length); 297 System.arraycopy(atoms2, 0, atoms, atoms1.length, atoms2.length); 298 299 return atoms; 300 } 301 302 /** 303 * Returns and array of all non-Hydrogen atoms in the given molecule, including all 304 * main chain HETATOM groups. Non main-chain HETATOM groups with fewer than minSizeHetAtomToInclude 305 * non-Hydrogen atoms are not included. 306 * @param m 307 * @param minSizeHetAtomToInclude HETATOM groups (non main-chain) with fewer number of 308 * non-Hydrogen atoms are not included 309 * @return 310 */ 311 private static Atom[] getAllNonHAtomArray(Atom[] m, int minSizeHetAtomToInclude) { 312 List<Atom> atoms = new ArrayList<>(); 313 314 for (Atom a:m){ 315 316 if (a.getElement()==Element.H) continue; 317 318 Group g = a.getGroup(); 319 if (g.getType().equals(GroupType.HETATM) && 320 !isInChain(g) && 321 getSizeNoH(g)<minSizeHetAtomToInclude) { 322 continue; 323 } 324 325 atoms.add(a); 326 327 } 328 return atoms.toArray(new Atom[0]); 329 } 330 331 /** 332 * Calculates the number of non-Hydrogen atoms in the given group 333 * @param g 334 * @return 335 */ 336 private static int getSizeNoH(Group g) { 337 int size = 0; 338 for (Atom a:g.getAtoms()) { 339 if (a.getElement()!=Element.H) 340 size++; 341 } 342 return size; 343 } 344 345 /** 346 * Returns true if the given group is part of the main chain, i.e. if it is 347 * a peptide-linked group or a nucleotide 348 * @param g 349 * @return 350 */ 351 private static boolean isInChain(Group g) { 352 ChemComp chemComp = g.getChemComp(); 353 354 if (chemComp==null) { 355 logger.warn("Can't determine PolymerType for group "+g.getResidueNumber()+" ("+g.getPDBName()+"). Will consider it as non-nucleotide/non-protein type."); 356 return false; 357 } 358 359 PolymerType polyType = chemComp.getPolymerType(); 360 for (PolymerType protOnlyType: PolymerType.PROTEIN_ONLY) { 361 if (polyType==protOnlyType) return true; 362 } 363 for (PolymerType protOnlyType: PolymerType.POLYNUCLEOTIDE_ONLY) { 364 if (polyType==protOnlyType) return true; 365 } 366 367 return false; 368 } 369 370 /** 371 * Tells whether the interface corresponds to one mediated by crystallographic symmetry, 372 * i.e. it is between symmetry-related molecules (with same chain identifier) 373 * @return 374 */ 375 public boolean isSymRelated() { 376 return moleculeIds.getFirst().equals(moleculeIds.getSecond()); 377 } 378 379 /** 380 * Returns true if the transformation applied to the second molecule of this interface 381 * has an infinite character (pure translation or screw rotation) 382 * and both molecules of the interface have the same asymmetric unit identifier (chain id): in such cases the 383 * interface would lead to infinite fiber-like (linear or helical) assemblies 384 * @return 385 */ 386 public boolean isInfinite() { 387 return ((isSymRelated() && transforms.getSecond().getTransformType().isInfinite())); 388 } 389 390 /** 391 * Returns true if the 2 molecules of this interface are the same entity (i.e. homomeric interface), false 392 * otherwise (i.e. heteromeric interface) 393 * @return true if homomeric or if either of the entities is unknonw (null Compounds), false otherwise 394 */ 395 public boolean isHomomeric() { 396 EntityInfo first = getParentChains().getFirst().getEntityInfo(); 397 EntityInfo second = getParentChains().getSecond().getEntityInfo(); 398 if (first==null || second==null) { 399 logger.warn("Some compound of interface {} is null, can't determine whether it is homo/heteromeric. Consider it homomeric", getId()); 400 return true; 401 } 402 return 403 first.getRepresentative().getId().equals(second.getRepresentative().getId()); 404 } 405 406 /** 407 * Gets a map of ResidueNumbers to GroupAsas for all groups of first chain. 408 * @return 409 */ 410 public Map<ResidueNumber, GroupAsa> getFirstGroupAsas() { 411 return groupAsas1; 412 } 413 414 /** 415 * Gets the GroupAsa for the corresponding residue number of first chain 416 * @param resNum 417 * @return 418 */ 419 public GroupAsa getFirstGroupAsa(ResidueNumber resNum) { 420 return groupAsas1.get(resNum); 421 } 422 423 public void setFirstGroupAsa(GroupAsa groupAsa) { 424 groupAsas1.put(groupAsa.getGroup().getResidueNumber(), groupAsa); 425 } 426 427 public void setFirstGroupAsas(Map<ResidueNumber, GroupAsa> firstGroupAsas) { 428 this.groupAsas1 = firstGroupAsas; 429 } 430 431 public void setSecondGroupAsas(Map<ResidueNumber, GroupAsa> secondGroupAsas) { 432 this.groupAsas2 = secondGroupAsas; 433 } 434 435 /** 436 * Gets a map of ResidueNumbers to GroupAsas for all groups of second chain. 437 * @return 438 */ 439 public Map<ResidueNumber, GroupAsa> getSecondGroupAsas() { 440 return groupAsas2; 441 } 442 443 public void setSecondGroupAsa(GroupAsa groupAsa) { 444 groupAsas2.put(groupAsa.getGroup().getResidueNumber(), groupAsa); 445 } 446 447 /** 448 * Gets the GroupAsa for the corresponding residue number of second chain 449 * @param resNum 450 * @return 451 */ 452 public GroupAsa getSecondGroupAsa(ResidueNumber resNum) { 453 return groupAsas2.get(resNum); 454 } 455 456 /** 457 * Returns the residues belonging to the interface core, defined as those residues at 458 * the interface (BSA>0) and for which the BSA/ASA ratio is above the given bsaToAsaCutoff 459 * @param bsaToAsaCutoff 460 * @param minAsaForSurface the minimum ASA to consider a residue on the surface 461 * @return 462 */ 463 public Pair<List<Group>> getCoreResidues(double bsaToAsaCutoff, double minAsaForSurface) { 464 465 List<Group> core1 = new ArrayList<>(); 466 List<Group> core2 = new ArrayList<>(); 467 468 for (GroupAsa groupAsa:groupAsas1.values()) { 469 470 if (groupAsa.getAsaU()>minAsaForSurface && groupAsa.getBsa()>0) { 471 if (groupAsa.getBsaToAsaRatio()<bsaToAsaCutoff) { 472 //rim1.add(groupAsa.getGroup()); 473 } else { 474 core1.add(groupAsa.getGroup()); 475 } 476 } 477 } 478 for (GroupAsa groupAsa:groupAsas2.values()) { 479 480 if (groupAsa.getAsaU()>minAsaForSurface && groupAsa.getBsa()>0) { 481 if (groupAsa.getBsaToAsaRatio()<bsaToAsaCutoff) { 482 //rim2.add(groupAsa.getGroup()); 483 } else { 484 core2.add(groupAsa.getGroup()); 485 } 486 } 487 } 488 489 return new Pair<>(core1, core2); 490 } 491 492 /** 493 * Returns the residues belonging to the interface rim, defined as those residues at 494 * the interface (BSA>0) and for which the BSA/ASA ratio is below the given bsaToAsaCutoff 495 * @param bsaToAsaCutoff 496 * @param minAsaForSurface the minimum ASA to consider a residue on the surface 497 * @return 498 */ 499 public Pair<List<Group>> getRimResidues(double bsaToAsaCutoff, double minAsaForSurface) { 500 501 List<Group> rim1 = new ArrayList<>(); 502 List<Group> rim2 = new ArrayList<>(); 503 504 for (GroupAsa groupAsa:groupAsas1.values()) { 505 506 if (groupAsa.getAsaU()>minAsaForSurface && groupAsa.getBsa()>0) { 507 if (groupAsa.getBsaToAsaRatio()<bsaToAsaCutoff) { 508 rim1.add(groupAsa.getGroup()); 509 } else { 510 //core1.add(groupAsa.getGroup()); 511 } 512 } 513 } 514 for (GroupAsa groupAsa:groupAsas2.values()) { 515 516 if (groupAsa.getAsaU()>minAsaForSurface && groupAsa.getBsa()>0) { 517 if (groupAsa.getBsaToAsaRatio()<bsaToAsaCutoff) { 518 rim2.add(groupAsa.getGroup()); 519 } else { 520 //core2.add(groupAsa.getGroup()); 521 } 522 } 523 } 524 525 return new Pair<List<Group>>(rim1, rim2); 526 } 527 528 /** 529 * Returns the residues belonging to the interface, i.e. the residues 530 * at the surface with BSA>0 531 * @param minAsaForSurface the minimum ASA to consider a residue on the surface 532 * @return 533 */ 534 public Pair<List<Group>> getInterfacingResidues(double minAsaForSurface) { 535 536 List<Group> interf1 = new ArrayList<>(); 537 List<Group> interf2 = new ArrayList<>(); 538 539 for (GroupAsa groupAsa:groupAsas1.values()) { 540 541 if (groupAsa.getAsaU()>minAsaForSurface && groupAsa.getBsa()>0) { 542 interf1.add(groupAsa.getGroup()); 543 } 544 } 545 for (GroupAsa groupAsa:groupAsas2.values()) { 546 547 if (groupAsa.getAsaU()>minAsaForSurface && groupAsa.getBsa()>0) { 548 interf2.add(groupAsa.getGroup()); 549 } 550 } 551 552 return new Pair<>(interf1, interf2); 553 } 554 555 /** 556 * Returns the residues belonging to the surface 557 * @param minAsaForSurface the minimum ASA to consider a residue on the surface 558 * @return 559 */ 560 public Pair<List<Group>> getSurfaceResidues(double minAsaForSurface) { 561 List<Group> surf1 = new ArrayList<>(); 562 List<Group> surf2 = new ArrayList<>(); 563 564 for (GroupAsa groupAsa:groupAsas1.values()) { 565 566 if (groupAsa.getAsaU()>minAsaForSurface) { 567 surf1.add(groupAsa.getGroup()); 568 } 569 } 570 for (GroupAsa groupAsa:groupAsas2.values()) { 571 572 if (groupAsa.getAsaU()>minAsaForSurface) { 573 surf2.add(groupAsa.getGroup()); 574 } 575 } 576 577 return new Pair<>(surf1, surf2); 578 } 579 580 public StructureInterfaceCluster getCluster() { 581 return cluster; 582 } 583 584 public void setCluster(StructureInterfaceCluster cluster) { 585 this.cluster = cluster; 586 } 587 588 /** 589 * Calculates the Jaccard contact set score (intersection over union) between this StructureInterface and 590 * the given one. The calculation assumes that both interfaces come from the same structure. The output 591 * will not necessarily make sense if the two interfaces come from different structures. 592 * The two sides of the given StructureInterface need to match this StructureInterface 593 * in the sense that they must come from the same Entity, i.e. 594 * their residue numbers need to align with 100% identity, except for unobserved 595 * density residues. The SEQRES indices obtained through {@link EntityInfo#getAlignedResIndex(Group, Chain)} are 596 * used to match residues, thus if no SEQRES is present or if {@link FileParsingParameters#setAlignSeqRes(boolean)} 597 * is not used, this calculation is not guaranteed to work properly. 598 * @param other the interface to be compared to this one 599 * @param invert if false the comparison will be done first-to-first and second-to-second, 600 * if true the match will be first-to-second and second-to-first 601 * @return the contact overlap score, range [0.0,1.0] 602 */ 603 public double getContactOverlapScore(StructureInterface other, boolean invert) { 604 605 Pair<Chain> thisChains = getParentChains(); 606 Pair<Chain> otherChains = other.getParentChains(); 607 608 if (thisChains.getFirst().getEntityInfo() == null || thisChains.getSecond().getEntityInfo() == null || 609 otherChains.getFirst().getEntityInfo() == null || otherChains.getSecond().getEntityInfo() == null ) { 610 // this happens in cases like 2uub 611 logger.warn("Found chains with null compounds while comparing interfaces {} and {}. Contact overlap score for them will be 0.", 612 this.getId(), other.getId()); 613 return 0; 614 } 615 616 Pair<EntityInfo> thisCompounds = new Pair<EntityInfo>(thisChains.getFirst().getEntityInfo(), thisChains.getSecond().getEntityInfo()); 617 Pair<EntityInfo> otherCompounds = new Pair<EntityInfo>(otherChains.getFirst().getEntityInfo(), otherChains.getSecond().getEntityInfo()); 618 619 if ( ( (thisCompounds.getFirst().getMolId() == otherCompounds.getFirst().getMolId()) && 620 (thisCompounds.getSecond().getMolId() == otherCompounds.getSecond().getMolId()) ) || 621 ( (thisCompounds.getFirst().getMolId() == otherCompounds.getSecond().getMolId()) && 622 (thisCompounds.getSecond().getMolId() == otherCompounds.getFirst().getMolId()) ) ) { 623 624 int common = 0; 625 GroupContactSet thisContacts = getGroupContacts(); 626 GroupContactSet otherContacts = other.getGroupContacts(); 627 628 for (GroupContact thisContact:thisContacts) { 629 630 ResidueIdentifier first; 631 ResidueIdentifier second; 632 633 if (!invert) { 634 first = new ResidueIdentifier(thisContact.getPair().getFirst()); 635 636 second = new ResidueIdentifier(thisContact.getPair().getSecond()); 637 } else { 638 first = new ResidueIdentifier(thisContact.getPair().getSecond()); 639 640 second = new ResidueIdentifier(thisContact.getPair().getFirst()); 641 } 642 643 if (otherContacts.hasContact(first,second)) { 644 common++; 645 } 646 } 647 return (2.0*common)/(thisContacts.size()+otherContacts.size()); 648 } else { 649 logger.debug("Chain pairs {},{} and {},{} belong to different compound pairs, contact overlap score will be 0 ", 650 thisChains.getFirst().getId(),thisChains.getSecond().getId(), 651 otherChains.getFirst().getId(),otherChains.getSecond().getId()); 652 return 0.0; 653 } 654 } 655 656 public GroupContactSet getGroupContacts() { 657 if (groupContacts==null) { 658 this.groupContacts = new GroupContactSet(contacts); 659 } 660 return this.groupContacts; 661 } 662 663 /** 664 * Tell whether the interface is isologous, i.e. it is formed 665 * by the same patches of same entity on both sides. 666 * 667 * @return true if isologous, false if heterologous 668 */ 669 public boolean isIsologous() { 670 double scoreInverse = this.getContactOverlapScore(this, true); 671 logger.debug("Interface {} contact overlap score with itself inverted: {}", 672 getId(), scoreInverse); 673 return (scoreInverse>SELF_SCORE_FOR_ISOLOGOUS); 674 } 675 676 /** 677 * Finds the parent chains by looking up the references of first atom of each side of this interface 678 * @return 679 */ 680 public Pair<Chain> getParentChains() { 681 Atom[] firstMol = this.molecules.getFirst(); 682 Atom[] secondMol = this.molecules.getSecond(); 683 if (firstMol.length==0 || secondMol.length==0) { 684 logger.warn("No atoms found in first or second molecule, can't get parent Chains"); 685 return null; 686 } 687 688 return new Pair<>(firstMol[0].getGroup().getChain(), secondMol[0].getGroup().getChain()); 689 } 690 691 /** 692 * Finds the parent entities by looking up the references of first atom of each side of this interface 693 * @return 694 */ 695 public Pair<EntityInfo> getParentCompounds() { 696 Pair<Chain> chains = getParentChains(); 697 if (chains == null) { 698 logger.warn("Could not find parents chains, compounds will be null"); 699 return null; 700 } 701 return new Pair<EntityInfo>(chains.getFirst().getEntityInfo(), chains.getSecond().getEntityInfo()); 702 } 703 704 private Structure getParentStructure() { 705 Atom[] firstMol = this.molecules.getFirst(); 706 if (firstMol.length==0) { 707 logger.warn("No atoms found in first molecule, can't get parent Structure"); 708 return null; 709 } 710 return firstMol[0].getGroup().getChain().getStructure(); 711 } 712 713 /** 714 * Return a String representing the 2 molecules of this interface in PDB format. 715 * If the molecule ids (i.e. chain ids) are the same for both molecules, then the second 716 * one will be replaced by the next letter in alphabet (or A for Z) 717 * @return the PDB-formatted string 718 */ 719 public String toPDB() { 720 721 String molecId1 = getMoleculeIds().getFirst(); 722 String molecId2 = getMoleculeIds().getSecond(); 723 724 if (molecId2.equals(molecId1)) { 725 // if both chains are named equally we want to still named them differently in the output pdb file 726 // so that molecular viewers can handle properly the 2 chains as separate entities 727 char letter = molecId1.charAt(0); 728 if (letter!='Z' && letter!='z') { 729 molecId2 = Character.toString((char)(letter+1)); // i.e. next letter in alphabet 730 } else { 731 molecId2 = Character.toString((char)(letter-25)); //i.e. 'A' or 'a' 732 } 733 } 734 735 StringBuilder sb = new StringBuilder(); 736 for (Atom atom:this.molecules.getFirst()) { 737 sb.append(FileConvert.toPDB(atom, molecId1)); 738 } 739 sb.append("TER"); 740 sb.append(System.getProperty("line.separator")); 741 for (Atom atom:this.molecules.getSecond()) { 742 sb.append(FileConvert.toPDB(atom,molecId2)); 743 } 744 sb.append("TER"); 745 sb.append(System.getProperty("line.separator")); 746 sb.append("END"); 747 sb.append(System.getProperty("line.separator")); 748 return sb.toString(); 749 } 750 751 /** 752 * Return a String representing the 2 molecules of this interface in mmCIF format. 753 * If the molecule ids (i.e. chain ids) are the same for both molecules, then the second 754 * one will be written as chainId_operatorId (with operatorId taken from {@link #getTransforms()} 755 * @return the mmCIF-formatted string 756 */ 757 public String toMMCIF() { 758 String molecId1 = getMoleculeIds().getFirst(); 759 String molecId2 = getMoleculeIds().getSecond(); 760 761 if (isSymRelated()) { 762 // if both chains are named equally we want to still named them differently in the output mmcif file 763 // so that molecular viewers can handle properly the 2 chains as separate entities 764 molecId2 = molecId2 + "_" + getTransforms().getSecond().getTransformId(); 765 } 766 767 MmCifBlockBuilder mmCifBlockBuilder = CifBuilder.enterFile(StandardSchemata.MMCIF) 768 .enterBlock("BioJava_interface_" + getId()); 769 770 // we reassign atom ids if sym related (otherwise atom ids would be duplicated and some molecular viewers can't cope with that) 771 int atomId = 1; 772 List<AbstractCifFileSupplier.WrappedAtom> wrappedAtoms = new ArrayList<>(); 773 for (Atom atom : this.molecules.getFirst()) { 774 if (isSymRelated()) { 775 wrappedAtoms.add(new AbstractCifFileSupplier.WrappedAtom(1, molecId1, molecId1, atom, atomId)); 776 } else { 777 wrappedAtoms.add(new AbstractCifFileSupplier.WrappedAtom(1, molecId1, molecId1, atom, atom.getPDBserial())); 778 } 779 atomId++; 780 } 781 for (Atom atom : this.molecules.getSecond()) { 782 if (isSymRelated()) { 783 wrappedAtoms.add(new AbstractCifFileSupplier.WrappedAtom(1, molecId2, molecId2, atom, atomId)); 784 } else { 785 wrappedAtoms.add(new AbstractCifFileSupplier.WrappedAtom(1, molecId2, molecId2, atom, atom.getPDBserial())); 786 } 787 atomId++; 788 } 789 790 Category atomSite = wrappedAtoms.stream().collect(AbstractCifFileSupplier.toAtomSite()); 791 mmCifBlockBuilder.addCategory(atomSite); 792 793 try { 794 return new String(CifIO.writeText(mmCifBlockBuilder.leaveBlock().leaveFile())); 795 } catch (IOException e) { 796 throw new UncheckedIOException(e); 797 } 798 } 799 800 @Override 801 public int compareTo(StructureInterface o) { 802 // this will sort descending on interface areas 803 return (Double.compare(o.totalArea,this.totalArea)); 804 } 805 806 @Override 807 public String toString() { 808 return String.format("StructureInterface %d (%s, %.0f A, <%s; %s>)", id, moleculeIds,totalArea,transforms.getFirst().toXYZString(),transforms.getSecond().toXYZString()); 809 } 810 811}