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

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

/**
	 * Calculates a structural alignment and returns alignment metrics.
	 * 
	 * @param alignmentAlgorithm name of the algorithm
	 * @param key unique identifier for protein chain pair
	 * @param points1 C-alpha positions of chain 1
	 * @param points2 C-alpha positions of chain 2
	 * @return
	 */
	public static List<Row> getAlignment(String alignmentAlgorithm, String key, Point3d[] points1, Point3d[] points2) {
		// create input for BioJava alignment method
		Atom[] ca1 = getCAAtoms(points1);
		Atom[] ca2 = getCAAtoms(points2);
		
		// calculate the alignment
		AFPChain afp = null;
		try {
			StructureAlignment algorithm  = StructureAlignmentFactory.getAlgorithm(alignmentAlgorithm);
			afp = algorithm.align(ca1,ca2);
			double tmScore = AFPChainScorer.getTMScore(afp, ca1, ca2);
			afp.setTMScore(tmScore);
		} catch (StructureException e) {
			e.printStackTrace();
			return Collections.emptyList();
		} 
		
		// TODO add alignments as arrays to results
//		int[][] alignment = afp.getAfpIndex();
//		for (int i = 0; i < alignment.length; i++) {
//			System.out.println(alignment[i][0] + " - " + alignment[i][1]);
//		}

		// record the alignment metrics
		Row row = RowFactory.create(key, afp.getOptLength(), afp.getCoverage1(), 
				afp.getCoverage2(), (float) afp.getTotalRmsdOpt(), (float) afp.getTMScore());

		return Collections.singletonList(row);
	}
 

public void nextStep( AFPChain afpChain,
		Atom[] ca1, Atom[] ca2) throws StructureException{


	if(nBestTrace>0) {
		checkBestTraces(afpChain,ca1,ca2);
	} else {
		noBestTrace();
	}

	convertAfpChain(afpChain, ca1, ca2);
	AFPAlignmentDisplay.getAlign(afpChain, ca1, ca2);
	double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2,false);
	afpChain.setTMScore(tmScore);
}
 

public static void main(String[] args) throws IOException, StructureException, Exception {
	Structure structure1 = StructureTools.getStructure("1w0p");
	Structure structure2 = StructureTools.getStructure("1w0p");
	String first = "alfdynatgdtefdspakqgwmqdntnngsgvltnadgmpawlvqgiggraqwtyslstnqhaqassfgwrmttemkvlsggmitnyyangtqrvlpiisldssgnlvvefegqtgrtvlatgtaateyhkfelvflpgsnpsasfyfdgklirdniqptaskQNMIVWGNGSSntdgvaayrdikfei------------------------------------------------------------------------------------------------------------------QGDVIf------------RGPDRIPSIVASsvTPGVVTAFAEKRVGGgdpgalsntNDIITRTSRDGGITWDTELNLTEQinvsdeFDFSDPRPIYDPs---SNTVLVSYARWPtdaaqngdrikpwmpNGIFYSVYDVASgnWQAPIDVTdqvkersfqiagwggselyrrntslnsqqdwqsnakirivdgaanqiqvadgsrkyvvtlsidesgglvanlngvsapiilqsehakvhsfhdyelqysalnhtttlfvdgqqittwagevsqenniqfgnadaqidgrlhvqkivltqqghnlvefdafylaqqtpevekdleklgwtkiktgntmslygNASVNPGpgHGITLtrqqnisgsqNGRLIYPAIVLdrfFLNVMSIYSDDGgsnwq-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------TGSTLpipfrwksssileTLEPSEADMVELQN--GDLLLTARLDFNQivngvny--SPRQQFLSKDGGITWSLLEANNANvfsnistgTVDASITRFEqsdgSHFLLFTNPQGnpagTNgr------------QNLGLWFSFDEG--VTWKGPIQ--LVNGasaysdiyqldsenaivivetdnsnmrilrmpitllkqklt";
	String second = "--------------------------------------------------------------------------------------------kirivdgaanqiqvadgsrkyvvtlsidesgglvanlngvsapiilqsehakvhsfhdyelqysalnhtttLFVDGQQITTWagevsqenniqfgnadaqidgrlhvqkivltqqghnlvefdafylaqqtpevekdleklgwtkiktgntmslygnasvnpgpghgitltrqqnisgsqngrliypaivldrfflnvmsiysddggsnwqTGSTLpipfrwksssileTLEPSEADMVEL--QNGDLLLTARLDFNQivngvny--SPRQQFLSKDGGITWSLLEANNANvfsnisTGTVDASITRFEqsdgSHFLLFTNPQGNpagtngr--------QNLGLWFSFDEG--VTWKGPIQlv---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------NGASAYS--DIYQLd---------SENAIVIVETD---NSNMRILRMPITllkqkltalfdynatgdtefdspakqgwmqdntnngsgvltnadgmpawlvqgiggraqwtyslstnqhaqassfgwrmttemkvlsggmitnyyangtqrvlpiisldssgnlvvefegqtgrtvlatgtaateyhkfelvflpgsnpsasfyfdgklirdniqptaskqnmivwgngssntdgvaayrdikfeiQGDVIf------------RGPDRIPSIVASSVtpGVVTAFAEKRVGGgdpgalsntNDIITRTSRDGGITWDTELNLTEQinvsdefdFSDPRPIYDPs---SNTVLVSYARW----PTdaaqngdrikpwmpNGIFYSVYDVASgnWQAPIDVTdqVKERsfqiagwggselyrrntslnsqqdwqsna------------";
	AFPChain afpChain = FastaAFPChainConverter.cpFastaToAfpChain(first, second, structure1, -393);
	Atom[] ca1 = StructureTools.getAtomCAArray(structure1);
	Atom[] ca2 = StructureTools.getAtomCAArray(structure2);
	String xml = AFPChainXMLConverter.toXML(afpChain);
	System.out.println(xml);
	double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2);
	afpChain.setTMScore(tmScore);
	System.out.println(AfpChainWriter.toScoresList(afpChain));
	StructureAlignmentDisplay.display(afpChain, ca1, ca2);
}
 

