Added new methods to the AlignmentTools class: a method to calculate the

lafita · lafita · commit 30cc0b46db42 · 2015-03-24T15:44:35.000+01:00
gaps in each block of an optimal alignment and a method to replace the
optimal alignment of an AFPChain.
diff --git a/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AlignmentTools.java b/biojava-structure/src/main/java/org/biojava/nbio/structure/align/util/AlignmentTools.java
@@ -672,6 +672,58 @@ public static AFPChain splitBlocksByTopology(AFPChain a, Atom[] ca1, Atom[] ca2)
 
 		return replaceOptAln(a, ca1, ca2, blocks.size(), newBlockLens, newOptAln);
 	}
+	
+	/**
+	 * It returns an AFPChain with a segmented optimal alignment, which means that it has <order of symmetry> blocks of
+	 * aligned subunits. The original AFPChain is not modified.
+	 * 
+	 * INPUT: the optimal alignment in a triple list (same format as optAln of AFPChain) and the Atom[] arrays.
+	 * OUTPUT: the optimal AFPChain alignment object divided into the subunits.
+	 */
+	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];
+		}
+		
+		/*//Print the sizes to check correctness
+		System.out.println("Number of subunits: "+newAlgn.length);
+		System.out.println("Subunit length: "+newAlgn[0][0].length);*/
+		
+		//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.cloneCAArray(ca2); // don't modify ca1 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;
+	}
 
 	/**
 	 * Takes an AFPChain and replaces the optimal alignment based on an alignment map
@@ -1031,4 +1083,46 @@ public static Map<Integer, Integer> fromConciseAlignmentString(String string) {
 		}
 		return map;
 	}
+	
+	/**
+	 * Method that calculates the number of gaps in each subunit block of an optimal AFP alignment.
+	 * 
+	 * INPUT: an optimal alignment in the format int[][][].
+	 * OUTPUT: an int[] array of <order> length containing the gaps in each block as int[block].
+	 */
+	public static int[] calculateBlockGap(int[][][] optAln){
+		
+		//Initialize the array to be returned
+		int [] blockGap = new int[optAln.length];
+		
+		//Loop for every block and look in both chains for non-contiguous residues.
+		for (int i=0; i<optAln.length; i++){
+			int gaps = 0; //the number of gaps in that block
+			int last1 = 0; //the last residue position in chain 1
+			int last2 = 0; //the last residue position in chain 2
+			//Loop for every position in the block
+			for (int j=0; j<optAln[i][0].length; j++){
+				//If the first position is evaluated initialize the last positions
+				if (j==0){
+					last1 = optAln[i][0][j];
+					last2 = optAln[i][1][j];
+				}
+				else{
+					//If one of the positions or both are not contiguous increment the number of gaps
+					if (optAln[i][0][j] > last1+1 || optAln[i][1][j] > last2+1){
+						gaps++;
+						last1 = optAln[i][0][j];
+						last2 = optAln[i][1][j];
+					}
+					//Otherwise just set the last position to the current one
+					else{
+						last1 = optAln[i][0][j];
+						last2 = optAln[i][1][j];
+					}
+				}
+			}
+			blockGap[i] = gaps;
+		}
+		return blockGap;
+	}
 }
