X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FAlignSeq.java;h=d04e897e201c49aaa69098d94611582a954baa14;hb=adca107ff292173bc844ad4a5c6c7eba958ad24a;hp=3ad31884642b1e734930db7aec7088641ab97f01;hpb=853624fb32058cccc544ae7d13af6ad4b0800b6c;p=jalview.git

diff --git a/src/jalview/analysis/AlignSeq.java b/src/jalview/analysis/AlignSeq.java
index 3ad3188..d04e897 100755
--- a/src/jalview/analysis/AlignSeq.java
+++ b/src/jalview/analysis/AlignSeq.java
@@ -20,13 +20,15 @@
  */
 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;
 import jalview.datamodel.Sequence;
 import jalview.datamodel.SequenceI;
-import jalview.schemes.ResidueProperties;
-import jalview.schemes.ScoreMatrix;
 import jalview.util.Comparison;
 import jalview.util.Format;
 import jalview.util.MapList;
@@ -53,17 +55,11 @@ public class AlignSeq
 
   private static final String NEWLINE = System.lineSeparator();
 
-  static String[] dna = { "A", "C", "G", "T", "-" };
+  float[][] score;
 
-  // "C", "T", "A", "G", "-"};
-  static String[] pep = { "A", "R", "N", "D", "C", "Q", "E", "G", "H", "I",
-      "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V", "B", "Z", "X", "-" };
+  float[][] E;
 
-  int[][] score;
-
-  int[][] E;
-
-  int[][] F;
+  float[][] F;
 
   int[][] traceback;
 
@@ -106,7 +102,7 @@ public class AlignSeq
   int count;
 
   /** DOCUMENT ME!! */
-  public int maxscore;
+  public float maxscore;
 
   float pid;
 
@@ -116,31 +112,26 @@ public class AlignSeq
 
   int gapExtend = 20;
 
-  int[][] lookup = ResidueProperties.getBLOSUM62();
+  float[][] lookup;
 
-  String[] intToStr = pep;
-
-  int defInt = 23;
+  int gapIndex = 23;
 
   StringBuffer output = new StringBuffer();
 
-  String type;
+  String type; // AlignSeq.PEP or AlignSeq.DNA
 
-  private int[] charToInt;
+  private ScoreMatrix scoreModel;
 
   /**
    * Creates a new AlignSeq object.
    * 
-   * @param s1
-   *          DOCUMENT ME!
-   * @param s2
-   *          DOCUMENT ME!
-   * @param type
-   *          DOCUMENT ME!
+   * @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)
   {
-    SeqInit(s1, s1.getSequenceAsString(), s2, s2.getSequenceAsString(),
+    seqInit(s1, s1.getSequenceAsString(), s2, s2.getSequenceAsString(),
             type);
   }
 
@@ -157,7 +148,7 @@ public class AlignSeq
   public AlignSeq(SequenceI s1, String string1, SequenceI s2,
           String string2, String type)
   {
-    SeqInit(s1, string1.toUpperCase(), s2, string2.toUpperCase(), type);
+    seqInit(s1, string1.toUpperCase(), s2, string2.toUpperCase(), type);
   }
 
   /**
@@ -165,7 +156,7 @@ public class AlignSeq
    * 
    * @return DOCUMENT ME!
    */
-  public int getMaxScore()
+  public float getMaxScore()
   {
     return maxscore;
   }
@@ -261,26 +252,6 @@ public class AlignSeq
   }
 
   /**
-   * DOCUMENT ME!
-   * 
-   * @return DOCUMENT ME!
-   */
-  public SequenceI getS1()
-  {
-    return s1;
-  }
-
-  /**
-   * DOCUMENT ME!
-   * 
-   * @return DOCUMENT ME!
-   */
-  public SequenceI getS2()
-  {
-    return s2;
-  }
-
-  /**
    * 
    * @return aligned instance of Seq 1
    */
@@ -322,36 +293,13 @@ public class AlignSeq
    * @param type
    *          DNA or PEPTIDE
    */
