org.biojava.nbio.structure.align.model.AFPChain Java Examples

@Test
public void testGetSymmetryOrder() throws IOException, StructureException, RefinerFailedException {
	// List of alignments to try, along with proper symmetry
	Map<String,Integer> orderMap = new HashMap<String,Integer>();
	orderMap.put("1itb.A",3); // b-trefoil, C3
	orderMap.put("1tim.A",2); // tim-barrel, C8
	//orderMap.put("d1p9ha_",-1); // not rotational symmetry
	orderMap.put("3HKE.A",2); // very questionable alignment
	orderMap.put("d1jlya1",3); // a very nice trefoil

	AtomCache cache = new AtomCache();

	for(String name : orderMap.keySet()) {
		CESymmParameters params = new CESymmParameters();
		params.setRefineMethod(RefineMethod.NOT_REFINED);
		Atom[] ca1 = cache.getAtoms(name);

		CeSymmResult result = CeSymm.analyzeLevel(ca1, params);
		AFPChain afpChain = result.getSelfAlignment();

		int order = new SequenceFunctionOrderDetector().calculateOrder(afpChain, ca1);

		assertEquals("Wrong order for "+name,orderMap.get(name).intValue(), order);
	}
}
 

/** Circular permutation specific code to be run after the standard CE alignment
 *
 * @param afpChain The finished alignment
 * @param ca1 CA atoms of the first protein
 * @param ca2m A duplicated copy of the second protein
 * @param calculator The CECalculator used to create afpChain
 * @param param Parameters
 * @throws StructureException
 */
public static AFPChain postProcessAlignment(AFPChain afpChain, Atom[] ca1, Atom[] ca2m,CECalculator calculator, CECPParameters param ) throws StructureException{

	// remove bottom half of the matrix
	Matrix doubledMatrix = afpChain.getDistanceMatrix();

	// the matrix can be null if the alignment is too short.
	if ( doubledMatrix != null ) {
		assert(doubledMatrix.getRowDimension() == ca1.length);
		assert(doubledMatrix.getColumnDimension() == ca2m.length);

		Matrix singleMatrix = doubledMatrix.getMatrix(0, ca1.length-1, 0, (ca2m.length/2)-1);
		assert(singleMatrix.getRowDimension() == ca1.length);
		assert(singleMatrix.getColumnDimension() == (ca2m.length/2));

		afpChain.setDistanceMatrix(singleMatrix);
	}
	// Check for circular permutations
	int alignLen = afpChain.getOptLength();
	if ( alignLen > 0) {
		afpChain = filterDuplicateAFPs(afpChain,calculator,ca1,ca2m,param);
	}
	return afpChain;
}
 

/**
 * Aligns ca1 with ca2 permuted by <i>cp</i> residues.
 * <p><strong>WARNING:</strong> Modifies ca2 during the permutation. Be sure
 * to make a copy before calling this method.
 *
 * @param ca1
 * @param ca2
 * @param param
 * @param cp
 * @return
 * @throws StructureException
 */
public AFPChain alignPermuted(Atom[] ca1, Atom[] ca2, Object param, int cp) throws StructureException {
	// initial permutation
	permuteArray(ca2,cp);

	// perform alignment
	AFPChain afpChain = super.align(ca1, ca2, param);

	// un-permute alignment
	permuteAFPChain(afpChain, -cp);

	if(afpChain.getName2() != null) {
		afpChain.setName2(afpChain.getName2()+" CP="+cp);
	}

	// Specify the permuted
	return afpChain;
}
 

/**
 * Creates a Map specifying the alignment as a mapping between residue indices
 * of protein 1 and residue indices of protein 2.
 *
 * <p>For example,<pre>
 * 1234
 * 5678</pre>
 * becomes<pre>
 * 1->5
 * 2->6
 * 3->7
 * 4->8</pre>
 *
 * @param afpChain An alignment
 * @return A mapping from aligned residues of protein 1 to their partners in protein 2.
 * @throws StructureException If afpChain is not one-to-one
 */
public static Map<Integer, Integer> alignmentAsMap(AFPChain afpChain) throws StructureException {
	Map<Integer,Integer> map = new HashMap<Integer,Integer>();

	if( afpChain.getAlnLength() < 1 ) {
		return map;
	}
	int[][][] optAln = afpChain.getOptAln();
	int[] optLen = afpChain.getOptLen();
	for(int block = 0; block < afpChain.getBlockNum(); block++) {
		for(int pos = 0; pos < optLen[block]; pos++) {
			int res1 = optAln[block][0][pos];
			int res2 = optAln[block][1][pos];
			if(map.containsKey(res1)) {
				throw new StructureException(String.format("Residue %d aligned to both %d and %d.", res1,map.get(res1),res2));
			}
			map.put(res1,res2);
		}
	}
	return map;
}
 

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

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

