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.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.util.ArrayList;
  43 import java.util.Arrays;
  44 import java.util.List;
  45 import java.util.StringTokenizer;
  46
  47 /**
  48  *
  49  *
  50  * @author $author$
  51  * @version $Revision$
  52  */
  53 public class AlignSeq
  54 {
  55   private static final int MAX_NAME_LENGTH = 30;
  56
  57   private static final int GAP_OPEN_COST = 120;
  58
  59   private static final int GAP_EXTEND_COST = 20;
  60
  61   private static final int GAP_INDEX = -1;
  62
  63   public static final String PEP = "pep";
  64
  65   public static final String DNA = "dna";
  66
  67   private static final String NEWLINE = System.lineSeparator();
  68
  69   float[][] score;
  70
  71   float[][] E;
  72
  73   float[][] F;
  74
  75   int[][] traceback; // todo is this actually used?
  76
  77   int[] seq1;
  78
  79   int[] seq2;
  80
  81   SequenceI s1;
  82
  83   SequenceI s2;
  84
  85   public String s1str;
  86
  87   public String s2str;
  88
  89   int maxi;
  90
  91   int maxj;
  92
  93   int[] aseq1;
  94
  95   int[] aseq2;
  96
  97   public String astr1 = "";
  98
  99   public String astr2 = "";
 100
 101   /** DOCUMENT ME!! */
 102   public int seq1start;
 103
 104   /** DOCUMENT ME!! */
 105   public int seq1end;
 106
 107   /** DOCUMENT ME!! */
 108   public int seq2start;
 109
 110   public int seq2end;
 111
 112   int count;
 113
 114   public float maxscore;
 115
 116   int prev = 0;
 117
 118   StringBuffer output = new StringBuffer();
 119
 120   String type; // AlignSeq.PEP or AlignSeq.DNA
 121
 122   private ScoreMatrix scoreMatrix;
 123
 124   /**
 125    * Creates a new AlignSeq object.
 126    *
 127    * @param s1
 128    *          first sequence for alignment
 129    * @param s2
 130    *          second sequence for alignment
 131    * @param type
 132    *          molecule type, either AlignSeq.PEP or AlignSeq.DNA
 133    */
 134   public AlignSeq(SequenceI s1, SequenceI s2, String type)
 135   {
 136     seqInit(s1, s1.getSequenceAsString(), s2, s2.getSequenceAsString(),
 137             type);
 138   }
 139
 140   /**
 141    * Creates a new AlignSeq object.
 142    *
 143    * @param s1
 144    *          DOCUMENT ME!
 145    * @param s2
 146    *          DOCUMENT ME!
 147    * @param type
 148    *          DOCUMENT ME!
 149    */
 150   public AlignSeq(SequenceI s1, String string1, SequenceI s2,
 151           String string2, String type)
 152   {
 153     seqInit(s1, string1.toUpperCase(Locale.ROOT), s2,
 154             string2.toUpperCase(Locale.ROOT), type);
 155   }
 156
 157   /**
 158    * DOCUMENT ME!
 159    *
 160    * @return DOCUMENT ME!
 161    */
 162   public float getMaxScore()
 163   {
 164     return maxscore;
 165   }
 166
 167   /**
 168    * DOCUMENT ME!
 169    *
 170    * @return DOCUMENT ME!
 171    */
 172   public int getSeq2Start()
 173   {
 174     return seq2start;
 175   }
 176
 177   /**
 178    * DOCUMENT ME!
 179    *
 180    * @return DOCUMENT ME!
 181    */
 182   public int getSeq2End()
 183   {
 184     return seq2end;
 185   }
 186
 187   /**
 188    * DOCUMENT ME!
 189    *
 190    * @return DOCUMENT ME!
 191    */
 192   public int getSeq1Start()
 193   {
 194     return seq1start;
 195   }
 196
 197   /**
 198    * DOCUMENT ME!
 199    *
 200    * @return DOCUMENT ME!
 201    */
 202   public int getSeq1End()
 203   {
 204     return seq1end;
 205   }
 206
 207   /**
 208    * DOCUMENT ME!
 209    *
 210    * @return DOCUMENT ME!
 211    */
 212   public String getOutput()
 213   {
 214     return output.toString();
 215   }
 216
 217   /**
 218    * DOCUMENT ME!
 219    *
 220    * @return DOCUMENT ME!
 221    */
 222   public String getAStr1()
 223   {
 224     return astr1;
 225   }
 226
 227   /**
 228    * DOCUMENT ME!
 229    *
 230    * @return DOCUMENT ME!
 231    */
 232   public String getAStr2()
 233   {
 234     return astr2;
 235   }
 236
 237   /**
 238    * DOCUMENT ME!
 239    *
 240    * @return DOCUMENT ME!
 241    */
 242   public int[] getASeq1()
 243   {
 244     return aseq1;
 245   }
 246
 247   /**
 248    * DOCUMENT ME!
 249    *
 250    * @return DOCUMENT ME!
 251    */
 252   public int[] getASeq2()
 253   {
 254     return aseq2;
 255   }
 256
 257   /**
 258    *
 259    * @return aligned instance of Seq 1
 260    */
 261   public SequenceI getAlignedSeq1()
 262   {
 263     SequenceI alSeq1 = new Sequence(s1.getName(), getAStr1());
 264     alSeq1.setStart(s1.getStart() + getSeq1Start() - 1);
 265     alSeq1.setEnd(s1.getStart() + getSeq1End() - 1);
 266     alSeq1.setDatasetSequence(
 267             s1.getDatasetSequence() == null ? s1 : s1.getDatasetSequence());
 268     return alSeq1;
 269   }
 270
 271   /**
 272    *
 273    * @return aligned instance of Seq 2
 274    */
 275   public SequenceI getAlignedSeq2()
 276   {
 277     SequenceI alSeq2 = new Sequence(s2.getName(), getAStr2());
 278     alSeq2.setStart(s2.getStart() + getSeq2Start() - 1);
 279     alSeq2.setEnd(s2.getStart() + getSeq2End() - 1);
 280     alSeq2.setDatasetSequence(
 281             s2.getDatasetSequence() == null ? s2 : s2.getDatasetSequence());
 282     return alSeq2;
 283   }
 284
 285   /**
 286    * Construct score matrix for sequences with standard DNA or PEPTIDE matrix
 287    *
 288    * @param s1
 289    *          - sequence 1
 290    * @param string1
 291    *          - string to use for s1
 292    * @param s2
 293    *          - sequence 2
 294    * @param string2
 295    *          - string to use for s2
 296    * @param type
 297    *          DNA or PEPTIDE
 298    */
 299   public void seqInit(SequenceI s1, String string1, SequenceI s2,
 300           String string2, String type)
 301   {
 302     this.s1 = s1;
 303     this.s2 = s2;
 304     setDefaultParams(type);
 305     seqInit(string1, string2);
 306   }
 307
 308   /**
 309    * construct score matrix for string1 and string2 (after removing any existing
 310    * gaps
 311    *
 312    * @param string1
 313    * @param string2
 314    */
 315   private void seqInit(String string1, String string2)
 316   {
 317     s1str = extractGaps(jalview.util.Comparison.GapChars, string1);
 318     s2str = extractGaps(jalview.util.Comparison.GapChars, string2);
 319
 320     if (s1str.length() == 0 || s2str.length() == 0)
 321     {
 322       output.append(
 323               "ALL GAPS: " + (s1str.length() == 0 ? s1.getName() : " ")
 324                       + (s2str.length() == 0 ? s2.getName() : ""));
 325       return;
 326     }
 327
 328     score = new float[s1str.length()][s2str.length()];
 329
 330     E = new float[s1str.length()][s2str.length()];
 331
 332     F = new float[s1str.length()][s2str.length()];
 333     traceback = new int[s1str.length()][s2str.length()];
 334
 335     seq1 = indexEncode(s1str);
 336
 337     seq2 = indexEncode(s2str);
 338   }
 339
 340   private void setDefaultParams(String moleculeType)
 341   {
 342     if (!PEP.equals(moleculeType) && !DNA.equals(moleculeType))
 343     {
 344       output.append("Wrong type = dna or pep only");
 345       throw new Error(MessageManager
 346               .formatMessage("error.unknown_type_dna_or_pep", new String[]
 347               { moleculeType }));
 348     }
 349
 350     type = moleculeType;
 351     scoreMatrix = ScoreModels.getInstance()
 352             .getDefaultModel(PEP.equals(type));
 353   }
 354
 355   /**
 356    * DOCUMENT ME!
 357    */
 358   public void traceAlignment()
 359   {
 360     // Find the maximum score along the rhs or bottom row
 361     float max = -Float.MAX_VALUE;
 362
 363     for (int i = 0; i < seq1.length; i++)
 364     {
 365       if (score[i][seq2.length - 1] > max)
 366       {
 367         max = score[i][seq2.length - 1];
 368         maxi = i;
 369         maxj = seq2.length - 1;
 370       }
 371     }
 372
 373     for (int j = 0; j < seq2.length; j++)
 374     {
 375       if (score[seq1.length - 1][j] > max)
 376       {
 377         max = score[seq1.length - 1][j];
 378         maxi = seq1.length - 1;
 379         maxj = j;
 380       }
 381     }
 382
 383     int i = maxi;
 384     int j = maxj;
 385     int trace;
 386     maxscore = score[i][j] / 10f;
 387
 388     seq1end = maxi + 1;
 389     seq2end = maxj + 1;
 390
 391     aseq1 = new int[seq1.length + seq2.length];
 392     aseq2 = new int[seq1.length + seq2.length];
 393
 394     StringBuilder sb1 = new StringBuilder(aseq1.length);
 395     StringBuilder sb2 = new StringBuilder(aseq2.length);
 396
 397     count = (seq1.length + seq2.length) - 1;
 398
 399     while (i > 0 && j > 0)
 400     {
 401       aseq1[count] = seq1[i];
 402       sb1.append(s1str.charAt(i));
 403       aseq2[count] = seq2[j];
 404       sb2.append(s2str.charAt(j));
 405
 406       trace = findTrace(i, j);
 407
 408       if (trace == 0)
 409       {
 410         i--;
 411         j--;
 412       }
 413       else if (trace == 1)
 414       {
 415         j--;
 416         aseq1[count] = GAP_INDEX;
 417         sb1.replace(sb1.length() - 1, sb1.length(), "-");
 418       }
 419       else if (trace == -1)
 420       {
 421         i--;
 422         aseq2[count] = GAP_INDEX;
 423         sb2.replace(sb2.length() - 1, sb2.length(), "-");
 424       }
 425
 426       count--;
 427     }
 428
 429     seq1start = i + 1;
 430     seq2start = j + 1;
 431
 432     if (aseq1[count] != GAP_INDEX)
 433     {
 434       aseq1[count] = seq1[i];
 435       sb1.append(s1str.charAt(i));
 436     }
 437
 438     if (aseq2[count] != GAP_INDEX)
 439     {
 440       aseq2[count] = seq2[j];
 441       sb2.append(s2str.charAt(j));
 442     }
 443
 444     /*
 445      * we built the character strings backwards, so now
 446      * reverse them to convert to sequence strings
 447      */
 448     astr1 = sb1.reverse().toString();
 449     astr2 = sb2.reverse().toString();
 450   }
 451
 452   /**
 453    * DOCUMENT ME!
 454    */
 455   public void printAlignment(PrintStream os)
 456   {
 457     // TODO: Use original sequence characters rather than re-translated
 458     // characters in output
 459     // Find the biggest id length for formatting purposes
 460     String s1id = getAlignedSeq1().getDisplayId(true);
 461     String s2id = getAlignedSeq2().getDisplayId(true);
 462     int nameLength = Math.max(s1id.length(), s2id.length());
 463     if (nameLength > MAX_NAME_LENGTH)
 464     {
 465       int truncateBy = nameLength - MAX_NAME_LENGTH;
 466       nameLength = MAX_NAME_LENGTH;
 467       // JAL-527 - truncate the sequence ids
 468       if (s1id.length() > nameLength)
 469       {
 470         int slashPos = s1id.lastIndexOf('/');
 471         s1id = s1id.substring(0, slashPos - truncateBy)
 472                 + s1id.substring(slashPos);
 473       }
 474       if (s2id.length() > nameLength)
 475       {
 476         int slashPos = s2id.lastIndexOf('/');
 477         s2id = s2id.substring(0, slashPos - truncateBy)
 478                 + s2id.substring(slashPos);
 479       }
 480     }
 481     int len = 72 - nameLength - 1;
 482     int nochunks = ((aseq1.length - count) / len)
 483             + ((aseq1.length - count) % len > 0 ? 1 : 0);
 484     float pid = 0f;
 485
 486     output.append("Score = ").append(score[maxi][maxj]).append(NEWLINE);
 487     output.append("Length of alignment = ")
 488             .append(String.valueOf(aseq1.length - count)).append(NEWLINE);
 489     output.append("Sequence ");
 490     Format nameFormat = new Format("%" + nameLength + "s");
 491     output.append(nameFormat.form(s1id));
 492     output.append(" (Sequence length = ")
 493             .append(String.valueOf(s1str.length())).append(")")
 494             .append(NEWLINE);
 495     output.append("Sequence ");
 496     output.append(nameFormat.form(s2id));
 497     output.append(" (Sequence length = ")
 498             .append(String.valueOf(s2str.length())).append(")")
 499             .append(NEWLINE).append(NEWLINE);
 500
 501     ScoreMatrix pam250 = ScoreModels.getInstance().getPam250();
 502
 503     for (int j = 0; j < nochunks; j++)
 504     {
 505       // Print the first aligned sequence
 506       output.append(nameFormat.form(s1id)).append(" ");
 507
 508       for (int i = 0; i < len; i++)
 509       {
 510         if ((i + (j * len)) < astr1.length())
 511         {
 512           output.append(astr1.charAt(i + (j * len)));
 513         }
 514       }
 515
 516       output.append(NEWLINE);
 517       output.append(nameFormat.form(" ")).append(" ");
 518
 519       /*
 520        * Print out the match symbols:
 521        * | for exact match (ignoring case)
 522        * . if PAM250 score is positive
 523        * else a space
 524        */
 525       for (int i = 0; i < len; i++)
 526       {
 527         if ((i + (j * len)) < astr1.length())
 528         {
 529           char c1 = astr1.charAt(i + (j * len));
 530           char c2 = astr2.charAt(i + (j * len));
 531           boolean sameChar = Comparison.isSameResidue(c1, c2, false);
 532           if (sameChar && !Comparison.isGap(c1))
 533           {
 534             pid++;
 535             output.append("|");
 536           }
 537           else if (PEP.equals(type))
 538           {
 539             if (pam250.getPairwiseScore(c1, c2) > 0)
 540             {
 541               output.append(".");
 542             }
 543             else
 544             {
 545               output.append(" ");
 546             }
 547           }
 548           else
 549           {
 550             output.append(" ");
 551           }
 552         }
 553       }
 554
 555       // Now print the second aligned sequence
 556       output = output.append(NEWLINE);
 557       output = output.append(nameFormat.form(s2id)).append(" ");
 558
 559       for (int i = 0; i < len; i++)
 560       {
 561         if ((i + (j * len)) < astr2.length())
 562         {
 563           output.append(astr2.charAt(i + (j * len)));
 564         }
 565       }
 566
 567       output.append(NEWLINE).append(NEWLINE);
 568     }
 569
 570     pid = pid / (aseq1.length - count) * 100;
 571     output.append(new Format("Percentage ID = %3.2f\n").form(pid));
 572     output.append(NEWLINE);
 573     try
 574     {
 575       os.print(output.toString());
 576     } catch (Exception ex)
 577     {
 578     }
 579   }
 580
 581   /**
 582    * DOCUMENT ME!
 583    *
 584    * @param i
 585    *          DOCUMENT ME!
 586    * @param j
 587    *          DOCUMENT ME!
 588    *
 589    * @return DOCUMENT ME!
 590    */
 591   public int findTrace(int i, int j)
 592   {
 593     int t = 0;
 594     float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 595             s2str.charAt(j));
 596     float max = score[i - 1][j - 1] + (pairwiseScore * 10);
 597
 598     if (F[i][j] > max)
 599     {
 600       max = F[i][j];
 601       t = -1;
 602     }
 603     else if (F[i][j] == max)
 604     {
 605       if (prev == -1)
 606       {
 607         max = F[i][j];
 608         t = -1;
 609       }
 610     }
 611
 612     if (E[i][j] >= max)
 613     {
 614       max = E[i][j];
 615       t = 1;
 616     }
 617     else if (E[i][j] == max)
 618     {
 619       if (prev == 1)
 620       {
 621         max = E[i][j];
 622         t = 1;
 623       }
 624     }
 625
 626     prev = t;
 627
 628     return t;
 629   }
 630
 631   /**
 632    * DOCUMENT ME!
 633    */
 634   public void calcScoreMatrix()
 635   {
 636     int n = seq1.length;
 637     int m = seq2.length;
 638
 639     // top left hand element
 640     score[0][0] = scoreMatrix.getPairwiseScore(s1str.charAt(0),
 641             s2str.charAt(0)) * 10;
 642     E[0][0] = -GAP_EXTEND_COST;
 643     F[0][0] = 0;
 644
 645     // Calculate the top row first
 646     for (int j = 1; j < m; j++)
 647     {
 648       // What should these values be? 0 maybe
 649       E[0][j] = max(score[0][j - 1] - GAP_OPEN_COST,
 650               E[0][j - 1] - GAP_EXTEND_COST);
 651       F[0][j] = -GAP_EXTEND_COST;
 652
 653       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(0),
 654               s2str.charAt(j));
 655       score[0][j] = max(pairwiseScore * 10, -GAP_OPEN_COST,
 656               -GAP_EXTEND_COST);
 657
 658       traceback[0][j] = 1;
 659     }
 660
 661     // Now do the left hand column
 662     for (int i = 1; i < n; i++)
 663     {
 664       E[i][0] = -GAP_OPEN_COST;
 665       F[i][0] = max(score[i - 1][0] - GAP_OPEN_COST,
 666               F[i - 1][0] - GAP_EXTEND_COST);
 667
 668       float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 669               s2str.charAt(0));
 670       score[i][0] = max(pairwiseScore * 10, E[i][0], F[i][0]);
 671       traceback[i][0] = -1;
 672     }
 673
 674     // Now do all the other rows
 675     for (int i = 1; i < n; i++)
 676     {
 677       for (int j = 1; j < m; j++)
 678       {
 679         E[i][j] = max(score[i][j - 1] - GAP_OPEN_COST,
 680                 E[i][j - 1] - GAP_EXTEND_COST);
 681         F[i][j] = max(score[i - 1][j] - GAP_OPEN_COST,
 682                 F[i - 1][j] - GAP_EXTEND_COST);
 683
 684         float pairwiseScore = scoreMatrix.getPairwiseScore(s1str.charAt(i),
 685                 s2str.charAt(j));
 686         score[i][j] = max(score[i - 1][j - 1] + (pairwiseScore * 10),
 687                 E[i][j], F[i][j]);
 688         traceback[i][j] = findTrace(i, j);
 689       }
 690     }
 691   }
 692
 693   /**
 694    * Returns the given sequence with all of the given gap characters removed.
 695    *
 696    * @param gapChars
 697    *          a string of characters to be treated as gaps
 698    * @param seq
 699    *          the input sequence
 700    *
 701    * @return
 702    */
 703   public static String extractGaps(String gapChars, String seq)
 704   {
 705     if (gapChars == null || seq == null)
 706     {
 707       return null;
 708     }
 709     StringTokenizer str = new StringTokenizer(seq, gapChars);
 710     StringBuilder newString = new StringBuilder(seq.length());
 711
 712     while (str.hasMoreTokens())
 713     {
 714       newString.append(str.nextToken());
 715     }
 716
 717     return newString.toString();
 718   }
 719
 720   /**
 721    * DOCUMENT ME!
 722    *
 723    * @param f1
 724    *          DOCUMENT ME!
 725    * @param f2
 726    *          DOCUMENT ME!
 727    * @param f3
 728    *          DOCUMENT ME!
 729    *
 730    * @return DOCUMENT ME!
 731    */
 732   private static float max(float f1, float f2, float f3)
 733   {
 734     float max = f1;
 735
 736     if (f2 > f1)
 737     {
 738       max = f2;
 739     }
 740
 741     if (f3 > max)
 742     {
 743       max = f3;
 744     }
 745
 746     return max;
 747   }
 748
 749   /**
 750    * DOCUMENT ME!
 751    *
 752    * @param f1
 753    *          DOCUMENT ME!
 754    * @param f2
 755    *          DOCUMENT ME!
 756    *
 757    * @return DOCUMENT ME!
 758    */
 759   private static float max(float f1, float f2)
 760   {
 761     float max = f1;
 762
 763     if (f2 > f1)
 764     {
 765       max = f2;
 766     }
 767
 768     return max;
 769   }
 770
 771   /**
 772    * Converts the character string to an array of integers which are the
 773    * corresponding indices to the characters in the score matrix
 774    *
 775    * @param s
 776    *
 777    * @return
 778    */
 779   int[] indexEncode(String s)
 780   {
 781     int[] encoded = new int[s.length()];
 782
 783     for (int i = 0; i < s.length(); i++)
 784     {
 785       char c = s.charAt(i);
 786       encoded[i] = scoreMatrix.getMatrixIndex(c);
 787     }
 788
 789     return encoded;
 790   }
 791
 792   /**
 793    * DOCUMENT ME!
 794    *
 795    * @param g
 796    *          DOCUMENT ME!
 797    * @param mat
 798    *          DOCUMENT ME!
 799    * @param n
 800    *          DOCUMENT ME!
 801    * @param m
 802    *          DOCUMENT ME!
 803    * @param psize
 804    *          DOCUMENT ME!
 805    */
 806   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
 807           int psize)
 808   {
 809     // TODO method doesn't seem to be referenced anywhere delete??
 810     int max = -1000;
 811     int min = 1000;
 812
 813     for (int i = 0; i < n; i++)
 814     {
 815       for (int j = 0; j < m; j++)
 816       {
 817         if (mat[i][j] >= max)
 818         {
 819           max = mat[i][j];
 820         }
 821
 822         if (mat[i][j] <= min)
 823         {
 824           min = mat[i][j];
 825         }
 826       }
 827     }
 828
 829     System.out.println(max + " " + min);
 830
 831     for (int i = 0; i < n; i++)
 832     {
 833       for (int j = 0; j < m; j++)
 834       {
 835         int x = psize * i;
 836         int y = psize * j;
 837
 838         // System.out.println(mat[i][j]);
 839         float score = (float) (mat[i][j] - min) / (float) (max - min);
 840         g.setColor(new Color(score, 0, 0));
 841         g.fillRect(x, y, psize, psize);
 842
 843         // System.out.println(x + " " + y + " " + score);
 844       }
 845     }
 846   }
 847
 848   /**
 849    * Compute a globally optimal needleman and wunsch alignment between two
 850    * sequences
 851    *
 852    * @param s1
 853    * @param s2
 854    * @param type
 855    *          AlignSeq.DNA or AlignSeq.PEP
 856    */
 857   public static AlignSeq doGlobalNWAlignment(SequenceI s1, SequenceI s2,
 858           String type)
 859   {
 860     AlignSeq as = new AlignSeq(s1, s2, type);
 861
 862     as.calcScoreMatrix();
 863     as.traceAlignment();
 864     return as;
 865   }
 866
 867   /**
 868    *
 869    * @return mapping from positions in S1 to corresponding positions in S2
 870    */
 871   public jalview.datamodel.Mapping getMappingFromS1(boolean allowmismatch)
 872   {
 873     ArrayList<Integer> as1 = new ArrayList<Integer>(),
 874             as2 = new ArrayList<Integer>();
 875     int pdbpos = s2.getStart() + getSeq2Start() - 2;
 876     int alignpos = s1.getStart() + getSeq1Start() - 2;
 877     int lp2 = pdbpos - 3, lp1 = alignpos - 3;
 878     boolean lastmatch = false;
 879     // and now trace the alignment onto the atom set.
 880     for (int i = 0; i < astr1.length(); i++)
 881     {
 882       char c1 = astr1.charAt(i), c2 = astr2.charAt(i);
 883       if (c1 != '-')
 884       {
 885         alignpos++;
 886       }
 887
 888       if (c2 != '-')
 889       {
 890         pdbpos++;
 891       }
 892
 893       if (allowmismatch || c1 == c2)
 894       {
 895         // extend mapping interval
 896         if (lp1 + 1 != alignpos || lp2 + 1 != pdbpos)
 897         {
 898           as1.add(Integer.valueOf(alignpos));
 899           as2.add(Integer.valueOf(pdbpos));
 900         }
 901         lastmatch = true;
 902         lp1 = alignpos;
 903         lp2 = pdbpos;
 904       }
 905       else
 906       {
 907         // extend mapping interval
 908         if (lastmatch)
 909         {
 910           as1.add(Integer.valueOf(lp1));
 911           as2.add(Integer.valueOf(lp2));
 912         }
 913         lastmatch = false;
 914       }
 915     }
 916     // construct range pairs
 917
 918     int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)],
 919             mapseq2 = new int[as2.size() + (lastmatch ? 1 : 0)];
 920     int i = 0;
 921     for (Integer ip : as1)
 922     {
 923       mapseq1[i++] = ip;
 924     }
 925     ;
 926     i = 0;
 927     for (Integer ip : as2)
 928     {
 929       mapseq2[i++] = ip;
 930     }
 931     ;
 932     if (lastmatch)
 933     {
 934       mapseq1[mapseq1.length - 1] = alignpos;
 935       mapseq2[mapseq2.length - 1] = pdbpos;
 936     }
 937     MapList map = new MapList(mapseq1, mapseq2, 1, 1);
 938
 939     jalview.datamodel.Mapping mapping = new Mapping(map);
 940     mapping.setTo(s2);
 941     return mapping;
 942   }
 943
 944   /**
 945    * matches ochains against al and populates seqs with the best match between
 946    * each ochain and the set in al
 947    *
 948    * @param ochains
 949    * @param al
 950    * @param dnaOrProtein
 951    * @param removeOldAnnots
 952    *          when true, old annotation is cleared before new annotation
 953    *          transferred
 954    * @return List<List<SequenceI> originals, List<SequenceI> replacement,
 955    *         List<AlignSeq> alignment between each>
 956    */
 957   public static List<List<? extends Object>> replaceMatchingSeqsWith(
 958           List<SequenceI> seqs, List<AlignmentAnnotation> annotations,
 959           List<SequenceI> ochains, AlignmentI al, String dnaOrProtein,
 960           boolean removeOldAnnots)
 961   {
 962     List<SequenceI> orig = new ArrayList<SequenceI>(),
 963             repl = new ArrayList<SequenceI>();
 964     List<AlignSeq> aligs = new ArrayList<AlignSeq>();
 965     if (al != null && al.getHeight() > 0)
 966     {
 967       ArrayList<SequenceI> matches = new ArrayList<SequenceI>();
 968       ArrayList<AlignSeq> aligns = new ArrayList<AlignSeq>();
 969
 970       for (SequenceI sq : ochains)
 971       {
 972         SequenceI bestm = null;
 973         AlignSeq bestaseq = null;
 974         float bestscore = 0;
 975         for (SequenceI msq : al.getSequences())
 976         {
 977           AlignSeq aseq = doGlobalNWAlignment(msq, sq, dnaOrProtein);
 978           if (bestm == null || aseq.getMaxScore() > bestscore)
 979           {
 980             bestscore = aseq.getMaxScore();
 981             bestaseq = aseq;
 982             bestm = msq;
 983           }
 984         }
 985         // System.out.println("Best Score for " + (matches.size() + 1) + " :"
 986         // + bestscore);
 987         matches.add(bestm);
 988         aligns.add(bestaseq);
 989         al.deleteSequence(bestm);
 990       }
 991       for (int p = 0, pSize = seqs.size(); p < pSize; p++)
 992       {
 993         SequenceI sq, sp = seqs.get(p);
 994         int q;
 995         if ((q = ochains.indexOf(sp)) > -1)
 996         {
 997           seqs.set(p, sq = matches.get(q));
 998           orig.add(sp);
 999           repl.add(sq);
1000           sq.setName(sp.getName());
1001           sq.setDescription(sp.getDescription());
1002           Mapping sp2sq;
1003           sq.transferAnnotation(sp,
1004                   sp2sq = aligns.get(q).getMappingFromS1(false));
1005           aligs.add(aligns.get(q));
1006           int inspos = -1;
1007           for (int ap = 0; ap < annotations.size();)
1008           {
1009             if (annotations.get(ap).sequenceRef == sp)
1010             {
1011               if (inspos == -1)
1012               {
1013                 inspos = ap;
1014               }
1015               if (removeOldAnnots)
1016               {
1017                 annotations.remove(ap);
1018               }
1019               else
1020               {
1021                 AlignmentAnnotation alan = annotations.remove(ap);
1022                 alan.liftOver(sq, sp2sq);
1023                 alan.setSequenceRef(sq);
1024                 sq.addAlignmentAnnotation(alan);
1025               }
1026             }
1027             else
1028             {
1029               ap++;
1030             }
1031           }
1032           if (sq.getAnnotation() != null && sq.getAnnotation().length > 0)
1033           {
1034             annotations.addAll(inspos == -1 ? annotations.size() : inspos,
1035                     Arrays.asList(sq.getAnnotation()));
1036           }
1037         }
1038       }
1039     }
1040     return Arrays.asList(orig, repl, aligs);
1041   }
1042
1043   /**
1044    * compute the PID vector used by the redundancy filter.
1045    *
1046    * @param originalSequences
1047    *          - sequences in alignment that are to filtered
1048    * @param omitHidden
1049    *          - null or strings to be analysed (typically, visible portion of
1050    *          each sequence in alignment)
1051    * @param start
1052    *          - first column in window for calculation
1053    * @param end
1054    *          - last column in window for calculation
1055    * @param ungapped
1056    *          - if true then use ungapped sequence to compute PID
1057    * @return vector containing maximum PID for i-th sequence and any sequences
1058    *         longer than that seuqence
1059    */
1060   public static float[] computeRedundancyMatrix(
1061           SequenceI[] originalSequences, String[] omitHidden, int start,
1062           int end, boolean ungapped)
1063   {
1064     int height = originalSequences.length;
1065     float[] redundancy = new float[height];
1066     int[] lngth = new int[height];
1067     for (int i = 0; i < height; i++)
1068     {
1069       redundancy[i] = 0f;
1070       lngth[i] = -1;
1071     }
1072
1073     // long start = System.currentTimeMillis();
1074
1075     SimilarityParams pidParams = new SimilarityParams(true, true, true,
1076             true);
1077     float pid;
1078     String seqi, seqj;
1079     for (int i = 0; i < height; i++)
1080     {
1081
1082       for (int j = 0; j < i; j++)
1083       {
1084         if (i == j)
1085         {
1086           continue;
1087         }
1088
1089         if (omitHidden == null)
1090         {
1091           seqi = originalSequences[i].getSequenceAsString(start, end);
1092           seqj = originalSequences[j].getSequenceAsString(start, end);
1093         }
1094         else
1095         {
1096           seqi = omitHidden[i];
1097           seqj = omitHidden[j];
1098         }
1099         if (lngth[i] == -1)
1100         {
1101           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
1102           lngth[i] = ug.length();
1103           if (ungapped)
1104           {
1105             seqi = ug;
1106           }
1107         }
1108         if (lngth[j] == -1)
1109         {
1110           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
1111           lngth[j] = ug.length();
1112           if (ungapped)
1113           {
1114             seqj = ug;
1115           }
1116         }
1117         pid = (float) PIDModel.computePID(seqi, seqj, pidParams);
1118
1119         // use real sequence length rather than string length
1120         if (lngth[j] < lngth[i])
1121         {
1122           redundancy[j] = Math.max(pid, redundancy[j]);
1123         }
1124         else
1125         {
1126           redundancy[i] = Math.max(pid, redundancy[i]);
1127         }
1128
1129       }
1130     }
1131     return redundancy;
1132   }
1133 }