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