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

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

/**
 * Sometimes it's convenient to store an alignment using java collections,
 * where <tt>blocks.get(blockNum).get(0).get(pos)</tt> specifies the aligned
 * residue at position <i>pos</i> of block <i>blockNum</i> of the first
 * protein.
 *
 * This method takes such a collection and stores it into the afpChain's
 * {@link AFPChain#setOptAln(int[][][]) optAln}, setting the associated
 * length variables as well.
 *
 * @param afpChain
 * @param blocks
 */
private static void assignOptAln(AFPChain afpChain, List<List<List<Integer>>> blocks)
{

	int[][][] optAln = new int[blocks.size()][][];
	int[] optLen = new int[blocks.size()];
	int optLength = 0;
	int numBlocks = blocks.size();

	for(int block = 0; block < numBlocks; block++) {
		// block should be 2xN rectangular
		assert(blocks.get(block).size() == 2);
		assert( blocks.get(block).get(0).size() == blocks.get(block).get(1).size());

		optLen[block] = blocks.get(block).get(0).size();
		optLength+=optLen[block];

		optAln[block] = new int[][] {
				new int[optLen[block]],
				new int[optLen[block]]
		};
		for(int pos = 0; pos < optLen[block]; pos++) {
			optAln[block][0][pos] = blocks.get(block).get(0).get(pos);
			optAln[block][1][pos] = blocks.get(block).get(1).get(pos);
		}
	}


	afpChain.setBlockNum(numBlocks);
	afpChain.setOptAln(optAln);
	afpChain.setOptLen(optLen);
	afpChain.setOptLength(optLength);

	// TODO I don't know what these do. Should they be set?
	//afpChain.setBlockSize(blockSize);
	//afpChain.setBlockResList(blockResList);
	//afpChain.setChainLen(chainLen);

}
 

/**
 * Creates a minimal AFPChain from the specified alignment and proteins
 * @param dupAlign
 * @param ca1
 * @param ca2
 * @return
 */
private AFPChain makeDummyAFPChain(int[][][] dupAlign, Atom[] ca1,Atom[] ca2) {
	AFPChain afp = new AFPChain(AFPChain.UNKNOWN_ALGORITHM);
	afp.setOptAln(dupAlign);
	afp.setOptLength(dupAlign[0][1].length);
	afp.setCa1Length(ca1.length);
	afp.setCa2Length(ca2.length);
	afp.setBlockNum(1);
	afp.setOptLen(new int[] {dupAlign[0][1].length});
	return afp;
}
 

/**
 * in some special cases, there is no maginificent twists in the
final chaining result; however, their rmsd (original and after
optimizing) are very large. Therefore, a post-process is taken
to split the blocks further at the ralative bad connections (
with relative high distance variation)
to be tested:
  split or not according to optimized or initial chaining???
 */

