src/jalview/analysis/AlignSeq.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.analysis;
  22
  23 import java.util.Locale;
  24
  25 import jalview.analysis.scoremodels.PIDModel;
  26 import jalview.analysis.scoremodels.ScoreMatrix;
  27 import jalview.analysis.scoremodels.ScoreModels;
  28 import jalview.analysis.scoremodels.SimilarityParams;
  29 import jalview.datamodel.AlignmentAnnotation;
  30 import jalview.datamodel.AlignmentI;
  31 import jalview.datamodel.Mapping;
  32 import jalview.datamodel.Sequence;
  33 import jalview.datamodel.SequenceI;
  34 import jalview.math.MiscMath;
  35 import jalview.util.Comparison;
  36 import jalview.util.Format;
  37 import jalview.util.MapList;
  38 import jalview.util.MessageManager;
  39
  40 import java.awt.Color;
  41 import java.awt.Graphics;
  42 import java.io.PrintStream;
  43 import java.lang.IllegalArgumentException;
  44 import java.util.ArrayList;
  45 import java.util.Arrays;
  46 import java.util.HashMap;
  47 import java.util.List;
  48 import java.util.StringTokenizer;
  49
  50 /**
  51  *
  52  *
  53  * @author $author$
  54  * @version $Revision$
  55  */
  56 public class AlignSeq
  57 {
  58   private static final int MAX_NAME_LENGTH = 30;
  59
  60   //&!
  61   //private static final int GAP_OPEN_COST = 120;
  62   private static final int GAP_OPEN_COST = 100;
  63
  64   //private static final int GAP_EXTEND_COST = 20;
  65   private static final int GAP_EXTEND_COST = 5;
  66
  67   private static final int GAP_INDEX = -1;
  68
  69   public static final String PEP = "pep";
  70
  71   public static final String DNA = "dna";
  72
  73   private static final String NEWLINE = System.lineSeparator();
  74
  75   float[][] score;
  76
  77   float alignmentScore;
  78
  79   float[][] E;
  80
  81   float[][] F;
  82
  83   int[][] traceback; // todo is this actually used?
  84
  85   int[] seq1;
  86
  87   int[] seq2;
  88
  89   SequenceI s1;
  90
  91   SequenceI s2;
  92
  93   public String s1str;
  94
  95   public String s2str;
  96
  97   int maxi;
  98
  99   int maxj;
 100
 101   int[] aseq1;
 102
 103   int[] aseq2;
 104
 105   public String astr1 = "";
 106
 107   public String astr2 = "";
 108
 109   public String indelfreeAstr1 = "";
 110
 111   public String indelfreeAstr2 = "";
 112
 113   /** DOCUMENT ME!! */
 114   public int seq1start;
 115
 116   /** DOCUMENT ME!! */
 117   public int seq1end;
 118
 119   /** DOCUMENT ME!! */
 120   public int seq2start;
 121
 122   public int seq2end;
 123
 124   int count;
 125
 126   public float maxscore;
 127
 128   public float meanScore;       //needed for PaSiMap
 129
 130   public int hypotheticMaxScore;        // needed for PaSiMap
 131
 132   int prev = 0;
 133
 134   StringBuffer output = new StringBuffer();
 135
 136   String type; // AlignSeq.PEP or AlignSeq.DNA
 137
 138   private ScoreMatrix scoreMatrix;
 139
 140   /**
 141    * Creates a new AlignSeq object.
 142    *
 143    * @param s1
 144    *          first sequence for alignment
 145    * @param s2
 146    *          second sequence for alignment
 147    * @param type
 148    *          molecule type, either AlignSeq.PEP or AlignSeq.DNA
 149    */
 150   public AlignSeq(SequenceI s1, SequenceI s2, String type)
 151   {
 152     seqInit(s1, s1.getSequenceAsString(), s2, s2.getSequenceAsString(),
 153             type);
 154   }
 155
 156   /**
 157    * Creates a new AlignSeq object.
 158    *
 159    * @param s1
 160    *          DOCUMENT ME!
 161    * @param s2
 162    *          DOCUMENT ME!
 163    * @param type
 164    *          DOCUMENT ME!
 165    */
 166   public AlignSeq(SequenceI s1, String string1, SequenceI s2,
 167           String string2, String type)
 168   {
 169     seqInit(s1, string1.toUpperCase(Locale.ROOT), s2,
 170             string2.toUpperCase(Locale.ROOT), type);
 171   }
 172
 173   /**
 174    * DOCUMENT ME!
 175    *
 176    * @return DOCUMENT ME!
 177    */
 178   public float getMaxScore()
 179   {
 180     return maxscore;
 181   }
 182
 183   /**
 184   * returns the overall score of the alignment
 185   *
 186   * @return
 187   */
 188   public float getAlignmentScore()
 189   {
 190     return alignmentScore;
 191   }
 192
 193   /**
 194    * DOCUMENT ME!
 195    *
 196    * @return DOCUMENT ME!
 197    */
 198   public int getSeq2Start()
 199   {
 200     return seq2start;
 201   }
 202
 203   /**
 204    * DOCUMENT ME!
 205    *
 206    * @return DOCUMENT ME!
 207    */
 208   public int getSeq2End()
 209   {
 210     return seq2end;
 211   }
 212
 213   /**
 214    * DOCUMENT ME!
 215    *
 216    * @return DOCUMENT ME!
 217    */
 218   public int getSeq1Start()
 219   {
 220     return seq1start;
 221   }
 222
 223   /**
 224    * DOCUMENT ME!
 225    *
 226    * @return DOCUMENT ME!
 227    */
 228   public int getSeq1End()
 229   {
 230     return seq1end;
 231   }
 232
 233   /**
 234    * DOCUMENT ME!
 235    *
 236    * @return DOCUMENT ME!
 237    */
 238   public String getOutput()
 239   {
 240     return output.toString();
 241   }
 242
 243   /**
 244    * DOCUMENT ME!
 245    *
 246    * @return DOCUMENT ME!
 247    */
 248   public String getAStr1()
 249   {
 250     return astr1;
 251   }
 252
 253   /**
 254    * DOCUMENT ME!
 255    *
 256    * @return DOCUMENT ME!
 257    */
 258   public String getAStr2()
 259   {
 260     return astr2;
 261   }
 262
 263   /**
 264    * DOCUMENT ME!
 265    *
 266    * @return DOCUMENT ME!
 267    */
 268   public int[] getASeq1()
 269   {
 270     return aseq1;
 271   }
 272
 273   /**
 274    * DOCUMENT ME!
 275    *
 276    * @return DOCUMENT ME!
 277    */
 278   public int[] getASeq2()
 279   {
 280     return aseq2;
 281   }
 282
 283   /**
 284    *
 285    * @return aligned instance of Seq 1
 286    */
 287   public SequenceI getAlignedSeq1()
 288   {
 289     SequenceI alSeq1 = new Sequence(s1.getName(), getAStr1());
 290     alSeq1.setStart(s1.getStart() + getSeq1Start() - 1);
 291     alSeq1.setEnd(s1.getStart() + getSeq1End() - 1);
 292     alSeq1.setDatasetSequence(
 293             s1.getDatasetSequence() == null ? s1 : s1.getDatasetSequence());
 294     return alSeq1;
 295   }
 296
 297   /**
 298    *
 299    * @return aligned instance of Seq 2
 300    */
 301   public SequenceI getAlignedSeq2()
 302   {
 303     SequenceI alSeq2 = new Sequence(s2.getName(), getAStr2());
 304     alSeq2.setStart(s2.getStart() + getSeq2Start() - 1);
 305     alSeq2.setEnd(s2.getStart() + getSeq2End() - 1);
 306     alSeq2.setDatasetSequence(
 307             s2.getDatasetSequence() == null ? s2 : s2.getDatasetSequence());
 308     return alSeq2;
 309   }
 310
 311   /**
 312    * Construct score matrix for sequences with standard DNA or PEPTIDE matrix
 313    *
 314    * @param s1
 315    *          - sequence 1
 316    * @param string1
 317    *          - string to use for s1
 318    * @param s2
 319    *          - sequence 2
 320    * @param string2
 321    *          - string to use for s2
 322    * @param type
 323    *          DNA or PEPTIDE
 324    */
 325   public void seqInit(SequenceI s1, String string1, SequenceI s2,
 326           String string2, String type)
 327   {
 328     this.s1 = s1;
 329     this.s2 = s2;
 330     setDefaultParams(type);
 331     seqInit(string1, string2);
 332   }
 333
 334   /**
 335    * construct score matrix for string1 and string2 (after removing any existing
 336    * gaps
 337    *
 338    * @param string1
 339    * @param string2
 340    */
 341   private void seqInit(String string1, String string2)
 342   {
 343     s1str = extractGaps(jalview.util.Comparison.GapChars, string1);
 344     s2str = extractGaps(jalview.util.Comparison.GapChars, string2);
 345
 346     if (s1str.length() == 0 || s2str.length() == 0)
 347     {
 348       output.append(
 349               "ALL GAPS: " + (s1str.length() == 0 ? s1.getName() : " ")
 350                       + (s2str.length() == 0 ? s2.getName() : ""));
 351       return;
 352     }
 353
 354     score = new float[s1str.length()][s2str.length()];
 355
 356     E = new float[s1str.length()][s2str.length()];
 357
 358     F = new float[s1str.length()][s2str.length()];
 359     traceback = new int[s1str.length()][s2str.length()];
 360
 361     seq1 = indexEncode(s1str);
 362
 363     seq2 = indexEncode(s2str);
 364   }
 365
 366   private void setDefaultParams(String moleculeType)
 367   {
 368     if (!PEP.equals(moleculeType) && !DNA.equals(moleculeType))
 369     {
 370       output.append("Wrong type = dna or pep only");
 371       throw new Error(MessageManager
 372               .formatMessage("error.unknown_type_dna_or_pep", new String[]
 373               { moleculeType }));
 374     }
 375
 376     type = moleculeType;
 377     scoreMatrix = ScoreModels.getInstance()
 378             .getDefaultModel(PEP.equals(type));
 379   }
 380
 381   /**
 382    * DOCUMENT ME!
 383    */
 384   public void traceAlignment()
 385   {
 386     // Find the maximum score along the rhs or bottom row
 387     float max = -Float.MAX_VALUE;
 388
 389     for (int i = 0; i < seq1.length; i++)
 390     {
 391       if (score[i][seq2.length - 1] > max)
 392       {
 393         max = score[i][seq2.length - 1];
 394         maxi = i;
 395         maxj = seq2.length - 1;
 396       }
 397     }
 398
 399     for (int j = 0; j < seq2.length; j++)
 400     {
 401       if (score[seq1.length - 1][j] > max)
 402       {
 403         max = score[seq1.length - 1][j];
 404         maxi = seq1.length - 1;
 405         maxj = j;
 406       }
 407     }
 408
 409     int i = maxi;
 410     int j = maxj;
 411     int trace;
 412     maxscore = score[i][j] / 10f;
 413
 414
 415     aseq1 = new int[seq1.length + seq2.length];
 416     aseq2 = new int[seq1.length + seq2.length];
 417
 418     StringBuilder sb1 = new StringBuilder(aseq1.length);
 419     StringBuilder sb2 = new StringBuilder(aseq2.length);
 420
 421     count = (seq1.length + seq2.length) - 1;
 422
 423
 424     while (i > 0 && j > 0)
 425     {
 426       aseq1[count] = seq1[i];
 427       sb1.append(s1str.charAt(i));
 428       aseq2[count] = seq2[j];
 429       sb2.append(s2str.charAt(j));
 430
 431       trace = findTrace(i, j);
 432
 433       if (trace == 0)
 434       {
 435         i--;
 436         j--;
 437       }
 438       else if (trace == 1)
 439       {
 440         j--;
 441         aseq1[count] = GAP_INDEX;
 442         sb1.replace(sb1.length() - 1, sb1.length(), "-");
 443       }
 444       else if (trace == -1)
 445       {
 446         i--;
 447         aseq2[count] = GAP_INDEX;
 448         sb2.replace(sb2.length() - 1, sb2.length(), "-");
 449       }
 450
 451       count--;
 452     }
 453
 454     seq1start = i + 1;
 455     seq2start = j + 1;
 456
 457     if (aseq1[count] != GAP_INDEX)
 458     {
 459       aseq1[count] = seq1[i];
 460       sb1.append(s1str.charAt(i));
 461     }
 462
 463     if (aseq2[count] != GAP_INDEX)
 464     {
 465       aseq2[count] = seq2[j];
 466       sb2.append(s2str.charAt(j));
 467     }
 468
 469
 470     /*
 471      * we built the character strings backwards, so now
 472      * reverse them to convert to sequence strings
 473      */
 474     astr1 = sb1.reverse().toString();
 475     astr2 = sb2.reverse().toString();
 476   }
 477
 478   /**
 479    * DOCUMENT ME!
 480    */
 481   public void traceAlignmentWithEndGaps()
 482   {
 483     // Find the maximum score along the rhs or bottom row
 484     float max = -Float.MAX_VALUE;
 485
 486     for (int i = 0; i < seq1.length; i++)
 487     {
 488       if (score[i][seq2.length - 1] > max)
 489       {
 490         max = score[i][seq2.length - 1];
 491         maxi = i;
 492         maxj = seq2.length - 1;
 493       }
 494     }
 495
 496     for (int j = 0; j < seq2.length; j++)
 497     {
 498       if (score[seq1.length - 1][j] > max)
 499       {
 500         max = score[seq1.length - 1][j];
 501         maxi = seq1.length - 1;
 502         maxj = j;
 503       }
 504     }
 505
 506     int i = maxi;
 507     int j = maxj;
 508     int trace;
 509     maxscore = score[i][j] / 10f;
 510
 511     //prepare trailing gaps
 512     while ((i < seq1.length - 1) || (j < seq2.length - 1))
 513     {
 514       i++;
 515       j++;
 516     }
 517     seq1end = i + 1;
 518     seq2end = j + 1;
 519
 520
 521     aseq1 = new int[seq1.length + seq2.length];
 522     aseq2 = new int[seq1.length + seq2.length];
 523
 524     StringBuilder sb1 = new StringBuilder(aseq1.length);
 525     StringBuilder sb2 = new StringBuilder(aseq2.length);
 526
 527     count = (seq1.length + seq2.length) - 1;
 528
 529     //get trailing gaps
 530     while ((i >= seq1.length) || (j >= seq2.length))
 531     {
 532       if (i >= seq1.length)
 533       {
 534         aseq1[count] = GAP_INDEX;
 535         sb1.append("-");
 536         aseq2[count] = seq2[j];
 537         sb2.append(s2str.charAt(j));
 538       } else if (j >= seq2.length) {
 539         aseq1[count] = seq1[i];
 540         sb1.append(s1str.charAt(i));
 541         aseq2[count] = GAP_INDEX;
 542         sb2.append("-");
 543       }
 544       i--;
 545       j--;
 546     }
 547
 548
 549     while (i > 0 && j > 0)
 550     {
 551       aseq1[count] = seq1[i];
 552       sb1.append(s1str.charAt(i));
 553       aseq2[count] = seq2[j];
 554       sb2.append(s2str.charAt(j));
 555
 556       trace = findTrace(i, j);
 557
 558       if (trace == 0)
 559       {
 560         i--;
 561         j--;
 562       }
 563       else if (trace == 1)
 564       {
 565         j--;
 566         aseq1[count] = GAP_INDEX;
 567         sb1.replace(sb1.length() - 1, sb1.length(), "-");
 568       }
 569       else if (trace == -1)
 570       {
 571         i--;
 572         aseq2[count] = GAP_INDEX;
 573         sb2.replace(sb2.length() - 1, sb2.length(), "-");
 574       }
 575
 576       count--;
 577     }
 578
 579     seq1start = i + 1;
 580     seq2start = j + 1;
 581
 582     if (aseq1[count] != GAP_INDEX)
 583     {
 584       aseq1[count] = seq1[i];
 585       sb1.append(s1str.charAt(i));
 586     }
 587
 588     if (aseq2[count] != GAP_INDEX)
 589     {
 590       aseq2[count] = seq2[j];
 591       sb2.append(s2str.charAt(j));
 592     }
 593
 594     //get initial gaps
 595     while (j > 0 || i > 0)
 596     {
 597       if (j > 0)
 598       {
 599         j--;
 600         sb1.append("-");
 601         sb2.append(s2str.charAt(j));
 602       } else if (i > 0) {
 603         i--;
 604         sb1.append(s1str.charAt(i));
 605         sb2.append("-");
 606       }
 607     }
 608
 609     /*
 610      * we built the character strings backwards, so now
 611      * reverse them to convert to sequence strings
 612      */
 613     astr1 = sb1.reverse().toString();
 614     astr2 = sb2.reverse().toString();
 615   }
 616
 617   /**
 618    * DOCUMENT ME!
 619    */
 620   public void printAlignment(PrintStream os)
 621   {
 622     // TODO: Use original sequence characters rather than re-translated
 623     // characters in output
 624     // Find the biggest id length for formatting purposes
 625     String s1id = getAlignedSeq1().getDisplayId(true);
 626     String s2id = getAlignedSeq2().getDisplayId(true);
 627     int nameLength = Math.max(s1id.length(), s2id.length());
 628     if (nameLength > MAX_NAME_LENGTH)
 629     {
 630       int truncateBy = nameLength - MAX_NAME_LENGTH;
 631       nameLength = MAX_NAME_LENGTH;
 632       // JAL-527 - truncate the sequence ids
 633       if (s1id.length() > nameLength)
 634       {
 635         int slashPos = s1id.lastIndexOf('/');
 636         s1id = s1id.substring(0, slashPos - truncateBy)
 637                 + s1id.substring(slashPos);
 638       }
 639       if (s2id.length() > nameLength)
 640       {
 641         int slashPos = s2id.lastIndexOf('/');
 642         s2id = s2id.substring(0, slashPos - truncateBy)
 643                 + s2id.substring(slashPos);
 644       }
 645     }
 646     int len = 72 - nameLength - 1;
 647     int nochunks = ((aseq1.length - count) / len)
 648             + ((aseq1.length - count) % len > 0 ? 1 : 0);
 649     float pid = 0f;
 650
 651     output.append("Score = ").append(score[maxi][maxj]).append(NEWLINE);
 652     output.append("Length of alignment = ")
 653             .append(String.valueOf(aseq1.length - count)).append(NEWLINE);
 654     output.append("Sequence ");
 655     Format nameFormat = new Format("%" + nameLength + "s");
 656     output.append(nameFormat.form(s1id));
 657     output.append(" (Sequence length = ")
 658             .append(String.valueOf(s1str.length())).append(")")
 659             .append(NEWLINE);
 660     output.append("Sequence ");
 661     output.append(nameFormat.form(s2id));
 662     output.append(" (Sequence length = ")
 663             .append(String.valueOf(s2str.length())).append(")")
 664             .append(NEWLINE).append(NEWLINE);
 665
 666     ScoreMatrix pam250 = ScoreModels.getInstance().getPam250();
 667
 668     for (int j = 0; j < nochunks; j++)
 669     {
 670       // Print the first aligned sequence
 671       output.append(nameFormat.form(s1id)).append(" ");
 672
 673       for (int i = 0; i < len; i++)
 674       {
 675         if ((i + (j * len)) < astr1.length())
 676         {
 677           output.append(astr1.charAt(i + (j * len)));
 678         }
 679       }
 680
 681       output.append(NEWLINE);
 682       output.append(nameFormat.form(" ")).append(" ");
 683
 684       /*
 685        * Print out the match symbols:
 686        * | for exact match (ignoring case)
 687        * . if PAM250 score is positive
 688        * else a space
 689        */
 690       for (int i = 0; i < len; i++)
 691       {
 692         if ((i + (j * len)) < astr1.length())
 693         {
 694           char c1 = astr1.charAt(i + (j * len));
 695           char c2 = astr2.charAt(i + (j * len));
 696           boolean sameChar = Comparison.isSameResidue(c1, c2, false);
 697           if (sameChar && !Comparison.isGap(c1))
 698           {
 699             pid++;
 700             output.append("|");
 701           }
 702           else if (PEP.equals(type))
 703           {
 704             if (pam250.getPairwiseScore(c1, c2) > 0)
 705             {
 706               output.append(".");
 707             }
 708             else
 709             {
 710               output.append(" ");
 711             }
 712           }
 713           else
 714           {
 715             output.append(" ");
 716           }
 717         }
 718       }
 719
 720       // Now print the second aligned sequence
 721       output = output.append(NEWLINE);
 722       output = output.append(nameFormat.form(s2id)).append(" ");
 723
 724       for (int i = 0; i < len; i++)
 725       {
 726         if ((i + (j * len)) < astr2.length())
 727         {
 728           output.append(astr2.charAt(i + (j * len)));
 729         }
 730       }
 731
 732       output.append(NEWLINE).append(NEWLINE);
 733     }
 734
 735     pid = pid / (aseq1.length - count) * 100;
 736     output.append(new Format("Percentage ID = %3.2f\n").form(pid));
 737     output.append(NEWLINE);
 738     try
 739     {
 740       os.print(output.toString());
 741     } catch (Exception ex)
 742     {
 743     }
 744   }
 745
 746   /**
 747    * DOCUMENT ME!
 748    *
 749    * @param i
 750    *          DOCUMENT ME!
 751    * @param j
 752    *          DOCUMENT ME!
 753    *
 754    * @return DOCUMENT ME!
 755    */
 756   public int findTrace(int i, int j)
 757   {
 758     int t = 0;
 759     float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 760             s2str.charAt(j));
 761     float max = score[i - 1][j - 1] + (pairwiseScore * 10);
 762
 763     if (F[i][j] > max)
 764     {
 765       max = F[i][j];
 766       t = -1;
 767     }
 768     else if (F[i][j] == max)
 769     {
 770       if (prev == -1)
 771       {
 772         max = F[i][j];
 773         t = -1;
 774       }
 775     }
 776
 777     if (E[i][j] >= max)
 778     {
 779       max = E[i][j];
 780       t = 1;
 781     }
 782     else if (E[i][j] == max)
 783     {
 784       if (prev == 1)
 785       {
 786         max = E[i][j];
 787         t = 1;
 788       }
 789     }
 790
 791     prev = t;
 792
 793     return t;
 794   }
 795
 796   /**
 797    * DOCUMENT ME!
 798    */
 799   public void calcScoreMatrix()
 800   {
 801     int n = seq1.length;
 802     int m = seq2.length;
 803
 804     // top left hand element
 805     score[0][0] = scoreMatrix.getPairwiseScore(s1str.charAt(0),
 806             s2str.charAt(0)) * 10;
 807     E[0][0] = -GAP_EXTEND_COST;
 808     F[0][0] = 0;
 809
 810     // Calculate the top row first
 811     for (int j = 1; j < m; j++)
 812     {
 813       // What should these values be? 0 maybe
 814       E[0][j] = max(score[0][j - 1] - GAP_OPEN_COST,
 815               E[0][j - 1] - GAP_EXTEND_COST);
 816       F[0][j] = -GAP_EXTEND_COST;
 817
 818       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(0),
 819               s2str.charAt(j));
 820       score[0][j] = max(pairwiseScore * 10, -GAP_OPEN_COST,
 821               -GAP_EXTEND_COST);
 822
 823       traceback[0][j] = 1;
 824     }
 825
 826     // Now do the left hand column
 827     for (int i = 1; i < n; i++)
 828     {
 829       E[i][0] = -GAP_OPEN_COST;
 830       F[i][0] = max(score[i - 1][0] - GAP_OPEN_COST,
 831               F[i - 1][0] - GAP_EXTEND_COST);
 832
 833       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 834               s2str.charAt(0));
 835       score[i][0] = max(pairwiseScore * 10, E[i][0], F[i][0]);
 836       traceback[i][0] = -1;
 837     }
 838
 839     // Now do all the other rows
 840     for (int i = 1; i < n; i++)
 841     {
 842       for (int j = 1; j < m; j++)
 843       {
 844         E[i][j] = max(score[i][j - 1] - GAP_OPEN_COST,
 845                 E[i][j - 1] - GAP_EXTEND_COST);
 846         F[i][j] = max(score[i - 1][j] - GAP_OPEN_COST,
 847                 F[i - 1][j] - GAP_EXTEND_COST);
 848
 849         float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 850                 s2str.charAt(j));
 851         score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore * 10),
 852                 E[i][j], F[i][j]);
 853         traceback[i][j] = findTrace(i, j);
 854       }
 855     }
 856   }
 857
 858   /**
 859    * Returns the given sequence with all of the given gap characters removed.
 860    *
 861    * @param gapChars
 862    *          a string of characters to be treated as gaps
 863    * @param seq
 864    *          the input sequence
 865    *
 866    * @return
 867    */
 868   public static String extractGaps(String gapChars, String seq)
 869   {
 870     if (gapChars == null || seq == null)
 871     {
 872       return null;
 873     }
 874     StringTokenizer str = new StringTokenizer(seq, gapChars);
 875     StringBuilder newString = new StringBuilder(seq.length());
 876
 877     while (str.hasMoreTokens())
 878     {
 879       newString.append(str.nextToken());
 880     }
 881
 882     return newString.toString();
 883   }
 884
 885   /**
 886    * DOCUMENT ME!
 887    *
 888    * @param f1
 889    *          DOCUMENT ME!
 890    * @param f2
 891    *          DOCUMENT ME!
 892    * @param f3
 893    *          DOCUMENT ME!
 894    *
 895    * @return DOCUMENT ME!
 896    */
 897   private static float max(float f1, float f2, float f3)
 898   {
 899     float max = f1;
 900
 901     if (f2 > f1)
 902     {
 903       max = f2;
 904     }
 905
 906     if (f3 > max)
 907     {
 908       max = f3;
 909     }
 910
 911     return max;
 912   }
 913
 914   /**
 915    * DOCUMENT ME!
 916    *
 917    * @param f1
 918    *          DOCUMENT ME!
 919    * @param f2
 920    *          DOCUMENT ME!
 921    *
 922    * @return DOCUMENT ME!
 923    */
 924   private static float max(float f1, float f2)
 925   {
 926     float max = f1;
 927
 928     if (f2 > f1)
 929     {
 930       max = f2;
 931     }
 932
 933     return max;
 934   }
 935
 936   /**
 937    * Converts the character string to an array of integers which are the
 938    * corresponding indices to the characters in the score matrix
 939    *
 940    * @param s
 941    *
 942    * @return
 943    */
 944   int[] indexEncode(String s)
 945   {
 946     int[] encoded = new int[s.length()];
 947
 948     for (int i = 0; i < s.length(); i++)
 949     {
 950       char c = s.charAt(i);
 951       encoded[i] = scoreMatrix.getMatrixIndex(c);
 952     }
 953
 954     return encoded;
 955   }
 956
 957   /**
 958    * DOCUMENT ME!
 959    *
 960    * @param g
 961    *          DOCUMENT ME!
 962    * @param mat
 963    *          DOCUMENT ME!
 964    * @param n
 965    *          DOCUMENT ME!
 966    * @param m
 967    *          DOCUMENT ME!
 968    * @param psize
 969    *          DOCUMENT ME!
 970    */
 971   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
 972           int psize)
 973   {
 974     // TODO method doesn't seem to be referenced anywhere delete??
 975     int max = -1000;
 976     int min = 1000;
 977
 978     for (int i = 0; i < n; i++)
 979     {
 980       for (int j = 0; j < m; j++)
 981       {
 982         if (mat[i][j] >= max)
 983         {
 984           max = mat[i][j];
 985         }
 986
 987         if (mat[i][j] <= min)
 988         {
 989           min = mat[i][j];
 990         }
 991       }
 992     }
 993
 994     System.out.println(max + " " + min);
 995
 996     for (int i = 0; i < n; i++)
 997     {
 998       for (int j = 0; j < m; j++)
 999       {
1000         int x = psize * i;
1001         int y = psize * j;
1002
1003         // System.out.println(mat[i][j]);
1004         float score = (float) (mat[i][j] - min) / (float) (max - min);
1005         g.setColor(new Color(score, 0, 0));
1006         g.fillRect(x, y, psize, psize);
1007
1008         // System.out.println(x + " " + y + " " + score);
1009       }
1010     }
1011   }
1012
1013   /**
1014    * Compute a globally optimal needleman and wunsch alignment between two
1015    * sequences
1016    *
1017    * @param s1
1018    * @param s2
1019    * @param type
1020    *          AlignSeq.DNA or AlignSeq.PEP
1021    */
1022   public static AlignSeq doGlobalNWAlignment(SequenceI s1, SequenceI s2,
1023           String type)
1024   {
1025     AlignSeq as = new AlignSeq(s1, s2, type);
1026
1027     as.calcScoreMatrix();
1028     as.traceAlignment();
1029     return as;
1030   }
1031
1032   /**
1033    *
1034    * @return mapping from positions in S1 to corresponding positions in S2
1035    */
1036   public jalview.datamodel.Mapping getMappingFromS1(boolean allowmismatch)
1037   {
1038     ArrayList<Integer> as1 = new ArrayList<Integer>(),
1039             as2 = new ArrayList<Integer>();
1040     int pdbpos = s2.getStart() + getSeq2Start() - 2;
1041     int alignpos = s1.getStart() + getSeq1Start() - 2;
1042     int lp2 = pdbpos - 3, lp1 = alignpos - 3;
1043     boolean lastmatch = false;
1044     // and now trace the alignment onto the atom set.
1045     for (int i = 0; i < astr1.length(); i++)
1046     {
1047       char c1 = astr1.charAt(i), c2 = astr2.charAt(i);
1048       if (c1 != '-')
1049       {
1050         alignpos++;
1051       }
1052
1053       if (c2 != '-')
1054       {
1055         pdbpos++;
1056       }
1057
1058       // ignore case differences
1059       if (allowmismatch || (c1 == c2) || (Math.abs(c2-c1)==('a'-'A')))
1060       {
1061         // extend mapping interval
1062         if (lp1 + 1 != alignpos || lp2 + 1 != pdbpos)
1063         {
1064           as1.add(Integer.valueOf(alignpos));
1065           as2.add(Integer.valueOf(pdbpos));
1066         }
1067         lastmatch = true;
1068         lp1 = alignpos;
1069         lp2 = pdbpos;
1070       }
1071       else
1072       {
1073         // extend mapping interval
1074         if (lastmatch)
1075         {
1076           as1.add(Integer.valueOf(lp1));
1077           as2.add(Integer.valueOf(lp2));
1078         }
1079         lastmatch = false;
1080       }
1081     }
1082     // construct range pairs
1083
1084     int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)],
1085             mapseq2 = new int[as2.size() + (lastmatch ? 1 : 0)];
1086     int i = 0;
1087     for (Integer ip : as1)
1088     {
1089       mapseq1[i++] = ip;
1090     }
1091     ;
1092     i = 0;
1093     for (Integer ip : as2)
1094     {
1095       mapseq2[i++] = ip;
1096     }
1097     ;
1098     if (lastmatch)
1099     {
1100       mapseq1[mapseq1.length - 1] = alignpos;
1101       mapseq2[mapseq2.length - 1] = pdbpos;
1102     }
1103     MapList map = new MapList(mapseq1, mapseq2, 1, 1);
1104
1105     jalview.datamodel.Mapping mapping = new Mapping(map);
1106     mapping.setTo(s2);
1107     return mapping;
1108   }
1109
1110   /**
1111    * matches ochains against al and populates seqs with the best match between
1112    * each ochain and the set in al
1113    *
1114    * @param ochains
1115    * @param al
1116    * @param dnaOrProtein
1117    * @param removeOldAnnots
1118    *          when true, old annotation is cleared before new annotation
1119    *          transferred
1120    * @return List<List<SequenceI> originals, List<SequenceI> replacement,
1121    *         List<AlignSeq> alignment between each>
1122    */
1123   public static List<List<? extends Object>> replaceMatchingSeqsWith(
1124           List<SequenceI> seqs, List<AlignmentAnnotation> annotations,
1125           List<SequenceI> ochains, AlignmentI al, String dnaOrProtein,
1126           boolean removeOldAnnots)
1127   {
1128     List<SequenceI> orig = new ArrayList<SequenceI>(),
1129             repl = new ArrayList<SequenceI>();
1130     List<AlignSeq> aligs = new ArrayList<AlignSeq>();
1131     if (al != null && al.getHeight() > 0)
1132     {
1133       ArrayList<SequenceI> matches = new ArrayList<SequenceI>();
1134       ArrayList<AlignSeq> aligns = new ArrayList<AlignSeq>();
1135
1136       for (SequenceI sq : ochains)
1137       {
1138         SequenceI bestm = null;
1139         AlignSeq bestaseq = null;
1140         float bestscore = 0;
1141         for (SequenceI msq : al.getSequences())
1142         {
1143           AlignSeq aseq = doGlobalNWAlignment(msq, sq, dnaOrProtein);
1144           if (bestm == null || aseq.getMaxScore() > bestscore)
1145           {
1146             bestscore = aseq.getMaxScore();
1147             bestaseq = aseq;
1148             bestm = msq;
1149           }
1150         }
1151         // System.out.println("Best Score for " + (matches.size() + 1) + " :"
1152         // + bestscore);
1153         matches.add(bestm);
1154         aligns.add(bestaseq);
1155         al.deleteSequence(bestm);
1156       }
1157       for (int p = 0, pSize = seqs.size(); p < pSize; p++)
1158       {
1159         SequenceI sq, sp = seqs.get(p);
1160         int q;
1161         if ((q = ochains.indexOf(sp)) > -1)
1162         {
1163           seqs.set(p, sq = matches.get(q));
1164           orig.add(sp);
1165           repl.add(sq);
1166           sq.setName(sp.getName());
1167           sq.setDescription(sp.getDescription());
1168           Mapping sp2sq;
1169           sq.transferAnnotation(sp,
1170                   sp2sq = aligns.get(q).getMappingFromS1(false));
1171           aligs.add(aligns.get(q));
1172           int inspos = -1;
1173           for (int ap = 0; ap < annotations.size();)
1174           {
1175             if (annotations.get(ap).sequenceRef == sp)
1176             {
1177               if (inspos == -1)
1178               {
1179                 inspos = ap;
1180               }
1181               if (removeOldAnnots)
1182               {
1183                 annotations.remove(ap);
1184               }
1185               else
1186               {
1187                 AlignmentAnnotation alan = annotations.remove(ap);
1188                 alan.liftOver(sq, sp2sq);
1189                 alan.setSequenceRef(sq);
1190                 sq.addAlignmentAnnotation(alan);
1191               }
1192             }
1193             else
1194             {
1195               ap++;
1196             }
1197           }
1198           if (sq.getAnnotation() != null && sq.getAnnotation().length > 0)
1199           {
1200             annotations.addAll(inspos == -1 ? annotations.size() : inspos,
1201                     Arrays.asList(sq.getAnnotation()));
1202           }
1203         }
1204       }
1205     }
1206     return Arrays.asList(orig, repl, aligs);
1207   }
1208
1209   /**
1210    * compute the PID vector used by the redundancy filter.
1211    *
1212    * @param originalSequences
1213    *          - sequences in alignment that are to filtered
1214    * @param omitHidden
1215    *          - null or strings to be analysed (typically, visible portion of
1216    *          each sequence in alignment)
1217    * @param start
1218    *          - first column in window for calculation
1219    * @param end
1220    *          - last column in window for calculation
1221    * @param ungapped
1222    *          - if true then use ungapped sequence to compute PID
1223    * @return vector containing maximum PID for i-th sequence and any sequences
1224    *         longer than that seuqence
1225    */
1226   public static float[] computeRedundancyMatrix(
1227           SequenceI[] originalSequences, String[] omitHidden, int start,
1228           int end, boolean ungapped)
1229   {
1230     int height = originalSequences.length;
1231     float[] redundancy = new float[height];
1232     int[] lngth = new int[height];
1233     for (int i = 0; i < height; i++)
1234     {
1235       redundancy[i] = 0f;
1236       lngth[i] = -1;
1237     }
1238
1239     // long start = System.currentTimeMillis();
1240
1241     SimilarityParams pidParams = new SimilarityParams(true, true, true,
1242             true);
1243     float pid;
1244     String seqi, seqj;
1245     for (int i = 0; i < height; i++)
1246     {
1247
1248       for (int j = 0; j < i; j++)
1249       {
1250         if (i == j)
1251         {
1252           continue;
1253         }
1254
1255         if (omitHidden == null)
1256         {
1257           seqi = originalSequences[i].getSequenceAsString(start, end);
1258           seqj = originalSequences[j].getSequenceAsString(start, end);
1259         }
1260         else
1261         {
1262           seqi = omitHidden[i];
1263           seqj = omitHidden[j];
1264         }
1265         if (lngth[i] == -1)
1266         {
1267           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
1268           lngth[i] = ug.length();
1269           if (ungapped)
1270           {
1271             seqi = ug;
1272           }
1273         }
1274         if (lngth[j] == -1)
1275         {
1276           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
1277           lngth[j] = ug.length();
1278           if (ungapped)
1279           {
1280             seqj = ug;
1281           }
1282         }
1283         pid = (float) PIDModel.computePID(seqi, seqj, pidParams);
1284
1285         // use real sequence length rather than string length
1286         if (lngth[j] < lngth[i])
1287         {
1288           redundancy[j] = Math.max(pid, redundancy[j]);
1289         }
1290         else
1291         {
1292           redundancy[i] = Math.max(pid, redundancy[i]);
1293         }
1294
1295       }
1296     }
1297     return redundancy;
1298   }
1299
1300   /**
1301   * calculate the mean score of the alignment
1302   * mean score is equal to the score of an alignmenet of two sequences with randomly shuffled AA sequence composited of the same AA as the two original sequences
1303   *
1304   */
1305   public void meanScore()
1306   {
1307     //int length = (indelfreeAstr1.length() > indelfreeAstr2.length()) ? indelfreeAstr1.length() : indelfreeAstr2.length();
1308     int length = indelfreeAstr1.length();       //both have the same length
1309     //create HashMap for counting residues in each sequence
1310     HashMap<Character, Integer> seq1ResCount = new HashMap<Character, Integer>();
1311     HashMap<Character, Integer> seq2ResCount = new HashMap<Character, Integer>();
1312
1313     // for both sequences (String indelfreeAstr1 or 2) create a key for the residue and add 1 each time its encountered
1314     for (char residue: indelfreeAstr1.toCharArray())
1315     {
1316       seq1ResCount.putIfAbsent(residue, 0);
1317       seq1ResCount.replace(residue, seq1ResCount.get(residue) + 1);
1318     }
1319     for (char residue: indelfreeAstr2.toCharArray())
1320     {
1321       seq2ResCount.putIfAbsent(residue, 0);
1322       seq2ResCount.replace(residue, seq2ResCount.get(residue) + 1);
1323     }
1324
1325     // meanscore = for each residue pair get the number of appearance and add (countA * countB * pairwiseScore(AB))
1326     // divide the meanscore by the sequence length afterwards
1327     float _meanscore = 0;
1328     for (char resA : seq1ResCount.keySet())
1329     {
1330       for (char resB : seq2ResCount.keySet())
1331       {
1332         int countA = seq1ResCount.get(resA);
1333         int countB = seq2ResCount.get(resB);
1334
1335         float scoreAB = scoreMatrix.getPairwiseScore(resA, resB);
1336
1337         _meanscore += countA *  countB * scoreAB;
1338       }
1339     }
1340     _meanscore /= length;
1341     this.meanScore = _meanscore;
1342   }
1343
1344   public float getMeanScore()
1345   {
1346     return this.meanScore;
1347   }
1348
1349   /**
1350   * calculate the hypothetic max score using the self-alignment of the sequences
1351   */
1352   public void hypotheticMaxScore()
1353   {
1354     int _hmsA = 0;
1355     int _hmsB = 0;
1356     for (char residue: indelfreeAstr1.toCharArray())
1357     {
1358       _hmsA += scoreMatrix.getPairwiseScore(residue, residue);
1359     }
1360     for (char residue: indelfreeAstr2.toCharArray())
1361     {
1362       _hmsB += scoreMatrix.getPairwiseScore(residue, residue);
1363     }
1364     this.hypotheticMaxScore = (_hmsA < _hmsB) ? _hmsA : _hmsB;  // take the lower self alignment
1365
1366   }
1367
1368   public int getHypotheticMaxScore()
1369   {
1370     return this.hypotheticMaxScore;
1371   }
1372
1373   /**
1374   * create strings based of astr1 and astr2 but without gaps
1375   */
1376   public void getIndelfreeAstr()
1377   {
1378     int n = astr1.length();     // both have the same length
1379     for (int i = 0; i < n; i++)
1380     {
1381       if (Character.isLetter(astr1.charAt(i)) && Character.isLetter(astr2.charAt(i)))   // if both sequences dont have a gap -> add to indelfreeAstr
1382       {
1383         this.indelfreeAstr1 += astr1.charAt(i);
1384         this.indelfreeAstr2 += astr2.charAt(i);
1385       }
1386     }
1387   }
1388
1389   /**
1390   * calculates the overall score of the alignment
1391   * preprescore = sum of all scores - all penalties
1392   * if preprescore < 1 ~ alignmentScore = Float.NaN     >
1393   * alignmentScore = ((preprescore - meanScore) / (hypotheticMaxScore - meanScore)) * coverage
1394   */
1395   public void scoreAlignment() throws RuntimeException
1396   {
1397
1398     getIndelfreeAstr();
1399     meanScore();
1400     hypotheticMaxScore();
1401     // cannot calculate score because denominator would be zero
1402     if (this.hypotheticMaxScore == this.meanScore)
1403     {
1404       throw new IllegalArgumentException(String.format("hypotheticMaxScore (%8.2f) == meanScore (%8.2f) - division by 0", hypotheticMaxScore, meanScore));
1405     }
1406     //int n = (astr1.length() > astr2.length()) ? astr1.length() : astr2.length();
1407     int n = indelfreeAstr1.length();
1408
1409     float score = 0;
1410     boolean aGapOpen = false;
1411     boolean bGapOpen = false;
1412     for (int i = 0; i < n; i++)
1413     {
1414       char char1 = indelfreeAstr1.charAt(i);
1415       char char2 = indelfreeAstr2.charAt(i);
1416       boolean aIsLetter = Character.isLetter(char1);
1417       boolean bIsLetter = Character.isLetter(char2);
1418       if (aIsLetter && bIsLetter)       // if pair -> get score
1419       {
1420         score += scoreMatrix.getPairwiseScore(char1, char2);
1421       } else if (!aIsLetter && !bIsLetter) {    // both are gap -> skip
1422       } else if ((!aIsLetter && aGapOpen) || (!bIsLetter && bGapOpen)) {        // one side gapopen -> score - gap_extend
1423         score -= GAP_EXTEND_COST;
1424       } else {          // no gap open -> score - gap_open
1425         score -= GAP_OPEN_COST;
1426       }
1427       // adjust GapOpen status in both sequences
1428       aGapOpen = (!aIsLetter) ? true : false;
1429       bGapOpen = (!bIsLetter) ? true : false;
1430     }
1431
1432     float preprescore = score;  // if this score < 1 --> alignment score = Float.NaN
1433     score = (score - this.meanScore) / (this.hypotheticMaxScore - this.meanScore);
1434     int[] _max = MiscMath.findMax(new int[]{astr1.replace("-","").length(), astr2.replace("-","").length()});   // {index of max, max}
1435     float coverage = (float) n / (float) _max[1];       // indelfreeAstr length / longest sequence length
1436     float prescore = score;     // only debug
1437     score *= coverage;
1438
1439     System.out.println(String.format("prepre-score: %f, pre-score: %f, longlength: %d\nscore: %1.16f, mean: %f, max: %d", preprescore, prescore, _max[1], score, this.meanScore, this.hypotheticMaxScore));
1440     float minScore = 1f;
1441     this.alignmentScore = (preprescore < minScore) ? Float.NaN : score;
1442   }
1443 }