fabfbff6d15915951891ca72b8a91f5740752724
[jalview.git] / src / jalview / util / Comparison.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.util;
22
23 import jalview.datamodel.SequenceI;
24
25 import java.util.ArrayList;
26 import java.util.List;
27
28 /**
29  * Assorted methods for analysing or comparing sequences.
30  */
31 public class Comparison
32 {
33   private static final int EIGHTY_FIVE = 85;
34
35   private static final int TO_UPPER_CASE = 'a' - 'A';
36
37   public static final char GAP_SPACE = ' ';
38
39   public static final char GAP_DOT = '.';
40
41   public static final char GAP_DASH = '-';
42
43   public static final String GapChars = new String(
44           new char[]
45           { GAP_SPACE, GAP_DOT, GAP_DASH });
46
47   /**
48    * DOCUMENT ME!
49    * 
50    * @param ii
51    *          DOCUMENT ME!
52    * @param jj
53    *          DOCUMENT ME!
54    * 
55    * @return DOCUMENT ME!
56    */
57   public static final float compare(SequenceI ii, SequenceI jj)
58   {
59     return Comparison.compare(ii, jj, 0, ii.getLength() - 1);
60   }
61
62   /**
63    * this was supposed to be an ungapped pid calculation
64    * 
65    * @param ii
66    *          SequenceI
67    * @param jj
68    *          SequenceI
69    * @param start
70    *          int
71    * @param end
72    *          int
73    * @return float
74    */
75   public static float compare(SequenceI ii, SequenceI jj, int start,
76           int end)
77   {
78     String si = ii.getSequenceAsString();
79     String sj = jj.getSequenceAsString();
80
81     int ilen = si.length() - 1;
82     int jlen = sj.length() - 1;
83
84     while (Comparison.isGap(si.charAt(start + ilen)))
85     {
86       ilen--;
87     }
88
89     while (Comparison.isGap(sj.charAt(start + jlen)))
90     {
91       jlen--;
92     }
93
94     int count = 0;
95     int match = 0;
96     float pid = -1;
97
98     if (ilen > jlen)
99     {
100       for (int j = 0; j < jlen; j++)
101       {
102         if (si.substring(start + j, start + j + 1)
103                 .equals(sj.substring(start + j, start + j + 1)))
104         {
105           match++;
106         }
107
108         count++;
109       }
110
111       pid = (float) match / (float) ilen * 100;
112     }
113     else
114     {
115       for (int j = 0; j < jlen; j++)
116       {
117         if (si.substring(start + j, start + j + 1)
118                 .equals(sj.substring(start + j, start + j + 1)))
119         {
120           match++;
121         }
122
123         count++;
124       }
125
126       pid = (float) match / (float) jlen * 100;
127     }
128
129     return pid;
130   }
131
132   /**
133    * this is a gapped PID calculation
134    * 
135    * @param s1
136    *          SequenceI
137    * @param s2
138    *          SequenceI
139    * @return float
140    * @deprecated use PIDModel.computePID()
141    */
142   @Deprecated
143   public final static float PID(String seq1, String seq2)
144   {
145     return PID(seq1, seq2, 0, seq1.length());
146   }
147
148   static final int caseShift = 'a' - 'A';
149
150   // Another pid with region specification
151   /**
152    * @deprecated use PIDModel.computePID()
153    */
154   @Deprecated
155   public final static float PID(String seq1, String seq2, int start,
156           int end)
157   {
158     return PID(seq1, seq2, start, end, true, false);
159   }
160
161   /**
162    * Calculate percent identity for a pair of sequences over a particular range,
163    * with different options for ignoring gaps.
164    * 
165    * @param seq1
166    * @param seq2
167    * @param start
168    *          - position in seqs
169    * @param end
170    *          - position in seqs
171    * @param wcGaps
172    *          - if true - gaps match any character, if false, do not match
173    *          anything
174    * @param ungappedOnly
175    *          - if true - only count PID over ungapped columns
176    * @return
177    * @deprecated use PIDModel.computePID()
178    */
179   @Deprecated
180   public final static float PID(String seq1, String seq2, int start,
181           int end, boolean wcGaps, boolean ungappedOnly)
182   {
183     int s1len = seq1.length();
184     int s2len = seq2.length();
185
186     int len = Math.min(s1len, s2len);
187
188     if (end < len)
189     {
190       len = end;
191     }
192
193     if (len < start)
194     {
195       start = len - 1; // we just use a single residue for the difference
196     }
197
198     int elen = len - start, bad = 0;
199     char chr1;
200     char chr2;
201     boolean agap;
202     for (int i = start; i < len; i++)
203     {
204       chr1 = seq1.charAt(i);
205
206       chr2 = seq2.charAt(i);
207       agap = isGap(chr1) || isGap(chr2);
208       if ('a' <= chr1 && chr1 <= 'z')
209       {
210         // TO UPPERCASE !!!
211         // Faster than toUpperCase
212         chr1 -= caseShift;
213       }
214       if ('a' <= chr2 && chr2 <= 'z')
215       {
216         // TO UPPERCASE !!!
217         // Faster than toUpperCase
218         chr2 -= caseShift;
219       }
220
221       if (chr1 != chr2)
222       {
223         if (agap)
224         {
225           if (ungappedOnly)
226           {
227             elen--;
228           }
229           else if (!wcGaps)
230           {
231             bad++;
232           }
233         }
234         else
235         {
236           bad++;
237         }
238       }
239
240     }
241     if (elen < 1)
242     {
243       return 0f;
244     }
245     return ((float) 100 * (elen - bad)) / elen;
246   }
247
248   /**
249    * Answers true if the supplied character is a recognised gap character, else
250    * false. Currently hard-coded to recognise '-', '-' or ' ' (hyphen / dot /
251    * space).
252    * 
253    * @param c
254    * 
255    * @return
256    */
257   public static final boolean isGap(char c)
258   {
259     return (c == GAP_DASH || c == GAP_DOT || c == GAP_SPACE);
260   }
261
262   /**
263    * Overloaded method signature to test whether a single sequence is nucleotide
264    * (that is, more than 85% CGTAUNX)
265    * 
266    * @param seq
267    * @return
268    */
269   public static final boolean isNucleotide(SequenceI seq)
270   {
271     if (seq == null)
272     {
273       return false;
274     }
275     long ntCount = 0;
276     long aaCount = 0;
277     long nCount = 0;
278
279     int len = seq.getLength();
280     for (int i = 0; i < len; i++)
281     {
282       char c = seq.getCharAt(i);
283       if (isNucleotide(c) || isX(c))
284       {
285         ntCount++;
286       }
287       else if (!isGap(c))
288       {
289         aaCount++;
290         if (isN(c))
291         {
292           nCount++;
293         }
294       }
295     }
296     /*
297      * Check for nucleotide count > 85% of total count (in a form that evades
298      * int / float conversion or divide by zero).
299      */
300     if ((ntCount + nCount) * 100 > EIGHTY_FIVE * (ntCount + aaCount))
301     {
302       return ntCount > 0; // all N is considered protein. Could use a threshold
303                           // here too
304     }
305     else
306     {
307       return false;
308     }
309   }
310
311   /**
312    * Answers true if more than 85% of the sequence residues (ignoring gaps) are
313    * A, G, C, T or U, else false. This is just a heuristic guess and may give a
314    * wrong answer (as AGCT are also amino acid codes).
315    * 
316    * @param seqs
317    * @return
318    */
319   public static final boolean isNucleotide(SequenceI[] seqs)
320   {
321     if (seqs == null)
322     {
323       return false;
324     }
325     // true if we have seen a nucleotide sequence
326     boolean na = false;
327     for (SequenceI seq : seqs)
328     {
329       if (seq == null)
330       {
331         continue;
332       }
333       na = true;
334       // TODO could possibly make an informed guess just from the first sequence
335       // to save a lengthy calculation
336       if (seq.isProtein())
337       {
338         // if even one looks like protein, the alignment is protein
339         return false;
340       }
341     }
342     return na;
343   }
344
345   /**
346    * Answers true if the character is one of aAcCgGtTuU
347    * 
348    * @param c
349    * @return
350    */
351   public static boolean isNucleotide(char c)
352   {
353     if ('a' <= c && c <= 'z')
354     {
355       c -= TO_UPPER_CASE;
356     }
357     switch (c)
358     {
359     case 'A':
360     case 'C':
361     case 'G':
362     case 'T':
363     case 'U':
364       return true;
365     }
366     return false;
367   }
368
369   public static boolean isN(char c)
370   {
371     switch (c)
372     {
373     case 'N':
374     case 'n':
375       return true;
376     }
377     return false;
378   }
379
380   public static boolean isX(char c)
381   {
382     switch (c)
383     {
384     case 'X':
385     case 'x':
386       return true;
387     }
388     return false;
389   }
390
391   /**
392    * Answers true if every character in the string is one of aAcCgGtTuU, or
393    * (optionally) a gap character (dot, dash, space), else false
394    * 
395    * @param s
396    * @param allowGaps
397    * @return
398    */
399   public static boolean isNucleotideSequence(String s, boolean allowGaps)
400   {
401     if (s == null)
402     {
403       return false;
404     }
405     for (int i = 0; i < s.length(); i++)
406     {
407       char c = s.charAt(i);
408       if (!isNucleotide(c))
409       {
410         if (!allowGaps || !isGap(c))
411         {
412           return false;
413         }
414       }
415     }
416     return true;
417   }
418
419   /**
420    * Convenience overload of isNucleotide
421    * 
422    * @param seqs
423    * @return
424    */
425   public static boolean isNucleotide(SequenceI[][] seqs)
426   {
427     if (seqs == null)
428     {
429       return false;
430     }
431     List<SequenceI> flattened = new ArrayList<SequenceI>();
432     for (SequenceI[] ss : seqs)
433     {
434       for (SequenceI s : ss)
435       {
436         flattened.add(s);
437       }
438     }
439     final SequenceI[] oneDArray = flattened
440             .toArray(new SequenceI[flattened.size()]);
441     return isNucleotide(oneDArray);
442   }
443
444   /**
445    * Compares two residues either case sensitively or case insensitively
446    * depending on the caseSensitive flag
447    * 
448    * @param c1
449    *          first char
450    * @param c2
451    *          second char to compare with
452    * @param caseSensitive
453    *          if true comparison will be case sensitive otherwise its not
454    * @return
455    */
456   public static boolean isSameResidue(char c1, char c2,
457           boolean caseSensitive)
458   {
459     if (caseSensitive)
460     {
461       return (c1 == c2);
462     }
463     else
464     {
465       return Character.toUpperCase(c1) == Character.toUpperCase(c2);
466     }
467   }
468 }