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