public void alignSuperfamily(){
	// download SCOP if required and load into memory
	ScopDatabase scop = ScopFactory.getSCOP();
	List<ScopDescription> superfams = scop.getByCategory(ScopCategory.Superfamily);

	System.out.println("Total nr. of superfamilies:" + superfams.size());


	// configure where to load PDB files from and
	// what information to load
	AtomCache cache = new AtomCache();
	FileParsingParameters fileparams = new FileParsingParameters() ;
	fileparams.setAlignSeqRes(false);
	fileparams.setParseSecStruc(false);
	cache.setFileParsingParams(fileparams);

	// get tge first superfamily
	ScopDescription superfam1 = superfams.get(0);
	System.out.println("First superfamily: " + superfam1);

	ScopNode node = scop.getScopNode(superfam1.getSunID());
	System.out.println("scopNode for first superfamily:" + node);

	List<ScopDomain> doms4superfam1 = scop.getScopDomainsBySunid(superfam1.getSunID());
	ScopDomain dom1 = doms4superfam1.get(0);

	// align the first domain against all others members of this superfamily
	for ( int i = 1 ; i < doms4superfam1.size() ; i ++){

		ScopDomain dom2 = doms4superfam1.get(i);

		try {
			Structure s1 = cache.getStructureForDomain(dom1);
			Structure s2 = cache.getStructureForDomain(dom2);

			Atom[] ca1 = StructureTools.getAtomCAArray(s1);
			Atom[] ca2 = StructureTools.getAtomCAArray(s2);
			StructureAlignment ce = StructureAlignmentFactory.getAlgorithm(CeMain.algorithmName);
			AFPChain afpChain = ce.align(ca1, ca2);

			//System.out.println(afpChain.toCE(ca1, ca2));

			//StructureAlignmentDisplay.display(afpChain, ca1, ca2);

			System.out.println(dom1.getScopId() + " vs. " + dom2.getScopId()+ " :" + afpChain.getProbability());
			double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2);
			afpChain.setTMScore(tmScore);
			System.out.println(AfpChainWriter.toScoresList(afpChain));

		} catch (Exception e){
			e.printStackTrace();
		}
	}

}
 

