Java Code Examples for org.biojava.nbio.structure.align.model.AFPChain#getOptLength()

/**
 * Refines a CE-Symm alignment so that it is perfectly symmetric.
 *
 * The resulting alignment will have a one-to-one correspondance between
 * aligned residues of each symmetric part.
 *
 * @param afpChain Input alignment from CE-Symm
 * @param k Symmetry order. This can be guessed by {@link CeSymm#getSymmetryOrder(AFPChain)}
 * @return The refined alignment
 * @throws StructureException
 * @throws RefinerFailedException
 */
public static AFPChain refineSymmetry(AFPChain afpChain, Atom[] ca1, Atom[] ca2, int k)
		throws StructureException, RefinerFailedException {

	// Transform alignment to a Map
	Map<Integer, Integer> alignment = AlignmentTools.alignmentAsMap(afpChain);

	// Refine the alignment Map
	Map<Integer, Integer> refined = refineSymmetry(alignment, k);
	if (refined.size() < 1)
		throw new RefinerFailedException("Refiner returned empty alignment");

	//Substitute and partition the alignment
	try {
		AFPChain refinedAFP = AlignmentTools.replaceOptAln(afpChain, ca1, ca2, refined);
		refinedAFP = partitionAFPchain(refinedAFP, ca1, ca2, k);
		if (refinedAFP.getOptLength() < 1)
			throw new IndexOutOfBoundsException("Refined alignment is empty.");
		return refinedAFP;
	} catch (IndexOutOfBoundsException e){
		// This Exception is thrown when the refined alignment is not consistent
		throw new RefinerFailedException("Refiner failure: non-consistent result", e);
	}
}
 

/** Circular permutation specific code to be run after the standard CE alignment
 *
 * @param afpChain The finished alignment
 * @param ca1 CA atoms of the first protein
 * @param ca2m A duplicated copy of the second protein
 * @param calculator The CECalculator used to create afpChain
 * @param param Parameters
 * @throws StructureException
 */
public static AFPChain postProcessAlignment(AFPChain afpChain, Atom[] ca1, Atom[] ca2m,CECalculator calculator, CECPParameters param ) throws StructureException{

	// remove bottom half of the matrix
	Matrix doubledMatrix = afpChain.getDistanceMatrix();

	// the matrix can be null if the alignment is too short.
	if ( doubledMatrix != null ) {
		assert(doubledMatrix.getRowDimension() == ca1.length);
		assert(doubledMatrix.getColumnDimension() == ca2m.length);

		Matrix singleMatrix = doubledMatrix.getMatrix(0, ca1.length-1, 0, (ca2m.length/2)-1);
		assert(singleMatrix.getRowDimension() == ca1.length);
		assert(singleMatrix.getColumnDimension() == (ca2m.length/2));

		afpChain.setDistanceMatrix(singleMatrix);
	}
	// Check for circular permutations
	int alignLen = afpChain.getOptLength();
	if ( alignLen > 0) {
		afpChain = filterDuplicateAFPs(afpChain,calculator,ca1,ca2m,param);
	}
	return afpChain;
}
 

/**
 *  Partitions an afpChain alignment into order blocks of aligned residues.
 *
 * @param afpChain
 * @param ca1
 * @param ca2
 * @param order
 * @return
 * @throws StructureException
 */
private static AFPChain partitionAFPchain(AFPChain afpChain,
		Atom[] ca1, Atom[] ca2, int order) throws StructureException{

	int[][][] newAlgn = new int[order][][];
	int repeatLen = afpChain.getOptLength()/order;

	//Extract all the residues considered in the first chain of the alignment
	List<Integer> alignedRes = new ArrayList<Integer>();
	for (int su=0; su<afpChain.getBlockNum(); su++){
		for (int i=0; i<afpChain.getOptLen()[su]; i++){
			alignedRes.add(afpChain.getOptAln()[su][0][i]);
		}
	}

	//Build the new alignment
	for (int su=0; su<order; su++){
		newAlgn[su] = new int[2][];
		newAlgn[su][0] = new  int[repeatLen];
		newAlgn[su][1] = new  int[repeatLen];
		for (int i=0; i<repeatLen; i++){
			newAlgn[su][0][i] = alignedRes.get(repeatLen*su+i);
			newAlgn[su][1][i] = alignedRes.get(
					(repeatLen*(su+1)+i)%alignedRes.size());
		}
	}

	return AlignmentTools.replaceOptAln(newAlgn, afpChain, ca1, ca2);
}
 

public double calSigAll(FatCatParameters params, AFPChain afpChain){

		int twist = params.getMaxTra();
		int sparse = params.getSparse();
		int len1 = afpChain.getCa1Length();
		int len2 = afpChain.getCa2Length();
		double score = afpChain.getAlignScore();
		double rmsd = afpChain.getTotalRmsdOpt();
		int optLen = afpChain.getOptLength();
		int r = afpChain.getBlockNum() -1;

		return calSigAll(twist,sparse,len1, len2,score,rmsd,optLen,r);

	}
 

public double calNS(FatCatParameters params, AFPChain afpChain){
	int len1 = afpChain.getCa1Length();
	int len2 = afpChain.getCa2Length();
	double score =afpChain.getAlignScore();
	double rmsd = afpChain.getTotalRmsdOpt();
	int optLen  = afpChain.getOptLength();
	int r = afpChain.getBlockNum() -1;
	return calNS(len1, len2,score,rmsd,optLen,r);
}
 