private static void splitBlock(FatCatParameters params, AFPChain afpChain, Atom[] ca1, Atom[] ca2)
{
	if ( debug)
		System.err.println("AFPPostProcessor: splitBlock");
	int     i, a, bk, cut;
	double  maxs, maxt;
	int blockNum = afpChain.getBlockNum();
	int maxTra = params.getMaxTra();
	double badRmsd = params.getBadRmsd();

	int     blockNum0 = blockNum;

	double[] blockRmsd = afpChain.getBlockRmsd();
	int[] blockSize = afpChain.getBlockSize();
	int[] block2Afp = afpChain.getBlock2Afp();
	double[] afpChainTwiList = afpChain.getAfpChainTwiList();

	bk = 0;
	while(blockNum < maxTra + 1)    {
		maxs = 0;
		for(i = 0; i < blockNum; i ++)   {
			if(blockRmsd[i] > maxs && blockSize[i] > 2) { //according to the optimized alignment
				maxs = blockRmsd[i];
				bk = i;
			} //!(Note: optRmsd, not blockRmsd, according to the optimized alignment
		}
		if(maxs < badRmsd)      break;
		maxt = 0;
		cut = 0;
		for(i = 1; i < blockSize[bk]; i ++)     {
			a = i + block2Afp[bk];
			if(afpChainTwiList[a] > maxt)   {
				maxt = afpChainTwiList[a];
				cut = i;

			}
		}
		if(debug)
			System.out.println(String.format("block %d original size %d rmsd %.3f maxt %.2f cut at %d\n", bk, blockSize[bk], maxs, maxt, cut));
		for(i = blockNum - 1; i > bk; i --)     {
			block2Afp[i + 1] = block2Afp[i];
			blockSize[i + 1] = blockSize[i];
			blockRmsd[i + 1] = blockRmsd[i];
		} //update block information
		block2Afp[bk + 1] = cut + block2Afp[bk];
		blockSize[bk + 1] = blockSize[bk] - cut;
		blockSize[bk] = cut;

		if(debug)
			System.out.println(String.format("  split into %d and %d sizes\n", blockSize[bk], blockSize[bk + 1]));


		int[] afpChainList = afpChain.getAfpChainList();
		//int[] subrange1    = getSubrange(afpChainList, block2Afp[bk + 1] );
		blockRmsd[bk + 1]  = AFPChainer.calAfpRmsd(blockSize[bk + 1],  afpChainList, block2Afp[bk + 1] , afpChain, ca1, ca2);

		//int[] subrange2    = getSubrange(afpChainList, block2Afp[bk] );
		blockRmsd[bk]      = AFPChainer.calAfpRmsd(blockSize[bk],      afpChainList, block2Afp[bk], afpChain, ca1, ca2);

		//split a block at the biggest position
		blockNum ++;
		afpChain.setAfpChainList(afpChainList);
	}
	if(blockNum - blockNum0 > 0)    {
		if(debug)
			System.out.println(String.format("Split %d times:\n", blockNum - blockNum0));
		for(i = 0; i < blockNum; i ++)  {
			if(debug)
				System.out.println(String.format("  block %d size %d from %d rmsd %.3f\n", i, blockSize[i], block2Afp[i], blockRmsd[i]));
		}
	}


	afpChain.setBlockNum(blockNum);
	afpChain.setBlockSize(blockSize);
	afpChain.setBlockRmsd(blockRmsd);
	afpChain.setBlock2Afp(block2Afp);


}
 

/**
 * remove the artifical small rigid-body superimpose in the middle
 clust the similar superimpositions (caused by the small flexible
 region, which is detected as a seperate rigid superimposing region by adding
 two twists before and after it(artifically!)
 one possible solution: allowing long enough loops in the chaining process,
 which however increase the calculation complexity
 */
private static void deleteBlock(FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[] ca2)
{
	int blockNum = afpChain.getBlockNum();
	List<AFP> afpSet = afpChain.getAfpSet();

	int[] afpChainList =afpChain.getAfpChainList();



	int[] block2Afp = afpChain.getBlock2Afp();
	int[] blockSize = afpChain.getBlockSize();

	double[] blockRmsd = afpChain.getBlockRmsd();

	int fragLen = params.getFragLen();

	//remove those blocks (both in terminals and in the middle) with only a AFP
	//but still keep those small blocks spaning large regions
	if(blockNum <= 1)       return;
	int     blockNumOld = blockNum;
	int     i, j, b1, b2, e1, e2, len;
	e1 = e2 = 0;
	for(i = 0; i < blockNum; i ++) {
		b1 = e1;
		b2 = e2;
		if(i < blockNum - 1)    {
			e1 = afpSet.get(afpChainList[block2Afp[i + 1]]).getP1();
			e2 = afpSet.get(afpChainList[block2Afp[i + 1]]).getP2();
		}
		else    {
			e1 = ca1.length;
			e2 = ca2.length;
		}
		if(blockSize[i] > 1)    continue;
		len = (e1 - b1) < (e2 - b2)?(e1 - b1):(e2 - b2);
		//if(i == blockNum - 1) blockNum --;
		if(len < 2 * fragLen)   {
			for(j = i; j < blockNum - 1; j ++)      {
				blockRmsd[j] = blockRmsd[j + 1];
				blockSize[j] = blockSize[j + 1];
				block2Afp[j] = block2Afp[j + 1];
			}
			blockNum --;
			i --;
		} //delete a block
	}
	if(blockNumOld > blockNum)
		if(debug)
			System.out.println(
					String.format("Delete %d small blocks\n", blockNumOld - blockNum)
			);


	if (debug)
		System.err.println("deleteBlock: end blockNum:"+ blockNum);
	afpChain.setBlock2Afp(block2Afp);
	afpChain.setBlockSize(blockSize);
	afpChain.setAfpChainList(afpChainList);
	afpChain.setBlockNum(blockNum);
	afpChain.setBlockRmsd(blockRmsd);
}
 

