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