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

/**
 * Tests that {@link AlignmentTools#updateSuperposition(AFPChain, Atom[], Atom[])} calculates the correct RMSD and TM-score for an AFPChain of 1 block.
 * TODO: Write a test with 2 blocks
 */
@Test
public void testUpdateSuperposition() throws IOException, StructureException {
	Structure s = StructureTools.getStructure("31BI");
	Atom[] ca1 = StructureTools.getRepresentativeAtomArray(s);
	Atom[] ca2 = StructureTools.getRepresentativeAtomArray(s);
	StringBuilder sb = new StringBuilder();
	BufferedReader br = new BufferedReader(new FileReader("src/test/resources/align/31BI_symm_align.xml"));
	String line = "";
	while ((line = br.readLine()) != null) {
		sb.append(line);
	}
	br.close();
	AFPChain afpChain = AFPChainXMLParser.fromXML(sb.toString(), ca1, ca2);
	afpChain.setTMScore(-1);
	afpChain.setTotalRmsdOpt(-1);
	AlignmentTools.updateSuperposition(afpChain, ca1, ca2);
	Assert.assertEquals("TM-score is wrong", 0.62779, afpChain.getTMScore(), 0.001);
	Assert.assertEquals("RMSD is wrong", 2.50569, afpChain.getTotalRmsdOpt(), 0.001);
}
 

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

}