-  public void SeqInit(SequenceI s1, String string1, SequenceI s2,
+  public void seqInit(SequenceI s1, String string1, SequenceI s2,
           String string2, String type)
   {
     this.s1 = s1;
     this.s2 = s2;
     setDefaultParams(type);
-    SeqInit(string1, string2);
-  }
-
-  /**
-   * Construct score matrix for sequences with custom substitution matrix
-   * 
-   * @param s1
-   *          - sequence 1
-   * @param string1
-   *          - string to use for s1
-   * @param s2
-   *          - sequence 2
-   * @param string2
-   *          - string to use for s2
-   * @param scoreMatrix
-   *          - substitution matrix to use for alignment
-   */
-  public void SeqInit(SequenceI s1, String string1, SequenceI s2,
-          String string2, ScoreMatrix scoreMatrix)
-  {
-    this.s1 = s1;
-    this.s2 = s2;
-    setType(scoreMatrix.isDNA() ? AlignSeq.DNA : AlignSeq.PEP);
-    lookup = scoreMatrix.getMatrix();
+    seqInit(string1, string2);
   }
 
   /**
@@ -361,7 +309,7 @@ public class AlignSeq
    * @param string1
    * @param string2
    */
-  private void SeqInit(String string1, String string2)
+  private void seqInit(String string1, String string2)
   {
     s1str = extractGaps(jalview.util.Comparison.GapChars, string1);
     s2str = extractGaps(jalview.util.Comparison.GapChars, string2);
@@ -374,84 +322,37 @@ public class AlignSeq
       return;
     }
 
-    // System.out.println("lookuip " + rt.freeMemory() + " "+ rt.totalMemory());
     seq1 = new int[s1str.length()];
 
-    // System.out.println("seq1 " + rt.freeMemory() +" " + rt.totalMemory());
     seq2 = new int[s2str.length()];
 
-    // System.out.println("seq2 " + rt.freeMemory() + " " + rt.totalMemory());
-    score = new int[s1str.length()][s2str.length()];
+    score = new float[s1str.length()][s2str.length()];
 
-    // System.out.println("score " + rt.freeMemory() + " " + rt.totalMemory());
-    E = new int[s1str.length()][s2str.length()];
+    E = new float[s1str.length()][s2str.length()];
 
-    // System.out.println("E " + rt.freeMemory() + " " + rt.totalMemory());
-    F = new int[s1str.length()][s2str.length()];
+    F = new float[s1str.length()][s2str.length()];
     traceback = new int[s1str.length()][s2str.length()];
 
-    // System.out.println("F " + rt.freeMemory() + " " + rt.totalMemory());
-    seq1 = stringToInt(s1str, type);
-
-    // System.out.println("seq1 " + rt.freeMemory() + " " + rt.totalMemory());
-    seq2 = stringToInt(s2str, type);
-
-    // System.out.println("Seq2 " + rt.freeMemory() + " " + rt.totalMemory());
-    // long tstart = System.currentTimeMillis();
-    // calcScoreMatrix();
-    // long tend = System.currentTimeMillis();
-    // System.out.println("Time take to calculate score matrix = " +
-    // (tend-tstart) + " ms");
-    // printScoreMatrix(score);
-    // System.out.println();
-    // printScoreMatrix(traceback);
-    // System.out.println();
-    // printScoreMatrix(E);
-    // System.out.println();
-    // /printScoreMatrix(F);
-    // System.out.println();
-    // tstart = System.currentTimeMillis();
-    // traceAlignment();
-    // tend = System.currentTimeMillis();
-    // System.out.println("Time take to traceback alignment = " + (tend-tstart)
-    // + " ms");
-  }
-
-  private void setDefaultParams(String type)
-  {
-    setType(type);
+    seq1 = indexEncode(s1str);
 
-    if (type.equals(AlignSeq.PEP))
-    {
-      lookup = ResidueProperties.getDefaultPeptideMatrix();
-    }
-    else if (type.equals(AlignSeq.DNA))
-    {
-      lookup = ResidueProperties.getDefaultDnaMatrix();
-    }
+    seq2 = indexEncode(s2str);
   }
 
-  private void setType(String type2)
+  private void setDefaultParams(String moleculeType)
   {
-    this.type = type2;
-    if (type.equals(AlignSeq.PEP))
-    {
-      intToStr = pep;
-      charToInt = ResidueProperties.aaIndex;
-      defInt = ResidueProperties.maxProteinIndex;
-    }
-    else if (type.equals(AlignSeq.DNA))
-    {
-      intToStr = dna;
-      charToInt = ResidueProperties.nucleotideIndex;
-      defInt = ResidueProperties.maxNucleotideIndex;
-    }
-    else
+    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[] { type2 }));
+              "error.unknown_type_dna_or_pep",
+              new String[] { moleculeType }));
     }
