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   //&! not / 10
 385   public void traceAlignment()
 386   {
 387     // Find the maximum score along the rhs or bottom row
 388     float max = -Float.MAX_VALUE;
 389
 390     for (int i = 0; i < seq1.length; i++)
 391     {
 392       if (score[i][seq2.length - 1] > max)
 393       {
 394         max = score[i][seq2.length - 1];
 395         maxi = i;
 396         maxj = seq2.length - 1;
 397       }
 398     }
 399
 400     for (int j = 0; j < seq2.length; j++)
 401     {
 402       if (score[seq1.length - 1][j] > max)
 403       {
 404         max = score[seq1.length - 1][j];
 405         maxi = seq1.length - 1;
 406         maxj = j;
 407       }
 408     }
 409
 410     int i = maxi;
 411     int j = maxj;
 412     int trace;
 413     maxscore = score[i][j] / 10f;
 414     //maxscore = score[i][j];
 415
 416     //&! get trailing gaps
 417     while ((i < seq1.length - 1) || (j < seq2.length - 1))
 418     {
 419       i++;
 420       j++;
 421     }
 422     seq1end = i + 1;
 423     seq2end = j + 1;
 424
 425
 426     aseq1 = new int[seq1.length + seq2.length];
 427     aseq2 = new int[seq1.length + seq2.length];
 428
 429     StringBuilder sb1 = new StringBuilder(aseq1.length);
 430     StringBuilder sb2 = new StringBuilder(aseq2.length);
 431
 432     count = (seq1.length + seq2.length) - 1;
 433
 434     //&! get trailing gaps
 435     while ((i >= seq1.length) || (j >= seq2.length))
 436     {
 437       if (i >= seq1.length)
 438       {
 439         aseq1[count] = GAP_INDEX;
 440         sb1.append("-");
 441         aseq2[count] = seq2[j];
 442         sb2.append(s2str.charAt(j));
 443       } else if (j >= seq2.length) {
 444         aseq1[count] = seq1[i];
 445         sb1.append(s1str.charAt(i));
 446         aseq2[count] = GAP_INDEX;
 447         sb2.append("-");
 448       }
 449       i--;
 450       j--;
 451     }
 452
 453
 454     while (i > 0 && j > 0)
 455     {
 456       aseq1[count] = seq1[i];
 457       sb1.append(s1str.charAt(i));
 458       aseq2[count] = seq2[j];
 459       sb2.append(s2str.charAt(j));
 460
 461       trace = findTrace(i, j);
 462
 463       if (trace == 0)
 464       {
 465         i--;
 466         j--;
 467       }
 468       else if (trace == 1)
 469       {
 470         j--;
 471         aseq1[count] = GAP_INDEX;
 472         sb1.replace(sb1.length() - 1, sb1.length(), "-");
 473       }
 474       else if (trace == -1)
 475       {
 476         i--;
 477         aseq2[count] = GAP_INDEX;
 478         sb2.replace(sb2.length() - 1, sb2.length(), "-");
 479       }
 480
 481       count--;
 482     }
 483
 484     seq1start = i + 1;
 485     seq2start = j + 1;
 486
 487     if (aseq1[count] != GAP_INDEX)
 488     {
 489       aseq1[count] = seq1[i];
 490       sb1.append(s1str.charAt(i));
 491     }
 492
 493     if (aseq2[count] != GAP_INDEX)
 494     {
 495       aseq2[count] = seq2[j];
 496       sb2.append(s2str.charAt(j));
 497     }
 498
 499     //&! get initial gaps
 500     while (j > 0 || i > 0)
 501     {
 502       if (j > 0)
 503       {
 504         sb1.append("-");
 505         sb2.append(s2str.charAt(j));
 506         j--;
 507       } else if (i > 0) {
 508         sb1.append(s1str.charAt(i));
 509         sb2.append("-");
 510         i--;
 511       }
 512     }
 513
 514     /*
 515      * we built the character strings backwards, so now
 516      * reverse them to convert to sequence strings
 517      */
 518     astr1 = sb1.reverse().toString();
 519     astr2 = sb2.reverse().toString();
 520   }
 521
 522   /**
 523    * DOCUMENT ME!
 524    */
 525   public void printAlignment(PrintStream os)
 526   {
 527     // TODO: Use original sequence characters rather than re-translated
 528     // characters in output
 529     // Find the biggest id length for formatting purposes
 530     String s1id = getAlignedSeq1().getDisplayId(true);
 531     String s2id = getAlignedSeq2().getDisplayId(true);
 532     int nameLength = Math.max(s1id.length(), s2id.length());
 533     if (nameLength > MAX_NAME_LENGTH)
 534     {
 535       int truncateBy = nameLength - MAX_NAME_LENGTH;
 536       nameLength = MAX_NAME_LENGTH;
 537       // JAL-527 - truncate the sequence ids
 538       if (s1id.length() > nameLength)
 539       {
 540         int slashPos = s1id.lastIndexOf('/');
 541         s1id = s1id.substring(0, slashPos - truncateBy)
 542                 + s1id.substring(slashPos);
 543       }
 544       if (s2id.length() > nameLength)
 545       {
 546         int slashPos = s2id.lastIndexOf('/');
 547         s2id = s2id.substring(0, slashPos - truncateBy)
 548                 + s2id.substring(slashPos);
 549       }
 550     }
 551     int len = 72 - nameLength - 1;
 552     int nochunks = ((aseq1.length - count) / len)
 553             + ((aseq1.length - count) % len > 0 ? 1 : 0);
 554     float pid = 0f;
 555
 556     output.append("Score = ").append(score[maxi][maxj]).append(NEWLINE);
 557     output.append("Length of alignment = ")
 558             .append(String.valueOf(aseq1.length - count)).append(NEWLINE);
 559     output.append("Sequence ");
 560     Format nameFormat = new Format("%" + nameLength + "s");
 561     output.append(nameFormat.form(s1id));
 562     output.append(" (Sequence length = ")
 563             .append(String.valueOf(s1str.length())).append(")")
 564             .append(NEWLINE);
 565     output.append("Sequence ");
 566     output.append(nameFormat.form(s2id));
 567     output.append(" (Sequence length = ")
 568             .append(String.valueOf(s2str.length())).append(")")
 569             .append(NEWLINE).append(NEWLINE);
 570
 571     ScoreMatrix pam250 = ScoreModels.getInstance().getPam250();
 572
 573     for (int j = 0; j < nochunks; j++)
 574     {
 575       // Print the first aligned sequence
 576       output.append(nameFormat.form(s1id)).append(" ");
 577
 578       for (int i = 0; i < len; i++)
 579       {
 580         if ((i + (j * len)) < astr1.length())
 581         {
 582           output.append(astr1.charAt(i + (j * len)));
 583         }
 584       }
 585
 586       output.append(NEWLINE);
 587       output.append(nameFormat.form(" ")).append(" ");
 588
 589       /*
 590        * Print out the match symbols:
 591        * | for exact match (ignoring case)
 592        * . if PAM250 score is positive
 593        * else a space
 594        */
 595       for (int i = 0; i < len; i++)
 596       {
 597         if ((i + (j * len)) < astr1.length())
 598         {
 599           char c1 = astr1.charAt(i + (j * len));
 600           char c2 = astr2.charAt(i + (j * len));
 601           boolean sameChar = Comparison.isSameResidue(c1, c2, false);
 602           if (sameChar && !Comparison.isGap(c1))
 603           {
 604             pid++;
 605             output.append("|");
 606           }
 607           else if (PEP.equals(type))
 608           {
 609             if (pam250.getPairwiseScore(c1, c2) > 0)
 610             {
 611               output.append(".");
 612             }
 613             else
 614             {
 615               output.append(" ");
 616             }
 617           }
 618           else
 619           {
 620             output.append(" ");
 621           }
 622         }
 623       }
 624
 625       // Now print the second aligned sequence
 626       output = output.append(NEWLINE);
 627       output = output.append(nameFormat.form(s2id)).append(" ");
 628
 629       for (int i = 0; i < len; i++)
 630       {
 631         if ((i + (j * len)) < astr2.length())
 632         {
 633           output.append(astr2.charAt(i + (j * len)));
 634         }
 635       }
 636
 637       output.append(NEWLINE).append(NEWLINE);
 638     }
 639
 640     pid = pid / (aseq1.length - count) * 100;
 641     output.append(new Format("Percentage ID = %3.2f\n").form(pid));
 642     output.append(NEWLINE);
 643     try
 644     {
 645       os.print(output.toString());
 646     } catch (Exception ex)
 647     {
 648     }
 649   }
 650
 651   /**
 652    * DOCUMENT ME!
 653    *
 654    * @param i
 655    *          DOCUMENT ME!
 656    * @param j
 657    *          DOCUMENT ME!
 658    *
 659    * @return DOCUMENT ME!
 660    */
 661   //&! not * 10
 662   public int findTrace(int i, int j)
 663   {
 664     int t = 0;
 665     float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 666             s2str.charAt(j));
 667     float max = score[i - 1][j - 1] + (pairwiseScore * 10);
 668     //float max = score[i - 1][j - 1] + (pairwiseScore);
 669
 670     if (F[i][j] > max)
 671     {
 672       max = F[i][j];
 673       t = -1;
 674     }
 675     else if (F[i][j] == max)
 676     {
 677       if (prev == -1)
 678       {
 679         max = F[i][j];
 680         t = -1;
 681       }
 682     }
 683
 684     if (E[i][j] >= max)
 685     {
 686       max = E[i][j];
 687       t = 1;
 688     }
 689     else if (E[i][j] == max)
 690     {
 691       if (prev == 1)
 692       {
 693         max = E[i][j];
 694         t = 1;
 695       }
 696     }
 697
 698     prev = t;
 699
 700     return t;
 701   }
 702
 703   /**
 704    * DOCUMENT ME!
 705    */
 706   //&! not * 10
 707   public void calcScoreMatrix()
 708   {
 709     int n = seq1.length;
 710     int m = seq2.length;
 711
 712     // top left hand element
 713     score[0][0] = scoreMatrix.getPairwiseScore(s1str.charAt(0),
 714             //s2str.charAt(0));
 715             s2str.charAt(0)) * 10;
 716     E[0][0] = -GAP_EXTEND_COST;
 717     F[0][0] = 0;
 718
 719     // Calculate the top row first
 720     for (int j = 1; j < m; j++)
 721     {
 722       // What should these values be? 0 maybe
 723       E[0][j] = max(score[0][j - 1] - GAP_OPEN_COST,
 724               E[0][j - 1] - GAP_EXTEND_COST);
 725       F[0][j] = -GAP_EXTEND_COST;
 726
 727       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(0),
 728               s2str.charAt(j));
 729       score[0][j] = max(pairwiseScore * 10, -GAP_OPEN_COST,
 730       //score[0][j] = max(pairwiseScore, -GAP_OPEN_COST,
 731               -GAP_EXTEND_COST);
 732
 733       traceback[0][j] = 1;
 734     }
 735
 736     // Now do the left hand column
 737     for (int i = 1; i < n; i++)
 738     {
 739       E[i][0] = -GAP_OPEN_COST;
 740       F[i][0] = max(score[i - 1][0] - GAP_OPEN_COST,
 741               F[i - 1][0] - GAP_EXTEND_COST);
 742
 743       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 744               s2str.charAt(0));
 745       score[i][0] = max(pairwiseScore * 10, E[i][0], F[i][0]);
 746       //score[i][0] = max(pairwiseScore, E[i][0], F[i][0]);
 747       traceback[i][0] = -1;
 748     }
 749
 750     // Now do all the other rows
 751     for (int i = 1; i < n; i++)
 752     {
 753       for (int j = 1; j < m; j++)
 754       {
 755         E[i][j] = max(score[i][j - 1] - GAP_OPEN_COST,
 756                 E[i][j - 1] - GAP_EXTEND_COST);
 757         F[i][j] = max(score[i - 1][j] - GAP_OPEN_COST,
 758                 F[i - 1][j] - GAP_EXTEND_COST);
 759
 760         float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 761                 s2str.charAt(j));
 762         score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore * 10),
 763         //score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore),
 764                 E[i][j], F[i][j]);
 765         traceback[i][j] = findTrace(i, j);
 766       }
 767     }
 768   }
 769
 770   /**
 771    * Returns the given sequence with all of the given gap characters removed.
 772    *
 773    * @param gapChars
 774    *          a string of characters to be treated as gaps
 775    * @param seq
 776    *          the input sequence
 777    *
 778    * @return
 779    */
 780   public static String extractGaps(String gapChars, String seq)
 781   {
 782     if (gapChars == null || seq == null)
 783     {
 784       return null;
 785     }
 786     StringTokenizer str = new StringTokenizer(seq, gapChars);
 787     StringBuilder newString = new StringBuilder(seq.length());
 788
 789     while (str.hasMoreTokens())
 790     {
 791       newString.append(str.nextToken());
 792     }
 793
 794     return newString.toString();
 795   }
 796
 797   /**
 798    * DOCUMENT ME!
 799    *
 800    * @param f1
 801    *          DOCUMENT ME!
 802    * @param f2
 803    *          DOCUMENT ME!
 804    * @param f3
 805    *          DOCUMENT ME!
 806    *
 807    * @return DOCUMENT ME!
 808    */
 809   private static float max(float f1, float f2, float f3)
 810   {
 811     float max = f1;
 812
 813     if (f2 > f1)
 814     {
 815       max = f2;
 816     }
 817
 818     if (f3 > max)
 819     {
 820       max = f3;
 821     }
 822
 823     return max;
 824   }
 825
 826   /**
 827    * DOCUMENT ME!
 828    *
 829    * @param f1
 830    *          DOCUMENT ME!
 831    * @param f2
 832    *          DOCUMENT ME!
 833    *
 834    * @return DOCUMENT ME!
 835    */
 836   private static float max(float f1, float f2)
 837   {
 838     float max = f1;
 839
 840     if (f2 > f1)
 841     {
 842       max = f2;
 843     }
 844
 845     return max;
 846   }
 847
 848   /**
 849    * Converts the character string to an array of integers which are the
 850    * corresponding indices to the characters in the score matrix
 851    *
 852    * @param s
 853    *
 854    * @return
 855    */
 856   int[] indexEncode(String s)
 857   {
 858     int[] encoded = new int[s.length()];
 859
 860     for (int i = 0; i < s.length(); i++)
 861     {
 862       char c = s.charAt(i);
 863       encoded[i] = scoreMatrix.getMatrixIndex(c);
 864     }
 865
 866     return encoded;
 867   }
 868
 869   /**
 870    * DOCUMENT ME!
 871    *
 872    * @param g
 873    *          DOCUMENT ME!
 874    * @param mat
 875    *          DOCUMENT ME!
 876    * @param n
 877    *          DOCUMENT ME!
 878    * @param m
 879    *          DOCUMENT ME!
 880    * @param psize
 881    *          DOCUMENT ME!
 882    */
 883   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
 884           int psize)
 885   {
 886     // TODO method doesn't seem to be referenced anywhere delete??
 887     int max = -1000;
 888     int min = 1000;
 889
 890     for (int i = 0; i < n; i++)
 891     {
 892       for (int j = 0; j < m; j++)
 893       {
 894         if (mat[i][j] >= max)
 895         {
 896           max = mat[i][j];
 897         }
 898
 899         if (mat[i][j] <= min)
 900         {
 901           min = mat[i][j];
 902         }
 903       }
 904     }
 905
 906     System.out.println(max + " " + min);
 907
 908     for (int i = 0; i < n; i++)
 909     {
 910       for (int j = 0; j < m; j++)
 911       {
 912         int x = psize * i;
 913         int y = psize * j;
 914
 915         // System.out.println(mat[i][j]);
 916         float score = (float) (mat[i][j] - min) / (float) (max - min);
 917         g.setColor(new Color(score, 0, 0));
 918         g.fillRect(x, y, psize, psize);
 919
 920         // System.out.println(x + " " + y + " " + score);
 921       }
 922     }
 923   }
 924
 925   /**
 926    * Compute a globally optimal needleman and wunsch alignment between two
 927    * sequences
 928    *
 929    * @param s1
 930    * @param s2
 931    * @param type
 932    *          AlignSeq.DNA or AlignSeq.PEP
 933    */
 934   public static AlignSeq doGlobalNWAlignment(SequenceI s1, SequenceI s2,
 935           String type)
 936   {
 937     AlignSeq as = new AlignSeq(s1, s2, type);
 938
 939     as.calcScoreMatrix();
 940     as.traceAlignment();
 941     return as;
 942   }
 943
 944   /**
 945    *
 946    * @return mapping from positions in S1 to corresponding positions in S2
 947    */
 948   public jalview.datamodel.Mapping getMappingFromS1(boolean allowmismatch)
 949   {
 950     ArrayList<Integer> as1 = new ArrayList<Integer>(),
 951             as2 = new ArrayList<Integer>();
 952     int pdbpos = s2.getStart() + getSeq2Start() - 2;
 953     int alignpos = s1.getStart() + getSeq1Start() - 2;
 954     int lp2 = pdbpos - 3, lp1 = alignpos - 3;
 955     boolean lastmatch = false;
 956     // and now trace the alignment onto the atom set.
 957     for (int i = 0; i < astr1.length(); i++)
 958     {
 959       char c1 = astr1.charAt(i), c2 = astr2.charAt(i);
 960       if (c1 != '-')
 961       {
 962         alignpos++;
 963       }
 964
 965       if (c2 != '-')
 966       {
 967         pdbpos++;
 968       }
 969
 970       // ignore case differences
 971       if (allowmismatch || (c1 == c2) || (Math.abs(c2-c1)==('a'-'A')))
 972       {
 973         // extend mapping interval
 974         if (lp1 + 1 != alignpos || lp2 + 1 != pdbpos)
 975         {
 976           as1.add(Integer.valueOf(alignpos));
 977           as2.add(Integer.valueOf(pdbpos));
 978         }
 979         lastmatch = true;
 980         lp1 = alignpos;
 981         lp2 = pdbpos;
 982       }
 983       else
 984       {
 985         // extend mapping interval
 986         if (lastmatch)
 987         {
 988           as1.add(Integer.valueOf(lp1));
 989           as2.add(Integer.valueOf(lp2));
 990         }
 991         lastmatch = false;
 992       }
 993     }
 994     // construct range pairs
 995
 996     int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)],
 997             mapseq2 = new int[as2.size() + (lastmatch ? 1 : 0)];
 998     int i = 0;
 999     for (Integer ip : as1)
