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.size() + " ladders!"); 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 " + i + " to " + 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" + i + " and " + j 366 + ". DSSP assignment might differ."); 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 "+one.getResidueNumber()+" has no H"); 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: " + one.getResidueNumber() + 852 " " + one.getPDBName() + " " +two.getResidueNumber()+ 853 " " + two.getPDBName() + String.format(" O ("+ 854 O.getPDBserial()+")..N ("+ N.getPDBserial()+ 855 "):%4.1f | ho:%4.1f - hc:%4.1f + nc:%4.1f - no:%4.1f ", 856 dno,dho,dhc,dnc,dno)); 857 858 //there seems to be a contact! 859 if ( (dno < MINDIST) || (dhc < MINDIST) || 860 (dnc < MINDIST) || (dno < MINDIST)) { 861 return HBONDLOWENERGY; 862 } 863 864 double e1 = Q / dho - Q / dhc; 865 double e2 = Q / dnc - Q / dno; 866 867 double energy = e1 + e2; 868 869 logger.debug(String.format(" N (%d) O(%d): %4.1f : %4.2f ", 870 N.getPDBserial(),O.getPDBserial(), (float) dno, energy)); 871 872 //Avoid too strong energy 873 if (energy > HBONDLOWENERGY) return energy; 874 875 return HBONDLOWENERGY ; 876 } 877 878 /** 879 * Store Hbonds in the Groups. 880 * DSSP allows two HBonds per aminoacids to allow bifurcated bonds. 881 */ 882 private void trackHBondEnergy(int i, int j, double energy) { 883 884 if (groups[i].getPDBName().equals("PRO")) { 885 logger.debug("Ignore: PRO " + groups[i].getResidueNumber()); 886 return; 887 } 888 889 SecStrucState stateOne = getSecStrucState(i); 890 SecStrucState stateTwo = getSecStrucState(j); 891 892 double acc1e = stateOne.getAccept1().getEnergy(); 893 double acc2e = stateOne.getAccept2().getEnergy(); 894 895 double don1e = stateTwo.getDonor1().getEnergy(); 896 double don2e = stateTwo.getDonor2().getEnergy(); 897 898 //Acceptor: N-H-->O 899 if (energy < acc1e) { 900 logger.debug(energy +"<"+acc1e); 901 stateOne.setAccept2(stateOne.getAccept1()); 902 903 HBond bond = new HBond(); 904 bond.setEnergy(energy); 905 bond.setPartner(j); 906 907 stateOne.setAccept1(bond); 908 909 } else if ( energy < acc2e ) { 910 logger.debug(energy +"<"+acc2e); 911 912 HBond bond = new HBond(); 913 bond.setEnergy(energy); 914 bond.setPartner(j); 915 916 stateOne.setAccept2(bond); 917 } 918 919 920 //The other side of the bond: donor O-->N-H 921 if (energy < don1e) { 922 logger.debug(energy +"<"+don1e); 923 stateTwo.setDonor2(stateTwo.getDonor1()); 924 925 HBond bond = new HBond(); 926 bond.setEnergy(energy); 927 bond.setPartner(i); 928 929 stateTwo.setDonor1(bond); 930 931 } else if ( energy < don2e ) { 932 logger.debug(energy +"<"+don2e); 933 934 HBond bond = new HBond(); 935 bond.setEnergy(energy); 936 bond.setPartner(i); 937 938 stateTwo.setDonor2(bond); 939 } 940 } 941 942 /** 943 * Detect helical turn patterns. 944 */ 945 private void calculateTurns(){ 946 947 for (int i = 0 ; i< groups.length; i++){ 948 for (int turn = 3; turn <= 5; turn++) { 949 950 if (i+turn >= groups.length) continue; 951 952 //Check for H bond from NH(i+n) to CO(i) 953 if (isBonded(i, i+turn)) { 954 logger.debug("Turn at ("+i+","+(i+turn)+") turn "+turn); 955 getSecStrucState(i).setTurn('>', turn); 956 getSecStrucState(i+turn).setTurn('<', turn); 957 //Bracketed residues get the helix number 958 for (int j=i+1; j<i+turn; j++){ 959 Integer t = turn; 960 char helix = t.toString().charAt(0); 961 getSecStrucState(j).setTurn(helix, turn); 962 } 963 } 964 } 965 } 966 } 967 968 /** 969 * Test if two groups are forming an H-Bond. The bond tested is 970 * from the CO of group i to the NH of group j. Acceptor (i) and 971 * donor (j). The donor of i has to be j, and the acceptor of j 972 * has to be i. 973 * DSSP defines H-Bonds if the energy < -500 cal/mol. 974 * 975 * @param i group one 976 * @param j group two 977 * @return flag if the two are forming an Hbond 978 */ 979 private boolean isBonded(int i, int j) { 980 981 SecStrucState one = getSecStrucState(i); 982 SecStrucState two = getSecStrucState(j); 983 984 double don1e = one.getDonor1().getEnergy(); 985 double don2e = one.getDonor2().getEnergy(); 986 double acc1e = two.getAccept1().getEnergy(); 987 double acc2e = two.getAccept2().getEnergy(); 988 989 int don1p = one.getDonor1().getPartner(); 990 int don2p = one.getDonor2().getPartner(); 991 int acc1p = two.getAccept1().getPartner(); 992 int acc2p = two.getAccept2().getPartner(); 993 994 //Either donor from i is j, or accept from j is i 995 boolean hbond = (don1p == j && don1e < HBONDHIGHENERGY) || 996 (don2p == j && don2e < HBONDHIGHENERGY) || 997 (acc1p == i && acc1e < HBONDHIGHENERGY) || 998 (acc2p == i && acc2e < HBONDHIGHENERGY); 999 1000 if (hbond){ 1001 logger.debug("*** H-bond from CO of " + i + " to NH of " + j); 1002 return true; 1003 } 1004 return false ; 1005 } 1006 1007 /** 1008 * Use unit vectors NC and NCalpha Add them. Calc unit vector and 1009 * substract it from N. 1010 * C coordinates are from amino acid i-1 1011 * N, CA atoms from amino acid i 1012 * 1013 * @link http://openbioinformatics.blogspot.com/ 1014 * 2009/08/how-to-calculate-h-atoms-for-nitrogens.html 1015 */ 1016 @SuppressWarnings("unused") 1017 private static Atom calc_H(Atom C, Atom N, Atom CA) 1018 throws StructureException { 1019 1020 Atom nc = Calc.subtract(N,C); 1021 Atom nca = Calc.subtract(N,CA); 1022 1023 Atom u_nc = Calc.unitVector(nc) ; 1024 Atom u_nca = Calc.unitVector(nca); 1025 1026 Atom added = Calc.add(u_nc,u_nca); 1027 1028 Atom U = Calc.unitVector(added); 1029 1030 // according to Creighton distance N-H is 1.03 +/- 0.02A 1031 Atom H = Calc.add(N,U); 1032 1033 H.setName("H"); 1034 // this atom does not have a pdbserial number ... 1035 return H; 1036 1037 } 1038 1039 private static Atom calcSimple_H(Atom c, Atom o, Atom n) { 1040 1041 Atom h = Calc.subtract(c,o); 1042 double dist = Calc.getDistance(o,c); 1043 //System.out.println(dist); 1044 double x = n.getX() + h.getX() / dist; 1045 double y = n.getY() + h.getY() / dist; 1046 double z = n.getZ() + h.getZ() / dist; 1047 1048 h.setX(x); 1049 h.setY(y); 1050 h.setZ(z); 1051 1052 h.setName("H"); 1053 return h; 1054 } 1055 1056 private void buildHelices(){ 1057 1058 //Alpha-helix (i+4), 3-10-helix (i+3), Pi-helix (i+5) 1059 checkSetHelix(4, SecStrucType.helix4); 1060 checkSetHelix(3, SecStrucType.helix3); 1061 checkSetHelix(5, SecStrucType.helix5); 1062 1063 checkSetTurns(); 1064 } 1065 1066 private void checkSetTurns() { 1067 1068 SecStrucType type = SecStrucType.turn; 1069 1070 for (int idx = 0; idx < 3; idx++) { 1071 for (int i = 0; i < groups.length-1; i++) { 1072 1073 SecStrucState state = getSecStrucState(i); 1074 char[] turn = state.getTurn(); 1075 1076 //Any turn opening matters 1077 if (turn[idx] == '>' || turn[idx] == 'X') { 1078 //Mark following n residues as turn 1079 for (int k=1; k<idx+3; k++){ 1080 setSecStrucType(i+k, type); 1081 } 1082 if (!DSSP_HELICES) { 1083 setSecStrucType(i, type); 1084 setSecStrucType(i+idx+3, type); 1085 } 1086 } 1087 } 1088 } 1089 } 1090 1091 /** 1092 * A minimal helix is defined by two consecutive n-turns. 1093 * For example, a 4-helix, of minimal length 4 from residues 1094 * i to (i+3), requires turns (of type 4) at residues (i-1) and i. 1095 * <p> 1096 * Note that the orignal DSSP implementation does not assign 1097 * helix type to residue (i-1) and residue (i+n+1), although 1098 * they contain a helix turn. As they state in the original paper, 1099 * "the helices are one residue shorter than they would be according 1100 * to rule 6.3 of IUPAC-IUB". 1101 * 1102 * @param n 1103 * @param type 1104 */ 1105 private void checkSetHelix(int n, SecStrucType type){ 1106 1107 int idx = n - 3; 1108 logger.debug("Set helix " + type + " " + n + " " + idx); 1109 1110 for (int i = 1; i < groups.length-n; i++) { 1111 1112 SecStrucState state = getSecStrucState(i); 1113 SecStrucState previousState = getSecStrucState(i-1); 1114 1115 //Check that no other helix was assgined to this range 1116 if (state.getType().compareTo(type) < 0) continue; 1117 if (getSecStrucState(i+1).getType().compareTo(type) < 0) continue; 1118 1119 char turn = state.getTurn()[idx]; 1120 char pturn = previousState.getTurn()[idx]; 1121 1122 //Two consecutive n-turns present to define a n-helix 1123 if ((turn=='>' || turn=='X') && (pturn=='>' || pturn=='X')) { 1124 //Mark following n residues as turn 1125 for (int k=0; k<n; k++){ 1126 setSecStrucType(i+k, type); 1127 } 1128 if (!DSSP_HELICES) { 1129 setSecStrucType(i-1, type); 1130 setSecStrucType(i+n, type); 1131 } 1132 } 1133 } 1134 } 1135 1136 /** 1137 * Set the new type only if it has more preference than the 1138 * current residue SS type. 1139 * @param pos 1140 * @param type 1141 */ 1142 private void setSecStrucType(int pos, SecStrucType type){ 1143 SecStrucState ss = getSecStrucState(pos); 1144 if (type.compareTo(ss.getType()) < 0) ss.setType(type); 1145 } 1146 1147 private SecStrucState getSecStrucState(int pos){ 1148 Group g = groups[pos]; 1149 SecStrucState state = (SecStrucState) g.getProperty(Group.SEC_STRUC); 1150 return state; 1151 } 1152 1153}