private void postProcessAlignment(AFPChain afpChain) {
	if (mirrorSequence) {
		// reverse the optAln
		reverseOptAln(afpChain);

		// reverse the distance matrix
		afpChain.setDistanceMatrix(reverseMatrixCols(afpChain
				.getDistanceMatrix()));

		// reverse the ca2 matrix
		Matrix distMat2 = afpChain.getDisTable2();
		distMat2 = reverseMatrixRows(distMat2);
		distMat2 = reverseMatrixCols(distMat2);
		afpChain.setDisTable2(distMat2);
	}

}
 

/**
 * This method takes the all-to-all pairwise alignments Matrix (as a
 * double List of AFPChain) and calculates the structure with the
 * lowest average RMSD against all others.
 * The index of this structure is returned.
 *
 * @param alignments List double containing all-to-all pairwise alignments
 * @return int reference index
 */
private static int chooseReferenceRMSD(List<List<AFPChain>> afpAlignments){

	int size = afpAlignments.size();

	List<Double> RMSDs = new ArrayList<Double>();
	for (int i=0; i<afpAlignments.size(); i++){
		double rmsd=0.0;
		for (int j=0; j<size; j++){
			if (i!=j)
				rmsd += afpAlignments.get(i).get(j).getTotalRmsdOpt();
		}
		RMSDs.add(rmsd);
	}
	int reference = 0;
	for (int i=1; i<size; i++){
		if (RMSDs.get(i) < RMSDs.get(reference)) reference = i;
	}
	logger.info("Reference structure is "+reference);
	return reference;
}
 

/**
 * Converts a set of AFP alignments into a Graph of aligned residues, where
 * each vertex is a residue and each edge means the connection between the
 * two residues in one of the alignments.
 *
 * @param afps
 *            List of AFPChains
 * @param atoms
 *            Atom array of the symmetric structure
 * @param undirected
 *            if true, the graph is undirected
 *
 * @return adjacency List of aligned residues
 */
public static List<List<Integer>> buildSymmetryGraph(List<AFPChain> afps,
		Atom[] atoms, boolean undirected) {

	List<List<Integer>> graph = new ArrayList<List<Integer>>();

	for (int n = 0; n < atoms.length; n++) {
		graph.add(new ArrayList<Integer>());
	}

	for (int k = 0; k < afps.size(); k++) {
		for (int i = 0; i < afps.get(k).getOptAln().length; i++) {
			for (int j = 0; j < afps.get(k).getOptAln()[i][0].length; j++) {
				Integer res1 = afps.get(k).getOptAln()[i][0][j];
				Integer res2 = afps.get(k).getOptAln()[i][1][j];
				graph.get(res1).add(res2);
				if (undirected)
					graph.get(res2).add(res1);
			}
		}
	}
	return graph;
}
 

private static void printScores(AFPChain afpChain) {
	System.out.println("=====================");
	System.out.println("The main scores for the alignment:");

	System.out.println("EQR       :\t" + afpChain.getNrEQR() + "\t The number of residues on structurally equivalent positions.")  ;
	System.out.println("RMSD      :\t" + String.format("%.2f",afpChain.getTotalRmsdOpt() )+ "\t The RMSD of the alignment");
	System.out.println("Z-score   :\t" + afpChain.getProbability() + "\t The Z-score of the alignment (CE)");
	System.out.println("TM-score  :\t" + String.format("%.2f",afpChain.getTMScore()) + "\t The TM-score of the alignment.");
	System.out.println("");
	System.out.println("Other scores:");
	System.out.println("Identity  :\t" + String.format("%.2f",afpChain.getIdentity())   + "\t The percent of residues that are sequence-identical in the alignment.");
	System.out.println("Similarity:\t" + String.format("%.2f",afpChain.getSimilarity()) + "\t The percent of residues in the alignment that are sequence-similar.");
	System.out.println("Coverage1 :\t" + afpChain.getCoverage1() + " %\t Percent of protein 1 that is covered with the alignment.");
	System.out.println("Coverage2 :\t" + afpChain.getCoverage2() + " %\t Percent of protein 2 that is covered with the alignment.");
	int dab = afpChain.getCa1Length()+afpChain.getCa2Length() - 2 * afpChain.getNrEQR();
	System.out.println("Distance  :\t" + dab + "\t Distance between folds a,b ");
	double sab = 2 * afpChain.getNrEQR() / (double)( afpChain.getCa1Length() + afpChain.getCa2Length());
	System.out.println("Rel. Sim. :\t" + String.format("%.2f",sab) + "\t Relative similarity");



}
 

