org.biojava.nbio.structure.align.util.AlignmentTools Maven / Gradle / Ivy
Show all versions of biojava-structure Show documentation
/*
* 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.
* 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
* checkWithinBlocks to true 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; blockFor example,
* 1234
* 5678
* becomes
* 1->5
* 2->6
* 3->7
* 4->8
*
* @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 alignmentAsMap(AFPChain afpChain) throws StructureException {
Map map = new HashMap<>();
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
* @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 Map applyAlignment(Map alignmentMap, int k) {
return applyAlignment(alignmentMap, new IdentityMap(), 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
* @param
* @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 ) and the domain (type ).
* @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 Map applyAlignment(Map alignmentMap, Map 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<>(alignmentMap);
}
// Convert to lists to establish a fixed order
List preimage = new ArrayList<>(alignmentMap.keySet()); // currently unmodified
List image = new ArrayList<>(preimage);
for(int n=1;n imageMap = new HashMap<>(alignmentMap.size());
//TODO handle nulls consistently.
// assure that all the residues in the domain are valid keys
/*
for(int i=0;iX).
*
* 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 alignment,
final int maxSymmetry, final float minimumMetricChange) {
return getSymmetryOrder(alignment, new IdentityMap(), maxSymmetry, minimumMetricChange);
}
/**
* Tries to detect symmetry in an alignment.
*
* 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 alignment against a reference function
* identity, which should have simple topology. It then tries to
* determine the symmetry order of alignment relative to
* identity, up to a maximum order of maxSymmetry.
*
*
*
Details
* Considers the offset (in number of residues) which a residue moves
* after undergoing n alternating transforms by alignment and
* identity. If n corresponds to the intrinsic order of the alignment,
* this will be small. This algorithm tries increasing values of n
* and looks for abrupt decreases in the root mean squared offset.
* If none are found at n<=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 alignment, Map identity,
final int maxSymmetry, final float minimumMetricChange) {
List preimage = new ArrayList<>(alignment.keySet()); // currently unmodified
List image = new ArrayList<>(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;iUses {@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.
*
* Note that the selection of the identity alignment here is very
* 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 alignment = AlignmentTools.alignmentAsMap(afpChain);
// Now construct identity to map aligned residues in sequential order
Map 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.
*
* 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)}.
*
* - Perfect self-alignments will have the same pre-image and image,
* so will map X->X
* - Gaps and alignment errors will cause errors in the resulting map,
* but only locally. Errors do not propagate through the whole
* alignment.
*
*
* Example:
* A non sequential alignment, represented schematically as
*
* 12456789
* 78912345
* would result in a map
*
* 12456789
* 12345789
* @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 guessSequentialAlignment(
Map alignment, boolean inverseAlignment) {
Map identity = new HashMap<>();
SortedSet aligned1 = new TreeSet<>();
SortedSet aligned2 = new TreeSet<>();
for(Entry 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 it1 = aligned1.iterator();
Iterator 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().
*
* List>> aln = getOptAlnAsList(AFPChain afpChain);
* aln.get(blockNum).get(structureNum={0,1}).get(pos)
*
* @param afpChain
* @return
*/
public static List>> getOptAlnAsList(AFPChain afpChain) {
int[][][] optAln = afpChain.getOptAln();
int[] optLen = afpChain.getOptLen();
List>> blocks = new ArrayList<>(afpChain.getBlockNum());
for(int blockNum=0;blockNum align1 = new ArrayList<>(optLen[blockNum]);
List align2 = new ArrayList<>(optLen[blockNum]);
for(int pos=0;pos> block = new ArrayList<>(2);
block.add(align1);
block.add(align2);
blocks.add(block);
}
return blocks;
}
/**
* A Map can be viewed as a function from K to V. This class represents
* the identity function. Getting a value results in the value itself.
*
* 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
*/
public static class IdentityMap extends AbstractMap {
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> 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.
*
* Parameters are filled with defaults (often null) or sometimes
* calculated.
*
*
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 newBlkLen = new ArrayList<>();
boolean blockChanged = false;
for(int blk=0;blk blocks = new ArrayList<>( 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; iParameters are filled with defaults (often null) or sometimes
* calculated.
*
* 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 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 blockLens = new ArrayList<>(2);
int optLength = 0;
Integer lastRes = alignment.get(res1[0]);
int blkLen = lastRes==null?0:1;
for(int i=1;i 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.
*
* Note that more concise representations may be possible.
*
* Examples:
* 1>2>3>1
* 1>2>3>2 4>3
*
* @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 ) and the domain (type ).
* @return
*/
public static String toConciseAlignmentString(Map alignment, Map identity) {
// Clone input to prevent changes
Map alig = new HashMap<>(alignment);
// Generate inverse alignment
Map> inverse = new HashMap<>();
for(Entry e: alig.entrySet()) {
S val = identity.get(e.getValue());
if( inverse.containsKey(val) ) {
List l = inverse.get(val);
l.add(e.getKey());
} else {
List l = new ArrayList<>();
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 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 String toConciseAlignmentString(Map alignment) {
return toConciseAlignmentString(alignment, new IdentityMap());
}
/**
* @see #toConciseAlignmentString(Map, Map)
*/
public static Map fromConciseAlignmentString(String string) {
Map map = new HashMap<>();
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 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 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;b0 && 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 getAlignedModel(Atom[] ca){
List model = new ArrayList<>();
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();
Listmodel1 = getAlignedModel(ca1);
Listmodel2 = getAlignedModel(ca2);
s.addModel(model1);
s.addModel(model2);
return s;*/
Structure s = new StructureImpl();
Listmodel1 = getAlignedModel(ca1);
s.addModel(model1);
List 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 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 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);
}
}