/* * 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.bio.structure.align; import java.io.FileOutputStream; import java.io.IOException; 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 org.biojava.bio.structure.Atom; import org.biojava.bio.structure.AtomImpl; import org.biojava.bio.structure.Calc; import org.biojava.bio.structure.SVDSuperimposer; import org.biojava.bio.structure.Structure; import org.biojava.bio.structure.StructureException; import org.biojava.bio.structure.StructureTools; import org.biojava.bio.structure.align.helper.AlignTools; import org.biojava.bio.structure.align.helper.JointFragments; import org.biojava.bio.structure.align.pairwise.AltAligComparator; import org.biojava.bio.structure.align.pairwise.AlternativeAlignment; import org.biojava.bio.structure.align.pairwise.FragmentJoiner; import org.biojava.bio.structure.align.pairwise.FragmentPair; import org.biojava.bio.structure.gui.BiojavaJmol; import org.biojava.bio.structure.io.PDBFileParser; import org.biojava.bio.structure.io.PDBFileReader; import org.biojava.bio.structure.jama.Matrix; /** * Perform a pairwise protein structure superimposition. * *
* 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) . *
* * * ** 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}(); try { structure1 = pdbpars.parsePDBFile(inStream1) ; structure2 = pdbpars.parsePDBFile(inStream2); } catch ({@link IOException} e) { e.printStackTrace(); return; } // calculate structure superimposition for two complete structures {@link StructurePairAligner} aligner = new {@link StructurePairAligner}(); try { // 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) {@link BiojavaJmol} jmol = new {@link 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;"); } catch ({@link StructureException} e){ e.printStackTrace(); } } ** * * * @author Andreas Prlic * @author Peter Lackner * @since 1.4 * @version %I% %G% */ public class StructurePairAligner { AlternativeAlignment[] alts; Matrix distanceMatrix; StrucAligParameters params; FragmentPair[] fragPairs; boolean debug = false; public StructurePairAligner() { super(); params = StrucAligParameters.getDefaultParameters(); reset(); alts = new AlternativeAlignment[0]; distanceMatrix = new Matrix(0,0); } /** example usage of this class * * @param args */ public static void main(String[] args){ try { // 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 System.out.println("aligning " + pdb1 + " vs. " + 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 (int i=0 ; i< aligs.length; i ++){ AlternativeAlignment aa = aligs[i]; System.out.println(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); System.out.println("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 (! BiojavaJmol.jmolInClassPath()){ System.err.println("Could not find Jmol in classpath, please install first!"); return; } AlternativeAlignment aa1 =aligs[0]; // first get an artificial structure for the alignment Structure artificial = aa1.getAlignedStructure(s1, s2); // 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[] pdbs1 = aa1.getPDBresnum1(); for (String res : pdbs1 ){ jmol.evalString("select " + res + "/1 ; backbone 0.6; color white;"); } jmol.evalString("select */2; color blue; model 2;"); String[] pdbs2 = aa1.getPDBresnum2(); for (String res :pdbs2 ){ jmol.evalString("select " + res + "/2 ; backbone 0.6; color yellow;"); } // now show both models again. jmol.evalString("model 0;"); jmol.evalString("echo done."); } } catch (Exception e){ e.printStackTrace(); } } 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 */ public boolean isDebug() { return debug; } /** set the debug flag * * @param debug flag */ public void setDebug(boolean debug) { this.debug = 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); align(ca1,ca2,params); } 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(); // step 1 get all Diagonals of length X that are similar between both structures if ( debug ) { System.out.println(" length atoms1:" + ca1.length); System.out.println(" length atoms2:" + ca2.length); System.out.println("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) || ca2.length < (fragmentLength + 1)) { throw new StructureException("structure too short, can not align"); } int rows = ca1.length - fragmentLength + 1; int cols = ca2.length - fragmentLength + 1; //System.out.println("rows " + rows + " " + cols + // " ca1 l " + ca1.length + " ca2 l " + ca2.length); 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}}; //Matrix unitv = new Matrix(utmp); Atom unitvector = new AtomImpl(); unitvector.setCoords(utmp[0]); List