public String alignPair(String name1, String name2, String algorithmName)
	throws StructureException, IOException {

	// 	make sure each thread has an independent instance of the algorithm object ...

	StructureAlignment algorithm = getAlgorithm(algorithmName);

	// we are running with default parameters

	if ( verbose ) {
		logger.debug("aligning {} against {}", name1, name2);
	}

	long startTime = System.currentTimeMillis();

	if ( prevName1 == null)
		initMaster(name1);

	if ( ! prevName1.equals(name1) ) {
		// we need to reload the master
		initMaster(name1);
	}

	// get a copy of the atoms, but clone them, since they will be rotated...
	Atom[] ca2 =  cache.getAtoms(name2);

	AFPChain afpChain = algorithm.align(ca1, ca2);

	afpChain.setName1(name1);
	afpChain.setName2(name2);

	try {
		// add tmScore
		double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2);
		afpChain.setTMScore(tmScore);
	} catch (RuntimeException e){
		logger.error("ca1 size: {} ca2 length: {} {}  {}", ca1.length, ca2.length, afpChain.getName1(), afpChain.getName2(), e);

	}
	long endTime = System.currentTimeMillis();

	long calcTime = (endTime-startTime);
	if ( verbose ){
		boolean isCP = !AlignmentTools.isSequentialAlignment(afpChain, false);
		String msg = "finished alignment: " + name1 + " vs. " + name2 + " in " + (calcTime) / 1000.0 + " sec.";
		msg += " algo: " + algorithmName + " v:" + version + " " + afpChain;

		if ( isCP ) msg += "HAS A CIRCULAR PERMUTATION!!!";
		logger.debug(msg);
	}
	if (verbose){
		printMemory();
	}
	afpChain.setCalculationTime(calcTime);

	return AFPChainXMLConverter.toXML(afpChain, ca1, ca2);
}
 

private void align (StructureAlignment algorithm, String name1, String name2)
		throws StructureException, IOException{


	AtomCache cache = new AtomCache();
	Atom[] ca1 = cache.getAtoms(name1);
	Atom[] ca2 = cache.getAtoms(name2);


	AFPChain afpChain = algorithm.align(ca1,ca2);
	afpChain.setName1(name1);
	afpChain.setName2(name2);
	double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2);
	afpChain.setTMScore(tmScore);

	String xml = AFPChainXMLConverter.toXML(afpChain, ca1, ca2);

	AFPChain newC    = AFPChainXMLParser.fromXML(xml, ca1, ca2);
	//System.out.println(xml);
	//System.out.println(AFPChainXMLConverter.toXML(newC));
	AFPChain flipped = AFPChainFlipper.flipChain(newC);

	Assert.assertEquals(afpChain.getName1(), flipped.getName2());
	Assert.assertEquals(afpChain.getName2(), flipped.getName1());
	Assert.assertEquals(afpChain.getCa1Length(), flipped.getCa2Length());
	Assert.assertEquals(afpChain.getCa2Length(), flipped.getCa1Length());
	Assert.assertEquals(String.format("%.2f", afpChain.getTMScore()), String.format("%.2f", flipped.getTMScore()));
	Assert.assertTrue(afpChain.getTMScore() != -1);

	String xmlNew = AFPChainXMLConverter.toXML(flipped, ca2, ca1);
	//System.out.println(xmlNew);
	AFPChain backChain = AFPChainXMLParser.fromXML(xmlNew, ca2, ca1);
	AFPChain origFlip  = AFPChainFlipper.flipChain(backChain);

	Assert.assertNotNull("Got null, instead of an AFPChain object!", origFlip);

	Assert.assertNotNull("could not get nr. of eqr: ", afpChain.getNrEQR());
	Assert.assertNotNull("could not get nr. of eqr: ", origFlip.getNrEQR());

	Assert.assertTrue("The nr. of equivalent positions is not equal!", afpChain.getNrEQR() == origFlip.getNrEQR());

	Atom shift1 = afpChain.getBlockShiftVector()[0];
	Atom shift2 = origFlip.getBlockShiftVector()[0];

	Assert.assertTrue("The shift vectors are not similar!", Calc.getDistance(shift1, shift2) < 0.1);

	//assert the RMSD in the flipped alignment is small
	double rmsd1 = getRMSD(afpChain,ca1,ca2);
	double rmsd2 = getRMSD(flipped,ca2,ca1);
	//System.out.println("rmsd:" +rmsd1 + " " + rmsd2);
	Assert.assertTrue("The RMSD are vastly different!", Math.abs(rmsd1 - rmsd2) < 0.01);


	// this can;t work any more because there is minor after comma mismatches..
	//String xmlBack = AFPChainXMLConverter.toXML(origFlip);
	//if ( ! xmlBack.equals(xml)){
	//	printFirstMismatch(xmlBack, xml);
	//}
	//assertEquals("The alignment representations are not the same!" , xmlBack, xml);
	AFPAlignmentDisplay.getAlign(origFlip, ca1, ca2);
	String img1 = AfpChainWriter.toDBSearchResult(afpChain);
	String img2 = AfpChainWriter.toDBSearchResult(origFlip);
	Assert.assertEquals("The alignment images do not match!", img1, img2);

	//System.out.println(xml);
	//System.out.println(xmlNew);



}
 

