JAL-1159 Merging 2.8.1 patches and Anne Menards work from 2012 into development branch
[jalview.git] / src / jalview / analysis / AlignSeq.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8)
3  * Copyright (C) 2012 J Procter, AM Waterhouse, LM Lui, J Engelhardt, G Barton, M Clamp, S Searle
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 of the License, or (at your option) any later version.
10  *  
11  * Jalview is distributed in the hope that it will be useful, but 
12  * WITHOUT ANY WARRANTY; without even the implied warranty 
13  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
14  * PURPOSE.  See the GNU General Public License for more details.
15  * 
16  * You should have received a copy of the GNU General Public License along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
17  */
18 package jalview.analysis;
19 import java.util.*;
20
21 import java.awt.*;
22
23 import jalview.datamodel.*;
24 import jalview.schemes.*;
25 import jalview.util.*;
26
27 /**
28  * 
29  * 
30  * @author $author$
31  * @version $Revision$
32  */
33 public class AlignSeq
34 {
35   public static final String PEP = "pep";
36
37   public static final String DNA = "dna";
38
39   static String[] dna =
40   { "A", "C", "G", "T", "-" };
41
42   // "C", "T", "A", "G", "-"};
43   static String[] pep =
44   { "A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F",
45       "P", "S", "T", "W", "Y", "V", "B", "Z", "X", "-" };
46
47   int[][] score;
48
49   int[][] E;
50
51   int[][] F;
52
53   int[][] traceback;
54
55   int[] seq1;
56
57   int[] seq2;
58
59   SequenceI s1;
60
61   SequenceI s2;
62
63   public String s1str;
64
65   public String s2str;
66
67   int maxi;
68
69   int maxj;
70
71   int[] aseq1;
72
73   int[] aseq2;
74
75   public String astr1 = "";
76
77   public String astr2 = "";
78
79   /** DOCUMENT ME!! */
80   public int seq1start;
81
82   /** DOCUMENT ME!! */
83   public int seq1end;
84
85   /** DOCUMENT ME!! */
86   public int seq2start;
87
88   /** DOCUMENT ME!! */
89   public int seq2end;
90
91   int count;
92
93   /** DOCUMENT ME!! */
94   public int maxscore;
95
96   float pid;
97
98   int prev = 0;
99
100   int gapOpen = 120;
101
102   int gapExtend = 20;
103
104   int[][] lookup = ResidueProperties.getBLOSUM62();
105
106   String[] intToStr = pep;
107
108   int defInt = 23;
109
110   StringBuffer output = new StringBuffer();
111
112   String type;
113
114   private int[] charToInt;
115
116   /**
117    * Creates a new AlignSeq object.
118    * 
119    * @param s1
120    *          DOCUMENT ME!
121    * @param s2
122    *          DOCUMENT ME!
123    * @param type
124    *          DOCUMENT ME!
125    */
126   public AlignSeq(SequenceI s1, SequenceI s2, String type)
127   {
128     SeqInit(s1, s1.getSequenceAsString(), s2, s2.getSequenceAsString(),
129             type);
130   }
131
132   /**
133    * Creates a new AlignSeq object.
134    * 
135    * @param s1
136    *          DOCUMENT ME!
137    * @param s2
138    *          DOCUMENT ME!
139    * @param type
140    *          DOCUMENT ME!
141    */
142   public AlignSeq(SequenceI s1, String string1, SequenceI s2,
143           String string2, String type)
144   {
145     SeqInit(s1, string1.toUpperCase(), s2, string2.toUpperCase(), type);
146   }
147
148   /**
149    * DOCUMENT ME!
150    * 
151    * @return DOCUMENT ME!
152    */
153   public int getMaxScore()
154   {
155     return maxscore;
156   }
157
158   /**
159    * DOCUMENT ME!
160    * 
161    * @return DOCUMENT ME!
162    */
163   public int getSeq2Start()
164   {
165     return seq2start;
166   }
167
168   /**
169    * DOCUMENT ME!
170    * 
171    * @return DOCUMENT ME!
172    */
173   public int getSeq2End()
174   {
175     return seq2end;
176   }
177
178   /**
179    * DOCUMENT ME!
180    * 
181    * @return DOCUMENT ME!
182    */
183   public int getSeq1Start()
184   {
185     return seq1start;
186   }
187
188   /**
189    * DOCUMENT ME!
190    * 
191    * @return DOCUMENT ME!
192    */
193   public int getSeq1End()
194   {
195     return seq1end;
196   }
197
198   /**
199    * DOCUMENT ME!
200    * 
201    * @return DOCUMENT ME!
202    */
203   public String getOutput()
204   {
205     return output.toString();
206   }
207
208   /**
209    * DOCUMENT ME!
210    * 
211    * @return DOCUMENT ME!
212    */
213   public String getAStr1()
214   {
215     return astr1;
216   }
217
218   /**
219    * DOCUMENT ME!
220    * 
221    * @return DOCUMENT ME!
222    */
223   public String getAStr2()
224   {
225     return astr2;
226   }
227
228   /**
229    * DOCUMENT ME!
230    * 
231    * @return DOCUMENT ME!
232    */
233   public int[] getASeq1()
234   {
235     return aseq1;
236   }
237
238   /**
239    * DOCUMENT ME!
240    * 
241    * @return DOCUMENT ME!
242    */
243   public int[] getASeq2()
244   {
245     return aseq2;
246   }
247
248   /**
249    * DOCUMENT ME!
250    * 
251    * @return DOCUMENT ME!
252    */
253   public SequenceI getS1()
254   {
255     return s1;
256   }
257
258   /**
259    * DOCUMENT ME!
260    * 
261    * @return DOCUMENT ME!
262    */
263   public SequenceI getS2()
264   {
265     return s2;
266   }
267
268   /**
269    * Construct score matrix for sequences with standard DNA or PEPTIDE matrix
270    * @param s1 - sequence 1
271    * @param string1 - string to use for s1
272    * @param s2 - sequence 2
273    * @param string2 - string to use for s2
274    * @param type
275    *          DNA or PEPTIDE
276    */
277   public void SeqInit(SequenceI s1, String string1, SequenceI s2,
278           String string2, String type)
279   {
280     this.s1 = s1;
281     this.s2 = s2;
282     setDefaultParams(type);
283     SeqInit(string1, string2);
284   }
285
286   /**
287    * Construct score matrix for sequences with custom substitution matrix
288    * @param s1 - sequence 1
289    * @param string1 - string to use for s1
290    * @param s2 - sequence 2
291    * @param string2 - string to use for s2
292    * @param scoreMatrix - substitution matrix to use for alignment
293    */
294   public void SeqInit(SequenceI s1, String string1, SequenceI s2,
295           String string2, ScoreMatrix scoreMatrix)
296   {
297     this.s1 = s1;
298     this.s2 = s2;
299     setType(scoreMatrix.isDNA() ? AlignSeq.DNA : AlignSeq.PEP);
300     lookup = scoreMatrix.getMatrix();
301   }
302
303   /**
304    * construct score matrix for string1 and string2 (after removing any existing
305    * gaps
306    * 
307    * @param string1 
308    * @param string2
309    */
310   private void SeqInit(String string1, String string2)
311   {
312     s1str = extractGaps(jalview.util.Comparison.GapChars, string1);
313     s2str = extractGaps(jalview.util.Comparison.GapChars, string2);
314
315     if (s1str.length() == 0 || s2str.length() == 0)
316     {
317       output.append("ALL GAPS: "
318               + (s1str.length() == 0 ? s1.getName() : " ")
319               + (s2str.length() == 0 ? s2.getName() : ""));
320       return;
321     }
322
323     // System.out.println("lookuip " + rt.freeMemory() + " "+ rt.totalMemory());
324     seq1 = new int[s1str.length()];
325
326     // System.out.println("seq1 " + rt.freeMemory() +" " + rt.totalMemory());
327     seq2 = new int[s2str.length()];
328
329     // System.out.println("seq2 " + rt.freeMemory() + " " + rt.totalMemory());
330     score = new int[s1str.length()][s2str.length()];
331
332     // System.out.println("score " + rt.freeMemory() + " " + rt.totalMemory());
333     E = new int[s1str.length()][s2str.length()];
334
335     // System.out.println("E " + rt.freeMemory() + " " + rt.totalMemory());
336     F = new int[s1str.length()][s2str.length()];
337     traceback = new int[s1str.length()][s2str.length()];
338
339     // System.out.println("F " + rt.freeMemory() + " " + rt.totalMemory());
340     seq1 = stringToInt(s1str, type);
341
342     // System.out.println("seq1 " + rt.freeMemory() + " " + rt.totalMemory());
343     seq2 = stringToInt(s2str, type);
344
345     // System.out.println("Seq2 " + rt.freeMemory() + " " + rt.totalMemory());
346     // long tstart = System.currentTimeMillis();
347     // calcScoreMatrix();
348     // long tend = System.currentTimeMillis();
349     // System.out.println("Time take to calculate score matrix = " +
350     // (tend-tstart) + " ms");
351     // printScoreMatrix(score);
352     // System.out.println();
353     // printScoreMatrix(traceback);
354     // System.out.println();
355     // printScoreMatrix(E);
356     // System.out.println();
357     // /printScoreMatrix(F);
358     // System.out.println();
359     // tstart = System.currentTimeMillis();
360     // traceAlignment();
361     // tend = System.currentTimeMillis();
362     // System.out.println("Time take to traceback alignment = " + (tend-tstart)
363     // + " ms");
364   }
365
366   private void setDefaultParams(String type)
367   {
368     setType(type);
369
370     if (type.equals(AlignSeq.PEP))
371     {
372       lookup = ResidueProperties.getDefaultPeptideMatrix();
373     }
374     else if (type.equals(AlignSeq.DNA))
375     {
376       lookup = ResidueProperties.getDefaultDnaMatrix();
377     }
378   }
379
380   private void setType(String type2)
381   {
382     this.type = type2;
383     if (type.equals(AlignSeq.PEP))
384     {
385       intToStr = pep;
386       charToInt = ResidueProperties.aaIndex;
387       defInt = ResidueProperties.maxProteinIndex;
388     }
389     else if (type.equals(AlignSeq.DNA))
390     {
391       intToStr = dna;
392       charToInt = ResidueProperties.nucleotideIndex;
393       defInt = ResidueProperties.maxNucleotideIndex;
394     }
395     else
396     {
397       output.append("Wrong type = dna or pep only");
398       throw new Error("Unknown Type " + type2
399               + " - dna or pep are the only allowed values.");
400     }
401   }
402
403   /**
404    * DOCUMENT ME!
405    */
406   public void traceAlignment()
407   {
408     // Find the maximum score along the rhs or bottom row
409     int max = -9999;
410
411     for (int i = 0; i < seq1.length; i++)
412     {
413       if (score[i][seq2.length - 1] > max)
414       {
415         max = score[i][seq2.length - 1];
416         maxi = i;
417         maxj = seq2.length - 1;
418       }
419     }
420
421     for (int j = 0; j < seq2.length; j++)
422     {
423       if (score[seq1.length - 1][j] > max)
424       {
425         max = score[seq1.length - 1][j];
426         maxi = seq1.length - 1;
427         maxj = j;
428       }
429     }
430
431     // System.out.println(maxi + " " + maxj + " " + score[maxi][maxj]);
432     int i = maxi;
433     int j = maxj;
434     int trace;
435     maxscore = score[i][j] / 10;
436
437     seq1end = maxi + 1;
438     seq2end = maxj + 1;
439
440     aseq1 = new int[seq1.length + seq2.length];
441     aseq2 = new int[seq1.length + seq2.length];
442
443     count = (seq1.length + seq2.length) - 1;
444
445     while ((i > 0) && (j > 0))
446     {
447       if ((aseq1[count] != defInt) && (i >= 0))
448       {
449         aseq1[count] = seq1[i];
450         astr1 = s1str.charAt(i) + astr1;
451       }
452
453       if ((aseq2[count] != defInt) && (j > 0))
454       {
455         aseq2[count] = seq2[j];
456         astr2 = s2str.charAt(j) + astr2;
457       }
458
459       trace = findTrace(i, j);
460
461       if (trace == 0)
462       {
463         i--;
464         j--;
465       }
466       else if (trace == 1)
467       {
468         j--;
469         aseq1[count] = defInt;
470         astr1 = "-" + astr1.substring(1);
471       }
472       else if (trace == -1)
473       {
474         i--;
475         aseq2[count] = defInt;
476         astr2 = "-" + astr2.substring(1);
477       }
478
479       count--;
480     }
481
482     seq1start = i + 1;
483     seq2start = j + 1;
484
485     if (aseq1[count] != defInt)
486     {
487       aseq1[count] = seq1[i];
488       astr1 = s1str.charAt(i) + astr1;
489     }
490
491     if (aseq2[count] != defInt)
492     {
493       aseq2[count] = seq2[j];
494       astr2 = s2str.charAt(j) + astr2;
495     }
496   }
497
498   /**
499    * DOCUMENT ME!
500    */
501   public void printAlignment(java.io.PrintStream os)
502   {
503     // TODO: Use original sequence characters rather than re-translated
504     // characters in output
505     // Find the biggest id length for formatting purposes
506     String s1id = s1.getName(), s2id = s2.getName();
507     int maxid = s1.getName().length();
508     if (s2.getName().length() > maxid)
509     {
510       maxid = s2.getName().length();
511     }
512     if (maxid > 30)
513     {
514       maxid = 30;
515       // JAL-527 - truncate the sequence ids
516       if (s1.getName().length() > maxid)
517       {
518         s1id = s1.getName().substring(0, 30);
519       }
520       if (s2.getName().length() > maxid)
521       {
522         s2id = s2.getName().substring(0, 30);
523       }
524     }
525     int len = 72 - maxid - 1;
526     int nochunks = ((aseq1.length - count) / len) + 1;
527     pid = 0;
528
529     output.append("Score = " + score[maxi][maxj] + "\n");
530     output.append("Length of alignment = " + (aseq1.length - count) + "\n");
531     output.append("Sequence ");
532     output.append(new Format("%" + maxid + "s").form(s1.getName()));
533     output.append(" :  " + s1.getStart() + " - " + s1.getEnd()
534             + " (Sequence length = " + s1str.length() + ")\n");
535     output.append("Sequence ");
536     output.append(new Format("%" + maxid + "s").form(s2.getName()));
537     output.append(" :  " + s2.getStart() + " - " + s2.getEnd()
538             + " (Sequence length = " + s2str.length() + ")\n\n");
539
540     for (int j = 0; j < nochunks; j++)
541     {
542       // Print the first aligned sequence
543       output.append(new Format("%" + (maxid) + "s").form(s1id) + " ");
544
545       for (int i = 0; i < len; i++)
546       {
547         if ((i + (j * len)) < astr1.length())
548         {
549           output.append(astr1.charAt(i + (j * len)));
550         }
551       }
552
553       output.append("\n");
554       output.append(new Format("%" + (maxid) + "s").form(" ") + " ");
555
556       // Print out the matching chars
557       for (int i = 0; i < len; i++)
558       {
559         if ((i + (j * len)) < astr1.length())
560         {
561           if (astr1.charAt(i + (j * len)) == astr2.charAt(i + (j * len))
562                   && !jalview.util.Comparison.isGap(astr1.charAt(i
563                           + (j * len))))
564           {
565             pid++;
566             output.append("|");
567           }
568           else if (type.equals("pep"))
569           {
570             if (ResidueProperties.getPAM250(astr1.charAt(i + (j * len)),
571                     astr2.charAt(i + (j * len))) > 0)
572             {
573               output.append(".");
574             }
575             else
576             {
577               output.append(" ");
578             }
579           }
580           else
581           {
582             output.append(" ");
583           }
584         }
585       }
586
587       // Now print the second aligned sequence
588       output = output.append("\n");
589       output = output.append(new Format("%" + (maxid) + "s").form(s2id)
590               + " ");
591
592       for (int i = 0; i < len; i++)
593       {
594         if ((i + (j * len)) < astr2.length())
595         {
596           output.append(astr2.charAt(i + (j * len)));
597         }
598       }
599
600       output = output.append("\n\n");
601     }
602
603     pid = pid / (float) (aseq1.length - count) * 100;
604     output = output.append(new Format("Percentage ID = %2.2f\n\n")
605             .form(pid));
606
607     try
608     {
609       os.print(output.toString());
610     } catch (Exception ex)
611     {
612     }
613   }
614
615   /**
616    * DOCUMENT ME!
617    * 
618    * @param mat
619    *          DOCUMENT ME!
620    */
621   public void printScoreMatrix(int[][] mat)
622   {
623     int n = seq1.length;
624     int m = seq2.length;
625
626     for (int i = 0; i < n; i++)
627     {
628       // Print the top sequence
629       if (i == 0)
630       {
631         Format.print(System.out, "%8s", s2str.substring(0, 1));
632
633         for (int jj = 1; jj < m; jj++)
634         {
635           Format.print(System.out, "%5s", s2str.substring(jj, jj + 1));
636         }
637
638         System.out.println();
639       }
640
641       for (int j = 0; j < m; j++)
642       {
643         if (j == 0)
644         {
645           Format.print(System.out, "%3s", s1str.substring(i, i + 1));
646         }
647
648         Format.print(System.out, "%3d ", mat[i][j] / 10);
649       }
650
651       System.out.println();
652     }
653   }
654
655   /**
656    * DOCUMENT ME!
657    * 
658    * @param i
659    *          DOCUMENT ME!
660    * @param j
661    *          DOCUMENT ME!
662    * 
663    * @return DOCUMENT ME!
664    */
665   public int findTrace(int i, int j)
666   {
667     int t = 0;
668     int max = score[i - 1][j - 1] + (lookup[seq1[i]][seq2[j]] * 10);
669
670     if (F[i][j] > max)
671     {
672       max = F[i][j];
673       t = -1;
674     }
675     else if (F[i][j] == max)
676     {
677       if (prev == -1)
678       {
679         max = F[i][j];
680         t = -1;
681       }
682     }
683
684     if (E[i][j] >= max)
685     {
686       max = E[i][j];
687       t = 1;
688     }
689     else if (E[i][j] == max)
690     {
691       if (prev == 1)
692       {
693         max = E[i][j];
694         t = 1;
695       }
696     }
697
698     prev = t;
699
700     return t;
701   }
702
703   /**
704    * DOCUMENT ME!
705    */
706   public void calcScoreMatrix()
707   {
708     int n = seq1.length;
709     int m = seq2.length;
710
711     // top left hand element
712     score[0][0] = lookup[seq1[0]][seq2[0]] * 10;
713     E[0][0] = -gapExtend;
714     F[0][0] = 0;
715
716     // Calculate the top row first
717     for (int j = 1; j < m; j++)
718     {
719       // What should these values be? 0 maybe
720       E[0][j] = max(score[0][j - 1] - gapOpen, E[0][j - 1] - gapExtend);
721       F[0][j] = -gapExtend;
722
723       score[0][j] = max(lookup[seq1[0]][seq2[j]] * 10, -gapOpen, -gapExtend);
724
725       traceback[0][j] = 1;
726     }
727
728     // Now do the left hand column
729     for (int i = 1; i < n; i++)
730     {
731       E[i][0] = -gapOpen;
732       F[i][0] = max(score[i - 1][0] - gapOpen, F[i - 1][0] - gapExtend);
733
734       score[i][0] = max(lookup[seq1[i]][seq2[0]] * 10, E[i][0], F[i][0]);
735       traceback[i][0] = -1;
736     }
737
738     // Now do all the other rows
739     for (int i = 1; i < n; i++)
740     {
741       for (int j = 1; j < m; j++)
742       {
743         E[i][j] = max(score[i][j - 1] - gapOpen, E[i][j - 1] - gapExtend);
744         F[i][j] = max(score[i - 1][j] - gapOpen, F[i - 1][j] - gapExtend);
745
746         score[i][j] = max(score[i - 1][j - 1]
747                 + (lookup[seq1[i]][seq2[j]] * 10), E[i][j], F[i][j]);
748         traceback[i][j] = findTrace(i, j);
749       }
750     }
751   }
752
753   /**
754    * DOCUMENT ME!
755    * 
756    * @param gapChar
757    *          DOCUMENT ME!
758    * @param seq
759    *          DOCUMENT ME!
760    * 
761    * @return DOCUMENT ME!
762    */
763   public static String extractGaps(String gapChar, String seq)
764   {
765     StringTokenizer str = new StringTokenizer(seq, gapChar);
766     StringBuffer newString = new StringBuffer();
767
768     while (str.hasMoreTokens())
769     {
770       newString.append(str.nextToken());
771     }
772
773     return newString.toString();
774   }
775
776   /**
777    * DOCUMENT ME!
778    * 
779    * @param i1
780    *          DOCUMENT ME!
781    * @param i2
782    *          DOCUMENT ME!
783    * @param i3
784    *          DOCUMENT ME!
785    * 
786    * @return DOCUMENT ME!
787    */
788   public int max(int i1, int i2, int i3)
789   {
790     int max = i1;
791
792     if (i2 > i1)
793     {
794       max = i2;
795     }
796
797     if (i3 > max)
798     {
799       max = i3;
800     }
801
802     return max;
803   }
804
805   /**
806    * DOCUMENT ME!
807    * 
808    * @param i1
809    *          DOCUMENT ME!
810    * @param i2
811    *          DOCUMENT ME!
812    * 
813    * @return DOCUMENT ME!
814    */
815   public int max(int i1, int i2)
816   {
817     int max = i1;
818
819     if (i2 > i1)
820     {
821       max = i2;
822     }
823
824     return max;
825   }
826
827   /**
828    * DOCUMENT ME!
829    * 
830    * @param s
831    *          DOCUMENT ME!
832    * @param type
833    *          DOCUMENT ME!
834    * 
835    * @return DOCUMENT ME!
836    */
837   public int[] stringToInt(String s, String type)
838   {
839     int[] seq1 = new int[s.length()];
840
841     for (int i = 0; i < s.length(); i++)
842     {
843       // String ss = s.substring(i, i + 1).toUpperCase();
844       char c = s.charAt(i);
845       if ('a' <= c && c <= 'z')
846       {
847         // TO UPPERCASE !!!
848         c -= ('a' - 'A');
849       }
850
851       try
852       {
853         seq1[i] = charToInt[c]; // set accordingly from setType
854         if (seq1[i] < 0 || seq1[i] > defInt) // set from setType: 23 for
855                                              // peptides, or 4 for NA.
856         {
857           seq1[i] = defInt;
858         }
859
860       } catch (Exception e)
861       {
862         seq1[i] = defInt;
863       }
864     }
865
866     return seq1;
867   }
868
869   /**
870    * DOCUMENT ME!
871    * 
872    * @param g
873    *          DOCUMENT ME!
874    * @param mat
875    *          DOCUMENT ME!
876    * @param n
877    *          DOCUMENT ME!
878    * @param m
879    *          DOCUMENT ME!
880    * @param psize
881    *          DOCUMENT ME!
882    */
883   public static void displayMatrix(Graphics g, int[][] mat, int n, int m,
884           int psize)
885   {
886     int max = -1000;
887     int min = 1000;
888
889     for (int i = 0; i < n; i++)
890     {
891       for (int j = 0; j < m; j++)
892       {
893         if (mat[i][j] >= max)
894         {
895           max = mat[i][j];
896         }
897
898         if (mat[i][j] <= min)
899         {
900           min = mat[i][j];
901         }
902       }
903     }
904
905     System.out.println(max + " " + min);
906
907     for (int i = 0; i < n; i++)
908     {
909       for (int j = 0; j < m; j++)
910       {
911         int x = psize * i;
912         int y = psize * j;
913
914         // System.out.println(mat[i][j]);
915         float score = (float) (mat[i][j] - min) / (float) (max - min);
916         g.setColor(new Color(score, 0, 0));
917         g.fillRect(x, y, psize, psize);
918
919         // System.out.println(x + " " + y + " " + score);
920       }
921     }
922   }
923
924   /**
925    * Compute a globally optimal needleman and wunsch alignment between two
926    * sequences
927    * 
928    * @param s1
929    * @param s2
930    * @param type
931    *          AlignSeq.DNA or AlignSeq.PEP
932    */
933   public static AlignSeq doGlobalNWAlignment(SequenceI s1, SequenceI s2,
934           String type)
935   {
936     AlignSeq as = new AlignSeq(s1, s2, type);
937
938     as.calcScoreMatrix();
939     as.traceAlignment();
940     return as;
941   }
942
943   /**
944    * 
945    * @return mapping from positions in S1 to corresponding positions in S2
946    */
947   public jalview.datamodel.Mapping getMappingFromS1(boolean allowmismatch)
948   {
949     ArrayList<Integer> as1 = new ArrayList<Integer>(), as2 = new ArrayList<Integer>();
950     int pdbpos = s2.getStart() + getSeq2Start() - 2;
951     int alignpos = s1.getStart() + getSeq1Start() - 2;
952     int lp2 = pdbpos - 3, lp1 = alignpos - 3;
953     boolean lastmatch = false;
954     // and now trace the alignment onto the atom set.
955     for (int i = 0; i < astr1.length(); i++)
956     {
957       char c1 = astr1.charAt(i), c2 = astr2.charAt(i);
958       if (c1 != '-')
959       {
960         alignpos++;
961       }
962
963       if (c2 != '-')
964       {
965         pdbpos++;
966       }
967
968       if (allowmismatch || c1 == c2)
969       {
970         lastmatch = true;
971         // extend mapping interval.
972         if (lp1 + 1 != alignpos || lp2+1 !=pdbpos)
973         {
974           as1.add(Integer.valueOf(alignpos));
975           as2.add(Integer.valueOf(pdbpos));
976         }
977         lp1 = alignpos;
978         lp2 = pdbpos;
979       }
980       else
981       {
982         lastmatch = false;
983       }
984     }
985     // construct range pairs
986     int[] mapseq1 = new int[as1.size() + (lastmatch ? 1 : 0)], mapseq2 = new int[as2
987             .size() + (lastmatch ? 1 : 0)];
988     int i = 0;
989     for (Integer ip : as1)
990     {
991       mapseq1[i++] = ip;
992     }
993     ;
994     i = 0;
995     for (Integer ip : as2)
996     {
997       mapseq2[i++] = ip;
998     }
999     ;
1000     if (lastmatch)
1001     {
1002       mapseq1[mapseq1.length - 1] = alignpos;
1003       mapseq2[mapseq2.length - 1] = pdbpos;
1004     }
1005     MapList map = new MapList(mapseq1, mapseq2, 1, 1);
1006
1007     jalview.datamodel.Mapping mapping = new Mapping(map);
1008     mapping.setTo(s2);
1009     return mapping;
1010   }
1011
1012   /**
1013    * compute the PID vector used by the redundancy filter.
1014    * 
1015    * @param originalSequences
1016    *          - sequences in alignment that are to filtered
1017    * @param omitHidden
1018    *          - null or strings to be analysed (typically, visible portion of
1019    *          each sequence in alignment)
1020    * @param start
1021    *          - first column in window for calculation
1022    * @param end
1023    *          - last column in window for calculation
1024    * @param ungapped
1025    *          - if true then use ungapped sequence to compute PID
1026    * @return vector containing maximum PID for i-th sequence and any sequences
1027    *         longer than that seuqence
1028    */
1029   public static float[] computeRedundancyMatrix(
1030           SequenceI[] originalSequences, String[] omitHidden, int start,
1031           int end, boolean ungapped)
1032   {
1033     int height = originalSequences.length;
1034     float[] redundancy = new float[height];
1035     int[] lngth = new int[height];
1036     for (int i = 0; i < height; i++)
1037     {
1038       redundancy[i] = 0f;
1039       lngth[i] = -1;
1040     }
1041
1042     // long start = System.currentTimeMillis();
1043
1044     float pid;
1045     String seqi, seqj;
1046     for (int i = 0; i < height; i++)
1047     {
1048
1049       for (int j = 0; j < i; j++)
1050       {
1051         if (i == j)
1052         {
1053           continue;
1054         }
1055
1056         if (omitHidden == null)
1057         {
1058           seqi = originalSequences[i].getSequenceAsString(start, end);
1059           seqj = originalSequences[j].getSequenceAsString(start, end);
1060         }
1061         else
1062         {
1063           seqi = omitHidden[i];
1064           seqj = omitHidden[j];
1065         }
1066         if (lngth[i] == -1)
1067         {
1068           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqi);
1069           lngth[i] = ug.length();
1070           if (ungapped)
1071           {
1072             seqi = ug;
1073           }
1074         }
1075         if (lngth[j] == -1)
1076         {
1077           String ug = AlignSeq.extractGaps(Comparison.GapChars, seqj);
1078           lngth[j] = ug.length();
1079           if (ungapped)
1080           {
1081             seqj = ug;
1082           }
1083         }
1084         pid = Comparison.PID(seqi, seqj);
1085
1086         // use real sequence length rather than string length
1087         if (lngth[j] < lngth[i])
1088         {
1089           redundancy[j] = Math.max(pid, redundancy[j]);
1090         }
1091         else
1092         {
1093           redundancy[i] = Math.max(pid, redundancy[i]);
1094         }
1095
1096       }
1097     }
1098     return redundancy;
1099   }
1100 }