@Test
public void testCalculateOrderByHarmonicsFloating() throws IOException,
		StructureException, RefinerFailedException {
	String name;
	RotationOrderDetector detector = new RotationOrderDetector(8,
			HARMONICS_FLOATING);

	// Perform alignment to determine axis
	Atom[] ca1;
	AFPChain alignment;
	int order;

	name = "1MER.A";
	ca1 = StructureTools.getRepresentativeAtomArray(StructureTools
			.getStructure(name));
	alignment = CeSymm.analyze(ca1, params).getSelfAlignment();
	order = detector.calculateOrder(alignment, ca1);
	assertEquals(name, 2, order);

	name = "d1ijqa1";
	ca1 = StructureTools.getRepresentativeAtomArray(StructureTools
			.getStructure(name));
	alignment = CeSymm.analyze(ca1, params).getSelfAlignment();
	order = detector.calculateOrder(alignment, ca1);
	assertEquals(name, 6, order);

	name = "1TIM.A";
	ca1 = StructureTools.getRepresentativeAtomArray(StructureTools
			.getStructure(name));
	alignment = CeSymm.analyze(ca1, params).getSelfAlignment();
	order = detector.calculateOrder(alignment, ca1);
	assertEquals(name, 6, order);// tough case

}
 

/**
 * Return a list of pdb Strings corresponding to the aligned positions of the molecule.
 * Only supports a pairwise alignment with the AFPChain DS.
 *
 * @param aligPos
 * @param afpChain
 * @param ca
 */
public static final List<String> getPDBresnum(int aligPos, AFPChain afpChain, Atom[] ca){
	List<String> lst = new ArrayList<String>();
	if ( aligPos > 1) {
		System.err.println("multiple alignments not supported yet!");
		return lst;
	}

	int blockNum = afpChain.getBlockNum();
	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	if ( optLen == null)
		return lst;

	for(int bk = 0; bk < blockNum; bk ++)       {

		for ( int i=0;i< optLen[bk];i++){

			int pos = optAln[bk][aligPos][i];
			if ( pos < ca.length) {
				String pdbInfo = JmolTools.getPdbInfo(ca[pos]);
				//lst.add(ca1[pos].getParent().getPDBCode());
				lst.add(pdbInfo);
			}
		}

	}
	return lst;
}
 

/**
 * Takes a structure and sequence corresponding to an alignment between a structure or sequence and itself (or even a structure with a sequence), where the result has a circular permutation site
 * {@link cpSite} residues to the right.
 *
 * @param fastaFile A FASTA file containing exactly 2 sequences, the first unpermuted and the second permuted
 * @param cpSite
 *            The number of residues from the beginning of the sequence at which the circular permutation site occurs; can be positive or negative; values greater than the length of the sequence
 *            are acceptable
 * @throws IOException
 * @throws StructureException
 */
public static AFPChain cpFastaToAfpChain(File fastaFile, Structure structure, int cpSite) throws IOException, StructureException {
	InputStream inStream = new FileInputStream(fastaFile);
	SequenceCreatorInterface<AminoAcidCompound> creator = new CasePreservingProteinSequenceCreator(AminoAcidCompoundSet.getAminoAcidCompoundSet());
	SequenceHeaderParserInterface<ProteinSequence, AminoAcidCompound> headerParser = new GenericFastaHeaderParser<ProteinSequence, AminoAcidCompound>();
	FastaReader<ProteinSequence, AminoAcidCompound> fastaReader = new FastaReader<ProteinSequence, AminoAcidCompound>(inStream, headerParser, creator);
	LinkedHashMap<String, ProteinSequence> sequences = fastaReader.process();
	inStream.close();
	Iterator<ProteinSequence> iter = sequences.values().iterator();
	ProteinSequence first = iter.next();
	ProteinSequence second = iter.next();
	return cpFastaToAfpChain(first, second, structure, cpSite);
}
 

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

