/*
* 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.align.util;
import org.biojava.nbio.structure.*;
import org.biojava.nbio.structure.align.AFPTwister;
import org.biojava.nbio.structure.align.ce.CECalculator;
import org.biojava.nbio.structure.align.fatcat.FatCatFlexible;
import org.biojava.nbio.structure.align.fatcat.FatCatRigid;
import org.biojava.nbio.structure.align.model.AFPChain;
import org.biojava.nbio.structure.align.xml.AFPChainXMLParser;
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.IOException;
import java.io.Writer;
import java.util.*;
import java.util.Map.Entry;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import javax.vecmath.Matrix4d;
/**
* Methods for analyzing and manipulating AFPChains and for
* other pairwise alignment utilities. <p>
* Current methods: replace optimal alignment, create new AFPChain,
* format conversion, update superposition, etc.
*
* @author Spencer Bliven
* @author Aleix Lafita
*
*/
public class AlignmentTools {
private static final Logger logger = LoggerFactory.getLogger(AlignmentTools.class);
public static boolean debug = false;
/**
* Checks that the alignment given by afpChain is sequential. This means
* that the residue indices of both proteins increase monotonically as
* a function of the alignment position (ie both proteins are sorted).
*
* This will return false for circularly permuted alignments or other
* non-topological alignments. It will also return false for cases where
* the alignment itself is sequential but it is not stored in the afpChain
* in a sorted manner.
*
* Since algorithms which create non-sequential alignments split the
* alignment into multiple blocks, some computational time can be saved
* by only checking block boundaries for sequentiality. Setting
* <tt>checkWithinBlocks</tt> to <tt>true</tt> makes this function slower,
* but detects AFPChains with non-sequential blocks.
*
* Note that this method should give the same results as
* {@link AFPChain#isSequentialAlignment()}. However, the AFPChain version
* relies on the StructureAlignment algorithm correctly setting this
* parameter, which is sadly not always the case.
*
* @param afpChain An alignment
* @param checkWithinBlocks Indicates whether individual blocks should be
* checked for sequentiality
* @return True if the alignment is sequential.
*/
public static boolean isSequentialAlignment(AFPChain afpChain, boolean checkWithinBlocks) {
int[][][] optAln = afpChain.getOptAln();
int[] alnLen = afpChain.getOptLen();
int blocks = afpChain.getBlockNum();
if(blocks < 1) return true; //trivial case
if ( alnLen[0] < 1) return true;
// Check that blocks are sequential
if(checkWithinBlocks) {
for(int block = 0; block<blocks; block++) {
if(alnLen[block] < 1 ) continue; //skip empty blocks
int prevRes1 = optAln[block][0][0];
int prevRes2 = optAln[block][1][0];
for(int pos = 1; pos<alnLen[block]; pos++) {
int currRes1 = optAln[block][0][pos];
int currRes2 = optAln[block][1][pos];
if(currRes1 < prevRes1) {
return false;
}
if(currRes2 < prevRes2) {
return false;
}
prevRes1 = currRes1;
prevRes2 = currRes2;
}
}
}
// Check that blocks are sequential
int prevRes1 = optAln[0][0][alnLen[0]-1];
int prevRes2 = optAln[0][1][alnLen[0]-1];
for(int block = 1; block<blocks;block++) {
if(alnLen[block] < 1 ) continue; //skip empty blocks
if(optAln[block][0][0]<prevRes1) {
return false;
}
if(optAln[block][1][0]<prevRes2) {
return false;
}
prevRes1 = optAln[block][0][alnLen[block]-1];
prevRes2 = optAln[block][1][alnLen[block]-1];
}
return true;
}
/**
* Creates a Map specifying the alignment as a mapping between residue indices
* of protein 1 and residue indices of protein 2.
*
* <p>For example,<pre>
* 1234
* 5678</pre>
* becomes<pre>
* 1->5
* 2->6
* 3->7
* 4->8</pre>
*
* @param afpChain An alignment
* @return A mapping from aligned residues of protein 1 to their partners in protein 2.
* @throws StructureException If afpChain is not one-to-one
*/
public static Map<Integer, Integer> alignmentAsMap(AFPChain afpChain) throws StructureException {
Map<Integer,Integer> map = new HashMap<Integer,Integer>();
if( afpChain.getAlnLength() < 1 ) {
return map;
}
int[][][] optAln = afpChain.getOptAln();
int[] optLen = afpChain.getOptLen();
for(int block = 0; block < afpChain.getBlockNum(); block++) {
for(int pos = 0; pos < optLen[block]; pos++) {
int res1 = optAln[block][0][pos];
int res2 = optAln[block][1][pos];
if(map.containsKey(res1)) {
throw new StructureException(String.format("Residue %d aligned to both %d and %d.", res1,map.get(res1),res2));
}
map.put(res1,res2);
}
}
return map;
}
/**
* Applies an alignment k times. Eg if alignmentMap defines function f(x),
* this returns a function f^k(x)=f(f(...f(x)...)).
*
* @param <T>
* @param alignmentMap The input function, as a map (see {@link AlignmentTools#alignmentAsMap(AFPChain)})
* @param k The number of times to apply the alignment
* @return A new alignment. If the input function is not automorphic
* (one-to-one), then some inputs may map to null, indicating that the
* function is undefined for that input.
*/
public static <T> Map<T,T> applyAlignment(Map<T, T> alignmentMap, int k) {
return applyAlignment(alignmentMap, new IdentityMap<T>(), k);
}
/**
* Applies an alignment k times. Eg if alignmentMap defines function f(x),
* this returns a function f^k(x)=f(f(...f(x)...)).
*
* To allow for functions with different domains and codomains, the identity
* function allows converting back in a reasonable way. For instance, if
* alignmentMap represented an alignment between two proteins with different
* numbering schemes, the identity function could calculate the offset
* between residue numbers, eg I(x) = x-offset.
*
* When an identity function is provided, the returned function calculates
* f^k(x) = f(I( f(I( ... f(x) ... )) )).
*
* @param <S>
* @param <T>
* @param alignmentMap The input function, as a map (see {@link AlignmentTools#alignmentAsMap(AFPChain)})
* @param identity An identity-like function providing the isomorphism between
* the codomain of alignmentMap (of type <T>) and the domain (type <S>).
* @param k The number of times to apply the alignment
* @return A new alignment. If the input function is not automorphic
* (one-to-one), then some inputs may map to null, indicating that the
* function is undefined for that input.
*/
public static <S,T> Map<S,T> applyAlignment(Map<S, T> alignmentMap, Map<T,S> identity, int k) {
// This implementation simply applies the map k times.
// If k were large, it would be more efficient to do this recursively,
// (eg f^4 = (f^2)^2) but k will usually be small.
if(k<0) throw new IllegalArgumentException("k must be positive");
if(k==1) {
return new HashMap<S,T>(alignmentMap);
}
// Convert to lists to establish a fixed order
List<S> preimage = new ArrayList<S>(alignmentMap.keySet()); // currently unmodified
List<S> image = new ArrayList<S>(preimage);
for(int n=1;n<k;n++) {
// apply alignment
for(int i=0;i<image.size();i++) {
S pre = image.get(i);
T intermediate = (pre==null?null: alignmentMap.get(pre));
S post = (intermediate==null?null: identity.get(intermediate));
image.set(i, post);
}
}
Map<S, T> imageMap = new HashMap<S,T>(alignmentMap.size());
//TODO handle nulls consistently.
// assure that all the residues in the domain are valid keys
/*
for(int i=0;i<preimage.size();i++) {
S pre = preimage.get(i);
T intermediate = (pre==null?null: alignmentMap.get(pre));
S post = (intermediate==null?null: identity.get(intermediate));
imageMap.put(post, null);
}
*/
// now populate with actual values
for(int i=0;i<preimage.size();i++) {
S pre = preimage.get(i);
// image is currently f^k-1(x), so take the final step
S preK1 = image.get(i);
T postK = (preK1==null?null: alignmentMap.get(preK1));
imageMap.put(pre,postK);
}
return imageMap;
}
/**
* Helper for {@link #getSymmetryOrder(Map, Map, int, float)} with a true
* identity function (X->X).
*
* <p>This method should only be used in cases where the two proteins
* aligned have identical numbering, as for self-alignments. See
* {@link #getSymmetryOrder(AFPChain, int, float)} for a way to guess
* the sequential correspondence between two proteins.
*
* @param alignment
* @param maxSymmetry
* @param minimumMetricChange
* @return
*/
public static int getSymmetryOrder(Map<Integer, Integer> alignment,
final int maxSymmetry, final float minimumMetricChange) {
return getSymmetryOrder(alignment, new IdentityMap<Integer>(), maxSymmetry, minimumMetricChange);
}
/**
* Tries to detect symmetry in an alignment.
*
* <p>Conceptually, an alignment is a function f:A->B between two sets of
* integers. The function may have simple topology (meaning that if two
* elements of A are close, then their images in B will also be close), or
* may have more complex topology (such as a circular permutation). This
* function checks <i>alignment</i> against a reference function
* <i>identity</i>, which should have simple topology. It then tries to
* determine the symmetry order of <i>alignment</i> relative to
* <i>identity</i>, up to a maximum order of <i>maxSymmetry</i>.
*
*
* <p><strong>Details</strong><br/>
* Considers the offset (in number of residues) which a residue moves
* after undergoing <i>n</i> alternating transforms by alignment and
* identity. If <i>n</i> corresponds to the intrinsic order of the alignment,
* this will be small. This algorithm tries increasing values of <i>n</i>
* and looks for abrupt decreases in the root mean squared offset.
* If none are found at <i>n</i><=maxSymmetry, the alignment is reported as
* non-symmetric.
*
* @param alignment The alignment to test for symmetry
* @param identity An alignment with simple topology which approximates
* the sequential relationship between the two proteins. Should map in the
* reverse direction from alignment.
* @param maxSymmetry Maximum symmetry to consider. High values increase
* the calculation time and can lead to overfitting.
* @param minimumMetricChange Percent decrease in root mean squared offsets
* in order to declare symmetry. 0.4f seems to work well for CeSymm.
* @return The order of symmetry of alignment, or 1 if no order <=
* maxSymmetry is found.
*
* @see IdentityMap For a simple identity function
*/
public static int getSymmetryOrder(Map<Integer, Integer> alignment, Map<Integer,Integer> identity,
final int maxSymmetry, final float minimumMetricChange) {
List<Integer> preimage = new ArrayList<Integer>(alignment.keySet()); // currently unmodified
List<Integer> image = new ArrayList<Integer>(preimage);
int bestSymmetry = 1;
double bestMetric = Double.POSITIVE_INFINITY; //lower is better
boolean foundSymmetry = false;
if(debug) {
logger.trace("Symm\tPos\tDelta");
}
for(int n=1;n<=maxSymmetry;n++) {
int deltasSq = 0;
int numDeltas = 0;
// apply alignment
for(int i=0;i<image.size();i++) {
Integer pre = image.get(i);
Integer intermediate = (pre==null?null: alignment.get(pre));
Integer post = (intermediate==null?null: identity.get(intermediate));
image.set(i, post);
if(post != null) {
int delta = post-preimage.get(i);
deltasSq += delta*delta;
numDeltas++;
if(debug) {
logger.debug("%d\t%d\t%d\n",n,preimage.get(i),delta);
}
}
}
// Metrics: RMS compensates for the trend of smaller numDeltas with higher order
// Not normalizing by numDeltas favors smaller orders
double metric = Math.sqrt((double)deltasSq/numDeltas); // root mean squared distance
if(!foundSymmetry && metric < bestMetric * minimumMetricChange) {
// n = 1 is never the best symmetry
if(bestMetric < Double.POSITIVE_INFINITY) {
foundSymmetry = true;
}
bestSymmetry = n;
bestMetric = metric;
}
// When debugging need to loop over everything. Unneeded in production
if(!debug && foundSymmetry) {
break;
}
}
if(foundSymmetry) {
return bestSymmetry;
} else {
return 1;
}
}
/**
* Guesses the order of symmetry in an alignment
*
* <p>Uses {@link #getSymmetryOrder(Map alignment, Map identity, int, float)}
* to determine the the symmetry order. For the identity alignment, sorts
* the aligned residues of each protein sequentially, then defines the ith
* residues of each protein to be equivalent.
*
* <p>Note that the selection of the identity alignment here is <i>very</i>
* naive, and only works for proteins with very good coverage. Wherever
* possible, it is better to construct an identity function explicitly
* from a sequence alignment (or use an {@link IdentityMap} for internally
* symmetric proteins) and use {@link #getSymmetryOrder(Map, Map, int, float)}.
*/
public static int getSymmetryOrder(AFPChain afpChain, int maxSymmetry, float minimumMetricChange) throws StructureException {
// alignment comes from the afpChain alignment
Map<Integer,Integer> alignment = AlignmentTools.alignmentAsMap(afpChain);
// Now construct identity to map aligned residues in sequential order
Map<Integer, Integer> identity = guessSequentialAlignment(alignment, true);
return AlignmentTools.getSymmetryOrder(alignment,
identity,
maxSymmetry, minimumMetricChange);
}
/**
* Takes a potentially non-sequential alignment and guesses a sequential
* version of it. Residues from each structure are sorted sequentially and
* then compared directly.
*
* <p>The results of this method are consistent with what one might expect
* from an identity function, and are therefore useful with
* {@link #getSymmetryOrder(Map, Map identity, int, float)}.
* <ul>
* <li>Perfect self-alignments will have the same pre-image and image,
* so will map X->X</li>
* <li>Gaps and alignment errors will cause errors in the resulting map,
* but only locally. Errors do not propagate through the whole
* alignment.</li>
* </ul>
*
* <h4>Example:</h4>
* A non sequential alignment, represented schematically as
* <pre>
* 12456789
* 78912345</pre>
* would result in a map
* <pre>
* 12456789
* 12345789</pre>
* @param alignment The non-sequential input alignment
* @param inverseAlignment If false, map from structure1 to structure2. If
* true, generate the inverse of that map.
* @return A mapping from sequential residues of one protein to those of the other
* @throws IllegalArgumentException if the input alignment is not one-to-one.
*/
public static Map<Integer, Integer> guessSequentialAlignment(
Map<Integer,Integer> alignment, boolean inverseAlignment) {
Map<Integer,Integer> identity = new HashMap<Integer,Integer>();
SortedSet<Integer> aligned1 = new TreeSet<Integer>();
SortedSet<Integer> aligned2 = new TreeSet<Integer>();
for(Entry<Integer,Integer> pair : alignment.entrySet()) {
aligned1.add(pair.getKey());
if( !aligned2.add(pair.getValue()) )
throw new IllegalArgumentException("Alignment is not one-to-one for residue "+pair.getValue()+" of the second structure.");
}
Iterator<Integer> it1 = aligned1.iterator();
Iterator<Integer> it2 = aligned2.iterator();
while(it1.hasNext()) {
if(inverseAlignment) { // 2->1
identity.put(it2.next(),it1.next());
} else { // 1->2
identity.put(it1.next(),it2.next());
}
}
return identity;
}
/**
* Retrieves the optimum alignment from an AFPChain and returns it as a
* java collection. The result is indexed in the same way as
* {@link AFPChain#getOptAln()}, but has the correct size().
* <pre>
* List<List<List<Integer>>> aln = getOptAlnAsList(AFPChain afpChain);
* aln.get(blockNum).get(structureNum={0,1}).get(pos)</pre>
*
* @param afpChain
* @return
*/
public static List<List<List<Integer>>> getOptAlnAsList(AFPChain afpChain) {
int[][][] optAln = afpChain.getOptAln();
int[] optLen = afpChain.getOptLen();
List<List<List<Integer>>> blocks = new ArrayList<List<List<Integer>>>(afpChain.getBlockNum());
for(int blockNum=0;blockNum<afpChain.getBlockNum();blockNum++) {
//TODO could improve speed an memory by wrapping the arrays with
// an unmodifiable list, similar to Arrays.asList(...) but with the
// correct size parameter.
List<Integer> align1 = new ArrayList<Integer>(optLen[blockNum]);
List<Integer> align2 = new ArrayList<Integer>(optLen[blockNum]);
for(int pos=0;pos<optLen[blockNum];pos++) {
align1.add(optAln[blockNum][0][pos]);
align2.add(optAln[blockNum][1][pos]);
}
List<List<Integer>> block = new ArrayList<List<Integer>>(2);
block.add(align1);
block.add(align2);
blocks.add(block);
}
return blocks;
}
/**
* A Map<K,V> can be viewed as a function from K to V. This class represents
* the identity function. Getting a value results in the value itself.
*
* <p>The class is a bit inconsistent when representing its contents. On
* the one hand, containsKey(key) is true for all objects. However,
* attempting to iterate through the values returns an empty set.
*
* @author Spencer Bliven
*
* @param <K>
*/
public static class IdentityMap<K> extends AbstractMap<K,K> {
public IdentityMap() {}
/**
* @param key
* @return the key
* @throws ClassCastException if key is not of type K
*/
@SuppressWarnings("unchecked")
@Override
public K get(Object key) {
return (K)key;
}
/**
* Always returns the empty set
*/
@Override
public Set<java.util.Map.Entry<K, K>> entrySet() {
return Collections.emptySet();
}
@Override
public boolean containsKey(Object key) {
return true;
}
}
/**
* Fundamentally, an alignment is just a list of aligned residues in each
* protein. This method converts two lists of ResidueNumbers into an
* AFPChain.
*
* <p>Parameters are filled with defaults (often null) or sometimes
* calculated.
*
* <p>For a way to modify the alignment of an existing AFPChain, see
* {@link AlignmentTools#replaceOptAln(AFPChain, Atom[], Atom[], Map)}
* @param ca1 CA atoms of the first protein
* @param ca2 CA atoms of the second protein
* @param aligned1 A list of aligned residues from the first protein
* @param aligned2 A list of aligned residues from the second protein.
* Must be the same length as aligned1.
* @return An AFPChain representing the alignment. Many properties may be
* null or another default.
* @throws StructureException if an error occured during superposition
* @throws IllegalArgumentException if aligned1 and aligned2 have different
* lengths
* @see AlignmentTools#replaceOptAln(AFPChain, Atom[], Atom[], Map)
*/
public static AFPChain createAFPChain(Atom[] ca1, Atom[] ca2,
ResidueNumber[] aligned1, ResidueNumber[] aligned2 ) throws StructureException {
//input validation
int alnLen = aligned1.length;
if(alnLen != aligned2.length) {
throw new IllegalArgumentException("Alignment lengths are not equal");
}
AFPChain a = new AFPChain(AFPChain.UNKNOWN_ALGORITHM);
try {
a.setName1(ca1[0].getGroup().getChain().getStructure().getName());
if(ca2[0].getGroup().getChain().getStructure() != null) {
// common case for cloned ca2
a.setName2(ca2[0].getGroup().getChain().getStructure().getName());
}
} catch(Exception e) {
// One of the structures wasn't fully created. Ignore
}
a.setBlockNum(1);
a.setCa1Length(ca1.length);
a.setCa2Length(ca2.length);
a.setOptLength(alnLen);
a.setOptLen(new int[] {alnLen});
Matrix[] ms = new Matrix[a.getBlockNum()];
a.setBlockRotationMatrix(ms);
Atom[] blockShiftVector = new Atom[a.getBlockNum()];
a.setBlockShiftVector(blockShiftVector);
String[][][] pdbAln = new String[1][2][alnLen];
for(int i=0;i<alnLen;i++) {
pdbAln[0][0][i] = aligned1[i].getChainName()+":"+aligned1[i];
pdbAln[0][1][i] = aligned2[i].getChainName()+":"+aligned2[i];
}
a.setPdbAln(pdbAln);
// convert pdbAln to optAln, and fill in some other basic parameters
AFPChainXMLParser.rebuildAFPChain(a, ca1, ca2);
return a;
// Currently a single block. Split into several blocks by sequence if needed
// return AlignmentTools.splitBlocksByTopology(a,ca1,ca2);
}
/**
*
* @param a
* @param ca1
* @param ca2
* @return
* @throws StructureException if an error occurred during superposition
*/
public static AFPChain splitBlocksByTopology(AFPChain a, Atom[] ca1, Atom[] ca2) throws StructureException {
int[][][] optAln = a.getOptAln();
int blockNum = a.getBlockNum();
int[] optLen = a.getOptLen();
// Determine block lengths
// Split blocks if residue indices don't increase monotonically
List<Integer> newBlkLen = new ArrayList<Integer>();
boolean blockChanged = false;
for(int blk=0;blk<blockNum;blk++) {
int currLen=1;
for(int pos=1;pos<optLen[blk];pos++) {
if( optAln[blk][0][pos] <= optAln[blk][0][pos-1]
|| optAln[blk][1][pos] <= optAln[blk][1][pos-1] )
{
//start a new block
newBlkLen.add(currLen);
currLen = 0;
blockChanged = true;
}
currLen++;
}
if(optLen[blk] < 2 ) {
newBlkLen.add(optLen[blk]);
} else {
newBlkLen.add(currLen);
}
}
// Check if anything needs to be split
if( !blockChanged ) {
return a;
}
// Split blocks
List<int[][]> blocks = new ArrayList<int[][]>( newBlkLen.size() );
int oldBlk = 0;
int pos = 0;
for(int blkLen : newBlkLen) {
if( blkLen == optLen[oldBlk] ) {
assert(pos == 0); //should be the whole block
// Use the old block
blocks.add(optAln[oldBlk]);
} else {
int[][] newBlock = new int[2][blkLen];
assert( pos+blkLen <= optLen[oldBlk] ); // don't overrun block
for(int i=0; i<blkLen;i++) {
newBlock[0][i] = optAln[oldBlk][0][pos + i];
newBlock[1][i] = optAln[oldBlk][1][pos + i];
}
pos += blkLen;
blocks.add(newBlock);
if( pos == optLen[oldBlk] ) {
// Finished this oldBlk, start the next
oldBlk++;
pos = 0;
}
}
}
// Store new blocks
int[][][] newOptAln = blocks.toArray(new int[0][][]);
int[] newBlockLens = new int[newBlkLen.size()];
for(int i=0;i<newBlkLen.size();i++) {
newBlockLens[i] = newBlkLen.get(i);
}
return replaceOptAln(a, ca1, ca2, blocks.size(), newBlockLens, newOptAln);
}
/**
* It replaces an optimal alignment of an AFPChain and calculates all the new alignment scores and variables.
*/
public static AFPChain replaceOptAln(int[][][] newAlgn, AFPChain afpChain, Atom[] ca1, Atom[] ca2) throws StructureException {
//The order is the number of groups in the newAlgn
int order = newAlgn.length;
//Calculate the alignment length from all the subunits lengths
int[] optLens = new int[order];
for(int s=0;s<order;s++) {
optLens[s] = newAlgn[s][0].length;
}
int optLength = 0;
for(int s=0;s<order;s++) {
optLength += optLens[s];
}
//Create a copy of the original AFPChain and set everything needed for the structure update
AFPChain copyAFP = (AFPChain) afpChain.clone();
//Set the new parameters of the optimal alignment
copyAFP.setOptLength(optLength);
copyAFP.setOptLen(optLens);
copyAFP.setOptAln(newAlgn);
//Set the block information of the new alignment
copyAFP.setBlockNum(order);
copyAFP.setBlockSize(optLens);
copyAFP.setBlockResList(newAlgn);
copyAFP.setBlockResSize(optLens);
copyAFP.setBlockGap(calculateBlockGap(newAlgn));
//Recalculate properties: superposition, tm-score, etc
Atom[] ca2clone = StructureTools.cloneAtomArray(ca2); // don't modify ca2 positions
AlignmentTools.updateSuperposition(copyAFP, ca1, ca2clone);
//It re-does the sequence alignment strings from the OptAlgn information only
copyAFP.setAlnsymb(null);
AFPAlignmentDisplay.getAlign(copyAFP, ca1, ca2clone);
return copyAFP;
}
/**
* Takes an AFPChain and replaces the optimal alignment based on an alignment map
*
* <p>Parameters are filled with defaults (often null) or sometimes
* calculated.
*
* <p>For a way to create a new AFPChain, see
* {@link AlignmentTools#createAFPChain(Atom[], Atom[], ResidueNumber[], ResidueNumber[])}
*
* @param afpChain The alignment to be modified
* @param alignment The new alignment, as a Map
* @throws StructureException if an error occurred during superposition
* @see AlignmentTools#createAFPChain(Atom[], Atom[], ResidueNumber[], ResidueNumber[])
*/
public static AFPChain replaceOptAln(AFPChain afpChain, Atom[] ca1, Atom[] ca2,
Map<Integer, Integer> alignment) throws StructureException {
// Determine block lengths
// Sort ca1 indices, then start a new block whenever ca2 indices aren't
// increasing monotonically.
Integer[] res1 = alignment.keySet().toArray(new Integer[0]);
Arrays.sort(res1);
List<Integer> blockLens = new ArrayList<Integer>(2);
int optLength = 0;
Integer lastRes = alignment.get(res1[0]);
int blkLen = lastRes==null?0:1;
for(int i=1;i<res1.length;i++) {
Integer currRes = alignment.get(res1[i]); //res2 index
assert(currRes != null);// could be converted to if statement if assertion doesn't hold; just modify below as well.
if(lastRes<currRes) {
blkLen++;
} else {
// CP!
blockLens.add(blkLen);
optLength+=blkLen;
blkLen = 1;
}
lastRes = currRes;
}
blockLens.add(blkLen);
optLength+=blkLen;
// Create array structure for alignment
int[][][] optAln = new int[blockLens.size()][][];
int pos1 = 0; //index into res1
for(int blk=0;blk<blockLens.size();blk++) {
optAln[blk] = new int[2][];
blkLen = blockLens.get(blk);
optAln[blk][0] = new int[blkLen];
optAln[blk][1] = new int[blkLen];
int pos = 0; //index into optAln
while(pos<blkLen) {
optAln[blk][0][pos]=res1[pos1];
Integer currRes = alignment.get(res1[pos1]);
optAln[blk][1][pos]=currRes;
pos++;
pos1++;
}
}
assert(pos1 == optLength);
// Create length array
int[] optLens = new int[blockLens.size()];
for(int i=0;i<blockLens.size();i++) {
optLens[i] = blockLens.get(i);
}
return replaceOptAln(afpChain, ca1, ca2, blockLens.size(), optLens, optAln);
}
/**
* @param afpChain Input afpchain. UNMODIFIED
* @param ca1
* @param ca2
* @param optLens
* @param optAln
* @return A NEW AfpChain based off the input but with the optAln modified
* @throws StructureException if an error occured during superposition
*/
public static AFPChain replaceOptAln(AFPChain afpChain, Atom[] ca1, Atom[] ca2,
int blockNum, int[] optLens, int[][][] optAln) throws StructureException {
int optLength = 0;
for( int blk=0;blk<blockNum;blk++) {
optLength += optLens[blk];
}
//set everything
AFPChain refinedAFP = (AFPChain) afpChain.clone();
refinedAFP.setOptLength(optLength);
refinedAFP.setBlockSize(optLens);
refinedAFP.setOptLen(optLens);
refinedAFP.setOptAln(optAln);
refinedAFP.setBlockNum(blockNum);
//TODO recalculate properties: superposition, tm-score, etc
Atom[] ca2clone = StructureTools.cloneAtomArray(ca2); // don't modify ca2 positions
AlignmentTools.updateSuperposition(refinedAFP, ca1, ca2clone);
AFPAlignmentDisplay.getAlign(refinedAFP, ca1, ca2clone);
return refinedAFP;
}
/**
* After the alignment changes (optAln, optLen, blockNum, at a minimum),
* many other properties which depend on the superposition will be invalid.
*
* This method re-runs a rigid superposition over the whole alignment
* and repopulates the required properties, including RMSD (TotalRMSD) and
* TM-Score.
* @param afpChain
* @param ca1
* @param ca2 Second set of ca atoms. Will be modified based on the superposition
* @throws StructureException
* @see {@link CECalculator#calc_rmsd(Atom[], Atom[], int, boolean)}
* contains much of the same code, but stores results in a CECalculator
* instance rather than an AFPChain
*/
public static void updateSuperposition(AFPChain afpChain, Atom[] ca1,
Atom[] ca2) throws StructureException {
//Update ca information, because the atom array might also be changed
afpChain.setCa1Length(ca1.length);
afpChain.setCa2Length(ca2.length);
//We need this to get the correct superposition
int[] focusRes1 = afpChain.getFocusRes1();
int[] focusRes2 = afpChain.getFocusRes2();
if (focusRes1 == null) {
focusRes1 = new int[afpChain.getCa1Length()];
afpChain.setFocusRes1(focusRes1);
}
if (focusRes2 == null) {
focusRes2 = new int[afpChain.getCa2Length()];
afpChain.setFocusRes2(focusRes2);
}
if (afpChain.getNrEQR() == 0) return;
// create new arrays for the subset of atoms in the alignment.
Atom[] ca1aligned = new Atom[afpChain.getOptLength()];
Atom[] ca2aligned = new Atom[afpChain.getOptLength()];
fillAlignedAtomArrays(afpChain, ca1, ca2, ca1aligned, ca2aligned);
//Superimpose the two structures in correspondance to the new alignment
Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1aligned),
Calc.atomsToPoints(ca2aligned));
Matrix matrix = Matrices.getRotationJAMA(trans);
Atom shift = Calc.getTranslationVector(trans);
Matrix[] blockMxs = new Matrix[afpChain.getBlockNum()];
Arrays.fill(blockMxs, matrix);
afpChain.setBlockRotationMatrix(blockMxs);
Atom[] blockShifts = new Atom[afpChain.getBlockNum()];
Arrays.fill(blockShifts, shift);
afpChain.setBlockShiftVector(blockShifts);
for (Atom a : ca2aligned) {
Calc.rotate(a, matrix);
Calc.shift(a, shift);
}
//Calculate the RMSD and TM score for the new alignment
double rmsd = Calc.rmsd(ca1aligned, ca2aligned);
double tmScore = Calc.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length);
afpChain.setTotalRmsdOpt(rmsd);
afpChain.setTMScore(tmScore);
int[] blockLens = afpChain.getOptLen();
int[][][] optAln = afpChain.getOptAln();
//Calculate the RMSD and TM score for every block of the new alignment
double[] blockRMSD = new double[afpChain.getBlockNum()];
double[] blockScore = new double[afpChain.getBlockNum()];
for (int k=0; k<afpChain.getBlockNum(); k++){
//Create the atom arrays corresponding to the aligned residues in the block
Atom[] ca1block = new Atom[afpChain.getOptLen()[k]];
Atom[] ca2block = new Atom[afpChain.getOptLen()[k]];
int position=0;
for(int i=0;i<blockLens[k];i++) {
int pos1 = optAln[k][0][i];
int pos2 = optAln[k][1][i];
Atom a1 = ca1[pos1];
Atom a2 = (Atom) ca2[pos2].clone();
ca1block[position] = a1;
ca2block[position] = a2;
position++;
}
if (position != afpChain.getOptLen()[k]){
logger.warn("AFPChainScorer getTMScore: Problems reconstructing block alignment! nr of loaded atoms is " + position + " but should be " + afpChain.getOptLen()[k]);
// we need to resize the array, because we allocated too many atoms earlier on.
ca1block = (Atom[]) resizeArray(ca1block, position);
ca2block = (Atom[]) resizeArray(ca2block, position);
}
//Superimpose the two block structures
Matrix4d transb = SuperPositions.superpose(Calc.atomsToPoints(ca1block),
Calc.atomsToPoints(ca2block));
blockMxs[k] = Matrices.getRotationJAMA(trans);
blockShifts[k] = Calc.getTranslationVector(trans);
Calc.transform(ca2block, transb);
//Calculate the RMSD and TM score for the block
double rmsdb = Calc.rmsd(ca1block, ca2block);
double tmScoreb = Calc.getTMScore(ca1block, ca2block, ca1.length, ca2.length);
blockRMSD[k] = rmsdb;
blockScore[k] = tmScoreb;
}
afpChain.setOptRmsd(blockRMSD);
afpChain.setBlockRmsd(blockRMSD);
afpChain.setBlockScore(blockScore);
}
/**
* Reallocates an array with a new size, and copies the contents
* of the old array to the new array.
* @param oldArray the old array, to be reallocated.
* @param newSize the new array size.
* @return A new array with the same contents.
*/
public static Object resizeArray (Object oldArray, int newSize) {
int oldSize = java.lang.reflect.Array.getLength(oldArray);
@SuppressWarnings("rawtypes")
Class elementType = oldArray.getClass().getComponentType();
Object newArray = java.lang.reflect.Array.newInstance(
elementType,newSize);
int preserveLength = Math.min(oldSize,newSize);
if (preserveLength > 0)
System.arraycopy (oldArray,0,newArray,0,preserveLength);
return newArray;
}
/**
* Print an alignment map in a concise representation. Edges are given
* as two numbers separated by '>'. They are chained together where possible,
* or separated by spaces where disjoint or branched.
*
* <p>Note that more concise representations may be possible.</p>
*
* Examples:
* <li>1>2>3>1</li>
* <li>1>2>3>2 4>3</li>
*
* @param alignment The input function, as a map (see {@link AlignmentTools#alignmentAsMap(AFPChain)})
* @param identity An identity-like function providing the isomorphism between
* the codomain of alignment (of type <T>) and the domain (type <S>).
* @return
*/
public static <S,T> String toConciseAlignmentString(Map<S,T> alignment, Map<T,S> identity) {
// Clone input to prevent changes
Map<S,T> alig = new HashMap<S,T>(alignment);
// Generate inverse alignment
Map<S,List<S>> inverse = new HashMap<S,List<S>>();
for(Entry<S,T> e: alig.entrySet()) {
S val = identity.get(e.getValue());
if( inverse.containsKey(val) ) {
List<S> l = inverse.get(val);
l.add(e.getKey());
} else {
List<S> l = new ArrayList<S>();
l.add(e.getKey());
inverse.put(val,l);
}
}
StringBuilder str = new StringBuilder();
while(!alig.isEmpty()){
// Pick an edge and work upstream to a root or cycle
S seedNode = alig.keySet().iterator().next();
S node = seedNode;
if( inverse.containsKey(seedNode)) {
node = inverse.get(seedNode).iterator().next();
while( node != seedNode && inverse.containsKey(node)) {
node = inverse.get(node).iterator().next();
}
}
// Now work downstream, deleting edges as we go
seedNode = node;
str.append(node);
while(alig.containsKey(node)) {
S lastNode = node;
node = identity.get( alig.get(lastNode) );
// Output
str.append('>');
str.append(node);
// Remove edge
alig.remove(lastNode);
List<S> inv = inverse.get(node);
if(inv.size() > 1) {
inv.remove(node);
} else {
inverse.remove(node);
}
}
if(!alig.isEmpty()) {
str.append(' ');
}
}
return str.toString();
}
/**
* @see #toConciseAlignmentString(Map, Map)
*/
public static <T> String toConciseAlignmentString(Map<T, T> alignment) {
return toConciseAlignmentString(alignment, new IdentityMap<T>());
}
/**
* @see #toConciseAlignmentString(Map, Map)
*/
public static Map<Integer, Integer> fromConciseAlignmentString(String string) {
Map<Integer, Integer> map = new HashMap<Integer, Integer>();
boolean matches = true;
while (matches) {
Pattern pattern = Pattern.compile("(\\d+)>(\\d+)");
Matcher matcher = pattern.matcher(string);
matches = matcher.find();
if (matches) {
Integer from = Integer.parseInt(matcher.group(1));
Integer to = Integer.parseInt(matcher.group(2));
map.put(from, to);
string = string.substring(matcher.end(1) + 1);
}
}
return map;
}
/**
* Method that calculates the number of gaps in each subunit block of an optimal AFP alignment.
*
* INPUT: an optimal alignment in the format int[][][].
* OUTPUT: an int[] array of <order> length containing the gaps in each block as int[block].
*/
public static int[] calculateBlockGap(int[][][] optAln){
//Initialize the array to be returned
int [] blockGap = new int[optAln.length];
//Loop for every block and look in both chains for non-contiguous residues.
for (int i=0; i<optAln.length; i++){
int gaps = 0; //the number of gaps in that block
int last1 = 0; //the last residue position in chain 1
int last2 = 0; //the last residue position in chain 2
//Loop for every position in the block
for (int j=0; j<optAln[i][0].length; j++){
//If the first position is evaluated initialize the last positions
if (j==0){
last1 = optAln[i][0][j];
last2 = optAln[i][1][j];
}
else{
//If one of the positions or both are not contiguous increment the number of gaps
if (optAln[i][0][j] > last1+1 || optAln[i][1][j] > last2+1){
gaps++;
last1 = optAln[i][0][j];
last2 = optAln[i][1][j];
}
//Otherwise just set the last position to the current one
else{
last1 = optAln[i][0][j];
last2 = optAln[i][1][j];
}
}
}
blockGap[i] = gaps;
}
return blockGap;
}
/**
* Creates a simple interaction format (SIF) file for an alignment.
*
* The SIF file can be read by network software (eg Cytoscape) to analyze
* alignments as graphs.
*
* This function creates a graph with residues as nodes and two types of edges:
* 1. backbone edges, which connect adjacent residues in the aligned protein
* 2. alignment edges, which connect aligned residues
*
* @param out Stream to write to
* @param afpChain alignment to write
* @param ca1 First protein, used to generate node names
* @param ca2 Second protein, used to generate node names
* @param backboneInteraction Two-letter string used to identify backbone edges
* @param alignmentInteraction Two-letter string used to identify alignment edges
* @throws IOException
*/
public static void alignmentToSIF(Writer out,AFPChain afpChain,
Atom[] ca1,Atom[] ca2, String backboneInteraction,
String alignmentInteraction) throws IOException {
//out.write("Res1\tInteraction\tRes2\n");
String name1 = afpChain.getName1();
String name2 = afpChain.getName2();
if(name1==null) name1=""; else name1+=":";
if(name2==null) name2=""; else name2+=":";
// Print alignment edges
int nblocks = afpChain.getBlockNum();
int[] blockLen = afpChain.getOptLen();
int[][][] optAlign = afpChain.getOptAln();
for(int b=0;b<nblocks;b++) {
for(int r=0;r<blockLen[b];r++) {
int res1 = optAlign[b][0][r];
int res2 = optAlign[b][1][r];
ResidueNumber rn1 = ca1[res1].getGroup().getResidueNumber();
ResidueNumber rn2 = ca2[res2].getGroup().getResidueNumber();
String node1 = name1+rn1.getChainName()+rn1.toString();
String node2 = name2+rn2.getChainName()+rn2.toString();
out.write(String.format("%s\t%s\t%s\n",node1, alignmentInteraction, node2));
}
}
// Print first backbone edges
ResidueNumber rn = ca1[0].getGroup().getResidueNumber();
String last = name1+rn.getChainName()+rn.toString();
for(int i=1;i<ca1.length;i++) {
rn = ca1[i].getGroup().getResidueNumber();
String curr = name1+rn.getChainName()+rn.toString();
out.write(String.format("%s\t%s\t%s\n",last, backboneInteraction, curr));
last = curr;
}
// Print second backbone edges, if the proteins differ
// Do some quick checks for whether the proteins differ
// (Not perfect, but should detect major differences and CPs.)
if(!name1.equals(name2) ||
ca1.length!=ca2.length ||
(ca1.length>0 && ca1[0].getGroup()!=null && ca2[0].getGroup()!=null &&
!ca1[0].getGroup().getResidueNumber().equals(ca2[0].getGroup().getResidueNumber()) ) ) {
rn = ca2[0].getGroup().getResidueNumber();
last = name2+rn.getChainName()+rn.toString();
for(int i=1;i<ca2.length;i++) {
rn = ca2[i].getGroup().getResidueNumber();
String curr = name2+rn.getChainName()+rn.toString();
out.write(String.format("%s\t%s\t%s\n",last, backboneInteraction, curr));
last = curr;
}
}
}
/** get an artificial List of chains containing the Atoms and groups.
* Does NOT rotate anything.
* @param ca
* @return a list of Chains that is built up from the Atoms in the ca array
* @throws StructureException
*/
public static final List<Chain> getAlignedModel(Atom[] ca){
List<Chain> model = new ArrayList<Chain>();
for ( Atom a: ca){
Group g = a.getGroup();
Chain parentC = g.getChain();
Chain newChain = null;
for ( Chain c : model) {
if ( c.getId().equals(parentC.getId())){
newChain = c;
break;
}
}
if ( newChain == null){
newChain = new ChainImpl();
newChain.setId(parentC.getId());
model.add(newChain);
}
newChain.addGroup(g);
}
return model;
}
/** Get an artifical Structure containing both chains.
* Does NOT rotate anything
* @param ca1
* @param ca2
* @return a structure object containing two models, one for each set of Atoms.
* @throws StructureException
*/
public static final Structure getAlignedStructure(Atom[] ca1, Atom[] ca2) throws StructureException{
/* Previous implementation commented
Structure s = new StructureImpl();
List<Chain>model1 = getAlignedModel(ca1);
List<Chain>model2 = getAlignedModel(ca2);
s.addModel(model1);
s.addModel(model2);
return s;*/
Structure s = new StructureImpl();
List<Chain>model1 = getAlignedModel(ca1);
s.addModel(model1);
List<Chain> model2 = getAlignedModel(ca2);
s.addModel(model2);
return s;
}
/** Rotate the Atoms/Groups so they are aligned for the 3D visualisation
*
* @param afpChain
* @param ca1
* @param ca2
* @return an array of Groups that are transformed for 3D display
* @throws StructureException
*/
public static Group[] prepareGroupsForDisplay(AFPChain afpChain, Atom[] ca1, Atom[] ca2) throws StructureException{
if ( afpChain.getBlockRotationMatrix().length == 0 ) {
// probably the alignment is too short!
System.err.println("No rotation matrix found to rotate 2nd structure!");
afpChain.setBlockRotationMatrix(new Matrix[]{Matrix.identity(3, 3)});
afpChain.setBlockShiftVector(new Atom[]{new AtomImpl()});
}
// List of groups to be rotated according to the alignment
Group[] twistedGroups = new Group[ ca2.length];
//int blockNum = afpChain.getBlockNum();
int i = -1;
// List of groups from the structure not included in ca2 (e.g. ligands)
// Will be rotated according to first block
List<Group> hetatms2 = StructureTools.getUnalignedGroups(ca2);
if ( (afpChain.getAlgorithmName().equals(FatCatRigid.algorithmName) ) || (afpChain.getAlgorithmName().equals(FatCatFlexible.algorithmName) ) ){
for (Atom a: ca2){
i++;
twistedGroups[i]=a.getGroup();
}
twistedGroups = AFPTwister.twistOptimized(afpChain, ca1, ca2);
//} else if (( blockNum == 1 ) || (afpChain.getAlgorithmName().equals(CeCPMain.algorithmName))) {
} else {
Matrix m = afpChain.getBlockRotationMatrix()[ 0];
Atom shift = afpChain.getBlockShiftVector() [ 0 ];
shiftCA2(afpChain, ca2, m,shift, twistedGroups);
}
if ( afpChain.getBlockNum() > 0){
// Superimpose ligands relative to the first block
if( hetatms2.size() > 0 ) {
if ( afpChain.getBlockRotationMatrix().length > 0 ) {
Matrix m1 = afpChain.getBlockRotationMatrix()[0];
//m1.print(3,3);
Atom vector1 = afpChain.getBlockShiftVector()[0];
//System.out.println("shift vector:" + vector1);
for ( Group g : hetatms2){
Calc.rotate(g, m1);
Calc.shift(g,vector1);
}
}
}
}
return twistedGroups;
}
/** only shift CA positions.
*
*/
public static void shiftCA2(AFPChain afpChain, Atom[] ca2, Matrix m, Atom shift, Group[] twistedGroups) {
int i = -1;
for (Atom a: ca2){
i++;
Group g = a.getGroup();
Calc.rotate(g,m);
Calc.shift(g, shift);
if (g.hasAltLoc()){
for (Group alt: g.getAltLocs()){
for (Atom alta : alt.getAtoms()){
if ( g.getAtoms().contains(alta))
continue;
Calc.rotate(alta,m);
Calc.shift(alta,shift);
}
}
}
twistedGroups[i]=g;
}
}
/**
* Fill the aligned Atom arrays with the equivalent residues in the afpChain.
* @param afpChain
* @param ca1
* @param ca2
* @param ca1aligned
* @param ca2aligned
*/
public static void fillAlignedAtomArrays(AFPChain afpChain, Atom[] ca1,
Atom[] ca2, Atom[] ca1aligned, Atom[] ca2aligned) {
int pos=0;
int[] blockLens = afpChain.getOptLen();
int[][][] optAln = afpChain.getOptAln();
assert(afpChain.getBlockNum() <= optAln.length);
for (int block=0; block < afpChain.getBlockNum(); block++) {
for(int i=0;i<blockLens[block];i++) {
int pos1 = optAln[block][0][i];
int pos2 = optAln[block][1][i];
Atom a1 = ca1[pos1];
Atom a2 = (Atom) ca2[pos2].clone();
ca1aligned[pos] = a1;
ca2aligned[pos] = a2;
pos++;
}
}
// this can happen when we load an old XML serialization which did not support modern ChemComp representation of modified residues.
if (pos != afpChain.getOptLength()){
logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + afpChain.getOptLength());
// we need to resize the array, because we allocated too many atoms earlier on.
ca1aligned = (Atom[]) resizeArray(ca1aligned, pos);
ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
}
}
/**
* Find the alignment position with the highest atomic distance between the
* equivalent atomic positions of the arrays and remove it from the
* alignment.
*
* @param afpChain
* original alignment, will be modified
* @param ca1
* atom array, will not be modified
* @param ca2
* atom array, will not be modified
* @return the original alignment, with the alignment position at the
* highest distance removed
* @throws StructureException
*/
public static AFPChain deleteHighestDistanceColumn(AFPChain afpChain,
Atom[] ca1, Atom[] ca2) throws StructureException {
int[][][] optAln = afpChain.getOptAln();
int maxBlock = 0;
int maxPos = 0;
double maxDistance = Double.MIN_VALUE;
for (int b = 0; b < optAln.length; b++) {
for (int p = 0; p < optAln[b][0].length; p++) {
Atom ca2clone = ca2[optAln[b][1][p]];
Calc.rotate(ca2clone, afpChain.getBlockRotationMatrix()[b]);
Calc.shift(ca2clone, afpChain.getBlockShiftVector()[b]);
double distance = Calc.getDistance(ca1[optAln[b][0][p]],
ca2clone);
if (distance > maxDistance) {
maxBlock = b;
maxPos = p;
maxDistance = distance;
}
}
}
return deleteColumn(afpChain, ca1, ca2, maxBlock, maxPos);
}
/**
* Delete an alignment position from the original alignment object.
*
* @param afpChain
* original alignment, will be modified
* @param ca1
* atom array, will not be modified
* @param ca2
* atom array, will not be modified
* @param block
* block of the alignment position
* @param pos
* position index in the block
* @return the original alignment, with the alignment position removed
* @throws StructureException
*/
public static AFPChain deleteColumn(AFPChain afpChain, Atom[] ca1,
Atom[] ca2, int block, int pos) throws StructureException {
// Check validity of the inputs
if (afpChain.getBlockNum() <= block) {
throw new IndexOutOfBoundsException(String.format(
"Block index requested (%d) is higher than the total number of AFPChain blocks (%d).",
block, afpChain.getBlockNum()));
}
if (afpChain.getOptAln()[block][0].length <= pos) {
throw new IndexOutOfBoundsException(String.format(
"Position index requested (%d) is higher than the total number of aligned position in the AFPChain block (%d).",
block, afpChain.getBlockSize()[block]));
}
int[][][] optAln = afpChain.getOptAln();
int[] newPos0 = new int[optAln[block][0].length - 1];
int[] newPos1 = new int[optAln[block][1].length - 1];
int position = 0;
for (int p = 0; p < optAln[block][0].length; p++) {
if (p == pos)
continue;
newPos0[position] = optAln[block][0][p];
newPos1[position] = optAln[block][1][p];
position++;
}
optAln[block][0] = newPos0;
optAln[block][1] = newPos1;
return AlignmentTools.replaceOptAln(optAln, afpChain, ca1, ca2);
}
}