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