JAL-2759 Moved propagateInsertions out of HiddenColumns
[jalview.git] / 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\n").form(pid));
569     output.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<Integer>(),
866             as2 = new ArrayList<Integer>();
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     ;
918     i = 0;
919     for (Integer ip : as2)
920     {
921       mapseq2[i++] = ip;
922     }
923     ;
924     if (lastmatch)
925     {
926       mapseq1[mapseq1.length - 1] = alignpos;
927       mapseq2[mapseq2.length - 1] = pdbpos;
928     }
929     MapList map = new MapList(mapseq1, mapseq2, 1, 1);
930
931     jalview.datamodel.Mapping mapping = new Mapping(map);
932     mapping.setTo(s2);
933     return mapping;
934   }
935
936   /**
937    * matches ochains against al and populates seqs with the best match between
938    * each ochain and the set in al
939    * 
940    * @param ochains
941    * @param al
942    * @param dnaOrProtein
943    * @param removeOldAnnots
944    *          when true, old annotation is cleared before new annotation
945    *          transferred
946    * @return List<List<SequenceI> originals, List<SequenceI> replacement,
947    *         List<AlignSeq> alignment between each>
948    */
949   public static List<List<? extends Object>> replaceMatchingSeqsWith(
950           List<SequenceI> seqs, List<AlignmentAnnotation> annotations,
951           List<SequenceI> ochains, AlignmentI al, String dnaOrProtein,
952           boolean removeOldAnnots)
953   {
954     List<SequenceI> orig = new ArrayList<SequenceI>(),
955             repl = new ArrayList<SequenceI>();
956     List<AlignSeq> aligs = new ArrayList<AlignSeq>();
957     if (al != null && al.getHeight() > 0)
958     {
959       ArrayList<SequenceI> matches = new ArrayList<SequenceI>();
960       ArrayList<AlignSeq> aligns = new ArrayList<AlignSeq>();
961
962       for (SequenceI sq : ochains)
963       {
964         SequenceI bestm = null;
965         AlignSeq bestaseq = null;
966         float bestscore = 0;
967         for (SequenceI msq : al.getSequences())
968         {
969           AlignSeq aseq = doGlobalNWAlignment(msq, sq, dnaOrProtein);
970           if (bestm == null || aseq.getMaxScore() > bestscore)
971           {
972             bestscore = aseq.getMaxScore();
973             bestaseq = aseq;
974             bestm = msq;
975           }
976         }
977         // System.out.println("Best Score for " + (matches.size() + 1) + " :"
978         // + bestscore);
979         matches.add(bestm);
980         aligns.add(bestaseq);
981         al.deleteSequence(bestm);
982       }
983       for (int p = 0, pSize = seqs.size(); p < pSize; p++)
984       {
985         SequenceI sq, sp = seqs.get(p);
986         int q;
987         if ((q = ochains.indexOf(sp)) > -1)
988         {
989           seqs.set(p, sq = matches.get(q));
990           orig.add(sp);
991           repl.add(sq);
992           sq.setName(sp.getName());
993           sq.setDescription(sp.getDescription());
994           Mapping sp2sq;
995           sq.transferAnnotation(sp,
996                   sp2sq = aligns.get(q).getMappingFromS1(false));
997           aligs.add(aligns.get(q));
998           int inspos = -1;
999           for (int ap = 0; ap < annotations.size();)
1000           {
1001             if (annotations.get(ap).sequenceRef == sp)
1002             {
1003               if (inspos == -1)
1004               {
1005                 inspos = ap;
1006               }
1007               if (removeOldAnnots)
1008               {
1009                 annotations.remove(ap);
1010               }
1011               else
1012               {
1013                 AlignmentAnnotation alan = annotations.remove(ap);
1014                 alan.liftOver(sq, sp2sq);
1015                 alan.setSequenceRef(sq);
1016                 sq.addAlignmentAnnotation(alan);
1017               }
1018             }
1019             else
1020             {
1021               ap++;
1022             }
1023           }
1024           if (sq.getAnnotation() != null && sq.getAnnotation().length > 0)
1025           {
1026             annotations.addAll(inspos == -1 ? annotations.size() : inspos,
1027                     Arrays.asList(sq.getAnnotation()));
1028           }
1029         }
1030       }
1031     }
1032     return Arrays.asList(orig, repl, aligs);
1033   }
1034
1035   /**
1036    * compute the PID vector used by the redundancy filter.
1037    * 
1038    * @param originalSequences
1039    *          - sequences in alignment that are to filtered
1040    * @param omitHidden
1041    *          - null or strings to be analysed (typically, visible portion of
1042    *          each sequence in alignment)
1043    * @param start
1044    *          - first column in window for calculation
1045    * @param end
1046    *          - last column in window for calculation
1047    * @param ungapped
1048    *          - if true then use ungapped sequence to compute PID
1049    * @return vector containing maximum PID for i-th sequence and any sequences
1050    *         longer than that seuqence
1051    */
1052   public static float[] computeRedundancyMatrix(
1053           SequenceI[] originalSequences, String[] omitHidden, int start,
1054           int end, boolean ungapped)
1055   {
1056     int height = originalSequences.length;
1057     float[] redundancy = new float[height];
1058     int[] lngth = new int[height];
1059     for (int i = 0; i < height; i++)
1060     {
1061       redundancy[i] = 0f;
1062       lngth[i] = -1;
1063     }
1064
1065     // long start = System.currentTimeMillis();
1066
1067     SimilarityParams pidParams = new SimilarityParams(true, true, true,
1068             true);
1069     float pid;
1070     String seqi, seqj;
1071     for (int i = 0; i < height; i++)
1072     {
1073
1074       for (int j = 0; j < i; j++)
1075       {
1076         if (i == j)
1077         {
1078           continue;
1079         }
1080
1081         if (omitHidden == null)
1082         {
1083           seqi = originalSequences[i].getSequenceAsString(start, end);
1084           seqj = originalSequences[j].getSequenceAsString(start, end);
1085         }
1086         else
1087         {
1088           seqi = omitHidden[i];
1089           seqj = omitHidden[j];
1090         }
1091         if (lngth[i] == -1)
1092         {
1093           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
1094           lngth[i] = ug.length();
1095           if (ungapped)
1096           {
1097             seqi = ug;
1098           }
1099         }
1100         if (lngth[j] == -1)
1101         {
1102           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
1103           lngth[j] = ug.length();
1104           if (ungapped)
1105           {
1106             seqj = ug;
1107           }
1108         }
1109         pid = (float) PIDModel.computePID(seqi, seqj, pidParams);
1110
1111         // use real sequence length rather than string length
1112         if (lngth[j] < lngth[i])
1113         {
1114           redundancy[j] = Math.max(pid, redundancy[j]);
1115         }
1116         else
1117         {
1118           redundancy[i] = Math.max(pid, redundancy[i]);
1119         }
1120
1121       }
1122     }
1123     return redundancy;
1124   }
1125 }