src/jalview/analysis/AlignSeq.java

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