* 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