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