1000     {
1001       mapseq1[i++] = ip;
1002     }
1003     ;
1004     i = 0;
1005     for (Integer ip : as2)
1006     {
1007       mapseq2[i++] = ip;
1008     }
1009     ;
1010     if (lastmatch)
1011     {
1012       mapseq1[mapseq1.length - 1] = alignpos;
1013       mapseq2[mapseq2.length - 1] = pdbpos;
1014     }
1015     MapList map = new MapList(mapseq1, mapseq2, 1, 1);
1016
1017     jalview.datamodel.Mapping mapping = new Mapping(map);
1018     mapping.setTo(s2);
1019     return mapping;
1020   }
1021
1022   /**
1023    * matches ochains against al and populates seqs with the best match between
1024    * each ochain and the set in al
1025    *
1026    * @param ochains
1027    * @param al
1028    * @param dnaOrProtein
1029    * @param removeOldAnnots
1030    *          when true, old annotation is cleared before new annotation
1031    *          transferred
1032    * @return List<List<SequenceI> originals, List<SequenceI> replacement,
1033    *         List<AlignSeq> alignment between each>
1034    */
1035   public static List<List<? extends Object>> replaceMatchingSeqsWith(
1036           List<SequenceI> seqs, List<AlignmentAnnotation> annotations,
1037           List<SequenceI> ochains, AlignmentI al, String dnaOrProtein,
1038           boolean removeOldAnnots)
1039   {
1040     List<SequenceI> orig = new ArrayList<SequenceI>(),
1041             repl = new ArrayList<SequenceI>();
1042     List<AlignSeq> aligs = new ArrayList<AlignSeq>();
1043     if (al != null && al.getHeight() > 0)
1044     {
1045       ArrayList<SequenceI> matches = new ArrayList<SequenceI>();
1046       ArrayList<AlignSeq> aligns = new ArrayList<AlignSeq>();
1047
1048       for (SequenceI sq : ochains)
1049       {
1050         SequenceI bestm = null;
1051         AlignSeq bestaseq = null;
1052         float bestscore = 0;
1053         for (SequenceI msq : al.getSequences())
1054         {
1055           AlignSeq aseq = doGlobalNWAlignment(msq, sq, dnaOrProtein);
1056           if (bestm == null || aseq.getMaxScore() > bestscore)
1057           {
1058             bestscore = aseq.getMaxScore();
1059             bestaseq = aseq;
1060             bestm = msq;
1061           }
1062         }
1063         // System.out.println("Best Score for " + (matches.size() + 1) + " :"
1064         // + bestscore);
1065         matches.add(bestm);
1066         aligns.add(bestaseq);
1067         al.deleteSequence(bestm);
1068       }
1069       for (int p = 0, pSize = seqs.size(); p < pSize; p++)
1070       {
1071         SequenceI sq, sp = seqs.get(p);
1072         int q;
1073         if ((q = ochains.indexOf(sp)) > -1)
1074         {
1075           seqs.set(p, sq = matches.get(q));
1076           orig.add(sp);
1077           repl.add(sq);
1078           sq.setName(sp.getName());
1079           sq.setDescription(sp.getDescription());
1080           Mapping sp2sq;
1081           sq.transferAnnotation(sp,
1082                   sp2sq = aligns.get(q).getMappingFromS1(false));
1083           aligs.add(aligns.get(q));
1084           int inspos = -1;
1085           for (int ap = 0; ap < annotations.size();)
1086           {
1087             if (annotations.get(ap).sequenceRef == sp)
1088             {
1089               if (inspos == -1)
1090               {
1091                 inspos = ap;
1092               }
1093               if (removeOldAnnots)
1094               {
1095                 annotations.remove(ap);
1096               }
1097               else
1098               {
1099                 AlignmentAnnotation alan = annotations.remove(ap);
1100                 alan.liftOver(sq, sp2sq);
1101                 alan.setSequenceRef(sq);
1102                 sq.addAlignmentAnnotation(alan);
1103               }
1104             }
1105             else
1106             {
1107               ap++;
1108             }
1109           }
1110           if (sq.getAnnotation() != null && sq.getAnnotation().length > 0)
1111           {
1112             annotations.addAll(inspos == -1 ? annotations.size() : inspos,
1113                     Arrays.asList(sq.getAnnotation()));
1114           }
1115         }
1116       }
1117     }
1118     return Arrays.asList(orig, repl, aligs);
1119   }
1120
1121   /**
1122    * compute the PID vector used by the redundancy filter.
1123    *
1124    * @param originalSequences
1125    *          - sequences in alignment that are to filtered
1126    * @param omitHidden
1127    *          - null or strings to be analysed (typically, visible portion of
1128    *          each sequence in alignment)
1129    * @param start
1130    *          - first column in window for calculation
1131    * @param end
1132    *          - last column in window for calculation
1133    * @param ungapped
1134    *          - if true then use ungapped sequence to compute PID
1135    * @return vector containing maximum PID for i-th sequence and any sequences
1136    *         longer than that seuqence
1137    */
1138   public static float[] computeRedundancyMatrix(
1139           SequenceI[] originalSequences, String[] omitHidden, int start,
1140           int end, boolean ungapped)
1141   {
1142     int height = originalSequences.length;
1143     float[] redundancy = new float[height];
1144     int[] lngth = new int[height];
1145     for (int i = 0; i < height; i++)
1146     {
1147       redundancy[i] = 0f;
1148       lngth[i] = -1;
1149     }
1150
1151     // long start = System.currentTimeMillis();
1152
1153     SimilarityParams pidParams = new SimilarityParams(true, true, true,
1154             true);
1155     float pid;
1156     String seqi, seqj;
1157     for (int i = 0; i < height; i++)
1158     {
1159
1160       for (int j = 0; j < i; j++)
1161       {
1162         if (i == j)
1163         {
1164           continue;
1165         }
1166
1167         if (omitHidden == null)
1168         {
1169           seqi = originalSequences[i].getSequenceAsString(start, end);
1170           seqj = originalSequences[j].getSequenceAsString(start, end);
1171         }
1172         else
1173         {
1174           seqi = omitHidden[i];
1175           seqj = omitHidden[j];
1176         }
1177         if (lngth[i] == -1)
1178         {
1179           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
1180           lngth[i] = ug.length();
1181           if (ungapped)
1182           {
1183             seqi = ug;
1184           }
1185         }
1186         if (lngth[j] == -1)
1187         {
1188           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
1189           lngth[j] = ug.length();
1190           if (ungapped)
1191           {
1192             seqj = ug;
1193           }
1194         }
1195         pid = (float) PIDModel.computePID(seqi, seqj, pidParams);
1196
1197         // use real sequence length rather than string length
1198         if (lngth[j] < lngth[i])
1199         {
1200           redundancy[j] = Math.max(pid, redundancy[j]);
1201         }
1202         else
1203         {
1204           redundancy[i] = Math.max(pid, redundancy[i]);
1205         }
1206
1207       }
1208     }
1209     return redundancy;
1210   }
1211
1212   /**
1213   * calculate the mean score of the alignment
1214   * 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
1215   *
1216   */
1217   public void meanScore()
1218   {
1219     //int length = (indelfreeAstr1.length() > indelfreeAstr2.length()) ? indelfreeAstr1.length() : indelfreeAstr2.length();
1220     int length = indelfreeAstr1.length();       //both have the same length
1221     //create HashMap for counting residues in each sequence
1222     HashMap<Character, Integer> seq1ResCount = new HashMap<Character, Integer>();
1223     HashMap<Character, Integer> seq2ResCount = new HashMap<Character, Integer>();
1224
1225     // for both sequences (String indelfreeAstr1 or 2) create a key for the residue and add 1 each time its encountered
1226     for (char residue: indelfreeAstr1.toCharArray())
1227     {
1228       seq1ResCount.putIfAbsent(residue, 0);
1229       seq1ResCount.replace(residue, seq1ResCount.get(residue) + 1);
1230     }
1231     for (char residue: indelfreeAstr2.toCharArray())
1232     {
1233       seq2ResCount.putIfAbsent(residue, 0);
1234       seq2ResCount.replace(residue, seq2ResCount.get(residue) + 1);
1235     }
1236
1237     // meanscore = for each residue pair get the number of appearance and add (countA * countB * pairwiseScore(AB))
1238     // divide the meanscore by the sequence length afterwards
1239     float _meanscore = 0;
1240     for (char resA : seq1ResCount.keySet())
1241     {
1242       for (char resB : seq2ResCount.keySet())
1243       {
1244         int countA = seq1ResCount.get(resA);
1245         int countB = seq2ResCount.get(resB);
1246
1247         float scoreAB = scoreMatrix.getPairwiseScore(resA, resB);
1248
1249         _meanscore += countA *  countB * scoreAB;
1250       }
1251     }
1252     _meanscore /= length;
1253     this.meanScore = _meanscore;
1254   }
1255
1256   public float getMeanScore()
1257   {
1258     return this.meanScore;
1259   }
1260
1261   /**
1262   * calculate the hypothetic max score using the self-alignment of the sequences
1263   */
1264   public void hypotheticMaxScore()
1265   {
1266     int _hmsA = 0;
1267     int _hmsB = 0;
1268     for (char residue: indelfreeAstr1.toCharArray())
1269     {
1270       _hmsA += scoreMatrix.getPairwiseScore(residue, residue);
1271     }
1272     for (char residue: indelfreeAstr2.toCharArray())
1273     {
1274       _hmsB += scoreMatrix.getPairwiseScore(residue, residue);
1275     }
1276     this.hypotheticMaxScore = (_hmsA < _hmsB) ? _hmsA : _hmsB;  // take the lower self alignment
1277
1278   }
1279
1280   public int getHypotheticMaxScore()
1281   {
1282     return this.hypotheticMaxScore;
1283   }
1284
1285   /**
1286   * create strings based of astr1 and astr2 but without gaps
1287   */
1288   public void getIndelfreeAstr()
1289   {
1290     int n = astr1.length();     // both have the same length
1291     for (int i = 0; i < n; i++)
1292     {
1293       if (Character.isLetter(astr1.charAt(i)) && Character.isLetter(astr2.charAt(i)))   // if both sequences dont have a gap -> add to indelfreeAstr
1294       {
1295         this.indelfreeAstr1 += astr1.charAt(i);
1296         this.indelfreeAstr2 += astr2.charAt(i);
1297       }
1298     }
1299   }
1300
1301   /**
1302   * calculates the overall score of the alignment
1303   * preprescore = sum of all scores - all penalties
1304   * if preprescore < 1 ~ alignmentScore = Float.NaN     >
1305   * alignmentScore = ((preprescore - meanScore) / (hypotheticMaxScore - meanScore)) * coverage
1306   */
1307   public void scoreAlignment() throws RuntimeException
1308   {
1309
1310     getIndelfreeAstr();
1311     meanScore();
1312     hypotheticMaxScore();
1313     // cannot calculate score because denominator would be zero
1314     if (this.hypotheticMaxScore == this.meanScore)
1315     {
1316       throw new IllegalArgumentException(String.format("hypotheticMaxScore (%8.2f) == meanScore (%8.2f) - division by 0", hypotheticMaxScore, meanScore));
1317     }
1318     //int n = (astr1.length() > astr2.length()) ? astr1.length() : astr2.length();
1319     int n = indelfreeAstr1.length();
1320
1321     float score = 0;
1322     boolean aGapOpen = false;
1323     boolean bGapOpen = false;
1324     for (int i = 0; i < n; i++)
1325     {
1326       char char1 = indelfreeAstr1.charAt(i);
1327       char char2 = indelfreeAstr2.charAt(i);
1328       boolean aIsLetter = Character.isLetter(char1);
1329       boolean bIsLetter = Character.isLetter(char2);
1330       if (aIsLetter && bIsLetter)       // if pair -> get score
1331       {
1332         score += scoreMatrix.getPairwiseScore(char1, char2);
1333       } else if (!aIsLetter && !bIsLetter) {    // both are gap -> skip
1334       } else if ((!aIsLetter && aGapOpen) || (!bIsLetter && bGapOpen)) {        // one side gapopen -> score - gap_extend
1335         score -= GAP_EXTEND_COST;
1336       } else {          // no gap open -> score - gap_open
1337         score -= GAP_OPEN_COST;
1338       }
1339       // adjust GapOpen status in both sequences
1340       aGapOpen = (!aIsLetter) ? true : false;
1341       bGapOpen = (!bIsLetter) ? true : false;
1342     }
1343
1344     float preprescore = score;  // if this score < 1 --> alignment score = Float.NaN
1345     score = (score - this.meanScore) / (this.hypotheticMaxScore - this.meanScore);
1346     int[] _max = MiscMath.findMax(new int[]{astr1.replace("-","").length(), astr2.replace("-","").length()});   // {index of max, max}
1347     float coverage = (float) n / (float) _max[1];       // indelfreeAstr length / longest sequence length
1348     float prescore = score;     // only debug
1349     score *= coverage;
1350
1351     System.out.println(String.format("prepre-score: %f, pre-score: %f, longlength: %d\nscore: %f, mean: %f, max: %d", preprescore, prescore, _max[1], score, this.meanScore, this.hypotheticMaxScore));
1352     this.alignmentScore = (preprescore < 1) ? Float.NaN : score;
1353   }
1354 }