/**
 * 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);
	}

}
 

public  static double getTMScore(AFPChain align, Atom[] ca1, Atom[] ca2, boolean normalizeMin) throws StructureException
{
	if ( align.getNrEQR() == 0)
		return -1;


	// Create new arrays for the subset of atoms in the alignment.
	Atom[] ca1aligned = new Atom[align.getOptLength()];
	Atom[] ca2aligned = new Atom[align.getOptLength()];
	int pos=0;
	int[] blockLens = align.getOptLen();
	int[][][] optAln = align.getOptAln();
	assert(align.getBlockNum() <= optAln.length);

	for(int block=0;block< align.getBlockNum();block++) {

		if ( ! ( blockLens[block] <= optAln[block][0].length)) {
			logger.warn("AFPChainScorer getTMScore: errors reconstructing alignment block [" + block + "]. Length is " + blockLens[block] + " but should be <=" + optAln[block][0].length);
		}

		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 != align.getOptLength()){
		logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + align.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);
	}
	//Superimpose
	Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1aligned),
			Calc.atomsToPoints(ca2aligned));

	Calc.transform(ca2aligned, trans);

	return Calc.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length, normalizeMin);
}
 

@Override
public AFPChain align(Atom[] ca1, Atom[] ca2, Object parameters)
throws StructureException {

	if ( parameters == null) {
		throw new IllegalArgumentException("Got null instead of SmithWaterman3DParameters!");
	}
	if ( ! (parameters instanceof SmithWaterman3DParameters))
		throw new IllegalArgumentException("provided parameter object is not of type SmithWaterman3DParameters, but " + parameters.getClass().getName());

	params = (SmithWaterman3DParameters) parameters;
	AFPChain afpChain = new AFPChain(algorithmName);


	// covert input to sequences
	String seq1 = StructureTools.convertAtomsToSeq(ca1);
	String seq2 = StructureTools.convertAtomsToSeq(ca2);

	ProteinSequence s1 = null;
	ProteinSequence s2 = null;

	try {
		s1 = new ProteinSequence(seq1);
		s2 = new ProteinSequence(seq2);
	} catch (CompoundNotFoundException e){

		throw new StructureException(e.getMessage(),e);
	}

	// default blosum62
	SubstitutionMatrix<AminoAcidCompound> matrix = SubstitutionMatrixHelper.getBlosum65();

	GapPenalty penalty = new SimpleGapPenalty();

	penalty.setOpenPenalty(params.getGapOpen());
	penalty.setExtensionPenalty(params.getGapExtend());

	PairwiseSequenceAligner<ProteinSequence, AminoAcidCompound> smithWaterman =
		Alignments.getPairwiseAligner(s1, s2, PairwiseSequenceAlignerType.LOCAL, penalty, matrix);

	SequencePair<ProteinSequence, AminoAcidCompound> pair = smithWaterman.getPair();

	if (pair.getTarget().toString().isEmpty() || pair.getQuery().toString().isEmpty()) {
		throw new StructureException("Empty alignment for sequences "+s1+" and "+s2);
	}

	logger.debug("Smith-Waterman alignment is: "+pair.toString(100));

	// convert to a 3D alignment...
	afpChain = convert(ca1,ca2,pair, smithWaterman);

	// Perform an iterative dropping of the columns
	while (afpChain.getOptLength() > params.getMinLen()
			&& afpChain.getTotalRmsdOpt() > params.getMaxRmsd()) {
		afpChain = AlignmentTools.deleteHighestDistanceColumn(afpChain, ca1, ca2);
	}

	return afpChain;
}
 

/**
 * superimposing according to the optimized alignment
 *
 * @param afpChain
 * @param ca1
 * @param ca2
 * @return Group array twisted.
 * @throws StructureException
 */
public static Group[] twistOptimized(AFPChain afpChain, Atom[] ca1,
		Atom[] ca2) throws StructureException {

	Atom[] optTwistPdb = new Atom[ca2.length];

	int gPos = -1;
	for (Atom a : ca2) {
		gPos++;
		optTwistPdb[gPos] = a;
	}

	int blockNum = afpChain.getBlockNum();

	int b2 = 0;
	int e2 = 0;
	int focusResn = 0;
	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);
	}

	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	for (int bk = 0; bk < blockNum; bk++) {
		// THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING...
		// copies the atoms over to iniTwistPdb later on in modifyCod
		transformOrigPDB(optLen[bk], optAln[bk][0], optAln[bk][1], ca1,
				ca2, afpChain, bk);

		// transform pro2 according to comparison of pro1 and pro2 at give
		// residues
		if (bk > 0) {
			b2 = e2;
		}
		if (bk < blockNum - 1) { // bend at the middle of two consecutive
									// blocks
			e2 = optAln[bk][1][optLen[bk] - 1];
			e2 = (optAln[bk + 1][1][0] - e2) / 2 + e2;
		} else {
			e2 = ca2.length;
		}
		cloneAtomRange(optTwistPdb, ca2, b2, e2);
		for (int i = 0; i < optLen[bk]; i++) {
			focusRes1[focusResn] = optAln[bk][0][i];
			focusRes2[focusResn] = optAln[bk][1][i];
			focusResn++;
		}
	}
	int totalLenOpt = focusResn;
	logger.debug("calrmsdopt for {} residues", focusResn);
	double totalRmsdOpt = calCaRmsd(ca1, optTwistPdb, focusResn, focusRes1,
			focusRes2);
	logger.debug("got opt RMSD: {}", totalRmsdOpt);
	int optLength = afpChain.getOptLength();

	if (totalLenOpt != optLength) {
		logger.warn("Final alignment length is different {} {}",
				totalLenOpt, optLength);
	}
	logger.debug("final alignment length {}, rmsd {}", focusResn,
			totalRmsdOpt);

	afpChain.setTotalLenOpt(totalLenOpt);
	afpChain.setTotalRmsdOpt(totalRmsdOpt);

	return StructureTools.cloneGroups(optTwistPdb);

}