/**
 * Merge consecutive blocks with similar transformation
 */
private static  void mergeBlock(FatCatParameters params, AFPChain afpChain,Atom[] ca1,Atom[] ca2)
{

	int blockNum = afpChain.getBlockNum();
	double badRmsd = params.getBadRmsd();

	int[] block2Afp = afpChain.getBlock2Afp();
	int[] blockSize = afpChain.getBlockSize();

	double[] blockRmsd = afpChain.getBlockRmsd();

	int afpChainTwiNum = afpChain.getAfpChainTwiNum();

	//clustering the neighbor blocks if their transformations are similar
	int     i, j, b1, b2, minb1, minb2;
	double  minrmsd;
	int     merge = 0;
	int     blockNumOld = blockNum;
	double[][]  rmsdlist = null;
	if(blockNum > 1)        {
		rmsdlist = new double[blockNumOld][blockNumOld];
		for(b1 = 0; b1 < blockNum - 1; b1 ++)   {
			for(b2 = b1 + 1; b2 < blockNum; b2 ++)  {
				rmsdlist[b1][b2] = combineRmsd(b1, b2,afpChain,ca1,ca2);
			}
		}
	}
	minb1 = 0;
	while(blockNum > 1)     {
		minrmsd = 1000;
		for(i = 0; i < blockNum - 1; i ++)      {
			j = i + 1; //only consider neighbor blocks
			if(minrmsd > rmsdlist[i][j])    {
				minrmsd = rmsdlist[i][j];
				minb1 = i;
			}
		}
		minb2 = minb1 + 1; //merge those most similar blocks
		//maxrmsd = (blockRmsd[minb1] > blockRmsd[minb2])?blockRmsd[minb1]:blockRmsd[minb2];
		if(minrmsd < badRmsd)   {
			if(debug)
				System.out.println(String.format("merge block %d (rmsd %.3f) and %d (rmsd %.3f), total rmsd %.3f\n",
						minb1, blockRmsd[minb1], minb2, blockRmsd[minb2], minrmsd));
			blockSize[minb1] += blockSize[minb2];
			blockRmsd[minb1] = minrmsd;
			for(i = minb2; i < blockNum - 1; i ++)  {
				block2Afp[i] = block2Afp[i + 1];
				blockSize[i] = blockSize[i + 1];
				blockRmsd[i] = blockRmsd[i + 1];
			} //update block information
			afpChainTwiNum --;
			blockNum --;
			for(b1 = 0; b1 < blockNum - 1; b1 ++)   {
				for(b2 = b1 + 1; b2 < blockNum; b2 ++) {
					if(b1 == minb1 || b2 == minb1)  {
						rmsdlist[b1][b2] = combineRmsd(b1, b2, afpChain,ca1,ca2);
					}
					else if(b2 < minb1)     continue;
					else if(b1 < minb1)     {
						rmsdlist[b1][b2] = rmsdlist[b1][b2 + 1];
					}
					else    {
						rmsdlist[b1][b2] = rmsdlist[b1 + 1][b2 + 1];
					}
				}
			} //update the rmsdlist
			merge ++;
		} //merge two blocks
		else if(minrmsd >= 100) break;
		else    {
			rmsdlist[minb1][minb2] += 100;
		} //not merge, modify the rmsdlist so that this combination is not considered in next iteration
	}

	if(merge > 0)       {
		if(debug)
			System.out.println(String.format("Merge %d blocks, remaining %d blocks\n", merge, blockNum));
	}

	if (debug){
		System.err.println("AFPPostProcessor: mergeBlock end blocknum:" + blockNum);
	}
	afpChain.setBlock2Afp(block2Afp);
	afpChain.setBlockSize(blockSize);
	afpChain.setBlockNum(blockNum);
	afpChain.setBlockRmsd(blockRmsd);
	afpChain.setAfpChainTwiNum(afpChainTwiNum);
}
 

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

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

