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   float[][] lookup;
 116
 117   int gapIndex = 23;
 118
 119   StringBuffer output = new StringBuffer();
 120
 121   String type; // AlignSeq.PEP or AlignSeq.DNA
 122
 123   private ScoreMatrix scoreModel;
 124
 125   /**
 126    * Creates a new AlignSeq object.
 127    *
 128    * @param s1 first sequence for alignment
 129    * @param s2 second sequence for alignment
 130    * @param type 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(s1.getDatasetSequence() == null ? s1 : s1
 264             .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(s2.getDatasetSequence() == null ? s2 : s2
 278             .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("ALL GAPS: "
 320               + (s1str.length() == 0 ? s1.getName() : " ")
 321               + (s2str.length() == 0 ? s2.getName() : ""));
 322       return;
 323     }
 324
 325     seq1 = new int[s1str.length()];
 326
 327     seq2 = new int[s2str.length()];
 328
 329     score = new float[s1str.length()][s2str.length()];
 330
 331     E = new float[s1str.length()][s2str.length()];
 332
 333     F = new float[s1str.length()][s2str.length()];
 334     traceback = new int[s1str.length()][s2str.length()];
 335
 336     seq1 = indexEncode(s1str);
 337
 338     seq2 = indexEncode(s2str);
 339   }
 340
 341   private void setDefaultParams(String moleculeType)
 342   {
 343     if (!PEP.equals(moleculeType) && !DNA.equals(moleculeType))
 344     {
 345       output.append("Wrong type = dna or pep only");
 346       throw new Error(MessageManager.formatMessage(
 347               "error.unknown_type_dna_or_pep",
 348               new String[] { moleculeType }));
 349     }
 350
 351     type = moleculeType;
 352     scoreModel = ScoreModels.getInstance().getDefaultModel(
 353             PEP.equals(type));
 354     lookup = scoreModel.getMatrix();
 355     gapIndex = scoreModel.getMatrixIndex(' ');
 356   }
 357
 358   /**
 359    * DOCUMENT ME!
 360    */
 361   public void traceAlignment()
 362   {
 363     // Find the maximum score along the rhs or bottom row
 364     float max = -Float.MAX_VALUE;
 365
 366     for (int i = 0; i < seq1.length; i++)
 367     {
 368       if (score[i][seq2.length - 1] > max)
 369       {
 370         max = score[i][seq2.length - 1];
 371         maxi = i;
 372         maxj = seq2.length - 1;
 373       }
 374     }
 375
 376     for (int j = 0; j < seq2.length; j++)
 377     {
 378       if (score[seq1.length - 1][j] > max)
 379       {
 380         max = score[seq1.length - 1][j];
 381         maxi = seq1.length - 1;
 382         maxj = j;
 383       }
 384     }
 385
 386     // System.out.println(maxi + " " + maxj + " " + score[maxi][maxj]);
 387     int i = maxi;
 388     int j = maxj;
 389     int trace;
 390     maxscore = score[i][j] / 10;
 391
 392     seq1end = maxi + 1;
 393     seq2end = maxj + 1;
 394
 395     aseq1 = new int[seq1.length + seq2.length];
 396     aseq2 = new int[seq1.length + seq2.length];
 397
 398     StringBuilder sb1 = new StringBuilder(aseq1.length);
 399     StringBuilder sb2 = new StringBuilder(aseq2.length);
 400
 401     count = (seq1.length + seq2.length) - 1;
 402
 403     while (i > 0 && j > 0)
 404     {
 405       aseq1[count] = seq1[i];
 406       sb1.append(s1str.charAt(i));
 407       aseq2[count] = seq2[j];
 408       sb2.append(s2str.charAt(j));
 409
 410       trace = findTrace(i, j);
 411
 412       if (trace == 0)
 413       {
 414         i--;
 415         j--;
 416       }
 417       else if (trace == 1)
 418       {
 419         j--;
 420         aseq1[count] = gapIndex;
 421         sb1.replace(sb1.length() - 1, sb1.length(), "-");
 422       }
 423       else if (trace == -1)
 424       {
 425         i--;
 426         aseq2[count] = gapIndex;
 427         sb2.replace(sb2.length() - 1, sb2.length(), "-");
 428       }
 429
 430       count--;
 431     }
 432
 433     seq1start = i + 1;
 434     seq2start = j + 1;
 435
 436     if (aseq1[count] != gapIndex)
 437     {
 438       aseq1[count] = seq1[i];
 439       sb1.append(s1str.charAt(i));
 440     }
 441
 442     if (aseq2[count] != gapIndex)
 443     {
 444       aseq2[count] = seq2[j];
 445       sb2.append(s2str.charAt(j));
 446     }
 447
 448     /*
 449      * we built the character strings backwards, so now
 450      * reverse them to convert to sequence strings
 451      */
 452     astr1 = sb1.reverse().toString();
 453     astr2 = sb2.reverse().toString();
 454   }
 455
 456   /**
 457    * DOCUMENT ME!
 458    */
 459   public void printAlignment(java.io.PrintStream os)
 460   {
 461     // TODO: Use original sequence characters rather than re-translated
 462     // characters in output
 463     // Find the biggest id length for formatting purposes
 464     String s1id = s1.getName(), s2id = s2.getName();
 465     int maxid = s1.getName().length();
 466     if (s2.getName().length() > maxid)
 467     {
 468       maxid = s2.getName().length();
 469     }
 470     if (maxid > 30)
 471     {
 472       maxid = 30;
 473       // JAL-527 - truncate the sequence ids
 474       if (s1.getName().length() > maxid)
 475       {
 476         s1id = s1.getName().substring(0, 30);
 477       }
 478       if (s2.getName().length() > maxid)
 479       {
 480         s2id = s2.getName().substring(0, 30);
 481       }
 482     }
 483     int len = 72 - maxid - 1;
 484     int nochunks = ((aseq1.length - count) / len)
 485             + ((aseq1.length - count) % len > 0 ? 1 : 0);
 486     pid = 0;
 487
 488     output.append("Score = ").append(score[maxi][maxj]).append(NEWLINE);
 489     output.append("Length of alignment = ")
 490             .append(String.valueOf(aseq1.length - count)).append(NEWLINE);
 491     output.append("Sequence ");
 492     output.append(new Format("%" + maxid + "s").form(s1.getName()));
 493     output.append(" :  ").append(String.valueOf(s1.getStart()))
 494             .append(" - ").append(String.valueOf(s1.getEnd()));
 495     output.append(" (Sequence length = ")
 496             .append(String.valueOf(s1str.length())).append(")")
 497             .append(NEWLINE);
 498     output.append("Sequence ");
 499     output.append(new Format("%" + maxid + "s").form(s2.getName()));
 500     output.append(" :  ").append(String.valueOf(s2.getStart()))
 501             .append(" - ").append(String.valueOf(s2.getEnd()));
 502     output.append(" (Sequence length = ")
 503             .append(String.valueOf(s2str.length())).append(")")
 504             .append(NEWLINE).append(NEWLINE);
 505
 506     ScoreMatrix pam250 = ScoreModels.getInstance().getPam250();
 507
 508     for (int j = 0; j < nochunks; j++)
 509     {
 510       // Print the first aligned sequence
 511       output.append(new Format("%" + (maxid) + "s").form(s1id)).append(" ");
 512
 513       for (int i = 0; i < len; i++)
 514       {
 515         if ((i + (j * len)) < astr1.length())
 516         {
 517           output.append(astr1.charAt(i + (j * len)));
 518         }
 519       }
 520
 521       output.append(NEWLINE);
 522       output.append(new Format("%" + (maxid) + "s").form(" ")).append(" ");
 523
 524       /*
 525        * Print out the match symbols:
 526        * | for exact match (ignoring case)
 527        * . if PAM250 score is positive
 528        * else a space
 529        */
 530       for (int i = 0; i < len; i++)
 531       {
 532         if ((i + (j * len)) < astr1.length())
 533         {
 534           char c1 = astr1.charAt(i + (j * len));
 535           char c2 = astr2.charAt(i + (j * len));
 536           boolean sameChar = Comparison.isSameResidue(c1, c2, false);
 537           if (sameChar && !Comparison.isGap(c1))
 538           {
 539             pid++;
 540             output.append("|");
 541           }
 542           else if (type.equals("pep"))
 543           {
 544             if (pam250.getPairwiseScore(c1, c2) > 0)
 545             {
 546               output.append(".");
 547             }
 548             else
 549             {
 550               output.append(" ");
 551             }
 552           }
 553           else
 554           {
 555             output.append(" ");
 556           }
 557         }
 558       }
 559
 560       // Now print the second aligned sequence
 561       output = output.append(NEWLINE);
 562       output = output.append(new Format("%" + (maxid) + "s").form(s2id))
 563               .append(" ");
 564
 565       for (int i = 0; i < len; i++)
 566       {
 567         if ((i + (j * len)) < astr2.length())
 568         {
 569           output.append(astr2.charAt(i + (j * len)));
 570         }
 571       }
 572
 573       output.append(NEWLINE).append(NEWLINE);
 574     }
 575
 576     pid = pid / (aseq1.length - count) * 100;
 577     output = output.append(new Format("Percentage ID = %2.2f\n").form(pid));
 578     try
 579     {
 580       os.print(output.toString());
 581     } catch (Exception ex)
 582     {
 583     }
 584   }
 585
 586   /**
 587    * DOCUMENT ME!
 588    *
 589    * @param i
 590    *          DOCUMENT ME!
 591    * @param j
 592    *          DOCUMENT ME!
 593    *
 594    * @return DOCUMENT ME!
 595    */
 596   public int findTrace(int i, int j)
 597   {
 598     int t = 0;
 599     float max = score[i - 1][j - 1] + (lookup[seq1[i]][seq2[j]] * 10);
 600
 601     if (F[i][j] > max)
 602     {
 603       max = F[i][j];
 604       t = -1;
 605     }
 606     else if (F[i][j] == max)
 607     {
 608       if (prev == -1)
 609       {
 610         max = F[i][j];
 611         t = -1;
 612       }
 613     }
 614
 615     if (E[i][j] >= max)
 616     {
 617       max = E[i][j];
 618       t = 1;
 619     }
 620     else if (E[i][j] == max)
 621     {
 622       if (prev == 1)
 623       {
 624         max = E[i][j];
 625         t = 1;
 626       }
 627     }
 628
 629     prev = t;
 630
 631     return t;
 632   }
 633
 634   /**
 635    * DOCUMENT ME!
 636    */
 637   public void calcScoreMatrix()
 638   {
 639     int n = seq1.length;
 640     int m = seq2.length;
 641
 642     // top left hand element
 643     score[0][0] = lookup[seq1[0]][seq2[0]] * 10;
 644     E[0][0] = -gapExtend;
 645     F[0][0] = 0;
 646
 647     // Calculate the top row first
 648     for (int j = 1; j < m; j++)
 649     {
 650       // What should these values be? 0 maybe
 651       E[0][j] = max(score[0][j - 1] - gapOpen, E[0][j - 1] - gapExtend);
 652       F[0][j] = -gapExtend;
 653
 654       score[0][j] = max(lookup[seq1[0]][seq2[j]] * 10, -gapOpen, -gapExtend);
 655
 656       traceback[0][j] = 1;
 657     }
 658
 659     // Now do the left hand column
 660     for (int i = 1; i < n; i++)
 661     {
 662       E[i][0] = -gapOpen;
 663       F[i][0] = max(score[i - 1][0] - gapOpen, F[i - 1][0] - gapExtend);
 664
 665       score[i][0] = max(lookup[seq1[i]][seq2[0]] * 10, E[i][0], F[i][0]);
 666       traceback[i][0] = -1;
 667     }
 668
 669     // Now do all the other rows
 670     for (int i = 1; i < n; i++)
 671     {
 672       for (int j = 1; j < m; j++)
 673       {
 674         E[i][j] = max(score[i][j - 1] - gapOpen, E[i][j - 1] - gapExtend);
 675         F[i][j] = max(score[i - 1][j] - gapOpen, F[i - 1][j] - gapExtend);
 676
 677         score[i][j] = max(score[i - 1][j - 1]
 678                 + (lookup[seq1[i]][seq2[j]] * 10), E[i][j], F[i][j]);
 679         traceback[i][j] = findTrace(i, j);
 680       }
 681     }
 682   }
 683
 684   /**
 685    * Returns the given sequence with all of the given gap characters removed.
 686    *
 687    * @param gapChars
 688    *          a string of characters to be treated as gaps
 689    * @param seq
 690    *          the input sequence
 691    *
 692    * @return
 693    */
 694   public static String extractGaps(String gapChars, String seq)
 695   {
 696     if (gapChars == null || seq == null)
 697     {
 698       return null;
 699     }
 700     StringTokenizer str = new StringTokenizer(seq, gapChars);
 701     StringBuilder newString = new StringBuilder(seq.length());
 702
 703     while (str.hasMoreTokens())
 704     {
 705       newString.append(str.nextToken());
 706     }
 707
 708     return newString.toString();
 709   }
 710
 711   /**
 712    * DOCUMENT ME!
 713    *
 714    * @param f1
 715    *          DOCUMENT ME!
 716    * @param f2
 717    *          DOCUMENT ME!
 718    * @param f3
 719    *          DOCUMENT ME!
 720    *
 721    * @return DOCUMENT ME!
 722    */
 723   private static float max(float f1, float f2, float f3)
 724   {
 725     float max = f1;
 726
 727     if (f2 > f1)
 728     {
 729       max = f2;
 730     }
 731
 732     if (f3 > max)
 733     {
 734       max = f3;
 735     }
 736
 737     return max;
 738   }
 739
 740   /**
 741    * DOCUMENT ME!
 742    *
 743    * @param f1
 744    *          DOCUMENT ME!
 745    * @param f2
 746    *          DOCUMENT ME!
 747    *
 748    * @return DOCUMENT ME!
 749    */
 750   private static float max(float f1, float f2)
 751   {
 752     float max = f1;
 753
 754     if (f2 > f1)
 755     {
 756       max = f2;
 757     }
 758
 759     return max;
 760   }
 761
 762   /**
 763    * Converts the character string to an array of integers which are the
 764    * corresponding indices to the characters in the score matrix
 765    *
 766    * @param s
 767    *
 768    * @return
 769    */
 770   int[] indexEncode(String s)
 771   {
 772     int[] encoded = new int[s.length()];
 773
 774     for (int i = 0; i < s.length(); i++)
 775     {
 776       char c = s.charAt(i);
 777       encoded[i] = scoreModel.getMatrixIndex(c);
 778     }
 779
 780     return encoded;
 781   }
 782
 783   /**
 784    * DOCUMENT ME!
 785    *
 786    * @param g
 787    *          DOCUMENT ME!
 788    * @param mat
 789    *          DOCUMENT ME!
 790    * @param n
 791    *          DOCUMENT ME!
 792    * @param m
 793    *          DOCUMENT ME!
 794    * @param psize
 795    *          DOCUMENT ME!
 796    */
 797   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
 798           int psize)
 799   {
 800     // TODO method dosen't seem to be referenced anywhere delete??
 801     int max = -1000;
 802     int min = 1000;
 803
 804     for (int i = 0; i < n; i++)
 805     {
 806       for (int j = 0; j < m; j++)
 807       {
 808         if (mat[i][j] >= max)
 809         {
 810           max = mat[i][j];
 811         }
 812
 813         if (mat[i][j] <= min)
 814         {
 815           min = mat[i][j];
 816         }
 817       }
 818     }
 819
 820     System.out.println(max + " " + min);
 821
 822     for (int i = 0; i < n; i++)
 823     {
 824       for (int j = 0; j < m; j++)
 825       {
 826         int x = psize * i;
 827         int y = psize * j;
 828
 829         // System.out.println(mat[i][j]);
 830         float score = (float) (mat[i][j] - min) / (float) (max - min);
 831         g.setColor(new Color(score, 0, 0));
 832         g.fillRect(x, y, psize, psize);
 833
 834         // System.out.println(x + " " + y + " " + score);
 835       }
 836     }
 837   }
 838
 839   /**
 840    * Compute a globally optimal needleman and wunsch alignment between two
 841    * sequences
 842    *
 843    * @param s1
 844    * @param s2
 845    * @param type
 846    *          AlignSeq.DNA or AlignSeq.PEP
 847    */
 848   public static AlignSeq doGlobalNWAlignment(SequenceI s1, SequenceI s2,
 849           String type)
 850   {
 851     AlignSeq as = new AlignSeq(s1, s2, type);
 852
 853     as.calcScoreMatrix();
 854     as.traceAlignment();
 855     return as;
 856   }
 857
 858   /**
 859    *
 860    * @return mapping from positions in S1 to corresponding positions in S2
 861    */
 862   public jalview.datamodel.Mapping getMappingFromS1(boolean allowmismatch)
 863   {
 864     ArrayList<Integer> as1 = new ArrayList<Integer>(), as2 = new ArrayList<Integer>();
 865     int pdbpos = s2.getStart() + getSeq2Start() - 2;
 866     int alignpos = s1.getStart() + getSeq1Start() - 2;
 867     int lp2 = pdbpos - 3, lp1 = alignpos - 3;
 868     boolean lastmatch = false;
 869     // and now trace the alignment onto the atom set.
 870     for (int i = 0; i < astr1.length(); i++)
 871     {
 872       char c1 = astr1.charAt(i), c2 = astr2.charAt(i);
 873       if (c1 != '-')
 874       {
 875         alignpos++;
 876       }
 877
 878       if (c2 != '-')
 879       {
 880         pdbpos++;
 881       }
 882
 883       if (allowmismatch || c1 == c2)
 884       {
 885         // extend mapping interval
 886         if (lp1 + 1 != alignpos || lp2 + 1 != pdbpos)
 887         {
 888           as1.add(Integer.valueOf(alignpos));
 889           as2.add(Integer.valueOf(pdbpos));
 890         }
 891         lastmatch = true;
 892         lp1 = alignpos;
 893         lp2 = pdbpos;
 894       }
 895       else
 896       {
 897         // extend mapping interval
 898         if (lastmatch)
 899         {
 900           as1.add(Integer.valueOf(lp1));
 901           as2.add(Integer.valueOf(lp2));
 902         }
 903         lastmatch = false;
 904       }
 905     }
 906     // construct range pairs
 907
 908     int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)], mapseq2 = new int[as2
 909             .size() + (lastmatch ? 1 : 0)];
 910     int i = 0;
 911     for (Integer ip : as1)
 912     {
 913       mapseq1[i++] = ip;
 914     }
 915     ;
 916     i = 0;
 917     for (Integer ip : as2)
 918     {
 919       mapseq2[i++] = ip;
 920     }
 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<SequenceI>(), repl = new ArrayList<SequenceI>();
 953     List<AlignSeq> aligs = new ArrayList<AlignSeq>();
 954     if (al != null && al.getHeight() > 0)
 955     {
 956       ArrayList<SequenceI> matches = new ArrayList<SequenceI>();
 957       ArrayList<AlignSeq> aligns = new ArrayList<AlignSeq>();
 958
 959       for (SequenceI sq : ochains)
 960       {
 961         SequenceI bestm = null;
 962         AlignSeq bestaseq = null;
 963         float bestscore = 0;
 964         for (SequenceI msq : al.getSequences())
 965         {
 966           AlignSeq aseq = doGlobalNWAlignment(msq, sq, dnaOrProtein);
 967           if (bestm == null || aseq.getMaxScore() > bestscore)
 968           {
 969             bestscore = aseq.getMaxScore();
 970             bestaseq = aseq;
 971             bestm = msq;
 972           }
 973         }
 974         System.out.println("Best Score for " + (matches.size() + 1) + " :"
 975                 + bestscore);
 976         matches.add(bestm);
 977         aligns.add(bestaseq);
 978         al.deleteSequence(bestm);
 979       }
 980       for (int p = 0, pSize = seqs.size(); p < pSize; p++)
 981       {
 982         SequenceI sq, sp = seqs.get(p);
 983         int q;
 984         if ((q = ochains.indexOf(sp)) > -1)
 985         {
 986           seqs.set(p, sq = matches.get(q));
 987           orig.add(sp);
 988           repl.add(sq);
 989           sq.setName(sp.getName());
 990           sq.setDescription(sp.getDescription());
 991           Mapping sp2sq;
 992           sq.transferAnnotation(sp,
 993                   sp2sq = aligns.get(q).getMappingFromS1(false));
 994           aligs.add(aligns.get(q));
 995           int inspos = -1;
 996           for (int ap = 0; ap < annotations.size();)
 997           {
 998             if (annotations.get(ap).sequenceRef == sp)
 999             {
1000               if (inspos == -1)
1001               {
1002                 inspos = ap;
1003               }
1004               if (removeOldAnnots)
1005               {
1006                 annotations.remove(ap);
1007               }
1008               else
1009               {
1010                 AlignmentAnnotation alan = annotations.remove(ap);
1011                 alan.liftOver(sq, sp2sq);
1012                 alan.setSequenceRef(sq);
1013                 sq.addAlignmentAnnotation(alan);
1014               }
1015             }
1016             else
1017             {
1018               ap++;
1019             }
1020           }
1021           if (sq.getAnnotation() != null && sq.getAnnotation().length > 0)
1022           {
1023             annotations.addAll(inspos == -1 ? annotations.size() : inspos,
1024                     Arrays.asList(sq.getAnnotation()));
1025           }
1026         }
1027       }
1028     }
1029     return Arrays.asList(orig, repl, aligs);
1030   }
1031
1032   /**
1033    * compute the PID vector used by the redundancy filter.
1034    *
1035    * @param originalSequences
1036    *          - sequences in alignment that are to filtered
1037    * @param omitHidden
1038    *          - null or strings to be analysed (typically, visible portion of
1039    *          each sequence in alignment)
1040    * @param start
1041    *          - first column in window for calculation
1042    * @param end
1043    *          - last column in window for calculation
1044    * @param ungapped
1045    *          - if true then use ungapped sequence to compute PID
1046    * @return vector containing maximum PID for i-th sequence and any sequences
1047    *         longer than that seuqence
1048    */
1049   public static float[] computeRedundancyMatrix(
1050           SequenceI[] originalSequences, String[] omitHidden, int start,
1051           int end, boolean ungapped)
1052   {
1053     int height = originalSequences.length;
1054     float[] redundancy = new float[height];
1055     int[] lngth = new int[height];
1056     for (int i = 0; i < height; i++)
1057     {
1058       redundancy[i] = 0f;
1059       lngth[i] = -1;
1060     }
1061
1062     // long start = System.currentTimeMillis();
1063
1064     SimilarityParams pidParams = new SimilarityParams(true, true, true,
1065             true);
1066     float pid;
1067     String seqi, seqj;
1068     for (int i = 0; i < height; i++)
1069     {
1070
1071       for (int j = 0; j < i; j++)
1072       {
1073         if (i == j)
1074         {
1075           continue;
1076         }
1077
1078         if (omitHidden == null)
1079         {
1080           seqi = originalSequences[i].getSequenceAsString(start, end);
1081           seqj = originalSequences[j].getSequenceAsString(start, end);
1082         }
1083         else
1084         {
1085           seqi = omitHidden[i];
1086           seqj = omitHidden[j];
1087         }
1088         if (lngth[i] == -1)
1089         {
1090           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
1091           lngth[i] = ug.length();
1092           if (ungapped)
1093           {
1094             seqi = ug;
1095           }
1096         }
1097         if (lngth[j] == -1)
1098         {
1099           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
1100           lngth[j] = ug.length();
1101           if (ungapped)
1102           {
1103             seqj = ug;
1104           }
1105         }
1106         pid = (float) PIDModel.computePID(seqi, seqj, pidParams);
1107
1108         // use real sequence length rather than string length
1109         if (lngth[j] < lngth[i])
1110         {
1111           redundancy[j] = Math.max(pid, redundancy[j]);
1112         }
1113         else
1114         {
1115           redundancy[i] = Math.max(pid, redundancy[i]);
1116         }
1117
1118       }
1119     }
1120     return redundancy;
1121   }
1122 }