JAL-2617 unit tests for CDS + stop codon, CDS with start phase
[jalview.git] / src / jalview / analysis / AlignSeq.java
index ceca6d6..34a21e6 100755 (executable)
  */
 package jalview.analysis;
 
+import jalview.analysis.scoremodels.PIDModel;
 import jalview.analysis.scoremodels.ScoreMatrix;
 import jalview.analysis.scoremodels.ScoreModels;
+import jalview.analysis.scoremodels.SimilarityParams;
 import jalview.datamodel.AlignmentAnnotation;
 import jalview.datamodel.AlignmentI;
 import jalview.datamodel.Mapping;
@@ -110,22 +112,23 @@ public class AlignSeq
 
   int gapExtend = 20;
 
-  float[][] lookup;
-
-  int gapIndex = 23;
-
   StringBuffer output = new StringBuffer();
 
   String type; // AlignSeq.PEP or AlignSeq.DNA
 
-  private ScoreMatrix scoreModel;
+  private ScoreMatrix scoreMatrix;
+
+  private static final int GAP_INDEX = -1;
 
   /**
    * Creates a new AlignSeq object.
    * 
-   * @param s1 first sequence for alignment
-   * @param s2 second sequence for alignment
-   * @param type molecule type, either AlignSeq.PEP or AlignSeq.DNA
+   * @param s1
+   *          first sequence for alignment
+   * @param s2
+   *          second sequence for alignment
+   * @param type
+   *          molecule type, either AlignSeq.PEP or AlignSeq.DNA
    */
   public AlignSeq(SequenceI s1, SequenceI s2, String type)
   {
@@ -258,8 +261,8 @@ public class AlignSeq
     SequenceI alSeq1 = new Sequence(s1.getName(), getAStr1());
     alSeq1.setStart(s1.getStart() + getSeq1Start() - 1);
     alSeq1.setEnd(s1.getStart() + getSeq1End() - 1);
-    alSeq1.setDatasetSequence(s1.getDatasetSequence() == null ? s1 : s1
-            .getDatasetSequence());
+    alSeq1.setDatasetSequence(
+            s1.getDatasetSequence() == null ? s1 : s1.getDatasetSequence());
     return alSeq1;
   }
 
@@ -272,8 +275,8 @@ public class AlignSeq
     SequenceI alSeq2 = new Sequence(s2.getName(), getAStr2());
     alSeq2.setStart(s2.getStart() + getSeq2Start() - 1);
     alSeq2.setEnd(s2.getStart() + getSeq2End() - 1);
-    alSeq2.setDatasetSequence(s2.getDatasetSequence() == null ? s2 : s2
-            .getDatasetSequence());
+    alSeq2.setDatasetSequence(
+            s2.getDatasetSequence() == null ? s2 : s2.getDatasetSequence());
     return alSeq2;
   }
 
@@ -314,16 +317,12 @@ public class AlignSeq
 
     if (s1str.length() == 0 || s2str.length() == 0)
     {
-      output.append("ALL GAPS: "
-              + (s1str.length() == 0 ? s1.getName() : " ")
-              + (s2str.length() == 0 ? s2.getName() : ""));
+      output.append(
+              "ALL GAPS: " + (s1str.length() == 0 ? s1.getName() : " ")
+                      + (s2str.length() == 0 ? s2.getName() : ""));
       return;
     }
 
-    seq1 = new int[s1str.length()];
-
-    seq2 = new int[s2str.length()];
-
     score = new float[s1str.length()][s2str.length()];
 
     E = new float[s1str.length()][s2str.length()];
@@ -341,16 +340,14 @@ public class AlignSeq
     if (!PEP.equals(moleculeType) && !DNA.equals(moleculeType))
     {
       output.append("Wrong type = dna or pep only");
-      throw new Error(MessageManager.formatMessage(
-              "error.unknown_type_dna_or_pep",
-              new String[] { moleculeType }));
+      throw new Error(MessageManager
+              .formatMessage("error.unknown_type_dna_or_pep", new String[]
+              { moleculeType }));
     }
 
     type = moleculeType;
