/*
* 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/
*
* Created on May 21, 2006
*
*/
package org.biojava.nbio.structure.align;
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.ce.GuiWrapper;
import org.biojava.nbio.structure.align.helper.AlignTools;
import org.biojava.nbio.structure.align.helper.JointFragments;
import org.biojava.nbio.structure.align.pairwise.*;
import org.biojava.nbio.structure.geometry.Matrices;
import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.io.PDBFileParser;
import org.biojava.nbio.structure.io.PDBFileReader;
import org.biojava.nbio.structure.jama.Matrix;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.io.FileOutputStream;
import java.io.InputStream;
import java.io.PrintStream;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.List;
import javax.vecmath.Matrix4d;
/**
* Perform a pairwise protein structure superimposition.
*
* <p>
* The algorithm is a distance matrix based, rigid body protein structure
* superimposition. It is based on a variation of the PSC++ algorithm provided
* by Peter Lackner (Peter.Lackner@sbg.ac.at, personal communication) .
* </p>
*
*
*
* <h2>Example</h2>
*
* <pre>
* public void run(){
*
* // first load two example structures
* {@link InputStream} inStream1 = this.getClass().getResourceAsStream("/files/5pti.pdb");
* {@link InputStream} inStream2 = this.getClass().getResourceAsStream("/files/1tap.pdb");
*
* {@link Structure} structure1 = null;
* {@link Structure} structure2 = null;
*
* {@link PDBFileParser} pdbpars = new {@link PDBFileParser}();
* structure1 = pdbpars.parsePDBFile(inStream1) ;
* structure2 = pdbpars.parsePDBFile(inStream2);
*
*
* // calculate structure superimposition for two complete structures
* {@link StructurePairAligner} aligner = new {@link StructurePairAligner}();
*
*
* // align the full 2 structures with default parameters.
* // see StructurePairAligner for more options and how to align
* // any set of Atoms
* aligner.align(structure1,structure2);
*
* {@link AlternativeAlignment}[] aligs = aligner.getAlignments();
* {@link AlternativeAlignment} a = aligs[0];
* System.out.println(a);
*
* //display the alignment in Jmol
*
* // first get an artificial structure for the alignment
* {@link Structure} artificial = a.getAlignedStructure(structure1, structure2);
*
*
* // and then send it to Jmol (only will work if Jmol is in the Classpath)
*
* BiojavaJmol jmol = new BiojavaJmol();
* jmol.setTitle(artificial.getName());
* jmol.setStructure(artificial);
*
* // color the two structures
*
*
* jmol.evalString("select *; backbone 0.4; wireframe off; spacefill off; " +
* "select not protein and not solvent; spacefill on;");
* jmol.evalString("select *"+"/1 ; color red; model 1; ");
*
*
* // now color the equivalent residues ...
*
* String[] pdb1 = a.getPDBresnum1();
* for (String res : pdb1 ){
* jmol.evalString("select " + res + "/1 ; backbone 0.6; color white;");
* }
*
* jmol.evalString("select *"+"/2; color blue; model 2;");
* String[] pdb2 = a.getPDBresnum2();
* for (String res :pdb2 ){
* jmol.evalString("select " + res + "/2 ; backbone 0.6; color yellow;");
* }
*
*
* // now show both models again.
* jmol.evalString("model 0;");
*
* }
* </pre>
*
*
*
* @author Andreas Prlic
* @author Peter Lackner
* @since 1.4
* @version %I% %G%
*/
public class StructurePairAligner {
private final static Logger logger = LoggerFactory
.getLogger(StructurePairAligner.class);
AlternativeAlignment[] alts;
Matrix distanceMatrix;
StrucAligParameters params;
FragmentPair[] fragPairs;
List<AlignmentProgressListener> listeners = new ArrayList<AlignmentProgressListener>();
public StructurePairAligner() {
super();
params = StrucAligParameters.getDefaultParameters();
reset();
alts = new AlternativeAlignment[0];
distanceMatrix = new Matrix(0, 0);
}
public void addProgressListener(AlignmentProgressListener li) {
listeners.add(li);
}
public void clearListeners() {
listeners.clear();
}
/**
* example usage of this class
*
* @param args
*/
public static void main(String[] args) throws Exception {
// UPDATE THE FOLLOWING LINES TO MATCH YOUR SETUP
PDBFileReader pdbr = new PDBFileReader();
pdbr.setPath("/Users/andreas/WORK/PDB/");
// String pdb1 = "1crl";
// String pdb2 = "1ede";
String pdb1 = "1buz";
String pdb2 = "1ali";
String outputfile = "/tmp/alig_" + pdb1 + "_" + pdb2 + ".pdb";
// NO NEED TO DO CHANGE ANYTHING BELOW HERE...
StructurePairAligner sc = new StructurePairAligner();
// step1 : read molecules
logger.info("aligning {} vs. {}", pdb1, pdb2);
Structure s1 = pdbr.getStructureById(pdb1);
Structure s2 = pdbr.getStructureById(pdb2);
// step 2 : do the calculations
sc.align(s1, s2);
AlternativeAlignment[] aligs = sc.getAlignments();
// cluster similar results together
ClusterAltAligs.cluster(aligs);
// print the result:
// the AlternativeAlignment object gives also access to rotation
// matrices / shift vectors.
for (AlternativeAlignment aa : aligs) {
logger.info("Alternative Alignment: ", aa);
}
// convert AlternativeAlignemnt 1 to PDB file, so it can be opened with
// a viewer (e.g. Jmol, Rasmol)
if (aligs.length > 0) {
AlternativeAlignment aa1 = aligs[0];
String pdbstr = aa1.toPDB(s1, s2);
logger.info("writing alignment to {}", outputfile);
FileOutputStream out = new FileOutputStream(outputfile);
PrintStream p = new PrintStream(out);
p.println(pdbstr);
p.close();
out.close();
}
// display the alignment in Jmol
// only will work if Jmol is in the Classpath
if (aligs.length > 0) {
if (!GuiWrapper.isGuiModuleInstalled()) {
logger.error("Could not find structure-gui modules in classpath, please install first!");
}
}
}
private void reset() {
alts = new AlternativeAlignment[0];
distanceMatrix = new Matrix(0, 0);
fragPairs = new FragmentPair[0];
}
/**
* get the results of step 1 - the FragmentPairs used for seeding the
* alignment
*
* @return a FragmentPair[] array
*/
public FragmentPair[] getFragmentPairs() {
return fragPairs;
}
public void setFragmentPairs(FragmentPair[] fragPairs) {
this.fragPairs = fragPairs;
}
/**
* return the alternative alignments that can be found for the two
* structures
*
* @return AlternativeAlignment[] array
*/
public AlternativeAlignment[] getAlignments() {
return alts;
}
/**
* return the difference of distance matrix between the two structures
*
* @return a Matrix
*/
public Matrix getDistMat() {
return distanceMatrix;
}
/**
* get the parameters.
*
* @return the Parameters.
*/
public StrucAligParameters getParams() {
return params;
}
/**
* set the parameters to be used for the algorithm
*
* @param params
* the Parameter object
*/
public void setParams(StrucAligParameters params) {
this.params = params;
}
/**
* check if debug mode is set on
*
* @return debug flag
*/
@Deprecated
public boolean isDebug() {
return false;
}
/**
* set the debug flag
*
* @param debug
* flag
*/
@Deprecated
public void setDebug(boolean debug) {
}
/**
* Calculate the alignment between the two full structures with default
* parameters
*
* @param s1
* @param s2
* @throws StructureException
*/
public void align(Structure s1, Structure s2) throws StructureException {
align(s1, s2, params);
}
/**
* Calculate the alignment between the two full structures with user
* provided parameters
*
* @param s1
* @param s2
* @param params
* @throws StructureException
*/
public void align(Structure s1, Structure s2, StrucAligParameters params)
throws StructureException {
// step 1 convert the structures to Atom Arrays
Atom[] ca1 = getAlignmentAtoms(s1);
Atom[] ca2 = getAlignmentAtoms(s2);
notifyStartingAlignment(s1.getName(), ca1, s2.getName(), ca2);
align(ca1, ca2, params);
}
/**
* Align two chains from the structures. Uses default parameters.
*
* @param s1
* @param chainId1
* @param s2
* @param chainId2
*/
public void align(Structure s1, String chainId1, Structure s2,
String chainId2) throws StructureException {
align(s1, chainId1, s2, chainId2, params);
}
/**
* Aligns two chains from the structures using user provided parameters.
*
* @param s1
* @param chainId1
* @param s2
* @param chainId2
* @param params
* @throws StructureException
*/
public void align(Structure s1, String chainId1, Structure s2,
String chainId2, StrucAligParameters params)
throws StructureException {
reset();
this.params = params;
Chain c1 = s1.getPolyChainByPDB(chainId1);
Chain c2 = s2.getPolyChainByPDB(chainId2);
Structure s3 = new StructureImpl();
s3.addChain(c1);
Structure s4 = new StructureImpl();
s4.addChain(c2);
Atom[] ca1 = getAlignmentAtoms(s3);
Atom[] ca2 = getAlignmentAtoms(s4);
notifyStartingAlignment(s1.getName(), ca1, s2.getName(), ca2);
align(ca1, ca2, params);
}
/**
* Returns the atoms that are being used for the alignment. (E.g. Calpha
* only, etc.)
*
* @param s
* @return an array of Atoms objects
*/
public Atom[] getAlignmentAtoms(Structure s) {
String[] atomNames = params.getUsedAtomNames();
return StructureTools.getAtomArray(s, atomNames);
}
/**
* calculate the protein structure superimposition, between two sets of
* atoms.
*
*
*
* @param ca1
* set of Atoms of structure 1
* @param ca2
* set of Atoms of structure 2
* @param params
* the parameters to use for the alignment
* @throws StructureException
*/
public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
throws StructureException {
reset();
this.params = params;
long timeStart = System.currentTimeMillis();
// step 1 get all Diagonals of length X that are similar between both
// structures
logger.debug(" length atoms1:" + ca1.length);
logger.debug(" length atoms2:" + ca2.length);
logger.debug("step 1 - get fragments with similar intramolecular distances ");
int k = params.getDiagonalDistance();
int k2 = params.getDiagonalDistance2();
int fragmentLength = params.getFragmentLength();
if (ca1.length < (fragmentLength + 1)) {
throw new StructureException("structure 1 too short (" + ca1.length
+ "), can not align");
}
if (ca2.length < (fragmentLength + 1)) {
throw new StructureException("structure 2 too short (" + ca2.length
+ "), can not align");
}
int rows = ca1.length - fragmentLength + 1;
int cols = ca2.length - fragmentLength + 1;
distanceMatrix = new Matrix(rows, cols, 0.0);
double[] dist1 = AlignTools.getDiagonalAtK(ca1, k);
double[] dist2 = AlignTools.getDiagonalAtK(ca2, k);
double[] dist3 = new double[0];
double[] dist4 = new double[0];
if (k2 > 0) {
dist3 = AlignTools.getDiagonalAtK(ca1, k2);
dist4 = AlignTools.getDiagonalAtK(ca2, k2);
}
double[][] utmp = new double[][] { { 0, 0, 1 } };
Atom unitvector = new AtomImpl();
unitvector.setCoords(utmp[0]);
List<FragmentPair> fragments = new ArrayList<FragmentPair>();
for (int i = 0; i < rows; i++) {
Atom[] catmp1 = AlignTools.getFragment(ca1, i, fragmentLength);
Atom center1 = AlignTools.getCenter(ca1, i, fragmentLength);
for (int j = 0; j < cols; j++) {
double rdd1 = AlignTools.rms_dk_diag(dist1, dist2, i, j,
fragmentLength, k);
double rdd2 = 0;
if (k2 > 0)
rdd2 = AlignTools.rms_dk_diag(dist3, dist4, i, j,
fragmentLength, k2);
double rdd = rdd1 + rdd2;
distanceMatrix.set(i, j, rdd);
if (rdd < params.getFragmentMiniDistance()) {
FragmentPair f = new FragmentPair(fragmentLength, i, j);
Atom[] catmp2 = AlignTools.getFragment(ca2, j,
fragmentLength);
Atom center2 = AlignTools.getCenter(ca2, j,
fragmentLength);
f.setCenter1(center1);
f.setCenter2(center2);
Matrix4d t = SuperPositions.superpose(
Calc.atomsToPoints(catmp1),
Calc.atomsToPoints(catmp2));
Matrix rotmat = Matrices.getRotationJAMA(t);
f.setRot(rotmat);
Atom aunitv = (Atom) unitvector.clone();
Calc.rotate(aunitv, rotmat);
f.setUnitv(aunitv);
boolean doNotAdd = false;
if (params.reduceInitialFragments()) {
doNotAdd = FragmentJoiner.reduceFragments(
fragments, f, distanceMatrix);
}
if (doNotAdd)
continue;
fragments.add(f);
}
}
}
notifyFragmentListeners(fragments);
FragmentPair[] fp = fragments
.toArray(new FragmentPair[fragments.size()]);
setFragmentPairs(fp);
logger.debug(" got # fragment pairs: {}", fp.length);
logger.debug("step 2 - join fragments");
// step 2 combine them to possible models
FragmentJoiner joiner = new FragmentJoiner();
JointFragments[] frags;
if (params.isJoinFast()) {
// apply the quick alignment procedure.
// less quality in alignments, better for DB searches...
frags = joiner.approach_ap3(ca1, ca2, fp, params);
joiner.extendFragments(ca1, ca2, frags, params);
} else if (params.isJoinPlo()) {
// this approach by StrComPy (peter lackner):
frags = joiner.frag_pairwise_compat(fp, params.getAngleDiff(),
params.getFragCompat(), params.getMaxrefine());
} else {
// my first implementation
frags = joiner.approach_ap3(ca1, ca2, fp, params);
}
notifyJointFragments(frags);
logger.debug(" number joint fragments: ", frags.length);
logger.debug("step 3 - refine alignments");
List<AlternativeAlignment> aas = new ArrayList<AlternativeAlignment>();
for (int i = 0; i < frags.length; i++) {
JointFragments f = frags[i];
AlternativeAlignment a = new AlternativeAlignment();
a.apairs_from_idxlst(f);
a.setAltAligNumber(i + 1);
a.setDistanceMatrix(distanceMatrix);
try {
if (params.getMaxIter() > 0) {
a.refine(params, ca1, ca2);
} else {
a.finish(params, ca1, ca2);
}
} catch (StructureException e) {
logger.error("Refinement of fragment {} failed", i, e);
}
a.calcScores(ca1, ca2);
aas.add(a);
}
// sort the alternative alignments
Comparator<AlternativeAlignment> comp = new AltAligComparator();
Collections.sort(aas, comp);
Collections.reverse(aas);
alts = aas.toArray(new AlternativeAlignment[aas.size()]);
// do final numbering of alternative solutions
int aanbr = 0;
for (AlternativeAlignment a : alts) {
aanbr++;
a.setAltAligNumber(aanbr);
}
logger.debug("total calculation time: {} ms.",
(System.currentTimeMillis() - timeStart));
}
private void notifyStartingAlignment(String name1, Atom[] ca1,
String name2, Atom[] ca2) {
for (AlignmentProgressListener li : listeners) {
li.startingAlignment(name1, ca1, name2, ca2);
}
}
private void notifyFragmentListeners(List<FragmentPair> fragments) {
for (AlignmentProgressListener li : listeners) {
li.calculatedFragmentPairs(fragments);
}
}
private void notifyJointFragments(JointFragments[] fragments) {
for (AlignmentProgressListener li : listeners) {
li.jointFragments(fragments);
}
}
}