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 å as "strong, mostly covalent", 833 * 2.5-3.2 å as "moderate, mostly electrostatic", 834 * 3.2-4.0 å 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}