-    scoreModel = ScoreModels.getInstance().getDefaultModel(
-            PEP.equals(type));
-    lookup = scoreModel.getMatrix();
-    gapIndex = scoreModel.getMatrixIndex(' ');
+    scoreMatrix = ScoreModels.getInstance()
+            .getDefaultModel(PEP.equals(type));
   }
 
   /**
@@ -415,13 +412,13 @@ public class AlignSeq
       else if (trace == 1)
       {
         j--;
-        aseq1[count] = gapIndex;
+        aseq1[count] = GAP_INDEX;
         sb1.replace(sb1.length() - 1, sb1.length(), "-");
       }
       else if (trace == -1)
       {
         i--;
-        aseq2[count] = gapIndex;
+        aseq2[count] = GAP_INDEX;
         sb2.replace(sb2.length() - 1, sb2.length(), "-");
       }
 
@@ -431,13 +428,13 @@ public class AlignSeq
     seq1start = i + 1;
     seq2start = j + 1;
 
-    if (aseq1[count] != gapIndex)
+    if (aseq1[count] != GAP_INDEX)
     {
       aseq1[count] = seq1[i];
       sb1.append(s1str.charAt(i));
     }
 
-    if (aseq2[count] != gapIndex)
+    if (aseq2[count] != GAP_INDEX)
     {
       aseq2[count] = seq2[j];
       sb2.append(s2str.charAt(j));
@@ -594,7 +591,10 @@ public class AlignSeq
   public int findTrace(int i, int j)
   {
     int t = 0;
-    float max = score[i - 1][j - 1] + (lookup[seq1[i]][seq2[j]] * 10);
+    // float pairwiseScore = lookup[seq1[i]][seq2[j]];
+    float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
+            s2str.charAt(j));
+    float max = score[i - 1][j - 1] + (pairwiseScore * 10);
 
     if (F[i][j] > max)
     {
@@ -638,7 +638,8 @@ public class AlignSeq
     int m = seq2.length;
 
     // top left hand element
-    score[0][0] = lookup[seq1[0]][seq2[0]] * 10;
+    score[0][0] = scoreMatrix.getPairwiseScore(s1str.charAt(0),
+            s2str.charAt(0)) * 10;
     E[0][0] = -gapExtend;
     F[0][0] = 0;
 
@@ -649,7 +650,9 @@ public class AlignSeq
       E[0][j] = max(score[0][j - 1] - gapOpen, E[0][j - 1] - gapExtend);
       F[0][j] = -gapExtend;
 
-      score[0][j] = max(lookup[seq1[0]][seq2[j]] * 10, -gapOpen, -gapExtend);
+      float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(0),
+              s2str.charAt(j));
+      score[0][j] = max(pairwiseScore * 10, -gapOpen, -gapExtend);
 
       traceback[0][j] = 1;
     }
@@ -660,7 +663,9 @@ public class AlignSeq
       E[i][0] = -gapOpen;
       F[i][0] = max(score[i - 1][0] - gapOpen, F[i - 1][0] - gapExtend);
 
-      score[i][0] = max(lookup[seq1[i]][seq2[0]] * 10, E[i][0], F[i][0]);
+      float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
+              s2str.charAt(0));
+      score[i][0] = max(pairwiseScore * 10, E[i][0], F[i][0]);
       traceback[i][0] = -1;
     }
 
@@ -672,8 +677,10 @@ public class AlignSeq
         E[i][j] = max(score[i][j - 1] - gapOpen, E[i][j - 1] - gapExtend);
         F[i][j] = max(score[i - 1][j] - gapOpen, F[i - 1][j] - gapExtend);
 
-        score[i][j] = max(score[i - 1][j - 1]
-                + (lookup[seq1[i]][seq2[j]] * 10), E[i][j], F[i][j]);
+        float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
+                s2str.charAt(j));
+        score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore * 10),
+                E[i][j], F[i][j]);
         traceback[i][j] = findTrace(i, j);
       }
     }
@@ -772,7 +779,7 @@ public class AlignSeq
     for (int i = 0; i < s.length(); i++)
     {
       char c = s.charAt(i);
-      encoded[i] = scoreModel.getMatrixIndex(c);
+      encoded[i] = scoreMatrix.getMatrixIndex(c);
     }
 
     return encoded;
@@ -795,7 +802,7 @@ public class AlignSeq
   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
           int psize)
   {
-    // TODO method dosen't seem to be referenced anywhere delete??
+    // TODO method doesn't seem to be referenced anywhere delete??
     int max = -1000;
     int min = 1000;
 
@@ -859,7 +866,8 @@ public class AlignSeq
    */
   public jalview.datamodel.Mapping getMappingFromS1(boolean allowmismatch)
   {
-    ArrayList<Integer> as1 = new ArrayList<Integer>(), as2 = new ArrayList<Integer>();
+    ArrayList<Integer> as1 = new ArrayList<Integer>(),
+            as2 = new ArrayList<Integer>();
     int pdbpos = s2.getStart() + getSeq2Start() - 2;
     int alignpos = s1.getStart() + getSeq1Start() - 2;
     int lp2 = pdbpos - 3, lp1 = alignpos - 3;
@@ -903,8 +911,8 @@ public class AlignSeq
     }
     // construct range pairs
 
-    int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)], mapseq2 = new int[as2
-            .size() + (lastmatch ? 1 : 0)];
+    int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)],
+            mapseq2 = new int[as2.size() + (lastmatch ? 1 : 0)];
     int i = 0;
     for (Integer ip : as1)
     {
@@ -947,7 +955,8 @@ public class AlignSeq
           List<SequenceI> ochains, AlignmentI al, String dnaOrProtein,
           boolean removeOldAnnots)
   {
-    List<SequenceI> orig = new ArrayList<SequenceI>(), repl = new ArrayList<SequenceI>();
+    List<SequenceI> orig = new ArrayList<SequenceI>(),
+            repl = new ArrayList<SequenceI>();
     List<AlignSeq> aligs = new ArrayList<AlignSeq>();
     if (al != null && al.getHeight() > 0)
     {
@@ -969,8 +978,8 @@ public class AlignSeq
             bestm = msq;
           }
         }
-        System.out.println("Best Score for " + (matches.size() + 1) + " :"
-                + bestscore);
+        // System.out.println("Best Score for " + (matches.size() + 1) + " :"
+        // + bestscore);
         matches.add(bestm);
         aligns.add(bestaseq);
         al.deleteSequence(bestm);
@@ -1059,6 +1068,8 @@ public class AlignSeq
 
     // long start = System.currentTimeMillis();
 
+    SimilarityParams pidParams = new SimilarityParams(true, true, true,
+            true);
     float pid;
     String seqi, seqj;
     for (int i = 0; i < height; i++)
@@ -1099,7 +1110,7 @@ public class AlignSeq
             seqj = ug;
           }
         }
-        pid = Comparison.PID(seqi, seqj);
+        pid = (float) PIDModel.computePID(seqi, seqj, pidParams);
 
         // use real sequence length rather than string length
         if (lngth[j] < lngth[i])