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