@Test
public void testAFPconversion() throws Exception{

	//Fill an AFPChain with the general information
	AFPChain afp = new AFPChain("algorithm");
	afp.setName1("name1");
	afp.setName2("name2");
	afp.setVersion("1.0");
	afp.setCalculationTime(System.currentTimeMillis());
	//Generate a optimal alignment with 3 blocks and 5 residues per block
	int[][][] optAln = new int[3][][];
	for (int b=0; b<optAln.length; b++){
		int[][] block = new int[2][];
		for (int c=0; c<block.length; c++){
			int[] residues = {b+5,b+6,b+7,b+8,b+9};
			block[c] = residues;
		}
		optAln[b] = block;
	}
	afp.setOptAln(optAln);
	afp.setBlockNum(optAln.length);
	//Set the rotation matrix and shift to random numbers
	double[][] mat = {{0.13,1.5,0.84},{1.3,0.44,2.3},{1.0,1.2,2.03}};
	Matrix rot = new Matrix(mat);
	Atom shift = new AtomImpl();
	shift.setX(0.44);
	shift.setY(0.21);
	shift.setZ(0.89);
	Matrix[] blockRot = {rot,rot,rot};
	afp.setBlockRotationMatrix(blockRot);
	Atom[] blockShift = {shift,shift,shift};
	afp.setBlockShiftVector(blockShift);

	//Convert the AFPChain into a MultipleAlignment (without Atoms)
	MultipleAlignmentEnsemble ensemble =
			new MultipleAlignmentEnsembleImpl(afp,null,null,true);
	MultipleAlignment msa = ensemble.getMultipleAlignment(0);

	//Test for all the information
	assertEquals(afp.getName1(),ensemble.getStructureIdentifiers().get(0).getIdentifier());
	assertEquals(afp.getName2(), ensemble.getStructureIdentifiers().get(1).getIdentifier());
	assertEquals(afp.getAlgorithmName(), ensemble.getAlgorithmName());
	assertEquals(afp.getVersion(),ensemble.getVersion());
	assertTrue(ensemble.getCalculationTime().equals(
			afp.getCalculationTime()));
	assertEquals(afp.getBlockNum(), msa.getBlockSets().size());
	for (int b = 0; b<afp.getBlockNum(); b++){
		assertEquals(Calc.getTransformation(
				afp.getBlockRotationMatrix()[b],
				afp.getBlockShiftVector()[b]),
				msa.getBlockSet(b).getTransformations().get(1));
	}

	//Test for the scores
	assertEquals(msa.getScore(MultipleAlignmentScorer.CE_SCORE),
			(Double) afp.getAlignScore());
	assertEquals(msa.getScore(MultipleAlignmentScorer.AVGTM_SCORE),
			(Double) afp.getTMScore());
	assertEquals(msa.getScore(MultipleAlignmentScorer.RMSD),
			(Double) afp.getTotalRmsdOpt());


	//Test for the optimal alignment
	for (int b=0; b<3; b++){
		for (int c=0; c<2; c++){
			for (int res=0; res<5; res++){
				Integer afpRes = afp.getOptAln()[b][c][res];
				assertEquals(afpRes, msa.getBlock(b).
						getAlignRes().get(c).get(res));
			}
		}
	}
}