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.secstruc;
022
023import org.biojava.nbio.structure.*;
024import org.biojava.nbio.structure.contact.AtomContact;
025import org.biojava.nbio.structure.contact.AtomContactSet;
026import org.biojava.nbio.structure.contact.Grid;
027import org.biojava.nbio.structure.contact.Pair;
028import org.slf4j.Logger;
029import org.slf4j.LoggerFactory;
030
031import java.util.ArrayList;
032import java.util.Arrays;
033import java.util.Collections;
034import java.util.Comparator;
035import java.util.HashMap;
036import java.util.Iterator;
037import java.util.List;
038import java.util.Map;
039
040/**
041 * Calculate and assign the secondary structure (SS) to the
042 * Groups of a Structure object. This object also stores the result
043 * of the calculation.
044 * <p>
045 * The rules for SS calculation are the ones defined by DSSP:
046 * Kabsch,W. and Sander,C. (1983) Biopolymers 22, 2577-2637.
047 * Original DSSP article see at:
048 * <a href="http://www.cmbi.kun.nl/gv/dssp/dssp.pdf">dssp.pdf</a>.
049 * Some parts are also taken from: T.E.Creighton, Proteins -
050 * Structure and Molecular Properties, 2nd Edition, Freeman 1994.
051 *
052 * @author Andreas Prlic
053 * @author Aleix Lafita
054 * @autho Anthony Bradley
055 *
056 */
057public class SecStrucCalc {
058
059        /**
060         * DSSP assigns helices one residue shorter at each end, because the
061         * residues at (i-1) and (i+n+1) are not assigned helix type although
062         * they contain a consistent turn (H-bond). If this parameter
063         * is true, the helices will be the length of the original DSSP
064         * convention. If it is false, they will be two residue longer.
065         */
066        private static final boolean DSSP_HELICES = true;
067
068        private static final Logger logger =
069                        LoggerFactory.getLogger(SecStrucCalc.class);
070
071        /** min distance between two residues */
072        public static final double MINDIST = 0.5;
073
074        /** min distance of two CA atoms if H-bonds are allowed to form */
075        public static final double CA_MIN_DIST = 9.0;
076
077        /** max distance CA atoms in peptide bond (backbone discontinuity) */
078        public static final double MAX_PEPTIDE_BOND_LENGTH = 2.5;
079
080        /** Minimal H-bond energy in cal/mol */
081        public static final int HBONDLOWENERGY  = -9900;
082
083        /** higher limit for H-bond energy */
084        public static final double HBONDHIGHENERGY = -500.0;
085
086        /** constant for electrostatic energy
087         * <pre>
088         *      f  *  q1 *   q2  *  scale
089         * Q = -332 * 0.42 * 0.20 * 1000.0
090         *</pre>
091         *
092         * q1 and q2 are partial charges which are placed on the C,O
093         * (+q1,-q1) and N,H (-q2,+q2)
094         */
095        public static final double Q = -27888.0;
096
097        // Three lists
098        private SecStrucGroup[] groups;
099        private List<Ladder> ladders;
100        private List<BetaBridge> bridges;
101        private Atom[] atoms;
102        // Added by Anthony - to speed up intergroup calculations
103        private AtomContactSet contactSet;
104        private Map<ResidueNumber, Integer> indResMap;
105        public SecStrucCalc(){
106                ladders = new ArrayList<Ladder>();
107                bridges = new ArrayList<BetaBridge>();
108        }
109
110
111        /**
112         * Predicts the secondary structure of this Structure object,
113         * using a DSSP implementation.
114         *
115         * @param s Structure to predict the SS
116         * @param assign sets the SS information to the Groups of s
117         * @return a List of SS annotation objects
118         */
119        public List<SecStrucState> calculate(Structure s, boolean assign)
120                        throws StructureException {
121
122                List<SecStrucState> secstruc = new ArrayList<SecStrucState>();
123                for(int i=0; i<s.nrModels(); i++) {
124                        // Reinitialise the global vars
125                        ladders = new ArrayList<Ladder>();
126                        bridges = new ArrayList<BetaBridge>();
127                        groups = initGroupArray(s, i);
128                        // Initialise the contact set for this structure
129                        initContactSet();
130                        if (groups.length < 5) {
131                                // not enough groups to do anything
132                                throw new StructureException("Not enough backbone groups in the"
133                                                + " Structure to calculate the secondary structure ("
134                                                + groups.length+" given, minimum 5)" );
135                        }
136
137                        calculateHAtoms();
138                        calculateHBonds();
139                        calculateDihedralAngles();
140                        calculateTurns();
141                        buildHelices();
142                        detectBends();
143                        detectStrands();
144
145                        for (SecStrucGroup sg : groups){
146                                SecStrucState ss = (SecStrucState)
147                                                sg.getProperty(Group.SEC_STRUC);
148                                // Add to return list and assign to original if flag is true
149                                secstruc.add(ss);
150                                if (assign) sg.getOriginal().setProperty(Group.SEC_STRUC, ss);
151                        }
152                }
153                return secstruc;
154        }
155
156        /**
157         * Function to generate the contact sets
158         */
159        private void initContactSet() {
160
161                // Initialise an array of atoms
162                atoms = new Atom[groups.length];
163                // Remake this local var
164                indResMap = new HashMap<>();
165                for (int i=0 ; i < groups.length ; i++){
166                        SecStrucGroup one = groups[i];
167                        indResMap.put(one.getResidueNumber(), i);
168                        atoms[i] = one.getCA();
169                }
170                Grid grid = new Grid(CA_MIN_DIST);
171                if(atoms.length==0){
172                        contactSet = new AtomContactSet(CA_MIN_DIST);
173                }
174                else{
175                        grid.addAtoms(atoms);
176                        contactSet = grid.getAtomContacts();
177                }
178        }
179
180        /**
181         * Updated code to detect strands
182         */
183        private void detectStrands() {
184
185                //Find all the beta bridges of the structure
186                findBridges();
187                //Create Ladders
188                createLadders();
189
190                //Detect beta bulges between ladders
191                connectLadders();
192
193                //AND store SS assignments for Sheets, Strands and Bridges
194                updateSheets();
195        }
196
197
198        private void createLadders(){
199
200                for (BetaBridge b : bridges){
201                        boolean found = false;
202                        for (Ladder ladder : ladders){
203                                if (shouldExtendLadder(ladder, b)) {
204                                        found = true;
205                                        ladder.to++; //we go forward in this direction
206                                        switch(b.type){
207                                        case parallel:
208                                                ladder.lto++; //increment second strand
209                                                break;
210                                        case antiparallel:
211                                                ladder.lfrom--; //decrement second strand
212                                                break;
213                                        }
214                                        break;
215                                }
216                        }
217                        if (!found){
218                                //Create new ladder with a single Bridge
219                                Ladder l = new Ladder();
220                                l.from = b.partner1;
221                                l.to = b.partner1;
222                                l.lfrom = b.partner2;
223                                l.lto = b.partner2;
224                                l.btype = b.type;
225                                ladders.add(l);
226                        }
227                }
228        }
229
230
231        private void updateSheets() {
232
233                logger.debug(" got {} ladders!", ladders.size());
234
235                for (Ladder ladder : ladders){
236                        logger.debug(ladder.toString());
237
238                        for (int lcount = ladder.from; lcount <= ladder.to; lcount++) {
239
240                                SecStrucState state = getSecStrucState(lcount);
241                                SecStrucType stype = state.getType();
242
243                                int diff = ladder.from - lcount;
244                                int l2count = ladder.lfrom - diff ;
245
246                                SecStrucState state2 = getSecStrucState(l2count);
247                                SecStrucType stype2 = state2.getType();
248
249                                if ( ladder.from != ladder.to ) {
250                                        setSecStrucType(lcount, SecStrucType.extended);
251                                        setSecStrucType(l2count, SecStrucType.extended);
252                                }
253                                else {
254                                        if ( !stype.isHelixType() &&
255                                                        ( !stype.equals(SecStrucType.extended)))
256                                                setSecStrucType(lcount,SecStrucType.bridge);
257
258                                        if ( ! stype2.isHelixType() &&
259                                                        (! stype2.equals(SecStrucType.extended)))
260                                                setSecStrucType(l2count,SecStrucType.bridge);
261                                }
262                        }
263
264                        // Check if two ladders are connected. both sides are 'E'
265
266                        if (ladder.connectedTo == 0) continue;
267                        Ladder conladder = ladders.get(ladder.connectedTo);
268
269                        if (ladder.btype.equals(BridgeType.antiparallel)) {
270                                /* set one side */
271                                for (int lcount = ladder.from; lcount <= conladder.to;
272                                                lcount++) {
273                                        setSecStrucType(lcount, SecStrucType.extended);
274
275                                }
276                                /* set other side */
277                                for (int lcount = conladder.lto;
278                                                lcount <= ladder.lfrom;
279                                                lcount++) {
280                                        setSecStrucType(lcount, SecStrucType.extended);
281                                }
282
283                        } else {
284                                /* set one side */
285                                for ( int lcount = ladder.from;
286                                                lcount <= conladder.to;
287                                                lcount++) {
288
289                                        setSecStrucType(lcount, SecStrucType.extended);
290                                }
291                                /* set other side */
292                                for ( int lcount =  ladder.lfrom;
293                                                lcount <= conladder.lto;
294                                                lcount++) {
295
296                                        setSecStrucType(lcount, SecStrucType.extended);
297                                }
298                        }
299                }
300        }
301
302        private void connectLadders() {
303
304                for (int i = 0 ; i < ladders.size(); i++) {
305                        for ( int j = i ; j < ladders.size(); j++){
306                                Ladder l1 = ladders.get(i);
307                                Ladder l2 = ladders.get(j);
308                                if (hasBulge(l1,l2)) {
309                                        l1.connectedTo = j;
310                                        l2.connectedFrom = i;
311                                        logger.debug("Bulge from {} to {}", i, j);
312                                }
313                        }
314                }
315
316
317        }
318
319        /**
320         * For beta structures, we define explicitly: a bulge-linked
321         * ladder consists of two (perfect) ladder or bridges of the
322         * same type connected by at most one extra residue on one
323         * strand and at most four extra residues on the other strand,
324         * all residues in bulge-linked ladders are marked "E,"
325         * including the extra residues.
326         */
327        private boolean hasBulge(Ladder l1, Ladder l2) {
328
329                boolean bulge = ((l1.btype.equals(l2.btype)) &&
330                                (l2.from - l1.to < 6) &&
331                                (l1.to < l2.from) &&
332                                (l2.connectedTo == 0));
333
334                if (!bulge) return bulge;
335
336                switch(l1.btype){
337                case parallel:
338                        bulge = ( (l2.lfrom - l1.lto > 0) &&
339                                        ((( l2.lfrom -l1.lto < 6) &&
340                                                        (l2.from - l1.to < 3)) ||
341                                                        ( l2.lfrom - l1.lto <3)));
342
343                        break;
344
345                case antiparallel:
346                        bulge = ( (l1.lfrom - l2.lto > 0) &&
347                                        (((l1.lfrom -l2.lto < 6) &&
348                                                        ( l2.from - l1.to < 3)) ||
349                                                        (l1.lfrom - l2.lto < 3)));
350
351                        break;
352                }
353
354                return bulge;
355        }
356
357        private void registerBridge(int i, int j, BridgeType btype) {
358
359                BetaBridge bridge = new BetaBridge(i,j,btype);
360
361                boolean b1 = getSecStrucState(i).addBridge(bridge);
362                boolean b2 = getSecStrucState(j).addBridge(bridge);
363
364                if (!b1 && !b2)
365                        logger.warn("Ignoring Bridge between residues {} and {}. DSSP assignment might differ.", i, j);
366
367                bridges.add(bridge);
368        }
369
370        /**
371         * Conditions to extend a ladder with a given beta Bridge:
372         * <li>The bridge and ladder are of the same type.
373         * <li>The smallest bridge residue is sequential to the first
374         *              strand ladder.
375         * <li>The second bridge residue is either sequential (parallel)
376         *              or previous (antiparallel) to the second strand of the ladder
377         * </li>
378         * @param ladder the ladder candidate to extend
379         * @param b the beta bridge that would extend the ladder
380         * @return true if the bridge b extends the ladder
381         */
382        private boolean shouldExtendLadder(Ladder ladder, BetaBridge b) {
383
384                //Only extend if they are of the same type
385                boolean sameType = b.type.equals(ladder.btype);
386                if (!sameType) return false;
387
388                //Only extend if residue 1 is sequential to ladder strand
389                boolean sequential = (b.partner1 == ladder.to+1);
390                if (!sequential) return false;
391
392                switch(b.type){
393                case parallel:
394                        //Residue 2 should be sequential to second strand
395                        if (b.partner2 == ladder.lto+1) return true;
396                        break;
397                case antiparallel:
398                        //Residue 2 should be previous to second strand
399                        if (b.partner2 == ladder.lfrom-1) return true;
400                        break;
401                }
402                return false;
403        }
404
405
406
407
408        /**
409         * Two nonoverlapping stretches of three residues each, i-1,i,i+1 and
410         * j-1,j,j+1, form either a parallel or antiparallel bridge, depending on
411         * which of two basic patterns is matched. We assign a bridge between
412         * residues i and j if there are two H bonds characteristic of beta-
413         * structure; in particular:
414         * <p>
415         * Parallel Bridge(i,j) =: [Hbond(i-1,j) and Hbond(j,i+1)]
416         *                                                      or [Hbond(j-1,i) and Hbond(i,j+1)]
417         * <p>
418         * Antiparallel Bridge(i,j) =: [Hbond(i,j) and Hbond(j,i)]
419         *                                                              or [Hbond(i-1,j+1) and Hbond(j-1,i+1)]
420         *
421         * Optimised to use the contact set
422         */
423        private void findBridges() {
424                // Get the interator of contacts
425                Iterator<AtomContact> myIter = contactSet.iterator();
426                List<Pair<Integer>> outList = new ArrayList<Pair<Integer>>();
427
428                // Now iterate through this
429                while(myIter.hasNext()){
430                        // Get the next atom contact
431                        AtomContact ac = myIter.next();
432                        Group g1 = ac.getPair().getFirst().getGroup();
433                        Group g2 = ac.getPair().getSecond().getGroup();
434                        // Get the indices
435                        int i = indResMap.get(g1.getResidueNumber());
436                        int j = indResMap.get(g2.getResidueNumber());
437                        // If i>j switch them over
438                        if(i>j){
439                                // Switch them over
440                                int old = i;
441                                i = j;
442                                j = old;
443                        }
444                        // Only these
445                        if(j<i+3){
446                                continue;
447                        }
448                        // If it's the first
449                        if(i==0){
450                                continue;
451                        }
452                        if(j==0){
453                                continue;
454                        }
455                        // If it's the last
456                        if(i==groups.length-1){
457                                continue;
458                        }
459                        if(j==groups.length-1){
460                                continue;
461                        }
462
463                        Pair<Integer> thisPair = new Pair<Integer>(i,j);
464                        outList.add(thisPair);
465                }
466                //
467                Collections.sort(outList, new Comparator<Pair<Integer>>() {
468                        @Override
469                        public int compare(Pair<Integer> o1, Pair<Integer> o2) {
470                                if(o1.getFirst()<o2.getFirst()){
471                                        return -1;
472                                }
473                                else if(o1.getFirst()>o2.getFirst()){
474                                        return +1;
475                                }
476                                else{
477                                        if(o1.getSecond()<o2.getSecond()){
478                                                return -1;
479                                        }
480                                        else if(o1.getSecond()>o2.getSecond()){
481                                                return 1;
482                                        }
483                                        else{
484                                                return 0;
485                                        }
486                                }
487                        }
488                });
489
490
491
492                for(Pair<Integer> p: outList){
493                        int i = p.getFirst();
494                        int j = p.getSecond();
495                        BridgeType btype = null;
496                        // Now do the bonding
497                        if ((isBonded(i-1,j) && isBonded(j,i+1)) ||
498                                        (isBonded(j-1,i) && isBonded(i,j+1))) {
499                                btype = BridgeType.parallel;
500                        }
501                        else if ((isBonded(i,j) && isBonded(j,i)) ||
502                                        (isBonded(i-1,j+1) && (isBonded(j-1,i+1)))) {
503                                btype = BridgeType.antiparallel;
504                        }
505                        if (btype != null){
506                                registerBridge(i, j, btype);
507                        }
508                }
509
510
511        }
512
513        private void detectBends() {
514
515                for (int i = 2 ; i < groups.length-2 ;i++){
516
517                        //Check if all atoms form peptide bonds (backbone discontinuity)
518                        boolean bonded = true;
519                        for (int k=0; k<4; k++){
520                                int index = i+k-2;
521                                Atom C = groups[index].getC();
522                                Atom N = groups[index+1].getN();
523                                //Peptide bond C-N
524                                if (Calc.getDistance(C, N) > MAX_PEPTIDE_BOND_LENGTH){
525                                        bonded = false;
526                                        break;
527                                }
528                        }
529                        if (!bonded) continue;
530
531                        SecStrucGroup im2 = groups[i-2];
532                        SecStrucGroup g = groups[i];
533                        SecStrucGroup ip2 = groups[i+2];
534
535                        Atom caim2 = im2.getCA();
536                        Atom cag   = g.getCA();
537                        Atom caip2 = ip2.getCA();
538
539                        //Create vectors ( Ca i to Ca i-2 ) ; ( Ca i to CA i + 2 )
540                        Atom caminus2 = Calc.subtract(caim2,cag);
541                        Atom caplus2  = Calc.subtract(cag,caip2);
542
543                        double angle = Calc.angle(caminus2, caplus2);
544
545                        SecStrucState state = getSecStrucState(i);
546                        state.setKappa((float) angle);
547
548                        //Angles = 360 should be discarded
549                        if (angle > 70.0 && angle < 359.99) {
550                                setSecStrucType(i, SecStrucType.bend);
551                                state.setBend(true);
552                        }
553                }
554        }
555
556        private void calculateDihedralAngles()  {
557
558                // dihedral angles
559                // phi: C-N-CA-C
560                // psi: N-CA-C-N
561                // Chi1: N-CA-CB-CG, N-CA-CB-OG(SER),N-CA-CB-OG1(Thr),
562                // N-CA-CB-CG1(ILE/VAL), N-CA-CB-SG(CYS)
563                // Omega: CA-C-N-CA
564
565                for (int i=0 ; i < groups.length-1 ;  i++){
566
567                        SecStrucGroup a = groups[i];
568                        SecStrucGroup b = groups[i+1];
569
570                        Atom a_N   = a.getN();
571                        Atom a_CA  = a.getCA();
572                        Atom a_C  = a.getC();
573
574                        Atom b_N  = b.getN();
575                        Atom b_CA = b.getCA();
576                        Atom b_C  = b.getC();
577
578                        double phi = Calc.torsionAngle(a_C,b_N,b_CA,b_C);
579                        double psi = Calc.torsionAngle(a_N,a_CA,a_C,b_N);
580                        double omega = Calc.torsionAngle(a_CA,a_C,b_N,b_CA);
581
582                        SecStrucState state1 = (SecStrucState)
583                                        a.getProperty(Group.SEC_STRUC);
584                        SecStrucState state2 = (SecStrucState)
585                                        b.getProperty(Group.SEC_STRUC);
586
587                        state2.setPhi(phi);
588                        state1.setPsi(psi);
589                        state1.setOmega(omega);
590                }
591        }
592
593        @Override
594        public String toString() {
595                return printDSSP();
596        }
597
598        /**
599         * Generate a DSSP file format ouput String of this SS prediction.
600         * @return String in DSSP output file format
601         */
602        public String printDSSP() {
603
604                StringBuffer buf = new StringBuffer();
605                String nl = System.getProperty("line.separator");
606
607                //Header Line
608                buf.append("==== Secondary Structure Definition by BioJava"
609                                + " DSSP implementation, Version October 2015 ===="+nl);
610
611                //First line with column definition
612                buf.append("  #  RESIDUE AA STRUCTURE BP1 BP2  ACC     "
613                                + "N-H-->O    O-->H-N    N-H-->O    O-->H-N    "
614                                + "TCO  KAPPA ALPHA  PHI    PSI    "
615                                + "X-CA   Y-CA   Z-CA ");
616
617                for (int i =0 ; i < groups.length ;i++){
618                        buf.append(nl);
619                        SecStrucState ss = getSecStrucState(i);
620                        buf.append(ss.printDSSPline(i));
621                }
622
623                return buf.toString();
624        }
625
626        /**
627         * Generate a summary of this SS prediction with information about
628         * the three types of helix turns in different row sequences.
629         * <p>
630         * This is similar to the summary output of Jmol, and useful to visualize
631         * the helix patterns.
632         *
633         * @return String helix summary
634         */
635        public String printHelixSummary() {
636
637                StringBuffer g = new StringBuffer(); //3-10 helix
638                StringBuffer h = new StringBuffer(); //alpha helix
639                StringBuffer i = new StringBuffer(); //pi-helix
640                StringBuffer ss = new StringBuffer(); //SS summary
641                StringBuffer aa = new StringBuffer(); //AA one-letter
642                String nl = System.getProperty("line.separator");
643
644                g.append(       "3 turn: ");
645                h.append(       "4 turn: ");
646                i.append(       "5 turn: ");
647                ss.append(      "SS:     ");
648                aa.append(      "AA:     ");
649
650                for (int k = 0; k < groups.length; k++){
651
652                        SecStrucState state = getSecStrucState(k);
653                        g.append(state.getTurn()[0]);
654                        h.append(state.getTurn()[1]);
655                        i.append(state.getTurn()[2]);
656                        ss.append(state.getType());
657                        aa.append(StructureTools.get1LetterCode(groups[k].getPDBName()));
658                }
659
660                return g.toString()+nl+h.toString()+nl+
661                                i.toString()+nl+ss.toString()+nl+aa.toString();
662        }
663
664        /**
665         * Generate a FASTA sequence with the SS annotation letters in the
666         * aminoacid sequence order.
667         * @return String in FASTA sequence format
668         */
669        public String printFASTA() {
670
671                StringBuffer buf = new StringBuffer();
672                String nl = System.getProperty("line.separator");
673                buf.append(">"+groups[0].getChain().getStructure().getIdentifier()+nl);
674
675                for (int g = 0; g < groups.length; g++){
676                        buf.append(getSecStrucState(g).getType());
677                }
678                return buf.toString();
679        }
680
681        @Override
682        public int hashCode() {
683            final int prime = 31;
684            int result = 1;
685            result = prime * result + Arrays.hashCode(atoms);
686            return result;
687        }
688
689        @Override
690        public boolean equals(Object o){
691
692                if (!(o instanceof SecStrucCalc)) return false;
693                else {
694                        SecStrucCalc ss = (SecStrucCalc) o;
695                        if (groups.length != ss.groups.length) return false;
696
697                        for (int g=0; g<groups.length; g++){
698                                SecStrucInfo g1 = getSecStrucState(g);
699                                SecStrucInfo g2 = ss.getSecStrucState(g);
700                                if (!g1.equals(g2)) return false;
701                        }
702                        return true;
703                }
704        }
705
706        private static SecStrucGroup[] initGroupArray(Structure s, int modelId) {
707                List<SecStrucGroup> groupList = new ArrayList<SecStrucGroup>();
708                //
709                for ( Chain c : s.getChains(modelId)){
710
711                        for (Group g : c.getAtomGroups()){
712
713                                //We can also calc secstruc if it is a modified amino acid
714                                if ( g.hasAminoAtoms()) {
715
716                                        SecStrucGroup sg = new SecStrucGroup();
717                                        sg.setResidueNumber(g.getResidueNumber());
718                                        sg.setPDBFlag(true);
719                                        sg.setPDBName(g.getPDBName());
720                                        sg.setChain(g.getChain());
721
722                                        Atom N = g.getAtom(StructureTools.N_ATOM_NAME);
723                                        Atom CA =  g.getAtom(StructureTools.CA_ATOM_NAME);
724                                        Atom C = g.getAtom(StructureTools.C_ATOM_NAME);
725                                        Atom O =  g.getAtom(StructureTools.O_ATOM_NAME);
726                                        if ( N == null || CA == null || C == null || O == null)
727                                                continue;
728
729                                        sg.setN((Atom)   N.clone());
730                                        sg.setCA((Atom) CA.clone());
731                                        sg.setC((Atom)   C.clone());
732                                        sg.setO((Atom)  O.clone());
733                                        sg.setOriginal(g);
734
735                                        SecStrucState state = new SecStrucState(sg,
736                                                        SecStrucInfo.BIOJAVA_ASSIGNMENT,
737                                                        SecStrucType.coil);
738
739                                        sg.setProperty(Group.SEC_STRUC, state);
740                                        groupList.add(sg);
741                                }
742                        }
743
744                }
745                return groupList.toArray(new SecStrucGroup[groupList.size()]);
746        }
747
748        /**
749         * Calculate the coordinates of the H atoms. They are usually
750         * missing in the PDB files as only few experimental methods allow
751         * to resolve their location.
752         */
753        private void calculateHAtoms() throws StructureException {
754
755                for ( int i = 0 ; i < groups.length-1  ; i++) {
756
757                        SecStrucGroup a  = groups[i];
758                        SecStrucGroup b  = groups[i+1];
759
760                        if ( !b.hasAtom("H") ) {
761                                //Atom H = calc_H(a.getC(), b.getN(), b.getCA());
762                                Atom H = calcSimple_H(a.getC(), a.getO(), b.getN());
763                                b.setH(H);
764                        }
765                }
766        }
767
768
769
770        /**
771         * Calculate the HBonds between different groups.
772         * see Creighton page 147 f
773         * Modified to use only the contact map
774         */
775        private void calculateHBonds() {
776                /**
777                 * More efficient method for calculating C-Alpha pairs
778                 */
779                if (groups.length < 5) return;
780                Iterator<AtomContact> otu = contactSet.iterator();
781                while(otu.hasNext()){
782                        AtomContact ac = otu.next();
783                        Pair<Atom> pair = ac.getPair();
784                        Group g1 = pair.getFirst().getGroup();
785                        Group g2 = pair.getSecond().getGroup();
786                        // Now I need to get the index of the Group in the list groups
787                        int i = indResMap.get(g1.getResidueNumber());
788                        int j = indResMap.get(g2.getResidueNumber());
789                        // Now check this
790                        checkAddHBond(i,j);
791                        //"backwards" hbonds are not allowed
792                        if (j!=(i+1)) checkAddHBond(j,i);
793                }
794        }
795
796        private void checkAddHBond(int i, int j){
797
798                SecStrucGroup one = groups[i];
799
800                if (one.getPDBName().equals("PRO")){
801                        logger.debug("Ignore: PRO {}", one.getResidueNumber());
802                        return;
803                }
804                if (!one.hasAtom("H")) {
805                        logger.debug("Residue {} has no H",one.getResidueNumber());
806                        return;
807                }
808
809                SecStrucGroup two = groups[j];
810
811                double energy = 0;
812
813                try {
814                        energy = calculateHBondEnergy(one,two);
815                } catch (Exception e){
816                        logger.warn("Energy calculation failed", e);
817                        return;
818                }
819                logger.debug("Energy between positions ({},{}): ",i,j,energy);
820
821                trackHBondEnergy(i,j,energy);
822        }
823
824        /**
825         * Calculate HBond energy of two groups in cal/mol
826         * see Creighton page 147 f
827         * <p>
828         * Jeffrey, George A., An introduction to hydrogen bonding,
829         * Oxford University Press, 1997.
830         * categorizes hbonds with donor-acceptor distances of
831         * 2.2-2.5 &aring; as "strong, mostly covalent",
832         * 2.5-3.2 &aring; as "moderate, mostly electrostatic",
833         * 3.2-4.0 &aring; as "weak, electrostatic".
834         * Energies are given as 40-14, 15-4, and <4 kcal/mol respectively.
835         */
836        private static double calculateHBondEnergy(SecStrucGroup one,
837                        SecStrucGroup two) {
838
839                Atom N = one.getN();
840                Atom H = one.getH();
841
842                Atom O = two.getO();
843                Atom C = two.getC();
844
845                double dno = Calc.getDistance(O,N);
846                double dhc = Calc.getDistance(C,H);
847                double dho = Calc.getDistance(O,H);
848                double dnc = Calc.getDistance(C,N);
849
850                logger.debug(" cccc: {} {} {} {} O ({})..N ({}):{}  |  ho:{} - hc:{} + nc:{} - no:{}",
851                                one.getResidueNumber(),one.getPDBName(),two.getResidueNumber(),two.getPDBName(),
852                                O.getPDBserial(),N.getPDBserial(),dno,dho,dhc,dnc,dno);
853
854                //there seems to be a contact!
855                if ( (dno < MINDIST) || (dhc < MINDIST) ||
856                                (dnc < MINDIST) || (dno < MINDIST)) {
857                        return HBONDLOWENERGY;
858                }
859
860                double e1 = Q / dho - Q / dhc;
861                double e2 = Q / dnc - Q / dno;
862
863                double energy = e1 + e2;
864
865                logger.debug(" N ({}) O({}): {} : {} ",N.getPDBserial(),O.getPDBserial(),(float) dno,energy);
866
867                //Avoid too strong energy
868                if (energy > HBONDLOWENERGY) return energy;
869
870                return HBONDLOWENERGY ;
871        }
872
873        /**
874         * Store Hbonds in the Groups.
875         * DSSP allows two HBonds per aminoacids to allow bifurcated bonds.
876         */
877        private  void trackHBondEnergy(int i, int j, double energy) {
878
879                if (groups[i].getPDBName().equals("PRO")) {
880                        logger.debug("Ignore: PRO {}",groups[i].getResidueNumber());
881                        return;
882                }
883
884                SecStrucState stateOne = getSecStrucState(i);
885                SecStrucState stateTwo = getSecStrucState(j);
886
887                double acc1e = stateOne.getAccept1().getEnergy();
888                double acc2e = stateOne.getAccept2().getEnergy();
889
890                double don1e = stateTwo.getDonor1().getEnergy();
891                double don2e = stateTwo.getDonor2().getEnergy();
892
893                //Acceptor: N-H-->O
894                if (energy < acc1e) {
895                        logger.debug("{} < {}",energy,acc1e);
896                        stateOne.setAccept2(stateOne.getAccept1());
897
898                        HBond bond = new HBond();
899                        bond.setEnergy(energy);
900                        bond.setPartner(j);
901
902                        stateOne.setAccept1(bond);
903
904                } else if ( energy < acc2e ) {
905                        logger.debug("{} < {}",energy,acc2e);
906
907                        HBond bond = new HBond();
908                        bond.setEnergy(energy);
909                        bond.setPartner(j);
910
911                        stateOne.setAccept2(bond);
912                }
913
914
915                //The other side of the bond: donor O-->N-H
916                if (energy <  don1e) {
917                        logger.debug("{} < {}",energy,don1e);
918                        stateTwo.setDonor2(stateTwo.getDonor1());
919
920                        HBond bond = new HBond();
921                        bond.setEnergy(energy);
922                        bond.setPartner(i);
923
924                        stateTwo.setDonor1(bond);
925
926                } else if ( energy < don2e ) {
927                        logger.debug("{} < {}",energy,don2e);
928
929                        HBond bond = new HBond();
930                        bond.setEnergy(energy);
931                        bond.setPartner(i);
932
933                        stateTwo.setDonor2(bond);
934                }
935        }
936
937        /**
938         * Detect helical turn patterns.
939         */
940        private void calculateTurns(){
941
942                for (int i = 0 ; i< groups.length; i++){
943                        for (int turn = 3; turn <= 5; turn++) {
944
945                                if (i+turn >= groups.length) continue;
946
947                                //Check for H bond from NH(i+n) to CO(i)
948                                if (isBonded(i, i+turn)) {
949                                        logger.debug("Turn at ({},{}) turn {}",i,(i+turn),turn);
950                                        getSecStrucState(i).setTurn('>', turn);
951                                        getSecStrucState(i+turn).setTurn('<', turn);
952                                        //Bracketed residues get the helix number
953                                        for (int j=i+1; j<i+turn; j++){
954                                                Integer t = turn;
955                                                char helix = t.toString().charAt(0);
956                                                getSecStrucState(j).setTurn(helix, turn);
957                                        }
958                                }
959                        }
960                }
961        }
962
963        /**
964         * Test if two groups are forming an H-Bond. The bond tested is
965         * from the CO of group i to the NH of group j. Acceptor (i) and
966         * donor (j). The donor of i has to be j, and the acceptor of j
967         * has to be i.
968         * DSSP defines H-Bonds if the energy < -500 cal/mol.
969         *
970         * @param i group one
971         * @param j group two
972         * @return flag if the two are forming an Hbond
973         */
974        private boolean isBonded(int i, int j) {
975
976                SecStrucState one = getSecStrucState(i);
977                SecStrucState two = getSecStrucState(j);
978
979                double don1e = one.getDonor1().getEnergy();
980                double don2e = one.getDonor2().getEnergy();
981                double acc1e = two.getAccept1().getEnergy();
982                double acc2e = two.getAccept2().getEnergy();
983
984                int don1p = one.getDonor1().getPartner();
985                int don2p = one.getDonor2().getPartner();
986                int acc1p = two.getAccept1().getPartner();
987                int acc2p = two.getAccept2().getPartner();
988
989                //Either donor from i is j, or accept from j is i
990                boolean hbond = (don1p == j && don1e < HBONDHIGHENERGY) ||
991                                (don2p == j && don2e < HBONDHIGHENERGY) ||
992                                (acc1p == i && acc1e < HBONDHIGHENERGY) ||
993                                (acc2p == i && acc2e < HBONDHIGHENERGY);
994
995                if (hbond){
996                        logger.debug("*** H-bond from CO of {} to NH of {}", i, j);
997                        return true;
998                }
999                return false ;
1000        }
1001
1002        /**
1003         * Use unit vectors NC and NCalpha Add them. Calc unit vector and
1004         * substract it from N.
1005         * C coordinates are from amino acid i-1
1006         * N, CA atoms from amino acid i
1007         *
1008         * @link http://openbioinformatics.blogspot.com/
1009         *              2009/08/how-to-calculate-h-atoms-for-nitrogens.html
1010         */
1011        @SuppressWarnings("unused")
1012        private static Atom calc_H(Atom C, Atom N, Atom CA)
1013                        throws StructureException {
1014
1015                Atom nc  = Calc.subtract(N,C);
1016                Atom nca = Calc.subtract(N,CA);
1017
1018                Atom u_nc  = Calc.unitVector(nc)   ;
1019                Atom u_nca = Calc.unitVector(nca);
1020
1021                Atom added = Calc.add(u_nc,u_nca);
1022
1023                Atom U = Calc.unitVector(added);
1024
1025                // according to Creighton distance N-H is 1.03 +/- 0.02A
1026                Atom H = Calc.add(N,U);
1027
1028                H.setName("H");
1029                // this atom does not have a pdbserial number ...
1030                return H;
1031
1032        }
1033
1034        private static Atom calcSimple_H(Atom c, Atom o, Atom n)  {
1035
1036                Atom h = Calc.subtract(c,o);
1037                double dist = Calc.getDistance(o,c);
1038                //System.out.println(dist);
1039                double x = n.getX() + h.getX() / dist;
1040                double y = n.getY() + h.getY() / dist;
1041                double z = n.getZ() + h.getZ() / dist;
1042
1043                h.setX(x);
1044                h.setY(y);
1045                h.setZ(z);
1046
1047                h.setName("H");
1048                return h;
1049        }
1050
1051        private void buildHelices(){
1052
1053                //Alpha-helix (i+4), 3-10-helix (i+3), Pi-helix (i+5)
1054                checkSetHelix(4, SecStrucType.helix4);
1055                checkSetHelix(3, SecStrucType.helix3);
1056                checkSetHelix(5, SecStrucType.helix5);
1057
1058                checkSetTurns();
1059        }
1060
1061        private void checkSetTurns() {
1062
1063                SecStrucType type = SecStrucType.turn;
1064
1065                for (int idx = 0; idx < 3; idx++) {
1066                        for (int i = 0; i < groups.length-1; i++) {
1067
1068                                SecStrucState state = getSecStrucState(i);
1069                                char[] turn = state.getTurn();
1070
1071                                //Any turn opening matters
1072                                if (turn[idx] == '>' || turn[idx] == 'X') {
1073                                        //Mark following n residues as turn
1074                                        for (int k=1; k<idx+3; k++){
1075                                                setSecStrucType(i+k, type);
1076                                        }
1077                                        if (!DSSP_HELICES) {
1078                                                setSecStrucType(i, type);
1079                                                setSecStrucType(i+idx+3, type);
1080                                        }
1081                                }
1082                        }
1083                }
1084        }
1085
1086        /**
1087         * A minimal helix is defined by two consecutive n-turns.
1088         * For example, a 4-helix, of minimal length 4 from residues
1089         * i to (i+3), requires turns (of type 4) at residues (i-1) and i.
1090         * <p>
1091         * Note that the orignal DSSP implementation does not assign
1092         * helix type to residue (i-1) and residue (i+n+1), although
1093         * they contain a helix turn. As they state in the original paper,
1094         * "the helices are one residue shorter than they would be according
1095         * to rule 6.3 of IUPAC-IUB".
1096         *
1097         * @param n
1098         * @param type
1099         */
1100        private void checkSetHelix(int n, SecStrucType type){
1101
1102                int idx = n - 3;
1103                logger.debug("Set helix {} {} {}", type, n, idx);
1104
1105                for (int i = 1; i < groups.length-n; i++) {
1106
1107                        SecStrucState state = getSecStrucState(i);
1108                        SecStrucState previousState = getSecStrucState(i-1);
1109
1110                        //Check that no other helix was assgined to this range
1111                        if (state.getType().compareTo(type) < 0) continue;
1112                        if (getSecStrucState(i+1).getType().compareTo(type) < 0) continue;
1113
1114                        char turn = state.getTurn()[idx];
1115                        char pturn = previousState.getTurn()[idx];
1116
1117                        //Two consecutive n-turns present to define a n-helix
1118                        if ((turn=='>' || turn=='X') && (pturn=='>' || pturn=='X')) {
1119                                //Mark following n residues as turn
1120                                for (int k=0; k<n; k++){
1121                                        setSecStrucType(i+k, type);
1122                                }
1123                                if (!DSSP_HELICES) {
1124                                        setSecStrucType(i-1, type);
1125                                        setSecStrucType(i+n, type);
1126                                }
1127                        }
1128                }
1129        }
1130
1131        /**
1132         * Set the new type only if it has more preference than the
1133         * current residue SS type.
1134         * @param pos
1135         * @param type
1136         */
1137        private void setSecStrucType(int pos, SecStrucType type){
1138                SecStrucState ss = getSecStrucState(pos);
1139                if (type.compareTo(ss.getType()) < 0) ss.setType(type);
1140        }
1141
1142        private SecStrucState getSecStrucState(int pos){
1143                Group g = groups[pos];
1144                SecStrucState state = (SecStrucState) g.getProperty(Group.SEC_STRUC);
1145                return state;
1146        }
1147
1148}