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

/**
 * get the afp list and residue list for each block
 */

public static void blockInfo(AFPChain afpChain)
{
	int     i, j, k, a, n;

	int blockNum = afpChain.getBlockNum();

	int[] blockSize =afpChain.getBlockSize();
	int[] afpChainList = afpChain.getAfpChainList();
	int[] block2Afp = afpChain.getBlock2Afp();
	int[][][]blockResList = afpChain.getBlockResList();

	List<AFP>afpSet = afpChain.getAfpSet();
	int[] blockResSize = afpChain.getBlockResSize();

	for(i = 0; i < blockNum; i ++)  {
		n = 0;
		for(j = 0; j < blockSize[i]; j ++)      {
			//the index in afpChainList, not in the whole afp set
			a = afpChainList[block2Afp[i] + j];
			for(k = 0; k < afpSet.get(a).getFragLen(); k ++)     {
				blockResList[i][0][n] = afpSet.get(a).getP1() + k;
				blockResList[i][1][n] = afpSet.get(a).getP2() + k;
				n ++;
			}
		}
		blockResSize[i] = n;
	}

	afpChain.setBlockResSize(blockResSize);
	afpChain.setBlockSize(blockSize);
	afpChain.setAfpChainList(afpChainList);
	afpChain.setBlock2Afp(block2Afp);
	afpChain.setBlockResList(blockResList);
}
 

/**
	//return the root mean square of the distance matrix between the residues
	//from the segments that form the given AFP list
	//this value can be a measurement (2) for the connectivity of the AFPs
	//and its calculation is quicker than the measurement (1), rmsd
	//currently only deal with AFP pair
//
//           |-d1--|
//          |--d2---|
//         |---d3----|
	//-----------------------------------------------------------------------
	//this module is optimized
	 *
	 * @param afp1
	 * @param afp2
	 * @return
	 */
	private static double calAfpDis(int afp1, int afp2, FatCatParameters params, AFPChain afpChain)
	{

		List<AFP> afpSet = afpChain.getAfpSet();

		Matrix disTable1 = afpChain.getDisTable1();
		Matrix disTable2 = afpChain.getDisTable2();

		int fragLen = params.getFragLen();
		double afpDisCut = params.getAfpDisCut();
		double disCut = params.getDisCut();
		double fragLenSq = params.getFragLenSq();

		int     i, j, ai, bi, aj, bj;
		double  d;
		double  rms = 0;
		for(i = 0; i < fragLen; i ++)   {
			ai = afpSet.get(afp1).getP1() + i;
			bi = afpSet.get(afp1).getP2() + i;
			for(j = 0; j < fragLen; j ++)   {
				aj = afpSet.get(afp2).getP1() + j;
				bj = afpSet.get(afp2).getP2() + j;
				d = disTable1.get(aj,ai) - disTable2.get(bj,bi);
				rms += d * d;
				if(rms > afpDisCut)     { return (disCut); }
			}
		}
		return (Math.sqrt(rms / fragLenSq));
	}
 

/**
 * Set the list of equivalent residues in the two proteins given a list of
 * AFPs
 *
 * WARNING: changes the values for FocusRes1, focusRes2 and FocusResn in
 * afpChain!
 *
 * @param afpChain
 *            the AFPChain to store resuts
 * @param afpn
 *            nr of afp
 * @param afpPositions
 * @param listStart
 * @return nr of eq residues
 */

public static int afp2Res(AFPChain afpChain, int afpn, int[] afpPositions,
		int listStart) {
	int[] res1 = afpChain.getFocusRes1();
	int[] res2 = afpChain.getFocusRes2();
	int minLen = afpChain.getMinLen();

	int n = 0;

	List<AFP> afpSet = afpChain.getAfpSet();

	for (int i = listStart; i < listStart + afpn; i++) {
		int a = afpPositions[i];
		for (int j = 0; j < afpSet.get(a).getFragLen(); j++) {
			if (n >= minLen) {
				throw new RuntimeException(
						"Error: too many residues in AFPChainer.afp2Res!");
			}
			res1[n] = afpSet.get(a).getP1() + j;
			res2[n] = afpSet.get(a).getP2() + j;
			n++;
		}
	}

	afpChain.setFocusRes1(res1);
	afpChain.setFocusRes2(res2);
	afpChain.setFocusResn(n);

	if (n == 0) {
		logger.warn("n=0!!! + " + afpn + " " + listStart + " "
				+ afpPositions.length);
	}
	return n;
}
 

/**
 * to update the chaining score after block delete and merge processed
 * the blockScore value is important for significance evaluation
 */
