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

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

/**
 * Swaps the order of structures in an AFPChain
 * @param a
 * @return
 */
public AFPChain invertAlignment(AFPChain a) {
	String name1 = a.getName1();
	String name2 = a.getName2();
	a.setName1(name2);
	a.setName2(name1);

	int len1 = a.getCa1Length();
	a.setCa1Length( a.getCa2Length() );
	a.setCa2Length( len1 );

	int beg1 = a.getAlnbeg1();
	a.setAlnbeg1(a.getAlnbeg2());
	a.setAlnbeg2(beg1);

	char[] alnseq1 = a.getAlnseq1();
	a.setAlnseq1(a.getAlnseq2());
	a.setAlnseq2(alnseq1);

	Matrix distab1 = a.getDisTable1();
	a.setDisTable1(a.getDisTable2());
	a.setDisTable2(distab1);

	int[] focusRes1 = a.getFocusRes1();
	a.setFocusRes1(a.getFocusRes2());
	a.setFocusRes2(focusRes1);

	//What are aftIndex and befIndex used for? How are they indexed?
	//a.getAfpAftIndex()


	String[][][] pdbAln = a.getPdbAln();
	if( pdbAln != null) {
		for(int block = 0; block < a.getBlockNum(); block++) {
			String[] paln1 = pdbAln[block][0];
			pdbAln[block][0] = pdbAln[block][1];
			pdbAln[block][1] = paln1;
		}
	}

	int[][][] optAln = a.getOptAln();
	if( optAln != null ) {
		for(int block = 0; block < a.getBlockNum(); block++) {
			int[] aln1 = optAln[block][0];
			optAln[block][0] = optAln[block][1];
			optAln[block][1] = aln1;
		}
	}
	a.setOptAln(optAln); // triggers invalidate()

	Matrix distmat = a.getDistanceMatrix();
	if(distmat != null)
		a.setDistanceMatrix(distmat.transpose());


	// invert the rotation matrices
	Matrix[] blockRotMat = a.getBlockRotationMatrix();
	Atom[] shiftVec = a.getBlockShiftVector();
	if( blockRotMat != null) {
		for(int block = 0; block < a.getBlockNum(); block++) {
			if(blockRotMat[block] != null) {
				// if y=x*A+b, then x=y*inv(A)-b*inv(A)
				blockRotMat[block] = blockRotMat[block].inverse();

				Calc.rotate(shiftVec[block],blockRotMat[block]);
				shiftVec[block] = Calc.invert(shiftVec[block]);
			}
		}
	}

	return a;
}
 

public static final  void extractAFPChains(FatCatParameters params, AFPChain afpChain,Atom[] ca1,Atom[] ca2) throws StructureException {



		List<AFP> afpSet = new ArrayList<AFP>();
		afpChain.setAfpSet(afpSet);

		if ( debug )
			System.err.println("nr of atoms ca1: " + ca1.length + " ca2: " +  ca2.length);



		int     p1, p2;
		@SuppressWarnings("unused")
		int n0, n, n1, n2;
		double  filter1;
		double rmsd = 0;

		Matrix r = new Matrix(3,3);
		Atom   t = new AtomImpl();


		int sparse = params.getSparse();
		int maxTra = params.getMaxTra();
		int fragLen = params.getFragLen();
		double disFilter = params.getDisFilter();
		double rmsdCut = params.getRmsdCut();
		double badRmsd = params.getBadRmsd();
		double fragScore = params.getFragScore();

		int     add = sparse + 1; //if add > 1, use sparse sampling
		n0 = n = n1 = n2 = 0;

		int minLen = 0;

		int prot1Length = ca1.length;
		int prot2Length = ca2.length;

		if(prot1Length < prot2Length)
			minLen = prot1Length;
		else
			minLen = prot2Length;
		afpChain.setMinLen(minLen);

		afpChain.setBlockResList(new int[maxTra+1][2][minLen]);
		afpChain.setFocusRes1(new int[minLen]);
		afpChain.setFocusRes2(new int[minLen]);

		for(p1 = 0; p1 < prot1Length - fragLen; p1 += add )    {
			for(p2 = 0; p2 < prot2Length - fragLen; p2 += add)     {
				n0 ++;
				filter1 = getEnd2EndDistance(ca1, ca2, p1, p1 + fragLen - 1, p2, p2 + fragLen - 1);
				//difference bewteen end-to-end distances
				if(filter1 > disFilter) { n1 ++; continue; }
				boolean filter2 = filterTerminal(ca1,ca2, p1, p1 + fragLen - 1, p2, p2 + fragLen - 1, fragLen, minLen);
				if(filter2)     {
					n2 ++;
					continue;

				} //be cautious to use this filter !!

				// here FATCAT does a a jacobi transformation
				//rmsd = kearsay(fragLen, ca1[p1], ca2[p2], r, t);
				// we use the BioJava SVD instead...

				//
				rmsd = getRmsd(ca1,ca2,fragLen, p1,p2,r,t);

				if(rmsd < rmsdCut)      {
					AFP     afptmp = new AFP();
					afptmp.setP1(p1);
					afptmp.setP2(p2);
					afptmp.setFragLen(fragLen);
					afptmp.setRmsd(rmsd);
					afptmp.setM(r);
					afptmp.setT(t.getCoords());
					afptmp.setScore(scoreAfp(afptmp,badRmsd,fragScore));
					afpSet.add(afptmp);
					n ++;
				}
			}
		}

		int afpNum = afpSet.size();

		if(debug) {
			String msg = String.format("possible AFP-pairs %d, remain %d after filter 1 remove %d; filter 2 remove %d\n",
					n0, afpNum, n1, n2);
			System.err.println(msg);
		}


	}
 

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

}