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