org.biojava.nbio.structure.secstruc.SecStrucCalc Maven / Gradle / Ivy
Show all versions of biojava-structure Show documentation
/*
* BioJava development code
*
* This code may be freely distributed and modified under the
* terms of the GNU Lesser General Public Licence. This should
* be distributed with the code. If you do not have a copy,
* see:
*
* http://www.gnu.org/copyleft/lesser.html
*
* Copyright for this code is held jointly by the individual
* authors. These should be listed in @author doc comments.
*
* For more information on the BioJava project and its aims,
* or to join the biojava-l mailing list, visit the home page
* at:
*
* http://www.biojava.org/
*
*/
package org.biojava.nbio.structure.secstruc;
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.contact.AtomContact;
import org.biojava.nbio.structure.contact.AtomContactSet;
import org.biojava.nbio.structure.contact.Grid;
import org.biojava.nbio.structure.contact.Pair;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import java.util.Map;
/**
* Calculate and assign the secondary structure (SS) to the
* Groups of a Structure object. This object also stores the result
* of the calculation.
*
* The rules for SS calculation are the ones defined by DSSP:
* Kabsch,W. and Sander,C. (1983) Biopolymers 22, 2577-2637.
* Original DSSP article see at:
* dssp.pdf.
* Some parts are also taken from: T.E.Creighton, Proteins -
* Structure and Molecular Properties, 2nd Edition, Freeman 1994.
*
* @author Andreas Prlic
* @author Aleix Lafita
* @author Anthony Bradley
*
*/
public class SecStrucCalc {
/**
* DSSP assigns helices one residue shorter at each end, because the
* residues at (i-1) and (i+n+1) are not assigned helix type although
* they contain a consistent turn (H-bond). If this parameter
* is true, the helices will be the length of the original DSSP
* convention. If it is false, they will be two residue longer.
*/
private static final boolean DSSP_HELICES = true;
private static final Logger logger =
LoggerFactory.getLogger(SecStrucCalc.class);
/** min distance between two residues */
public static final double MINDIST = 0.5;
/** min distance of two CA atoms if H-bonds are allowed to form */
public static final double CA_MIN_DIST = 9.0;
/** max distance CA atoms in peptide bond (backbone discontinuity) */
public static final double MAX_PEPTIDE_BOND_LENGTH = 2.5;
/** Minimal H-bond energy in cal/mol */
public static final int HBONDLOWENERGY = -9900;
/** higher limit for H-bond energy */
public static final double HBONDHIGHENERGY = -500.0;
/**
* constant for electrostatic energy
*
* f * q1 * q2 * scale
* Q = -332 * 0.42 * 0.20 * 1000.0
*
*
* q1 and q2 are partial charges which are placed on the C,O
* (+q1,-q1) and N,H (-q2,+q2)
*/
public static final double Q = -27888.0;
// Three lists
private SecStrucGroup[] groups;
private List ladders;
private List bridges;
private Atom[] atoms;
// Added by Anthony - to speed up intergroup calculations
private AtomContactSet contactSet;
private Map indResMap;
public SecStrucCalc(){
ladders = new ArrayList<>();
bridges = new ArrayList<>();
}
/**
* Predicts the secondary structure of this Structure object,
* using a DSSP implementation.
*
* @param s Structure to predict the SS
* @param assign sets the SS information to the Groups of s
* @return a List of SS annotation objects
*/
public List calculate(Structure s, boolean assign)
throws StructureException {
List secstruc = new ArrayList<>();
for(int i=0; i();
bridges = new ArrayList<>();
groups = initGroupArray(s, i);
// Initialise the contact set for this structure
initContactSet();
if (groups.length < 5) {
// not enough groups to do anything
throw new StructureException("Not enough backbone groups in the"
+ " Structure to calculate the secondary structure ("
+ groups.length+" given, minimum 5)" );
}
calculateHAtoms();
calculateHBonds();
calculateDihedralAngles();
calculateTurns();
buildHelices();
detectBends();
detectStrands();
for (SecStrucGroup sg : groups){
SecStrucState ss = (SecStrucState)
sg.getProperty(Group.SEC_STRUC);
// Add to return list and assign to original if flag is true
secstruc.add(ss);
if (assign) sg.getOriginal().setProperty(Group.SEC_STRUC, ss);
}
}
return secstruc;
}
/**
* Function to generate the contact sets
*/
private void initContactSet() {
// Initialise an array of atoms
atoms = new Atom[groups.length];
// Remake this local var
indResMap = new HashMap<>();
for (int i=0 ; i < groups.length ; i++){
SecStrucGroup one = groups[i];
indResMap.put(one.getResidueNumber(), i);
atoms[i] = one.getCA();
}
Grid grid = new Grid(CA_MIN_DIST);
if(atoms.length==0){
contactSet = new AtomContactSet(CA_MIN_DIST);
}
else{
grid.addAtoms(atoms);
contactSet = grid.getAtomContacts();
}
}
/**
* Updated code to detect strands
*/
private void detectStrands() {
//Find all the beta bridges of the structure
findBridges();
//Create Ladders
createLadders();
//Detect beta bulges between ladders
connectLadders();
//AND store SS assignments for Sheets, Strands and Bridges
updateSheets();
}
private void createLadders(){
for (BetaBridge b : bridges){
boolean found = false;
for (Ladder ladder : ladders){
if (shouldExtendLadder(ladder, b)) {
found = true;
ladder.to++; //we go forward in this direction
switch(b.type){
case parallel:
ladder.lto++; //increment second strand
break;
case antiparallel:
ladder.lfrom--; //decrement second strand
break;
}
break;
}
}
if (!found){
//Create new ladder with a single Bridge
Ladder l = new Ladder();
l.from = b.partner1;
l.to = b.partner1;
l.lfrom = b.partner2;
l.lto = b.partner2;
l.btype = b.type;
ladders.add(l);
}
}
}
private void updateSheets() {
logger.debug(" got {} ladders!", ladders.size());
for (Ladder ladder : ladders){
logger.debug(ladder.toString());
for (int lcount = ladder.from; lcount <= ladder.to; lcount++) {
SecStrucState state = getSecStrucState(lcount);
SecStrucType stype = state.getType();
int diff = ladder.from - lcount;
int l2count = ladder.lfrom - diff ;
SecStrucState state2 = getSecStrucState(l2count);
SecStrucType stype2 = state2.getType();
if ( ladder.from != ladder.to ) {
setSecStrucType(lcount, SecStrucType.extended);
setSecStrucType(l2count, SecStrucType.extended);
}
else {
if ( !stype.isHelixType() &&
( !stype.equals(SecStrucType.extended)))
setSecStrucType(lcount,SecStrucType.bridge);
if ( ! stype2.isHelixType() &&
(! stype2.equals(SecStrucType.extended)))
setSecStrucType(l2count,SecStrucType.bridge);
}
}
// Check if two ladders are connected. both sides are 'E'
if (ladder.connectedTo == 0) continue;
Ladder conladder = ladders.get(ladder.connectedTo);
if (ladder.btype.equals(BridgeType.antiparallel)) {
/* set one side */
for (int lcount = ladder.from; lcount <= conladder.to;
lcount++) {
setSecStrucType(lcount, SecStrucType.extended);
}
/* set other side */
for (int lcount = conladder.lto;
lcount <= ladder.lfrom;
lcount++) {
setSecStrucType(lcount, SecStrucType.extended);
}
} else {
/* set one side */
for ( int lcount = ladder.from;
lcount <= conladder.to;
lcount++) {
setSecStrucType(lcount, SecStrucType.extended);
}
/* set other side */
for ( int lcount = ladder.lfrom;
lcount <= conladder.lto;
lcount++) {
setSecStrucType(lcount, SecStrucType.extended);
}
}
}
}
private void connectLadders() {
for (int i = 0 ; i < ladders.size(); i++) {
for ( int j = i ; j < ladders.size(); j++){
Ladder l1 = ladders.get(i);
Ladder l2 = ladders.get(j);
if (hasBulge(l1,l2)) {
l1.connectedTo = j;
l2.connectedFrom = i;
logger.debug("Bulge from {} to {}", i, j);
}
}
}
}
/**
* For beta structures, we define explicitly: a bulge-linked
* ladder consists of two (perfect) ladder or bridges of the
* same type connected by at most one extra residue on one
* strand and at most four extra residues on the other strand,
* all residues in bulge-linked ladders are marked "E,"
* including the extra residues.
*/
private boolean hasBulge(Ladder l1, Ladder l2) {
boolean bulge = ((l1.btype.equals(l2.btype)) &&
(l2.from - l1.to < 6) &&
(l1.to < l2.from) &&
(l2.connectedTo == 0));
if (!bulge) return bulge;
switch(l1.btype){
case parallel:
bulge = ( (l2.lfrom - l1.lto > 0) &&
((( l2.lfrom -l1.lto < 6) &&
(l2.from - l1.to < 3)) ||
( l2.lfrom - l1.lto <3)));
break;
case antiparallel:
bulge = ( (l1.lfrom - l2.lto > 0) &&
(((l1.lfrom -l2.lto < 6) &&
( l2.from - l1.to < 3)) ||
(l1.lfrom - l2.lto < 3)));
break;
}
return bulge;
}
private void registerBridge(int i, int j, BridgeType btype) {
BetaBridge bridge = new BetaBridge(i,j,btype);
boolean b1 = getSecStrucState(i).addBridge(bridge);
boolean b2 = getSecStrucState(j).addBridge(bridge);
if (!b1 && !b2)
logger.warn("Ignoring Bridge between residues {} and {}. DSSP assignment might differ.", i, j);
bridges.add(bridge);
}
/**
* Conditions to extend a ladder with a given beta Bridge:
* The bridge and ladder are of the same type.
* The smallest bridge residue is sequential to the first
* strand ladder.
* The second bridge residue is either sequential (parallel)
* or previous (antiparallel) to the second strand of the ladder
*
* @param ladder the ladder candidate to extend
* @param b the beta bridge that would extend the ladder
* @return true if the bridge b extends the ladder
*/
private boolean shouldExtendLadder(Ladder ladder, BetaBridge b) {
//Only extend if they are of the same type
boolean sameType = b.type.equals(ladder.btype);
if (!sameType) return false;
//Only extend if residue 1 is sequential to ladder strand
boolean sequential = (b.partner1 == ladder.to+1);
if (!sequential) return false;
switch(b.type){
case parallel:
//Residue 2 should be sequential to second strand
if (b.partner2 == ladder.lto+1) return true;
break;
case antiparallel:
//Residue 2 should be previous to second strand
if (b.partner2 == ladder.lfrom-1) return true;
break;
}
return false;
}
/**
* Two nonoverlapping stretches of three residues each, i-1,i,i+1 and
* j-1,j,j+1, form either a parallel or antiparallel bridge, depending on
* which of two basic patterns is matched. We assign a bridge between
* residues i and j if there are two H bonds characteristic of beta-
* structure; in particular:
*
* Parallel Bridge(i,j) =: [Hbond(i-1,j) and Hbond(j,i+1)]
* or [Hbond(j-1,i) and Hbond(i,j+1)]
*
* Antiparallel Bridge(i,j) =: [Hbond(i,j) and Hbond(j,i)]
* or [Hbond(i-1,j+1) and Hbond(j-1,i+1)]
*
* Optimised to use the contact set
*/
private void findBridges() {
// Get the interator of contacts
Iterator myIter = contactSet.iterator();
List> outList = new ArrayList<>();
// Now iterate through this
while(myIter.hasNext()){
// Get the next atom contact
AtomContact ac = myIter.next();
Group g1 = ac.getPair().getFirst().getGroup();
Group g2 = ac.getPair().getSecond().getGroup();
// Get the indices
int i = indResMap.get(g1.getResidueNumber());
int j = indResMap.get(g2.getResidueNumber());
// If i>j switch them over
if(i>j){
// Switch them over
int old = i;
i = j;
j = old;
}
// Only these
if(j thisPair = new Pair<>(i,j);
outList.add(thisPair);
}
//
Collections.sort(outList, new Comparator>() {
@Override
public int compare(Pair o1, Pair o2) {
if(o1.getFirst()o2.getFirst()){
return +1;
}
else{
if(o1.getSecond()o2.getSecond()){
return 1;
}
else{
return 0;
}
}
}
});
for(Pair p: outList){
int i = p.getFirst();
int j = p.getSecond();
BridgeType btype = null;
// Now do the bonding
if ((isBonded(i-1,j) && isBonded(j,i+1)) ||
(isBonded(j-1,i) && isBonded(i,j+1))) {
btype = BridgeType.parallel;
}
else if ((isBonded(i,j) && isBonded(j,i)) ||
(isBonded(i-1,j+1) && (isBonded(j-1,i+1)))) {
btype = BridgeType.antiparallel;
}
if (btype != null){
registerBridge(i, j, btype);
}
}
}
private void detectBends() {
for (int i = 2 ; i < groups.length-2 ;i++){
//Check if all atoms form peptide bonds (backbone discontinuity)
boolean bonded = true;
for (int k=0; k<4; k++){
int index = i+k-2;
Atom C = groups[index].getC();
Atom N = groups[index+1].getN();
//Peptide bond C-N
if (Calc.getDistance(C, N) > MAX_PEPTIDE_BOND_LENGTH){
bonded = false;
break;
}
}
if (!bonded) continue;
SecStrucGroup im2 = groups[i-2];
SecStrucGroup g = groups[i];
SecStrucGroup ip2 = groups[i+2];
Atom caim2 = im2.getCA();
Atom cag = g.getCA();
Atom caip2 = ip2.getCA();
//Create vectors ( Ca i to Ca i-2 ) ; ( Ca i to CA i + 2 )
Atom caminus2 = Calc.subtract(caim2,cag);
Atom caplus2 = Calc.subtract(cag,caip2);
double angle = Calc.angle(caminus2, caplus2);
SecStrucState state = getSecStrucState(i);
state.setKappa((float) angle);
//Angles = 360 should be discarded
if (angle > 70.0 && angle < 359.99) {
setSecStrucType(i, SecStrucType.bend);
state.setBend(true);
}
}
}
private void calculateDihedralAngles() {
// dihedral angles
// phi: C-N-CA-C
// psi: N-CA-C-N
// Chi1: N-CA-CB-CG, N-CA-CB-OG(SER),N-CA-CB-OG1(Thr),
// N-CA-CB-CG1(ILE/VAL), N-CA-CB-SG(CYS)
// Omega: CA-C-N-CA
for (int i=0 ; i < groups.length-1 ; i++){
SecStrucGroup a = groups[i];
SecStrucGroup b = groups[i+1];
Atom a_N = a.getN();
Atom a_CA = a.getCA();
Atom a_C = a.getC();
Atom b_N = b.getN();
Atom b_CA = b.getCA();
Atom b_C = b.getC();
double phi = Calc.torsionAngle(a_C,b_N,b_CA,b_C);
double psi = Calc.torsionAngle(a_N,a_CA,a_C,b_N);
double omega = Calc.torsionAngle(a_CA,a_C,b_N,b_CA);
SecStrucState state1 = (SecStrucState)
a.getProperty(Group.SEC_STRUC);
SecStrucState state2 = (SecStrucState)
b.getProperty(Group.SEC_STRUC);
state2.setPhi(phi);
state1.setPsi(psi);
state1.setOmega(omega);
}
}
@Override
public String toString() {
return printDSSP();
}
/**
* Generate a DSSP file format ouput String of this SS prediction.
* @return String in DSSP output file format
*/
public String printDSSP() {
StringBuffer buf = new StringBuffer();
String nl = System.getProperty("line.separator");
//Header Line
buf.append("==== Secondary Structure Definition by BioJava"
+ " DSSP implementation, Version October 2015 ===="+nl);
//First line with column definition
buf.append(" # RESIDUE AA STRUCTURE BP1 BP2 ACC "
+ "N-H-->O O-->H-N N-H-->O O-->H-N "
+ "TCO KAPPA ALPHA PHI PSI "
+ "X-CA Y-CA Z-CA ");
for (int i =0 ; i < groups.length ;i++){
buf.append(nl);
SecStrucState ss = getSecStrucState(i);
buf.append(ss.printDSSPline(i));
}
return buf.toString();
}
/**
* Generate a summary of this SS prediction with information about
* the three types of helix turns in different row sequences.
*
* This is similar to the summary output of Jmol, and useful to visualize
* the helix patterns.
*
* @return String helix summary
*/
public String printHelixSummary() {
StringBuffer g = new StringBuffer(); //3-10 helix
StringBuffer h = new StringBuffer(); //alpha helix
StringBuffer i = new StringBuffer(); //pi-helix
StringBuffer ss = new StringBuffer(); //SS summary
StringBuffer aa = new StringBuffer(); //AA one-letter
String nl = System.getProperty("line.separator");
g.append( "3 turn: ");
h.append( "4 turn: ");
i.append( "5 turn: ");
ss.append( "SS: ");
aa.append( "AA: ");
for (int k = 0; k < groups.length; k++){
SecStrucState state = getSecStrucState(k);
g.append(state.getTurn()[0]);
h.append(state.getTurn()[1]);
i.append(state.getTurn()[2]);
ss.append(state.getType());
aa.append(StructureTools.get1LetterCode(groups[k].getPDBName()));
}
return g.toString()+nl+h.toString()+nl+
i.toString()+nl+ss.toString()+nl+aa.toString();
}
/**
* Generate a FASTA sequence with the SS annotation letters in the
* aminoacid sequence order.
* @return String in FASTA sequence format
*/
public String printFASTA() {
StringBuffer buf = new StringBuffer();
String nl = System.getProperty("line.separator");
buf.append(">"+groups[0].getChain().getStructure().getIdentifier()+nl);
for (int g = 0; g < groups.length; g++){
buf.append(getSecStrucState(g).getType());
}
return buf.toString();
}
@Override
public int hashCode() {
final int prime = 31;
int result = 1;
result = prime * result + Arrays.hashCode(atoms);
return result;
}
@Override
public boolean equals(Object o){
if (!(o instanceof SecStrucCalc)) return false;
else {
SecStrucCalc ss = (SecStrucCalc) o;
if (groups.length != ss.groups.length) return false;
for (int g=0; g groupList = new ArrayList<>();
//
for ( Chain c : s.getChains(modelId)){
for (Group g : c.getAtomGroups()){
//We can also calc secstruc if it is a modified amino acid
if ( g.hasAminoAtoms()) {
SecStrucGroup sg = new SecStrucGroup();
sg.setResidueNumber(g.getResidueNumber());
sg.setPDBFlag(true);
sg.setPDBName(g.getPDBName());
sg.setChain(g.getChain());
Atom N = g.getAtom(StructureTools.N_ATOM_NAME);
Atom CA = g.getAtom(StructureTools.CA_ATOM_NAME);
Atom C = g.getAtom(StructureTools.C_ATOM_NAME);
Atom O = g.getAtom(StructureTools.O_ATOM_NAME);
if ( N == null || CA == null || C == null || O == null)
continue;
sg.setN((Atom) N.clone());
sg.setCA((Atom) CA.clone());
sg.setC((Atom) C.clone());
sg.setO((Atom) O.clone());
sg.setOriginal(g);
SecStrucState state = new SecStrucState(sg,
SecStrucInfo.BIOJAVA_ASSIGNMENT,
SecStrucType.coil);
sg.setProperty(Group.SEC_STRUC, state);
groupList.add(sg);
}
}
}
return groupList.toArray(new SecStrucGroup[groupList.size()]);
}
/**
* Calculate the coordinates of the H atoms. They are usually
* missing in the PDB files as only few experimental methods allow
* to resolve their location.
*/
private void calculateHAtoms() throws StructureException {
for ( int i = 0 ; i < groups.length-1 ; i++) {
SecStrucGroup a = groups[i];
SecStrucGroup b = groups[i+1];
if ( !b.hasAtom("H") ) {
//Atom H = calc_H(a.getC(), b.getN(), b.getCA());
Atom H = calcSimple_H(a.getC(), a.getO(), b.getN());
b.setH(H);
}
}
}
/**
* Calculate the HBonds between different groups.
* see Creighton page 147 f
* Modified to use only the contact map
*/
private void calculateHBonds() {
/**
* More efficient method for calculating C-Alpha pairs
*/
if (groups.length < 5) return;
Iterator otu = contactSet.iterator();
while(otu.hasNext()){
AtomContact ac = otu.next();
Pair pair = ac.getPair();
Group g1 = pair.getFirst().getGroup();
Group g2 = pair.getSecond().getGroup();
// Now I need to get the index of the Group in the list groups
int i = indResMap.get(g1.getResidueNumber());
int j = indResMap.get(g2.getResidueNumber());
// Now check this
checkAddHBond(i,j);
//"backwards" hbonds are not allowed
if (j!=(i+1)) checkAddHBond(j,i);
}
}
private void checkAddHBond(int i, int j){
SecStrucGroup one = groups[i];
if ("PRO".equals(one.getPDBName())){
logger.debug("Ignore: PRO {}", one.getResidueNumber());
return;
}
if (!one.hasAtom("H")) {
logger.debug("Residue {} has no H",one.getResidueNumber());
return;
}
SecStrucGroup two = groups[j];
double energy = 0;
try {
energy = calculateHBondEnergy(one,two);
} catch (Exception e){
logger.warn("Energy calculation failed", e);
return;
}
logger.debug("Energy between positions ({},{}): ",i,j,energy);
trackHBondEnergy(i,j,energy);
}
/**
* Calculate HBond energy of two groups in cal/mol
* see Creighton page 147 f
*
* Jeffrey, George A., An introduction to hydrogen bonding,
* Oxford University Press, 1997.
* categorizes hbonds with donor-acceptor distances of
* 2.2-2.5 å as "strong, mostly covalent",
* 2.5-3.2 å as "moderate, mostly electrostatic",
* 3.2-4.0 å as "weak, electrostatic".
* Energies are given as 40-14, 15-4, and <4 kcal/mol respectively.
*/
private static double calculateHBondEnergy(SecStrucGroup one,
SecStrucGroup two) {
Atom N = one.getN();
Atom H = one.getH();
Atom O = two.getO();
Atom C = two.getC();
double dno = Calc.getDistance(O,N);
double dhc = Calc.getDistance(C,H);
double dho = Calc.getDistance(O,H);
double dnc = Calc.getDistance(C,N);
logger.debug(" cccc: {} {} {} {} O ({})..N ({}):{} | ho:{} - hc:{} + nc:{} - no:{}",
one.getResidueNumber(),one.getPDBName(),two.getResidueNumber(),two.getPDBName(),
O.getPDBserial(),N.getPDBserial(),dno,dho,dhc,dnc,dno);
//there seems to be a contact!
if ( (dno < MINDIST) || (dhc < MINDIST) ||
(dnc < MINDIST) || (dno < MINDIST)) {
return HBONDLOWENERGY;
}
double e1 = Q / dho - Q / dhc;
double e2 = Q / dnc - Q / dno;
double energy = e1 + e2;
logger.debug(" N ({}) O({}): {} : {} ",N.getPDBserial(),O.getPDBserial(),(float) dno,energy);
//Avoid too strong energy
if (energy > HBONDLOWENERGY) return energy;
return HBONDLOWENERGY ;
}
/**
* Store Hbonds in the Groups.
* DSSP allows two HBonds per aminoacids to allow bifurcated bonds.
*/
private void trackHBondEnergy(int i, int j, double energy) {
if ("PRO".equals(groups[i].getPDBName())) {
logger.debug("Ignore: PRO {}",groups[i].getResidueNumber());
return;
}
SecStrucState stateOne = getSecStrucState(i);
SecStrucState stateTwo = getSecStrucState(j);
double acc1e = stateOne.getAccept1().getEnergy();
double acc2e = stateOne.getAccept2().getEnergy();
double don1e = stateTwo.getDonor1().getEnergy();
double don2e = stateTwo.getDonor2().getEnergy();
//Acceptor: N-H-->O
if (energy < acc1e) {
logger.debug("{} < {}",energy,acc1e);
stateOne.setAccept2(stateOne.getAccept1());
HBond bond = new HBond();
bond.setEnergy(energy);
bond.setPartner(j);
stateOne.setAccept1(bond);
} else if ( energy < acc2e ) {
logger.debug("{} < {}",energy,acc2e);
HBond bond = new HBond();
bond.setEnergy(energy);
bond.setPartner(j);
stateOne.setAccept2(bond);
}
//The other side of the bond: donor O-->N-H
if (energy < don1e) {
logger.debug("{} < {}",energy,don1e);
stateTwo.setDonor2(stateTwo.getDonor1());
HBond bond = new HBond();
bond.setEnergy(energy);
bond.setPartner(i);
stateTwo.setDonor1(bond);
} else if ( energy < don2e ) {
logger.debug("{} < {}",energy,don2e);
HBond bond = new HBond();
bond.setEnergy(energy);
bond.setPartner(i);
stateTwo.setDonor2(bond);
}
}
/**
* Detect helical turn patterns.
*/
private void calculateTurns(){
for (int i = 0 ; i< groups.length; i++){
for (int turn = 3; turn <= 5; turn++) {
if (i+turn >= groups.length) continue;
//Check for H bond from NH(i+n) to CO(i)
if (isBonded(i, i+turn)) {
logger.debug("Turn at ({},{}) turn {}",i,(i+turn),turn);
getSecStrucState(i).setTurn('>', turn);
getSecStrucState(i+turn).setTurn('<', turn);
//Bracketed residues get the helix number
for (int j=i+1; j' || turn[idx] == 'X') {
//Mark following n residues as turn
for (int k=1; k
* Note that the orignal DSSP implementation does not assign
* helix type to residue (i-1) and residue (i+n+1), although
* they contain a helix turn. As they state in the original paper,
* "the helices are one residue shorter than they would be according
* to rule 6.3 of IUPAC-IUB".
*
* @param n
* @param type
*/
private void checkSetHelix(int n, SecStrucType type){
int idx = n - 3;
logger.debug("Set helix {} {} {}", type, n, idx);
for (int i = 1; i < groups.length-n; i++) {
SecStrucState state = getSecStrucState(i);
SecStrucState previousState = getSecStrucState(i-1);
//Check that no other helix was assgined to this range
if (state.getType().compareTo(type) < 0) continue;
if (getSecStrucState(i+1).getType().compareTo(type) < 0) continue;
char turn = state.getTurn()[idx];
char pturn = previousState.getTurn()[idx];
//Two consecutive n-turns present to define a n-helix
if ((turn=='>' || turn=='X') && (pturn=='>' || pturn=='X')) {
//Mark following n residues as turn
for (int k=0; k