org.biojava.nbio.structure.align.pairwise.AlternativeAlignment Maven / Gradle / Ivy
/*
* 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 Jan 28, 2006
*
*/
package org.biojava.nbio.structure.align.pairwise;
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.StrucAligParameters;
import org.biojava.nbio.structure.align.helper.AligMatEl;
import org.biojava.nbio.structure.align.helper.IndexPair;
import org.biojava.nbio.structure.align.helper.JointFragments;
import org.biojava.nbio.structure.geometry.Matrices;
import org.biojava.nbio.structure.geometry.SuperPositions;
import org.biojava.nbio.structure.jama.Matrix;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.io.Serializable;
import java.text.DecimalFormat;
import java.util.List;
import javax.vecmath.Matrix4d;
/**
* Implements a class which handles one possible (alternative) solution.
*
* Alternative alignments arise from different seed
* alignments or seed FPairs. The AltAlg class contains methods
* for refinement (Dynamic Programming based) and filtering
* (i.e. removing probably wrongly matched APairs). In the refinement
* phase, different seed alignments can converge to the same solution.
*
* @author Andreas Prlic,
* @author Peter Lackner (original Python and C code)
* @since 3:04:26 PM
* @version %I% %G%
*/
public class AlternativeAlignment implements Serializable{
private static final long serialVersionUID = -6226717654562221241L;
private int[] idx1;
private int[] idx2;
private String[] pdbresnum1;
private String[] pdbresnum2;
//short[] alig1;
//short[] alig2;
private int nfrags;
private Atom center;
private Matrix rot;
private Atom tr;
// the scores...
private int gaps0;
private int eqr0;
private int rms0;
private int joined;
private int percId;
private int cluster;
private float score;
private IndexPair[] aligpath;
private int fromia;
private Matrix currentRotMatrix;
private Atom currentTranMatrix;
private double rms;
private Matrix distanceMatrix;
public static final Logger logger = LoggerFactory.getLogger(AlternativeAlignment.class);
public AlternativeAlignment() {
super();
nfrags = 0;
aligpath = new IndexPair[0];
//alig1 = new short[0];
//alig2 = new short[0];
idx1 = new int[0];
idx2 = new int[0];
center = new AtomImpl();
rot = null;
tr = new AtomImpl();
eqr0 = -99;
rms0 = 99;
joined = 0;
gaps0 = -99;
fromia = 0;
currentRotMatrix = new Matrix(0,0);
currentTranMatrix = new AtomImpl();
distanceMatrix = new Matrix(0,0);
}
/** print the idx positions of this alignment
*
* @return a String representation
*/
@Override
public String toString(){
DecimalFormat d2 = new DecimalFormat();
// the result can be localized. To change this and enforce UK local do...
//(DecimalFormat)NumberFormat.getInstance(java.util.Locale.UK);
d2.setMaximumIntegerDigits(3);
d2.setMinimumFractionDigits(2);
d2.setMaximumFractionDigits(2);
StringBuffer s = new StringBuffer();
s.append("#" + getAltAligNumber() +
" cluster:" + cluster +
" eqr:" + getEqr() +
" rmsd:" + d2.format(getRmsd()) +
" %id:" + getPercId() +
" gaps:" + getGaps() +
" score:" + d2.format(score) );
return s.toString();
}
/** get the number of the cluster this alignment belongs to
*
* @return an int giving the number of the cluster
*/
public int getCluster() {
return cluster;
}
/** set the number of the cluster this alignment belongs to.
* All alignments in a cluster are quite similar.
*
* @param cluster the number of the cluster
*/
public void setCluster(int cluster) {
this.cluster = cluster;
}
public double getRmsd() {
return rms;
}
/** the rms in the structurally equivalent regions
*
* @param rms
*/
public void setRms(double rms) {
this.rms = rms;
}
/** the alignment score
*
* @return the score of this alignment
*/
public float getScore() {
return score;
}
public void setScore(float score) {
this.score = score;
}
/** return the number of gaps in this alignment
*
* @return the number of Gaps
*/
public int getGaps(){
return gaps0;
}
/** returns the number of euqivalent residues in this alignment
*
* @return the number of equivalent residues
*/
public int getEqr(){
return eqr0 ;
}
/** the positions of the structure equivalent positions in atom set 1
*
* @return the array of the positions
*/
public int[] getIdx1(){
return idx1;
}
/** the positions of the structure equivalent atoms in atom set 2
*
* @return the array of the positions
*/
public int[] getIdx2(){
return idx2;
}
public int getPercId() {
return percId;
}
public void setPercId(int percId) {
this.percId = percId;
}
/** Set apairs according to a seed position.
*
* @param l
* @param i
* @param j
*/
public void apairs_from_seed(int l,int i, int j){
aligpath = new IndexPair[l];
idx1 = new int[l];
idx2 = new int[l];
for (int x = 0 ; x < l ; x++) {
idx1[x]=i+x;
idx2[x]=j+x;
aligpath[x] = new IndexPair((short)(i+x),(short)(j+x));
}
}
/** Set apairs according to a list of (i,j) tuples.
*
* @param jf a JoingFragment
*/
public void apairs_from_idxlst(JointFragments jf) {
List il = jf.getIdxlist();
//System.out.println("Alt Alig apairs_from_idxlst");
aligpath = new IndexPair[il.size()];
idx1 = new int[il.size()];
idx2 = new int[il.size()];
for (int i =0 ; i < il.size();i++) {
int[] p = il.get(i);
//System.out.print(" idx1 " + p[0] + " idx2 " + p[1]);
idx1[i] = p[0];
idx2[i] = p[1];
aligpath[i] = new IndexPair((short)p[0],(short)p[1]);
}
eqr0 = idx1.length;
//System.out.println("eqr " + eqr0);
gaps0 = count_gaps(idx1,idx2);
}
/** returns the sequential number of this alternative alignment
*
* @return the sequential number of this alternative alignment
*/
public int getAltAligNumber() {
return fromia;
}
public void setAltAligNumber(int fromia) {
this.fromia = fromia;
}
/** rotate and shift atoms with currentRotMatrix and current Tranmatrix
*
* @param ca
*/
private void rotateShiftAtoms(Atom[] ca){
for (int i = 0 ; i < ca.length; i++){
Atom c = ca[i];
Calc.rotate(c,currentRotMatrix);
Calc.shift(c,currentTranMatrix);
//System.out.println("after " + c);
ca[i] = c;
}
//System.out.println("after " + ca[0]);
}
public void finish(StrucAligParameters params,Atom[]ca1,Atom[]ca2) throws StructureException{
Atom[] ca3 = new Atom[ca2.length];
for ( int i = 0 ; i < ca2.length;i++){
ca3[i] = (Atom) ca2[i].clone();
}
// do the inital superpos...
super_pos_alig(ca1,ca3,idx1,idx2,true);
rotateShiftAtoms(ca3);
calcScores(ca1,ca2);
logger.debug("eqr " + eqr0 + " " + gaps0 + " " +idx1[0] + " " +idx1[1]);
getPdbRegions(ca1,ca2);
}
public static Matrix getDistanceMatrix(Atom[] ca1, Atom[]ca2){
int r = ca1.length;
int c = ca2.length;
Matrix out = new Matrix(r,c);
for (int i=0; i co) {
e.setValue(0);
} else {
double s = 2 * co / ( 1 + ( d/co) * (d/co)) - co;
e.setValue((int) Math.round(Gotoh.ALIGFACTOR * s));
}
aligmat[i+1][j+1] = e;
}
}
return al;
}
/*private Alignable getAlignableFromSimMatrix (Matrix sim,StrucAligParameters params){
//System.out.println("align_NPE");
float gapOpen = params.getGapOpen();
float gapExtension = params.getGapExtension();
int rows = sim.getRowDimension();
int cols = sim.getColumnDimension();
Alignable al = new StrCompAlignment(rows,cols);
al.setGapExtCol(gapExtension);
al.setGapExtRow(gapExtension);
al.setGapOpenCol(gapOpen);
al.setGapOpenRow(gapOpen);
//System.out.println("size of aligmat: " + rows+1 + " " + cols+1);
//AligMatEl[][] aligmat = new AligMatEl[rows+1][cols+1];
AligMatEl[][] aligmat = al.getAligMat();
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
int e=0;
//if ( ( i < rows) &&
// ( j < cols)) {
//TODO: the ALIGFACTOR calc should be hidden in Gotoh!!
e = (int)Math.round(Gotoh.ALIGFACTOR * sim.get(i,j));
//}
//System.out.println(e);
//AligMatEl am = new AligMatEl();
aligmat[i+1][j+1].setValue(e);
//.setValue(e);
//am.setTrack(new IndexPair((short)-99,(short)-99));
//igmat[i+1][j+1] = am;
}
}
//al.setAligMat(aligmat);
return al;
}*/
/** Refinement procedure based on superposition and dynamic programming.
* Performs an iterative refinement. Several methods apply such a procedure,
* e.g. CE or ProSup. Here we additionally test for circular permutation,
* which are in the same frame of superposition as the optimal alignment.
* This feature may be switched off by setting permsize to -1.
*
* @param params the parameters
* @param ca1 atoms of structure 1
* @param ca2 atoms of structure 2
* @throws StructureException
*/
public void refine(StrucAligParameters params,Atom[]ca1,Atom[]ca2) throws StructureException{
// System.out.println("refine Alternative Alignment #"+ getAltAligNumber()+" l1:" + ca1.length + " l2:" + ca2.length);
// for ( int i= 0 ; i < idx1.length;i++){
// System.out.println(idx1[i] + " " + idx2[i]);
// }
// avoid any operations on the original Atoms ...
Atom[] ca3 = new Atom[ca2.length];
for ( int i = 0 ; i < ca2.length;i++){
ca3[i] = (Atom) ca2[i].clone();
}
//Atom[] ca3 = (Atom[]) ca2.clone();
// do the inital superpos...
super_pos_alig(ca1,ca3,idx1,idx2,false);
//JmolDisplay jmol1 = new JmolDisplay(ca1,ca2,"after first alig " );
//jmol1.selectEquiv(idx1,idx2);
//int[] jmoltmp1 = idx1;
//int[] jmoltmp2 = idx2;
int lenalt =idx1.length;
int lenneu = aligpath.length;
//Matrix rmsalt = currentRotMatrix;
int ml = Math.max(ca1.length, ca3.length);
idx1 = new int[ml];
idx2 = new int[ml];
//int m1l = ca1.length;
//int m2l = ca3.length;
//Matrix sim = new Matrix(ca1.length,ca3.length,0);
int maxiter = params.getMaxIter();
for (int iter = 0 ; iter< maxiter; iter++){
float subscore = 0.0f;
rotateShiftAtoms(ca3);
// do the dynamic alignment...
Alignable ali = getInitalStrCompAlignment(ca1,ca3, params);
new Gotoh(ali);
this.score = ali.getScore();
//System.out.println("score " + score);
IndexPair[] path = ali.getPath();
int pathsize = ali.getPathSize();
// now for a superimposable permutation. First we need to check the size and
// position of the quadrant left out by the optimal path
int firsta = path[0].getRow();
int firstb = path[0].getCol();
int lasta = path[pathsize-1].getRow();
int lastb = path[pathsize-1].getCol();
int quada_beg, quada_end, quadb_beg,quadb_end;
if ( firsta == 0){
quada_beg = lasta +1;
quada_end = ca1.length -1;
quadb_beg = 0;
quadb_end = firstb-1;
} else {
quada_beg = 0;
quada_end = firsta-1;
quadb_beg = lastb+1;
quadb_end = ca3.length -1;
}
// if the unaligned quadrant is larger than permsize
int permsize = params.getPermutationSize();
if ( permsize > 0 &&
(quada_end - quada_beg >= permsize) &&
( quadb_end - quadb_beg >= permsize)) {
AligMatEl[][] aligmat = ali.getAligMat();
Matrix submat = new Matrix(quada_end-quada_beg, quadb_end - quadb_beg) ;
//then we copy the score values of the quadrant into a new matrix
//System.out.println(quada_beg + " " + quada_end + " " +quadb_beg+" " + quadb_end);
//System.out.println(sim.getRowDimension() + " " + sim.getColumnDimension());
Atom[] tmp1 = new Atom[quada_end - quada_beg ];
Atom[] tmp2 = new Atom[quadb_end - quadb_beg ];
for ( int s = quada_beg; s < quada_end; s++){
//System.out.println(s + " " +( quada_end-s));
tmp1[quada_end - s -1] = ca1[s];
for ( int t = quadb_beg; t < quadb_end; t++){
if ( s == quada_beg )
tmp2[quadb_end - t -1 ] = ca3[t];
//System.out.println(s+" " +t);
//double val = sim.get(s,t);
double val = aligmat[s][t].getValue();
submat.set(s-quada_beg,t-quadb_beg,val);
}
}
// and perform a dp alignment again. (Note, that we fix the superposition).
//Alignable subali = getAlignableFromSimMatrix(submat,params);
Alignable subali = getInitalStrCompAlignment(tmp1, tmp2, params);
subscore = subali.getScore();
//System.out.println("join score" + score + " " + subscore);
this.score = score+subscore;
IndexPair[] subpath = subali.getPath();
int subpathsize = subali.getPathSize();
for ( int p=0;p permsize){
IndexPair[] wholepath = new IndexPair[pathsize+subpathsize];
for ( int t=0; t < pathsize; t++){
wholepath[t] = path[t];
}
for ( int t=0 ; t < subpathsize; t++){
wholepath[t+pathsize] = subpath[t];
}
pathsize += subpathsize;
path = wholepath;
ali.setPath(path);
ali.setPathSize(pathsize);
}
}
int j =0 ;
//System.out.println("pathsize,path " + pathsize + " " + path.length);
for (int i=0; i < pathsize; i++){
int x = path[i].getRow();
int y = path[i].getCol();
//System.out.println(x+ " " + y);
double d = Calc.getDistance(ca1[x], ca3[y]);
//double d = dismat.get(x,y);
// now we apply the evaluation distance cutoff
if ( d < params.getEvalCutoff()){
//System.out.println("j:"+j+ " " + x+" " + y + " " + d );
idx1[j] = x;
idx2[j] = y;
j+=1;
}
}
lenneu = j;
//System.out.println("lenalt,neu:" + lenalt + " " + lenneu);
int[] tmpidx1 = new int[j];
int[] tmpidx2 = new int[j];
//String idx1str ="idx1: ";
//String idx2str ="idx2: ";
for (int i = 0 ; i 0)
// System.out.println("do new superimpos " + idx1.length + " " + idx2.length + " " + idx1[0] + " " + idx2[0]);
super_pos_alig(ca1,ca3,tmpidx1,tmpidx2,false);
this.aligpath = path;
if (lenneu == lenalt)
break;
}
idx1 = (int[]) FragmentJoiner.resizeArray(idx1,lenneu);
idx2 = (int[]) FragmentJoiner.resizeArray(idx2,lenneu);
// new ... get rms...
// now use the original atoms to get the rotmat relative to the original structure...
super_pos_alig(ca1,ca2,idx1,idx2,true);
eqr0 = idx1.length;
gaps0 = count_gaps(idx1,idx2);
getPdbRegions(ca1,ca2);
//System.out.println("eqr " + eqr0 + " aligpath len:"+aligpath.length+ " gaps:" + gaps0 + " score " + score);
}
private void getPdbRegions(Atom[] ca1, Atom[] ca2){
pdbresnum1 = new String[idx1.length];
pdbresnum2 = new String[idx2.length];
for (int i =0 ; i < idx1.length;i++){
Atom a1 = ca1[idx1[i]];
Atom a2 = ca2[idx2[i]];
Group p1 = a1.getGroup();
Group p2 = a2.getGroup();
Chain c1 = p1.getChain();
Chain c2 = p2.getChain();
String cid1 = c1.getId();
String cid2 = c2.getId();
String pdb1 = p1.getResidueNumber().toString();
String pdb2 = p2.getResidueNumber().toString();
if ( ! cid1.equals(" "))
pdb1 += ":" + cid1;
if ( ! cid2.equals(" "))
pdb2 += ":" + cid2;
pdbresnum1[i] = pdb1;
pdbresnum2[i] = pdb2;
}
}
public String[] getPDBresnum1() {
return pdbresnum1;
}
public void setPDBresnum1(String[] pdbresnum1) {
this.pdbresnum1 = pdbresnum1;
}
public String[] getPDBresnum2() {
return pdbresnum2;
}
public void setPDBresnum2(String[] pdbresnum2) {
this.pdbresnum2 = pdbresnum2;
}
/**
* Count the number of gaps in an alignment represented by idx1,idx2.
*
* @param i1
* @param i2
* @return the number of gaps in this alignment
*/
private int count_gaps(int[] i1, int[] i2){
int i0 = i1[0];
int j0 = i2[0];
int gaps = 0;
for (int i =1 ; i © 2015 - 2025 Weber Informatics LLC | Privacy Policy