package bgu.bio.algorithms.alignment;
import gnu.trove.list.array.TCharArrayList;
import bgu.bio.util.ScoringMatrix;
import bgu.bio.util.alphabet.AlphabetUtils;
public abstract class AbsSequenceAlignment implements SequenceAlignment {
/**
* The first input string
*/
protected TCharArrayList str1;
/**
* The second input String
*/
protected TCharArrayList str2;
/**
* The lengths of the input strings
*/
protected int length1, length2;
/**
* The score matrix. The true scores should be divided by the normalization
* factor.
*/
protected double[][] dpTable;
protected ScoringMatrix scoringMatrix;
/**
* The normalization factor. To get the true score, divide the integer score
* used in computation by the normalization factor.
*/
private double normFactor = 1.0;
/**
* Constants of directions. Multiple directions are stored by bits. The zero
* direction is the starting point.
*/
protected static final int DR_LEFT = 1; // 0001
protected static final int DR_UP = 2; // 0010
protected static final int DR_DIAG = 4; // 0100
protected static final int DR_ZERO = 8; // 1000
/**
* The directions pointing to the cells that give the maximum score at the
* current cell.
*/
protected int[][] prevCells;
protected AlphabetUtils alphabet;
public AbsSequenceAlignment(String str1, String str2,
AlphabetUtils alphabet, ScoringMatrix matrix) {
this(str1.length(), str2.length(), alphabet, matrix);
this.setSequences(str1, str2);
}
public AbsSequenceAlignment(int size1, int size2, AlphabetUtils alphabet,
ScoringMatrix matrix) {
this.scoringMatrix = matrix;
this.alphabet = alphabet;
this.str1 = new TCharArrayList();
this.str2 = new TCharArrayList();
length1 = size1;
length2 = size2;
dpTable = new double[length1 + 1][length2 + 1];
prevCells = new int[length1 + 1][length2 + 1];
}
@Override
public void setSequences(String str1, String str2) {
this.setSequences(str1.toCharArray(), str2.toCharArray());
}
@Override
public void setSequences(TCharArrayList str1, TCharArrayList str2) {
this.str1.resetQuick();
for (int i = 0; i < str1.size(); i++) {
this.str1.add(str1.get(i));
}
this.str2.resetQuick();
for (int i = 0; i < str2.size(); i++) {
this.str2.add(str2.get(i));
}
boolean shouldInitTables = str1.size() > dpTable.length
|| str2.size() > dpTable[0].length;
length1 = str1.size();
length2 = str2.size();
if (shouldInitTables) {
dpTable = new double[length1 + 1][length2 + 1];
prevCells = new int[length1 + 1][length2 + 1];
}
}
@Override
public void setSequences(char[] str1, char[] str2) {
setSequences(str1, str1.length, str2, str2.length);
}
@Override
public void setSequences(char[] str1, int size1, char[] str2, int size2) {
this.str1.resetQuick();
for (int i = 0; i < str1.length; i++) {
this.str1.add(str1[i]);
}
this.str2.resetQuick();
for (int i = 0; i < str2.length; i++) {
this.str2.add(str2[i]);
}
boolean init = size1 > dpTable.length - 1
|| size2 > dpTable[0].length - 1;
length1 = size1;
length2 = size2;
if (init) {
dpTable = new double[length1 + 1][length2 + 1];
prevCells = new int[length1 + 1][length2 + 1];
}
}
/**
* Compute the similarity score of substitution: use a substitution matrix
* if the cost model The position of the first character is 1. A position of
* 0 represents a gap.
*
* @param i
* Position of the character in str1
* @param j
* Position of the character in str2
* @return cost of substitution of the character in str1 by the one in str2
*/
protected double similarity(int i, int j) {
final char c1 = i == 0 ? this.alphabet.emptyLetter() : str1.get(i - 1);
final char c2 = j == 0 ? this.alphabet.emptyLetter() : str2.get(j - 1);
return this.scoringMatrix.score(c1, c2);
}
/**
* Build the score matrix using dynamic programming. Note: The indel scores
* must be negative. Otherwise, the part handling the first row and column
* has to be modified.
*/
@Override
public void buildMatrix() {
int i; // length of prefix substring of str1
int j; // length of prefix substring of str2
initTable();
// the rest of the matrix
for (i = 1; i <= length1; i++) {
for (j = 1; j <= length2; j++) {
computeCell(i, j);
}
}
}
protected abstract void initTable();
protected abstract void computeCell(int i, int j);
/**
* Get the maximum value in the score matrix.
*/
protected abstract double getMaxScore();
/**
* Get the alignment score between the two input strings.
*/
@Override
public double getAlignmentScore() {
return getMaxScore() / normFactor;
}
/**
* Output the local alignments ending in the (i, j) cell. aligned1 and
* aligned2 are suffixes of final aligned strings found in backtracking
* before calling this function. Note: the strings are replicated at each
* recursive call. Use buffers or stacks to improve efficiency.
*/
private String[] getAlignment(int i, int j, String aligned1, String aligned2) {
// we've reached the starting point, so print the alignments
if ((prevCells[i][j] & DR_ZERO) > 0) {
return new String[] { aligned1, aligned2 };
}
// find out which directions to backtrack
if ((prevCells[i][j] & DR_UP) > 0 && i != 0) {
return getAlignment(i - 1, j, str1.get(i - 1) + aligned1, "_"
+ aligned2);
} else if ((prevCells[i][j] & DR_LEFT) > 0 && j != 0) {
return getAlignment(i, j - 1, "_" + aligned1, str2.get(j - 1)
+ aligned2);
} else if ((prevCells[i][j] & DR_DIAG) > 0) {
return getAlignment(i - 1, j - 1, str1.get(i - 1) + aligned1,
str2.get(j - 1) + aligned2);
}
return null;
}
public String[] getAlignment() {
return getAlignment(length1, length2, "", "");
}
/**
* Output the local alignments ending in the (i, j) cell. aligned1 and
* aligned2 are suffixes of final aligned strings found in backtracking
* before calling this function. Note: the strings are replicated at each
* recursive call. Use buffers or stacks to improve efficiency.
*/
protected void printAlignments(int i, int j, String aligned1,
String aligned2) {
// we've reached the starting point, so print the alignments
if ((prevCells[i][j] & DR_ZERO) > 0) {
System.out.println(aligned1);
System.out.println(aligned2);
System.out.println("");
// Note: we could check other directions for longer alignments
// with the same score. we don't do it here.
return;
}
// find out which directions to backtrack
if ((prevCells[i][j] & DR_UP) > 0 && i != 0) {
printAlignments(i - 1, j, str1.get(i - 1) + aligned1, "_"
+ aligned2);
}
if ((prevCells[i][j] & DR_LEFT) > 0 && j != 0) {
printAlignments(i, j - 1, "_" + aligned1, str2.get(j - 1)
+ aligned2);
}
if ((prevCells[i][j] & DR_DIAG) > 0) {
printAlignments(i - 1, j - 1, str1.get(i - 1) + aligned1,
str2.get(j - 1) + aligned2);
}
}
public void printAlignments() {
printAlignments(length1, length2, "", "");
}
/**
* print the dynamic programming matrix
*/
public void printDPMatrix() {
System.out.print(" ");
for (int j = 1; j <= length2; j++)
System.out.print(" " + str2.get(j - 1));
System.out.println();
for (int i = 0; i <= length1; i++) {
if (i > 0)
System.out.print(str1.get(i - 1) + " ");
else
System.out.print(" ");
for (int j = 0; j <= length2; j++) {
System.out.print(dpTable[i][j] / normFactor + " ");
}
System.out.println();
}
}
/**
* print the dynamic programming matrix
*/
public void printDPMatrixInLatexMatrix() {
System.out
.println("\\begin{tikzpicture}[auto,transform shape,scale=1]\n"
+ "\\tikzset{node style ge/.style={circle,draw,font=\\scriptsize,minimum size=4mm,inner sep=0pt}}\n"
+ "\\tikzset{node style ne/.style={draw=none,font=\\scriptsize,minimum size=4mm,inner sep=0pt}}");
System.out
.println("\\matrix (A) [matrix of nodes,nodes = {node style ge},ampersand replacement=\\&,column sep = 5mm,row sep = 5mm] {");
// System.out.print("\\node[node style ne,opacity=0.0]{-};");
// for (int j=1; j<=length2;j++)
// System.out.print ("\\& \\node[node style ne]{"+str2[j-1]+"};");
// System.out.println("\\\\");
for (int i = 0; i <= length1; i++) {
// if (i>0)
// System.out.print ("\\node[node style ne]{"+str2[i-1]+"};");
System.out.print("\\node[fill=black]{" + 0 + "};");
for (int j = 1; j <= length2; j++) {
System.out.print("\\& \\node[fill=black]{" + 0 + "};"); // dpTable[i][j]
}
System.out.println("\\\\");
}
// end of matrix
System.out.println("};");
for (int i = 0; i <= length1; i++) {
for (int j = 0; j <= length2; j++) {
// right arrow
if (j != length2) {
String print = Math.round(scoringMatrix.score(
alphabet.emptyLetter(), str2.get(j))) == 0 ? ""
: ""
+ Math.round(scoringMatrix.score(
alphabet.emptyLetter(), str2.get(j)));
System.out.println("\\draw [draw,->,midway] (A-" + (i + 1)
+ "-" + (j + 1)
+ ".east) -- node[above,font=\\scriptsize]{"
+ print + "} (A-" + (i + 1) + "-" + (j + 2)
+ ".west);");
}
// down arrow
if (i != length1) {
String print = Math.round(scoringMatrix.score(str1.get(i),
alphabet.emptyLetter())) == 0 ? "" : ""
+ Math.round(scoringMatrix.score(str1.get(i),
alphabet.emptyLetter()));
System.out.println("\\draw [draw,->,midway] (A-" + (i + 1)
+ "-" + (j + 1)
+ ".south) -- node[right=0,font=\\scriptsize]{"
+ print + "} (A-" + (i + 2) + "-" + (j + 1)
+ ".north);");
}
if (i != length1 && j != length2) {
String print = Math.round(scoringMatrix.score(str1.get(i),
str2.get(j))) == 0 ? "" : ""
+ Math.round(scoringMatrix.score(str1.get(i),
str2.get(j)));
System.out
.println("\\draw [draw,->,midway] (A-"
+ (i + 1)
+ "-"
+ (j + 1)
+ ".south east) -- node[right=0.4,pos=0.25,font=\\scriptsize]{"
+ print + "} (A-" + (i + 2) + "-" + (j + 2)
+ ".north west);");
}
}
}
for (int i = 0; i < length2; i++) {
System.out.println("\\draw [midway] (A-1-" + (i + 1)
+ ".east) -- node[above=10,pos=0.5,font=\\scriptsize]{"
+ str2.get(i) + "} (A-1-" + (i + 2) + ".west);");
}
for (int i = 0; i < length1; i++) {
System.out.println("\\draw [midway] (A-" + (i + 1)
+ "-1.south) -- node[left=20,pos=0.5,font=\\scriptsize]{"
+ str1.get(i) + "} (A-" + (i + 2) + "-1.north);");
}
System.out.println("\\end{tikzpicture}");
}
@Override
public void setNormFactor(double factor) {
this.normFactor = factor;
}
}