public static void updateScore(FatCatParameters params, AFPChain afpChain)
{
	int     i, j, bknow, bkold, g1, g2;


	afpChain.setConn(0d);
	afpChain.setDVar(0d);

	int blockNum = afpChain.getBlockNum();
	int alignScoreUpdate = 0;
	double[] blockScore = afpChain.getBlockScore();
	int[] blockGap = afpChain.getBlockGap();
	int[] blockSize =afpChain.getBlockSize();
	int[] afpChainList = afpChain.getAfpChainList();
	List<AFP>afpSet = afpChain.getAfpSet();
	int[] block2Afp = afpChain.getBlock2Afp();

	double torsionPenalty = params.getTorsionPenalty();


	bkold = 0;
	for(i = 0; i < blockNum; i ++)  {
		blockScore[i] = 0;
		blockGap[i] = 0;
		for(j = 0; j < blockSize[i]; j ++)      {
			bknow = afpChainList[block2Afp[i] + j];
			if(j == 0)      {
				blockScore[i] = afpSet.get(bknow).getScore();
			}
			else    {
				AFPChainer.afpPairConn(bkold, bknow, params, afpChain); //note: j, i
				Double conn = afpChain.getConn();
				blockScore[i] += afpSet.get(bknow).getScore() + conn;
				g1 = afpSet.get(bknow).getP1() - afpSet.get(bkold).getP1() - afpSet.get(bkold).getFragLen();
				g2 = afpSet.get(bknow).getP2() - afpSet.get(bkold).getP2() - afpSet.get(bkold).getFragLen();
				blockGap[i] += (g1 > g2)?g1:g2;
			}
			bkold = bknow;
		}
		alignScoreUpdate += blockScore[i];
	}
	if(blockNum >= 2)       {
		alignScoreUpdate += (blockNum - 1) * torsionPenalty;
	}

	afpChain.setBlockGap(blockGap);
	afpChain.setAlignScoreUpdate(alignScoreUpdate);
	afpChain.setBlockScore(blockScore);
	afpChain.setBlockSize(blockSize);
	afpChain.setAfpChainList(afpChainList);
	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);
}
 

private  static int getCompatibleAfps(int afp, int[] list, FatCatParameters params, AFPChain afpChain){

		int     i, j, i1, j1, f, G, c, a1, a2, a3, b1, b2, b3, s1, s2;

		int fragLen = params.getFragLen();
		int maxGapFrag = params.getMaxGapFrag();
		int misCut = params.getMisCut();
		int maxTra = params.getMaxTra();
		List<AFP> afpSet = afpChain.getAfpSet();

		f = fragLen;
		G = maxGapFrag;
		c = misCut;

		i1 = afpSet.get(afp).getP1();
		j1 = afpSet.get(afp).getP2();
		a3 = i1 - f;
		a2 = a3 - c;
		a1 = i1 - G;
		a2 = a2>0?a2:0;
		a1 = a1>0?a1:0;

		b3 = j1 - f;
		b2 = b3 - c;
		b1 = j1 - G;
		b2 = (b2 > 0)?b2:0;
		b1 = (b1 > 0)?b1:0;

		int[][] afpAftIndex = afpChain.getAfpAftIndex();
		int[][] afpBefIndex = afpChain.getAfpBefIndex();
		int[] twi = afpChain.getTwi();



		int     n = 0;
		//compatible region 1-2, [a1,a3][b2,b3]
				for(i = a1; i <= a3; i ++)      {//i <= a3 instead of i < a3
					s1 = afpAftIndex[i][b2]; //note afpAftIndex, not afpIndex
					if(s1 < 0)      continue;//no AFP for the given i with j > b2
					s2 = afpBefIndex[i][b3]; //afps is sorted by j given a i,it's sparse matrix
					if(s2 < 0)      continue;//no AFP for the given i with j < b3
					for(j = s1; j <= s2; j ++)      { //j <= s2 instead of j < s2
						if(twi[j] <= maxTra)    {
							list[n ++] = j;
						}
					}
				}

				//compatible region 3  [a2,a3][b1,b2]
				for(i = a2; i <= a3; i ++)      {
					s1 = afpAftIndex[i][b1];
					if(s1 < 0)      continue;
					s2 = afpBefIndex[i][b2]; //afps is sorted by j given a i
					if(s2 < 0)      continue;
					//note j < s2, as the cases of j == s2 is alread considered in previous region
					for(j = s1; j < s2; j ++)       {
						if(twi[j] <= maxTra)    {
							list[n ++] = j;
						}
					}
				}

				return n;

	}
 

/**
//Key function: calculate the connectivity of AFP pairs
//no compatibility criteria is executed
//note: afp1 is previous to afp2 in terms of the position
	 //this module must be optimized
 *
 * @param afp1
 * @param afp2
 * @return flag if they are connected
 */
public static boolean afpPairConn(int afp1, int afp2,  FatCatParameters params, AFPChain afpChain)

