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