Implemented least-squares optimisation in ccAnalysis
[jalview.git] / src / jalview / analysis / AlignSeq.java
index 0731173..420aa78 100755 (executable)
@@ -58,11 +58,11 @@ public class AlignSeq
   private static final int MAX_NAME_LENGTH = 30;
 
   //&!
-  //private static final int GAP_OPEN_COST = 120;
-  private static final int GAP_OPEN_COST = 10;
+  private static final int GAP_OPEN_COST = 120;
+  //private static final int GAP_OPEN_COST = 100;
 
-  //private static final int GAP_EXTEND_COST = 20;
-  private static final float GAP_EXTEND_COST = 0.5f;
+  private static final int GAP_EXTEND_COST = 20;
+  //private static final int GAP_EXTEND_COST = 5;
 
   private static final int GAP_INDEX = -1;
 
@@ -127,7 +127,7 @@ public class AlignSeq
 
   public float meanScore;      //needed for PaSiMap
 
-  public float hypotheticMaxScore;     // needed for PaSiMap
+  public int hypotheticMaxScore;       // needed for PaSiMap
 
   int prev = 0;
 
@@ -410,8 +410,8 @@ public class AlignSeq
     int i = maxi;
     int j = maxj;
     int trace;
-    //maxscore = score[i][j] / 10f;
-    maxscore = score[i][j];
+    maxscore = score[i][j] / 10f;
+    //maxscore = score[i][j];
 
     //&! get trailing gaps
     while ((i < seq1.length - 1) || (j < seq2.length - 1))
@@ -664,8 +664,8 @@ public class AlignSeq
     int t = 0;
     float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
             s2str.charAt(j));
-    //float max = score[i - 1][j - 1] + (pairwiseScore * 10);
-    float max = score[i - 1][j - 1] + (pairwiseScore);
+    float max = score[i - 1][j - 1] + (pairwiseScore * 10);
+    //float max = score[i - 1][j - 1] + (pairwiseScore);
 
     if (F[i][j] > max)
     {
@@ -711,8 +711,8 @@ public class AlignSeq
 
     // top left hand element
     score[0][0] = scoreMatrix.getPairwiseScore(s1str.charAt(0),
-            s2str.charAt(0));
-            //s2str.charAt(0)) * 10;
+            //s2str.charAt(0));
+            s2str.charAt(0)) * 10;
     E[0][0] = -GAP_EXTEND_COST;
     F[0][0] = 0;
 
@@ -726,8 +726,8 @@ public class AlignSeq
 
       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(0),
               s2str.charAt(j));
-      //score[0][j] = max(pairwiseScore * 10, -GAP_OPEN_COST,
-      score[0][j] = max(pairwiseScore, -GAP_OPEN_COST,
+      score[0][j] = max(pairwiseScore * 10, -GAP_OPEN_COST,
+      //score[0][j] = max(pairwiseScore, -GAP_OPEN_COST,
               -GAP_EXTEND_COST);
 
       traceback[0][j] = 1;
@@ -742,8 +742,8 @@ public class AlignSeq
 
       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
               s2str.charAt(0));
-      //score[i][0] = max(pairwiseScore * 10, E[i][0], F[i][0]);
-      score[i][0] = max(pairwiseScore, E[i][0], F[i][0]);
+      score[i][0] = max(pairwiseScore * 10, E[i][0], F[i][0]);
+      //score[i][0] = max(pairwiseScore, E[i][0], F[i][0]);
       traceback[i][0] = -1;
     }
 
@@ -759,8 +759,8 @@ public class AlignSeq
 
         float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
                 s2str.charAt(j));
-        //score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore * 10),
-        score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore),
+        score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore * 10),
+        //score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore),
                 E[i][j], F[i][j]);
         traceback[i][j] = findTrace(i, j);
       }