{

	Double conn = afpChain.getConn();
	Double dvar = afpChain.getDVar();

	double misScore = params.getMisScore();
	double maxPenalty = params.getMaxPenalty();
	double disCut = params.getDisCut();
	double gapExtend = params.getGapExtend();
	double torsionPenalty = params.getTorsionPenalty();
	double disSmooth = params.getDisSmooth();

	List<AFP> afpSet = afpChain.getAfpSet();

	int     m = calcGap(afpSet.get(afp2),afpSet.get(afp1));
	int     g = calcMismatch(afpSet.get(afp2),afpSet.get(afp1));


	double  gp = misScore * m;      //on average, penalty for a mismatch is misScore, no modification on score
	if(g > 0)       {
		gp += gapExtend * g;
	}
	if(gp < maxPenalty)     gp = maxPenalty; //penalty cut-off
	//note: use < (smaller) instead of >, because maxPenalty is a negative number

	double  d;
	d = calAfpDis(afp1, afp2,params, afpChain);
	//note: the 'dis' value is numerically equivalent to the 'rms' with exceptions

	boolean     ch = false;
	double  tp = 0.0;
	if(d >= disCut) {
		tp = torsionPenalty;
		ch = true;
	} //use the variation of the distances between AFPs
	else  if(d > disCut - disSmooth)        {
		double  wt = Math.sqrt((d - disCut + disSmooth) / disSmooth);
		//using sqrt: penalty increase with dis more quicker than linear function
		tp = torsionPenalty * wt;
	}

	dvar = d;
	conn = tp + gp;

	afpChain.setConn(conn);
	afpChain.setDVar(dvar);
	return ch;
}
 

/**
 * run AFP chaining allowing up to maxTra flexible regions.
 * Input is original coordinates.
 *
 */
private static Group[] chainAfp(FatCatParameters params,AFPChain afpChain, Atom[] ca1, Atom[] ca2) throws StructureException{

	// we don;t want to rotate input atoms, do we?
	Atom[] ca2clone = StructureTools.cloneAtomArray(ca2);

	List<AFP> afpSet = afpChain.getAfpSet();

	if (debug)
		System.out.println("entering chainAfp");
	int afpNum = afpSet.size();

	if ( afpNum < 1)
		return new Group[0];

	long bgtime = System.currentTimeMillis();
	if(debug)    {
		System.out.println(String.format("total AFP %d\n", afpNum));
	}

	//run AFP chaining

	AFPChainer.doChainAfp(params,afpChain ,ca1,ca2);

	int afpChainLen = afpChain.getAfpChainLen();

	if(afpChainLen < 1)     {

		afpChain.setShortAlign(true);
		return new Group[0];
	} //very short alignment

	long chaintime = System.currentTimeMillis();
	if(debug)    {

		System.out.println("Afp chaining: time " + (chaintime-bgtime));
	}

	// do post processing

	AFPPostProcessor.postProcess(params, afpChain,ca1,ca2);

	// Optimize the final alignment

	AFPOptimizer.optimizeAln(params, afpChain,ca1,ca2);

	AFPOptimizer.blockInfo( afpChain);

	AFPOptimizer.updateScore(params,afpChain);

	AFPAlignmentDisplay.getAlign(afpChain,ca1,ca2);

	Group[] twistedPDB = AFPTwister.twistPDB(afpChain, ca1, ca2clone);

	SigEva sig =  new SigEva();
	double probability = sig.calSigAll(params, afpChain);
	afpChain.setProbability(probability);
	double normAlignScore = sig.calNS(params,afpChain);
	afpChain.setNormAlignScore(normAlignScore);
	double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2,false);
	afpChain.setTMScore(tmScore);

	/*

	SIGEVA  sig;
	probability = sig.calSigAll(maxTra, sparse, pro1Len, pro2Len, alignScore, totalRmsdOpt, optLength, blockNum - 1);
	normAlignScore = sig.calNS(pro1Len, pro2Len, alignScore, totalRmsdOpt, optLength, blockNum - 1);

	 */

	//if(maxTra == 0)       probability = sig.calSigRigid(pro1Len, pro2Len, alignScore, totalRmsdOpt, optLength);
	//else  probability = sig.calSigFlexi(pro1Len, pro2Len, alignScore, totalRmsdOpt, optLength, blockNum - 1);

	if(debug)    {

		long nowtime = System.currentTimeMillis();
		long diff = nowtime - chaintime;
		System.out.println("Alignment optimization: time "+ diff);

		System.out.println("score:      " + afpChain.getAlignScore());
		System.out.println("opt length: " + afpChain.getOptLength());
		System.out.println("opt rmsd:   "+ afpChain.getTotalRmsdOpt());

	}
	return twistedPDB;

}