@Test
public void testCpSymmetric1() throws IOException,StructureException, CompoundNotFoundException {
	//cat 2GG6-best.fasta |tr -d \\n|pbcopy
	String a = "-SSRPATAR-KSSGLSGTVRIPGDKSISHRSFMFGGLA-SGETRITGLLEG-EDvINTGKAMQAMGARIRKEGd---------TWIIDGVgngglLAPEAPLD---FGNAATGCRLTMGLVGvydFDSTFIGDASLtkrp---MGRVLNPLREMGVQVKSEDgdrLPVTLRGPK---TPT---PITYRVpMASAQVKSAVLLAGLNTPGITTVIEpi---MTRDHTEKMLQGFGANLTVEtdadGVRTIRLEgRGKLTGQVIDVPGDPSSTAFPLVAALLVpGSDVTILNVLMNpTR-TGLILTLQEMGADIEVINprlaggedvaDLRVRSS-----TLKGVTVPedrAPSMIDEYPILAVAAAFAEGATVMNGLEELrvkesdrLSAVANGLKLNGVDCDEGE---TSLVVRGRPdgkGLGNasgAAVAT-HLDHRIAMSFLVMGLVSENPVTVDDatmIATSFPEFMDLMAGLGAKIELS---";
	String b = "dGVRTIRLEgRGKLTGQVIDVPGDPSSTAFPLVAALLVpGSDVTILNVLMNpTR-TGLILTLQEMGADIEVINprlaggedvaDLRVRSS-----TLKGVTVPedrAPSMIDEYPILAVAAAfaeGATVMNGLEELrvkesdrLSAVANGLKLNGVDCDEGE---TSLVVRGRPdgkGLGnasGAAVAT-HLDHRIAMSFLVMGLVSENPVTVDDatmiaTSFPEFMDLMAGLGAKIELS----SSRPATAR-KSSGLSGTVRIPGDKSISHRSFMFGGLA-SGETRITGLLEG-EDvINTGKAMQAMGARIRKEGd---------TWIIDGVgngglLAPEAPLD---FGNAATGCRLTMGLVGVYDFDSTFIGDASLtkrp---MGRVLNPLREMGVQVKSEDgdrLPVTLRGPK---TPTP---ITYRVpMASAQVKSAVLLAGLNTPGITTVIE---PIMTRDHTEKMLQGFGANLTVEtda";
	Structure structure = StructureTools.getStructure("2GG6");
	AFPChain afpChain = FastaAFPChainConverter.cpFastaToAfpChain(a, b, structure, -230); // 215
	assertEquals("Wrong TM-score", 0.7701, afpChain.getTMScore(), 0.001);
	assertEquals("Wrong RMSD", 3.035, afpChain.getTotalRmsdOpt(), 0.001);
}
 

@Override
public AFPChain align(Atom[] ca1, Atom[] ca2) throws StructureException {

	AFPChain afpChain = alignRigid(ca1, ca2, params);
	afpChain.setAlgorithmName(algorithmName);
	afpChain.setVersion(VERSION+"");
	return afpChain;
}
 

private static void reverseOptAln(AFPChain afpChain) {
	int ca2len = afpChain.getCa2Length();

	int[][][] optAln = afpChain.getOptAln();
	int[] optLen = afpChain.getOptLen();

	for (int block = 0; block < afpChain.getBlockNum(); block++) {
		for (int pos = 0; pos < optLen[block]; pos++) {
			optAln[block][1][pos] = ca2len - 1 - optAln[block][1][pos];
		}
	}

	afpChain.setOptAln(optAln);
}
 

/**
 * Calculate the rotation axis for the first block of an AFPChain
 * @param afpChain
 * @throws StructureException
 * @throws NullPointerException if afpChain does not contain a valid rotation matrix and shift vector
 */
public RotationAxis(AFPChain afpChain) throws StructureException {
	if(afpChain.getAlnLength() < 1) {
		throw new StructureException("No aligned residues");
	}
	init(afpChain.getBlockRotationMatrix()[0],afpChain.getBlockShiftVector()[0]);
}
 

@Override
public AFPChain align(Atom[] ca1, Atom[] ca2, Object param)
throws StructureException {

	if ( ! (param instanceof FatCatParameters)){
		throw new IllegalArgumentException("FatCat algorithm needs FatCatParameters object as argument.");
	}

	params = (FatCatParameters) param;

	AFPChain afpChain= alignFlexible(ca1, ca2, params);
	afpChain.setAlgorithmName(algorithmName);
	afpChain.setVersion(VERSION+"");
	return afpChain;
}
 