@@ -1215,21 +1215,26 @@ public class AlignSeq
   */
   public void meanScore()
   {
-    int length = (indelfreeAstr1.length() > indelfreeAstr2.length()) ? indelfreeAstr1.length() : indelfreeAstr2.length();
+    //int length = (indelfreeAstr1.length() > indelfreeAstr2.length()) ? indelfreeAstr1.length() : indelfreeAstr2.length();
+    int length = indelfreeAstr1.length();      //both have the same length
+    //create HashMap for counting residues in each sequence
     HashMap<Character, Integer> seq1ResCount = new HashMap<Character, Integer>();
     HashMap<Character, Integer> seq2ResCount = new HashMap<Character, Integer>();
 
-    for (char residue: s1str.toCharArray())
+    // for both sequences (String indelfreeAstr1 or 2) create a key for the residue and add 1 each time its encountered
+    for (char residue: indelfreeAstr1.toCharArray())
     {
       seq1ResCount.putIfAbsent(residue, 0);
       seq1ResCount.replace(residue, seq1ResCount.get(residue) + 1);
     }
-    for (char residue: s2str.toCharArray())
+    for (char residue: indelfreeAstr2.toCharArray())
     {
       seq2ResCount.putIfAbsent(residue, 0);
       seq2ResCount.replace(residue, seq2ResCount.get(residue) + 1);
     }
 
+    // meanscore = for each residue pair get the number of appearance and add (countA * countB * pairwiseScore(AB))
+    // divide the meanscore by the sequence length afterwards
     float _meanscore = 0;
     for (char resA : seq1ResCount.keySet())
     {
@@ -1257,20 +1262,21 @@ public class AlignSeq
   */
   public void hypotheticMaxScore()
   {
-    float _hmsA = 0f;
-    float _hmsB = 0f;
-    for (char residue: s1str.toCharArray())
+    int _hmsA = 0;
+    int _hmsB = 0;
+    for (char residue: indelfreeAstr1.toCharArray())
     {
       _hmsA += scoreMatrix.getPairwiseScore(residue, residue);
     }
-    for (char residue: s2str.toCharArray())
+    for (char residue: indelfreeAstr2.toCharArray())
     {
       _hmsB += scoreMatrix.getPairwiseScore(residue, residue);
     }
-    this.hypotheticMaxScore = (_hmsA < _hmsB) ? _hmsA : _hmsB;
+    this.hypotheticMaxScore = (_hmsA < _hmsB) ? _hmsA : _hmsB; // take the lower self alignment
+
   }
 
-  public float getHypotheticMaxScore()
+  public int getHypotheticMaxScore()
   {
     return this.hypotheticMaxScore;
   }
@@ -1283,7 +1289,7 @@ public class AlignSeq
     int n = astr1.length();    // both have the same length
     for (int i = 0; i < n; i++)
     {
-      if (Character.isLetter(astr1.charAt(i)) && Character.isLetter(astr2.charAt(i)))
+      if (Character.isLetter(astr1.charAt(i)) && Character.isLetter(astr2.charAt(i)))  // if both sequences dont have a gap -> add to indelfreeAstr
       {
        this.indelfreeAstr1 += astr1.charAt(i);
        this.indelfreeAstr2 += astr2.charAt(i);
@@ -1293,11 +1299,13 @@ public class AlignSeq
 
   /**
   * calculates the overall score of the alignment
-  * alignmentScore = sum of all scores - all penalties
-  *
+  * preprescore = sum of all scores - all penalties
+  * if preprescore < 1 ~ alignmentScore = Float.NaN    >
+  * alignmentScore = ((preprescore - meanScore) / (hypotheticMaxScore - meanScore)) * coverage
   */
   public void scoreAlignment() throws RuntimeException
   {
+
     getIndelfreeAstr();
     meanScore();
     hypotheticMaxScore();
@@ -1306,17 +1314,16 @@ public class AlignSeq
     {
       throw new IllegalArgumentException(String.format("hypotheticMaxScore (%8.2f) == meanScore (%8.2f) - division by 0", hypotheticMaxScore, meanScore));
     }
-    int n = (astr1.length() > astr2.length()) ? astr1.length() : astr2.length();
+    //int n = (astr1.length() > astr2.length()) ? astr1.length() : astr2.length();
+    int n = indelfreeAstr1.length();
 
-    System.out.println(String.format("astr1: %c .. %c ; astr2: %c .. %c", astr1.toCharArray()[0], astr1.toCharArray()[astr1.length() - 1], astr2.toCharArray()[0], astr2.toCharArray()[astr2.length() - 1]));
-    
     float score = 0;
     boolean aGapOpen = false;
     boolean bGapOpen = false;
     for (int i = 0; i < n; i++)
     {
-      char char1 = astr1.charAt(i);
-      char char2 = astr2.charAt(i);
+      char char1 = indelfreeAstr1.charAt(i);
+      char char2 = indelfreeAstr2.charAt(i);
       boolean aIsLetter = Character.isLetter(char1);
       boolean bIsLetter = Character.isLetter(char2);
       if (aIsLetter && bIsLetter)      // if pair -> get score
@@ -1328,18 +1335,19 @@ public class AlignSeq
       } else {         // no gap open -> score - gap_open
        score -= GAP_OPEN_COST;
       }
+      // adjust GapOpen status in both sequences
       aGapOpen = (!aIsLetter) ? true : false;
       bGapOpen = (!bIsLetter) ? true : false;
     }
 
-    float preprescore = score;
+    float preprescore = score; // if this score < 1 --> alignment score = Float.NaN
     score = (score - this.meanScore) / (this.hypotheticMaxScore - this.meanScore);
     int[] _max = MiscMath.findMax(new int[]{astr1.replace("-","").length(), astr2.replace("-","").length()});  // {index of max, max}
-    float coverage = (float) indelfreeAstr1.length() / (float) _max[1];
-    float prescore = score;
+    float coverage = (float) n / (float) _max[1];      // indelfreeAstr length / longest sequence length
+    float prescore = score;    // only debug
     score *= coverage;
 
-    System.out.println(String.format("prepre-score: %f, pre-score: %f, longlength: %d\nscore: %f, mean: %f, max: %f", preprescore, prescore, indelfreeAstr1.length(), score, this.meanScore, this.hypotheticMaxScore));
+    System.out.println(String.format("prepre-score: %f, pre-score: %f, longlength: %d\nscore: %f, mean: %f, max: %d", preprescore, prescore, _max[1], score, this.meanScore, this.hypotheticMaxScore));
     this.alignmentScore = (preprescore < 1) ? Float.NaN : score;
   }
 }