+
+    type = moleculeType;
+    scoreModel = ScoreModels.getInstance().getDefaultModel(
+            PEP.equals(type));
+    lookup = scoreModel.getMatrix();
+    gapIndex = scoreModel.getGapIndex();
   }
 
   /**
@@ -460,7 +361,7 @@ public class AlignSeq
   public void traceAlignment()
   {
     // Find the maximum score along the rhs or bottom row
-    int max = -9999;
+    float max = -Float.MAX_VALUE;
 
     for (int i = 0; i < seq1.length; i++)
     {
@@ -494,21 +395,17 @@ public class AlignSeq
     aseq1 = new int[seq1.length + seq2.length];
     aseq2 = new int[seq1.length + seq2.length];
 
+    StringBuilder sb1 = new StringBuilder(aseq1.length);
+    StringBuilder sb2 = new StringBuilder(aseq2.length);
+
     count = (seq1.length + seq2.length) - 1;
 
-    while ((i > 0) && (j > 0))
+    while (i > 0 && j > 0)
     {
-      if ((aseq1[count] != defInt) && (i >= 0))
-      {
-        aseq1[count] = seq1[i];
-        astr1 = s1str.charAt(i) + astr1;
-      }
-
-      if ((aseq2[count] != defInt) && (j > 0))
-      {
-        aseq2[count] = seq2[j];
-        astr2 = s2str.charAt(j) + astr2;
-      }
+      aseq1[count] = seq1[i];
+      sb1.append(s1str.charAt(i));
+      aseq2[count] = seq2[j];
+      sb2.append(s2str.charAt(j));
 
       trace = findTrace(i, j);
 
@@ -520,14 +417,14 @@ public class AlignSeq
       else if (trace == 1)
       {
         j--;
-        aseq1[count] = defInt;
-        astr1 = "-" + astr1.substring(1);
+        aseq1[count] = gapIndex;
+        sb1.replace(sb1.length() - 1, sb1.length(), "-");
       }
       else if (trace == -1)
       {
         i--;
-        aseq2[count] = defInt;
-        astr2 = "-" + astr2.substring(1);
+        aseq2[count] = gapIndex;
+        sb2.replace(sb2.length() - 1, sb2.length(), "-");
       }
 
       count--;
@@ -536,17 +433,24 @@ public class AlignSeq
     seq1start = i + 1;
     seq2start = j + 1;
 
-    if (aseq1[count] != defInt)
+    if (aseq1[count] != gapIndex)
     {
       aseq1[count] = seq1[i];
-      astr1 = s1str.charAt(i) + astr1;
+      sb1.append(s1str.charAt(i));
     }
 
-    if (aseq2[count] != defInt)
+    if (aseq2[count] != gapIndex)
     {
       aseq2[count] = seq2[j];
-      astr2 = s2str.charAt(j) + astr2;
+      sb2.append(s2str.charAt(j));
     }
+
+    /*
+     * we built the character strings backwards, so now
+     * reverse them to convert to sequence strings
+     */
+    astr1 = sb1.reverse().toString();
+    astr2 = sb2.reverse().toString();
   }
 
   /**
@@ -599,6 +503,8 @@ public class AlignSeq
             .append(String.valueOf(s2str.length())).append(")")
             .append(NEWLINE).append(NEWLINE);
 
+    ScoreMatrix pam250 = ScoreModels.getInstance().getPam250();
+
     for (int j = 0; j < nochunks; j++)
     {
       // Print the first aligned sequence
@@ -615,22 +521,27 @@ public class AlignSeq
       output.append(NEWLINE);
       output.append(new Format("%" + (maxid) + "s").form(" ")).append(" ");
 
-      // Print out the matching chars
+      /*
+       * Print out the match symbols:
+       * | for exact match (ignoring case)
+       * . if PAM250 score is positive
+       * else a space
+       */
       for (int i = 0; i < len; i++)
       {
         if ((i + (j * len)) < astr1.length())
         {
-          if (astr1.charAt(i + (j * len)) == astr2.charAt(i + (j * len))
-                  && !jalview.util.Comparison.isGap(astr1.charAt(i
-                          + (j * len))))
+          char c1 = astr1.charAt(i + (j * len));
+          char c2 = astr2.charAt(i + (j * len));
+          boolean sameChar = Comparison.isSameResidue(c1, c2, false);
+          if (sameChar && !Comparison.isGap(c1))
           {
             pid++;
             output.append("|");
           }
           else if (type.equals("pep"))
           {
-            if (ResidueProperties.getPAM250(astr1.charAt(i + (j * len)),
-                    astr2.charAt(i + (j * len))) > 0)
+            if (pam250.getPairwiseScore(c1, c2) > 0)
             {
               output.append(".");
             }
@@ -675,46 +586,6 @@ public class AlignSeq
   /**
    * DOCUMENT ME!
    * 
-   * @param mat
-   *          DOCUMENT ME!
-   */
-  public void printScoreMatrix(int[][] mat)
-  {
-    int n = seq1.length;
-    int m = seq2.length;
-
-    for (int i = 0; i < n; i++)
-    {
-      // Print the top sequence
-      if (i == 0)
-      {
-        Format.print(System.out, "%8s", s2str.substring(0, 1));
-
-        for (int jj = 1; jj < m; jj++)
-        {
-          Format.print(System.out, "%5s", s2str.substring(jj, jj + 1));
-        }
-
-        System.out.println();
-      }
-
-      for (int j = 0; j < m; j++)
-      {
-        if (j == 0)
-        {
-          Format.print(System.out, "%3s", s1str.substring(i, i + 1));
-        }
-
-        Format.print(System.out, "%3d ", mat[i][j] / 10);
-      }
-
-      System.out.println();
-    }
-  }
-
-  /**
-   * DOCUMENT ME!
-   * 
    * @param i
    *          DOCUMENT ME!
    * @param j
@@ -725,7 +596,7 @@ public class AlignSeq
   public int findTrace(int i, int j)
   {
     int t = 0;
-    int max = score[i - 1][j - 1] + (lookup[seq1[i]][seq2[j]] * 10);
+    float max = score[i - 1][j - 1] + (lookup[seq1[i]][seq2[j]] * 10);
 
     if (F[i][j] > max)
     {
@@ -840,27 +711,27 @@ public class AlignSeq
   /**
    * DOCUMENT ME!
    * 
-   * @param i1
+   * @param f1
    *          DOCUMENT ME!
-   * @param i2
+   * @param f2
    *          DOCUMENT ME!
-   * @param i3
+   * @param f3
    *          DOCUMENT ME!
    * 
    * @return DOCUMENT ME!
    */
-  public int max(int i1, int i2, int i3)
+  private static float max(float f1, float f2, float f3)
   {
-    int max = i1;
+    float max = f1;
 
-    if (i2 > i1)
+    if (f2 > f1)
     {
-      max = i2;
+      max = f2;
     }
 
-    if (i3 > max)
+    if (f3 > max)
     {
-      max = i3;
+      max = f3;
     }
 
     return max;
@@ -869,65 +740,44 @@ public class AlignSeq
   /**
    * DOCUMENT ME!
    * 
-   * @param i1
+   * @param f1
    *          DOCUMENT ME!
-   * @param i2
+   * @param f2
    *          DOCUMENT ME!
    * 
    * @return DOCUMENT ME!
    */
-  public int max(int i1, int i2)
+  private static float max(float f1, float f2)
   {
-    int max = i1;
+    float max = f1;
 
-    if (i2 > i1)
+    if (f2 > f1)
     {
-      max = i2;
+      max = f2;
     }
 
     return max;
   }
 
   /**
-   * DOCUMENT ME!
+   * Converts the character string to an array of integers which are the
+   * corresponding indices to the characters in the score matrix
    * 
    * @param s
-   *          DOCUMENT ME!
-   * @param type
-   *          DOCUMENT ME!
    * 
-   * @return DOCUMENT ME!
+   * @return
    */
-  public int[] stringToInt(String s, String type)
+  int[] indexEncode(String s)
   {
-    int[] seq1 = new int[s.length()];
+    int[] encoded = new int[s.length()];
 
     for (int i = 0; i < s.length(); i++)
     {
-      // String ss = s.substring(i, i + 1).toUpperCase();
       char c = s.charAt(i);
-      if ('a' <= c && c <= 'z')
-      {
-        // TO UPPERCASE !!!
-        c -= ('a' - 'A');
-      }
-
-      try
-      {
-        seq1[i] = charToInt[c]; // set accordingly from setType
-        if (seq1[i] < 0 || seq1[i] > defInt) // set from setType: 23 for
-                                             // peptides, or 4 for NA.
-        {
-          seq1[i] = defInt;
-        }
-
-      } catch (Exception e)
-      {
-        seq1[i] = defInt;
-      }
+      encoded[i] = scoreModel.getMatrixIndex(c);
     }
 
-    return seq1;
+    return encoded;
   }
 
   /**
@@ -1110,7 +960,7 @@ public class AlignSeq
       {
         SequenceI bestm = null;
         AlignSeq bestaseq = null;
-        int bestscore = 0;
+        float bestscore = 0;
         for (SequenceI msq : al.getSequences())
         {
           AlignSeq aseq = doGlobalNWAlignment(msq, sq, dnaOrProtein);
@@ -1211,6 +1061,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++)
@@ -1251,7 +1103,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])