/**
 * Constructor using an afpChain and the atom arrays for pairwise
 * alignments. The AFPChain is converted into a MultipleAlignment.
 *
 * @param afpChain
 * @param ca1
 * @param ca2
 * @throws StructureException
 */
public MultipleAligPanel(AFPChain afpChain, Atom[] ca1, Atom[] ca2,
		AbstractAlignmentJmol jmol) throws StructureException {

	this();

	String algorithm = afpChain.getAlgorithmName();
	boolean flex = false;
	if (algorithm != null){
		if (algorithm.contains("flexible")) flex = true;
	}

	//Convert the apfChain into a MultipleAlignment object
	MultipleAlignmentEnsembleImpl ensemble =
			new MultipleAlignmentEnsembleImpl(afpChain, ca1, ca2, flex);
	this.multAln = ensemble.getMultipleAlignment(0);

	//Create the sequence alignment and the structure-sequence mapping.
	this.mapSeqToStruct = new ArrayList<Integer>();
	this.alnSeq = MultipleAlignmentTools.getSequenceAlignment(
			this.multAln, this.mapSeqToStruct);

	//Initialize other memeber variables of the panel
	this.size = multAln.size();
	this.length = alnSeq.get(0).length();

	coordManager = new MultipleAlignmentCoordManager(size, length);
	this.jmol = jmol;
}
 

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

private static void printXMLMatrixShift(PrettyXMLWriter xml,
		AFPChain afpChain, int blockNr)  throws IOException {

	Matrix[] ms     = afpChain.getBlockRotationMatrix();
	if ( ms == null || ms.length == 0)
		return;

	Matrix matrix = ms[blockNr];
	if ( matrix == null)
		return;
	xml.openTag("matrix");


	for (int x=0;x<3;x++){
		for (int y=0;y<3;y++){
			String key = "mat"+(x+1)+(y+1);
			xml.attribute(key,String.format("%.6f",matrix.get(x,y)));
		}
	}
	xml.closeTag("matrix");

	Atom[]   shifts = afpChain.getBlockShiftVector();
	Atom shift = shifts[blockNr];
	xml.openTag("shift");
	xml.attribute("x", String.format("%.3f",shift.getX()));
	xml.attribute("y", String.format("%.3f",shift.getY()));
	xml.attribute("z", String.format("%.3f",shift.getZ()));
	xml.closeTag("shift");

}
 

/** new utility method that checks that the order of the pair in the XML alignment is correct and flips the direction if needed
*
* @param xml
* @param name1
* @param name1
* @param ca1
* @param ca2
* @return
*/
public static AFPChain fromXML(String xml, String name1, String name2, Atom[] ca1, Atom[] ca2) throws IOException, StructureException{
	AFPChain[] afps = parseMultiXML( xml);
	if ( afps.length > 0 ) {

		AFPChain afpChain = afps[0];

		String n1 = afpChain.getName1();
		String n2 = afpChain.getName2();

		if ( n1 == null )
			n1 = "";
		if ( n2 == null)
			n2 = "";

		//System.out.println("from AFPCHAIN: " + n1 + " " + n2);
		if ( n1.equals(name2) && n2.equals(name1)){
			// flipped order
			//System.out.println("AfpChain in wrong order, flipping...");
			afpChain  = AFPChainFlipper.flipChain(afpChain);
		}
		rebuildAFPChain(afpChain, ca1, ca2);

		return afpChain;
	}
	return null;

}
 

/**
 * Find the alignment position with the highest atomic distance between the
 * equivalent atomic positions of the arrays and remove it from the
 * alignment.
 *
 * @param afpChain
 *            original alignment, will be modified
 * @param ca1
 *            atom array, will not be modified
 * @param ca2
 *            atom array, will not be modified
 * @return the original alignment, with the alignment position at the
 *         highest distance removed
 * @throws StructureException
 */