public static void main(String[] args){

		//String name1 = "4hhb.A";
		//String name2 = "4hhb.B";

		String name1 = "1cdg.A";
		String name2 = "1tim.B";



		AtomCache cache = new AtomCache();

		Structure structure1 = null;
		Structure structure2 = null;

		try {

		   StructureAlignment algorithm  = StructureAlignmentFactory.getAlgorithm(CeMain.algorithmName);

			structure1 = cache.getStructure(name1);
			structure2 = cache.getStructure(name2);

			Atom[] ca1 = StructureTools.getAtomCAArray(structure1);
			Atom[] ca2 = StructureTools.getAtomCAArray(structure2);

			// get default parameters
			CeParameters params = new CeParameters();

			// add more print
			params.setShowAFPRanges(true);

			// set the maximum gap size to unlimited
			params.setMaxGapSize(-1);

			AFPChain afpChain = algorithm.align(ca1,ca2,params);

			afpChain.setName1(name1);
			afpChain.setName2(name2);

			// flexible original results:
			System.out.println(afpChain.toFatcat(ca1,ca2));

			System.out.println(afpChain.toRotMat());
			//System.out.println(afpChain.toCE(ca1, ca2));

			System.out.println(AFPChainXMLConverter.toXML(afpChain,ca1,ca2));

			double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2);
			afpChain.setTMScore(tmScore);

			//System.out.println(AfpChainWriter.toWebSiteDisplay(afpChain, ca1, ca2));

			printScores(afpChain);
		} catch (Exception e) {
			e.printStackTrace();
			return;
		}
	}
 

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

}
 

public static void main(String[] args){
	//String name1 = "4hhb.A";
	//String name2 = "4hhb.B";

	String name1 = "d1pqsa_";
	String name2 = "d1poha_";

	//String name1 = "5AZQ.A";
	//String name2 = "4ODC.A";

	AtomCache cache = new AtomCache();

	DownloadChemCompProvider prov = new DownloadChemCompProvider();
	prov.setDownloadAll(true);

	ChemCompGroupFactory.setChemCompProvider(prov);

	Structure structure1 = null;
	Structure structure2 = null;

	try {

		StructureAlignment algorithm  = StructureAlignmentFactory.getAlgorithm(CeMain.algorithmName);

		structure1 = cache.getStructure(name1);
		structure2 = cache.getStructure(name2);

		Atom[] ca1 = StructureTools.getAtomCAArray(structure1);
		Atom[] ca2 = StructureTools.getAtomCAArray(structure2);

		// get default parameters
		CeParameters params = new CeParameters();

		// add more print
		params.setShowAFPRanges(true);

		// set the maximum gap size to unlimited
		params.setMaxGapSize(-1);

		AFPChain afpChain = algorithm.align(ca1,ca2,params);

		afpChain.setName1(name1);
		afpChain.setName2(name2);

		// show a nice summary print
		System.out.println(AfpChainWriter.toWebSiteDisplay(afpChain, ca1, ca2));

		// print rotation matrices
		System.out.println(afpChain.toRotMat());
		//System.out.println(afpChain.toCE(ca1, ca2));

		// print XML representation
		//System.out.println(AFPChainXMLConverter.toXML(afpChain,ca1,ca2));

		StructureAlignmentDisplay.display(afpChain, ca1, ca2);


		double tmScore = AFPChainScorer.getTMScore(afpChain, ca1, ca2);
		afpChain.setTMScore(tmScore);
		System.out.println(AfpChainWriter.toScoresList(afpChain));

	} catch (Exception e) {
		e.printStackTrace();
		return;
	}
}