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       if (allowmismatch || c1 == c2)
 971       {
 972         // extend mapping interval
 973         if (lp1 + 1 != alignpos || lp2 + 1 != pdbpos)
 974         {
 975           as1.add(Integer.valueOf(alignpos));
 976           as2.add(Integer.valueOf(pdbpos));
 977         }
 978         lastmatch = true;
 979         lp1 = alignpos;
 980         lp2 = pdbpos;
 981       }
 982       else
 983       {
 984         // extend mapping interval
 985         if (lastmatch)
 986         {
 987           as1.add(Integer.valueOf(lp1));
 988           as2.add(Integer.valueOf(lp2));
 989         }
 990         lastmatch = false;
 991       }
 992     }
 993     // construct range pairs
 994
 995     int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)],
 996             mapseq2 = new int[as2.size() + (lastmatch ? 1 : 0)];
 997     int i = 0;
 998     for (Integer ip : as1)
 999     {
1000       mapseq1[i++] = ip;
1001     }
1002     ;
1003     i = 0;
1004     for (Integer ip : as2)
1005     {
1006       mapseq2[i++] = ip;
1007     }
1008     ;
1009     if (lastmatch)
1010     {
1011       mapseq1[mapseq1.length - 1] = alignpos;
1012       mapseq2[mapseq2.length - 1] = pdbpos;
1013     }
1014     MapList map = new MapList(mapseq1, mapseq2, 1, 1);
1015
1016     jalview.datamodel.Mapping mapping = new Mapping(map);
1017     mapping.setTo(s2);
1018     return mapping;
1019   }
1020
1021   /**
1022    * matches ochains against al and populates seqs with the best match between
1023    * each ochain and the set in al
1024    *
1025    * @param ochains
1026    * @param al
1027    * @param dnaOrProtein
1028    * @param removeOldAnnots
1029    *          when true, old annotation is cleared before new annotation
1030    *          transferred
1031    * @return List<List<SequenceI> originals, List<SequenceI> replacement,
1032    *         List<AlignSeq> alignment between each>
1033    */
1034   public static List<List<? extends Object>> replaceMatchingSeqsWith(
1035           List<SequenceI> seqs, List<AlignmentAnnotation> annotations,
1036           List<SequenceI> ochains, AlignmentI al, String dnaOrProtein,
1037           boolean removeOldAnnots)
1038   {
1039     List<SequenceI> orig = new ArrayList<SequenceI>(),
1040             repl = new ArrayList<SequenceI>();
1041     List<AlignSeq> aligs = new ArrayList<AlignSeq>();
1042     if (al != null && al.getHeight() > 0)
1043     {
1044       ArrayList<SequenceI> matches = new ArrayList<SequenceI>();
1045       ArrayList<AlignSeq> aligns = new ArrayList<AlignSeq>();
1046
1047       for (SequenceI sq : ochains)
1048       {
1049         SequenceI bestm = null;
1050         AlignSeq bestaseq = null;
1051         float bestscore = 0;
1052         for (SequenceI msq : al.getSequences())
1053         {
1054           AlignSeq aseq = doGlobalNWAlignment(msq, sq, dnaOrProtein);
1055           if (bestm == null || aseq.getMaxScore() > bestscore)
1056           {
1057             bestscore = aseq.getMaxScore();
1058             bestaseq = aseq;
1059             bestm = msq;
1060           }
1061         }
1062         // System.out.println("Best Score for " + (matches.size() + 1) + " :"
1063         // + bestscore);
1064         matches.add(bestm);
1065         aligns.add(bestaseq);
1066         al.deleteSequence(bestm);
1067       }
1068       for (int p = 0, pSize = seqs.size(); p < pSize; p++)
1069       {
1070         SequenceI sq, sp = seqs.get(p);
1071         int q;
1072         if ((q = ochains.indexOf(sp)) > -1)
1073         {
1074           seqs.set(p, sq = matches.get(q));
1075           orig.add(sp);
1076           repl.add(sq);
1077           sq.setName(sp.getName());
1078           sq.setDescription(sp.getDescription());
1079           Mapping sp2sq;
1080           sq.transferAnnotation(sp,
1081                   sp2sq = aligns.get(q).getMappingFromS1(false));
1082           aligs.add(aligns.get(q));
1083           int inspos = -1;
1084           for (int ap = 0; ap < annotations.size();)
1085           {
1086             if (annotations.get(ap).sequenceRef == sp)
1087             {
1088               if (inspos == -1)
1089               {
1090                 inspos = ap;
1091               }
1092               if (removeOldAnnots)
1093               {
1094                 annotations.remove(ap);
1095               }
1096               else
1097               {
1098                 AlignmentAnnotation alan = annotations.remove(ap);
1099                 alan.liftOver(sq, sp2sq);
1100                 alan.setSequenceRef(sq);
1101                 sq.addAlignmentAnnotation(alan);
1102               }
1103             }
1104             else
1105             {
1106               ap++;
1107             }
1108           }
1109           if (sq.getAnnotation() != null && sq.getAnnotation().length > 0)
1110           {
1111             annotations.addAll(inspos == -1 ? annotations.size() : inspos,
1112                     Arrays.asList(sq.getAnnotation()));
1113           }
1114         }
1115       }
1116     }
1117     return Arrays.asList(orig, repl, aligs);
1118   }
1119
1120   /**
1121    * compute the PID vector used by the redundancy filter.
1122    *
1123    * @param originalSequences
1124    *          - sequences in alignment that are to filtered
1125    * @param omitHidden
1126    *          - null or strings to be analysed (typically, visible portion of
1127    *          each sequence in alignment)
1128    * @param start
1129    *          - first column in window for calculation
1130    * @param end
1131    *          - last column in window for calculation
1132    * @param ungapped
1133    *          - if true then use ungapped sequence to compute PID
1134    * @return vector containing maximum PID for i-th sequence and any sequences
1135    *         longer than that seuqence
1136    */
1137   public static float[] computeRedundancyMatrix(
1138           SequenceI[] originalSequences, String[] omitHidden, int start,
1139           int end, boolean ungapped)
1140   {
1141     int height = originalSequences.length;
1142     float[] redundancy = new float[height];
1143     int[] lngth = new int[height];
1144     for (int i = 0; i < height; i++)
1145     {
1146       redundancy[i] = 0f;
1147       lngth[i] = -1;
1148     }
1149
1150     // long start = System.currentTimeMillis();
1151
1152     SimilarityParams pidParams = new SimilarityParams(true, true, true,
1153             true);
1154     float pid;
1155     String seqi, seqj;
1156     for (int i = 0; i < height; i++)
1157     {
1158
1159       for (int j = 0; j < i; j++)
1160       {
1161         if (i == j)
1162         {
1163           continue;
1164         }
1165
1166         if (omitHidden == null)
1167         {
1168           seqi = originalSequences[i].getSequenceAsString(start, end);
1169           seqj = originalSequences[j].getSequenceAsString(start, end);
1170         }
1171         else
1172         {
1173           seqi = omitHidden[i];
1174           seqj = omitHidden[j];
1175         }
1176         if (lngth[i] == -1)
1177         {
1178           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
1179           lngth[i] = ug.length();
1180           if (ungapped)
1181           {
1182             seqi = ug;
1183           }
1184         }
1185         if (lngth[j] == -1)
1186         {
1187           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
1188           lngth[j] = ug.length();
1189           if (ungapped)
1190           {
1191             seqj = ug;
1192           }
1193         }
1194         pid = (float) PIDModel.computePID(seqi, seqj, pidParams);
1195
1196         // use real sequence length rather than string length
1197         if (lngth[j] < lngth[i])
1198         {
1199           redundancy[j] = Math.max(pid, redundancy[j]);
1200         }
1201         else
1202         {
1203           redundancy[i] = Math.max(pid, redundancy[i]);
1204         }
1205
1206       }
1207     }
1208     return redundancy;
1209   }
1210
1211   /**
1212   * calculate the mean score of the alignment
1213   * 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
1214   *
1215   */
1216   public void meanScore()
1217   {
1218     //int length = (indelfreeAstr1.length() > indelfreeAstr2.length()) ? indelfreeAstr1.length() : indelfreeAstr2.length();
1219     int length = indelfreeAstr1.length();       //both have the same length
1220     //create HashMap for counting residues in each sequence
1221     HashMap<Character, Integer> seq1ResCount = new HashMap<Character, Integer>();
1222     HashMap<Character, Integer> seq2ResCount = new HashMap<Character, Integer>();
1223
1224     // for both sequences (String indelfreeAstr1 or 2) create a key for the residue and add 1 each time its encountered
1225     for (char residue: indelfreeAstr1.toCharArray())
1226     {
1227       seq1ResCount.putIfAbsent(residue, 0);
1228       seq1ResCount.replace(residue, seq1ResCount.get(residue) + 1);
1229     }
1230     for (char residue: indelfreeAstr2.toCharArray())
1231     {
1232       seq2ResCount.putIfAbsent(residue, 0);
1233       seq2ResCount.replace(residue, seq2ResCount.get(residue) + 1);
1234     }
1235
1236     // meanscore = for each residue pair get the number of appearance and add (countA * countB * pairwiseScore(AB))
1237     // divide the meanscore by the sequence length afterwards
1238     float _meanscore = 0;
1239     for (char resA : seq1ResCount.keySet())
1240     {
1241       for (char resB : seq2ResCount.keySet())
1242       {
1243         int countA = seq1ResCount.get(resA);
1244         int countB = seq2ResCount.get(resB);
1245
1246         float scoreAB = scoreMatrix.getPairwiseScore(resA, resB);
1247
1248         _meanscore += countA *  countB * scoreAB;
1249       }
1250     }
1251     _meanscore /= length;
1252     this.meanScore = _meanscore;
1253   }
1254
1255   public float getMeanScore()
1256   {
1257     return this.meanScore;
1258   }
1259
1260   /**
1261   * calculate the hypothetic max score using the self-alignment of the sequences
1262   */
1263   public void hypotheticMaxScore()
1264   {
1265     int _hmsA = 0;
1266     int _hmsB = 0;
1267     for (char residue: indelfreeAstr1.toCharArray())
1268     {
1269       _hmsA += scoreMatrix.getPairwiseScore(residue, residue);
1270     }
1271     for (char residue: indelfreeAstr2.toCharArray())
1272     {
1273       _hmsB += scoreMatrix.getPairwiseScore(residue, residue);
1274     }
1275     this.hypotheticMaxScore = (_hmsA < _hmsB) ? _hmsA : _hmsB;  // take the lower self alignment
1276
1277   }
1278
1279   public int getHypotheticMaxScore()
1280   {
1281     return this.hypotheticMaxScore;
1282   }
1283
1284   /**
1285   * create strings based of astr1 and astr2 but without gaps
1286   */
1287   public void getIndelfreeAstr()
1288   {
1289     int n = astr1.length();     // both have the same length
1290     for (int i = 0; i < n; i++)
1291     {
1292       if (Character.isLetter(astr1.charAt(i)) && Character.isLetter(astr2.charAt(i)))   // if both sequences dont have a gap -> add to indelfreeAstr
1293       {
1294         this.indelfreeAstr1 += astr1.charAt(i);
1295         this.indelfreeAstr2 += astr2.charAt(i);
1296       }
1297     }
1298   }
1299
1300   /**
1301   * calculates the overall score of the alignment
1302   * preprescore = sum of all scores - all penalties
1303   * if preprescore < 1 ~ alignmentScore = Float.NaN     >
1304   * alignmentScore = ((preprescore - meanScore) / (hypotheticMaxScore - meanScore)) * coverage
1305   */
1306   public void scoreAlignment() throws RuntimeException
1307   {
1308
1309     getIndelfreeAstr();
1310     meanScore();
1311     hypotheticMaxScore();
1312     // cannot calculate score because denominator would be zero
1313     if (this.hypotheticMaxScore == this.meanScore)
1314     {
1315       throw new IllegalArgumentException(String.format("hypotheticMaxScore (%8.2f) == meanScore (%8.2f) - division by 0", hypotheticMaxScore, meanScore));
1316     }
1317     //int n = (astr1.length() > astr2.length()) ? astr1.length() : astr2.length();
1318     int n = indelfreeAstr1.length();
1319
1320     float score = 0;
1321     boolean aGapOpen = false;
1322     boolean bGapOpen = false;
1323     for (int i = 0; i < n; i++)
1324     {
1325       char char1 = indelfreeAstr1.charAt(i);
1326       char char2 = indelfreeAstr2.charAt(i);
1327       boolean aIsLetter = Character.isLetter(char1);
1328       boolean bIsLetter = Character.isLetter(char2);
1329       if (aIsLetter && bIsLetter)       // if pair -> get score
1330       {
1331         score += scoreMatrix.getPairwiseScore(char1, char2);
1332       } else if (!aIsLetter && !bIsLetter) {    // both are gap -> skip
1333       } else if ((!aIsLetter && aGapOpen) || (!bIsLetter && bGapOpen)) {        // one side gapopen -> score - gap_extend
1334         score -= GAP_EXTEND_COST;
1335       } else {          // no gap open -> score - gap_open
1336         score -= GAP_OPEN_COST;
1337       }
1338       // adjust GapOpen status in both sequences
1339       aGapOpen = (!aIsLetter) ? true : false;
1340       bGapOpen = (!bIsLetter) ? true : false;
1341     }
1342
1343     float preprescore = score;  // if this score < 1 --> alignment score = Float.NaN
1344     score = (score - this.meanScore) / (this.hypotheticMaxScore - this.meanScore);
1345     int[] _max = MiscMath.findMax(new int[]{astr1.replace("-","").length(), astr2.replace("-","").length()});   // {index of max, max}
1346     float coverage = (float) n / (float) _max[1];       // indelfreeAstr length / longest sequence length
1347     float prescore = score;     // only debug
1348     score *= coverage;
1349
1350     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));
1351     this.alignmentScore = (preprescore < 1) ? Float.NaN : score;
1352   }
1353 }