public static AFPChain deleteHighestDistanceColumn(AFPChain afpChain,
		Atom[] ca1, Atom[] ca2) throws StructureException {

	int[][][] optAln = afpChain.getOptAln();

	int maxBlock = 0;
	int maxPos = 0;
	double maxDistance = Double.MIN_VALUE;

	for (int b = 0; b < optAln.length; b++) {
		for (int p = 0; p < optAln[b][0].length; p++) {
			Atom ca2clone = ca2[optAln[b][1][p]];
			Calc.rotate(ca2clone, afpChain.getBlockRotationMatrix()[b]);
			Calc.shift(ca2clone, afpChain.getBlockShiftVector()[b]);

			double distance = Calc.getDistance(ca1[optAln[b][0][p]],
					ca2clone);
			if (distance > maxDistance) {
				maxBlock = b;
				maxPos = p;
				maxDistance = distance;
			}
		}
	}

	return deleteColumn(afpChain, ca1, ca2, maxBlock, maxPos);
}
 

@Test
public void testCalculateOrderByHarmonics() throws IOException,
		StructureException, RefinerFailedException {
	String name;
	RotationOrderDetector detector = new RotationOrderDetector(8, HARMONICS);

	// Perform alignment to determine axis
	Atom[] ca1;
	AFPChain alignment;
	int order;

	name = "1MER.A";
	ca1 = StructureTools.getRepresentativeAtomArray(StructureTools
			.getStructure(name));
	alignment = CeSymm.analyze(ca1, params).getSelfAlignment();
	order = detector.calculateOrder(alignment, ca1);
	assertEquals(name, 2, order);

	name = "d1ijqa1";
	ca1 = StructureTools.getRepresentativeAtomArray(StructureTools
			.getStructure(name));
	alignment = CeSymm.analyze(ca1, params).getSelfAlignment();
	order = detector.calculateOrder(alignment, ca1);
	assertEquals(name, 6, order);

	name = "1TIM.A";
	ca1 = StructureTools.getRepresentativeAtomArray(StructureTools
			.getStructure(name));
	alignment = CeSymm.analyze(ca1, params).getSelfAlignment();
	order = detector.calculateOrder(alignment, ca1);
	assertEquals(name, 1, order);// tough case
}
 

/**
 * Calculate the rotation angle for a structure
 * @param afpChain
 * @return The rotation angle, in radians
 * @throws StructureException If the alignment doesn't contain any blocks
 * @throws NullPointerException If the alignment doesn't have a rotation matrix set
 */
public static double getAngle(AFPChain afpChain) throws StructureException {
	if(afpChain.getBlockNum() < 1) {
		throw new StructureException("No aligned residues");
	}
	Matrix rotation = afpChain.getBlockRotationMatrix()[0];

	if(rotation == null) {
		throw new NullPointerException("AFPChain does not contain a rotation matrix");
	}
	return getAngle(rotation);
}
 

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

}
 

@Test
public void testCpAsymmetric() throws IOException, StructureException, CompoundNotFoundException {
	Structure structure = cache.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, structure, -393);
	assertEquals("Wrong TM-score", 0.2949, afpChain.getTMScore(), 0.001);
	assertEquals("Wrong RMSD", 3.605, afpChain.getTotalRmsdOpt(), 0.001);
}
 

/**
return the rmsd of two blocks
 */
private static double combineRmsd(int b1, int b2, AFPChain afpChain,Atom[] ca1,Atom[] ca2)
{
	int     i;
	int     afpn = 0;

	int[] afpChainList =afpChain.getAfpChainList();

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


	int[]   list = new int[blockSize[b1]+blockSize[b2]];
	for(i = block2Afp[b1]; i < block2Afp[b1] + blockSize[b1]; i ++) {
		list[afpn ++] = afpChainList[i];
	}
	for(i = block2Afp[b2]; i < block2Afp[b2] + blockSize[b2]; i ++) {
		list[afpn ++] = afpChainList[i];
	}
	double  rmsd = AFPChainer.calAfpRmsd(afpn, list,0, afpChain,ca1,ca2);

	afpChain.setBlock2Afp(block2Afp);
	afpChain.setBlockSize(blockSize);
	afpChain.setAfpChainList(afpChainList);

	return rmsd;
}
 

@SuppressWarnings({ "rawtypes", "unchecked" })
public static void showAlignmentImage(AFPChain afpChain, Atom[] ca1,
		Atom[] ca2, Object jmol)
				throws ClassNotFoundException, NoSuchMethodException, InvocationTargetException, IllegalAccessException{

	Class abstractAligJmolC = Class.forName(abstractAligJmol);

	Class c = Class.forName(displayAFP);
	Method show = c.getMethod("showAlignmentPanel", new Class[] {AFPChain.class, Atom[].class, Atom[].class, abstractAligJmolC});

	show.invoke(null,afpChain, ca1, ca2, jmol);
}