src/jalview/analysis/Conservation.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.ResidueCount.SymbolCounts;
  24 import jalview.datamodel.AlignmentAnnotation;
  25 import jalview.datamodel.Annotation;
  26 import jalview.datamodel.Sequence;
  27 import jalview.datamodel.SequenceI;
  28 import jalview.schemes.ResidueProperties;
  29 import jalview.util.Comparison;
  30
  31 import java.awt.Color;
  32 import java.util.List;
  33 import java.util.Map;
  34 import java.util.Map.Entry;
  35 import java.util.TreeMap;
  36 import java.util.Vector;
  37
  38 /**
  39  * Calculates conservation values for a given set of sequences
  40  *
  41  * @author $author$
  42  * @version $Revision$
  43  */
  44 public class Conservation
  45 {
  46   private static final int TOUPPERCASE = 'a' - 'A';
  47
  48   SequenceI[] sequences;
  49
  50   int start;
  51
  52   int end;
  53
  54   Vector<int[]> seqNums; // vector of int vectors where first is sequence
  55                          // checksum
  56
  57   int maxLength = 0; // used by quality calcs
  58
  59   boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
  60
  61   /*
  62    * a map per column with {property, conservation} where conservation value is
  63    * 1 (property is conserved), 0 (property is negatively conserved) or -1
  64    * (property is not conserved i.e. column has residues with and without it)
  65    */
  66   Map<String, Integer>[] total;
  67
  68   boolean canonicaliseAa = true; // if true then conservation calculation will
  69
  70   // map all symbols to canonical aa numbering
  71   // rather than consider conservation of that
  72   // symbol
  73
  74   /** Stores calculated quality values */
  75   private Vector<Double> quality;
  76
  77   /** Stores maximum and minimum values of quality values */
  78   private double[] qualityRange = new double[2];
  79
  80   private Sequence consSequence;
  81
  82   /*
  83    * percentage of residues in a column to qualify for counting conservation
  84    */
  85   private int threshold;
  86
  87   private String name = "";
  88
  89   private int[][] cons2;
  90
  91   private String[] consSymbs;
  92
  93   /**
  94    * Creates a new Conservation object.
  95    *
  96    * @param name
  97    *          Name of conservation
  98    * @param threshold
  99    *          to count the residues in residueHash(). commonly used value is 3
 100    * @param sequences
 101    *          sequences to be used in calculation
 102    * @param start
 103    *          start residue position
 104    * @param end
 105    *          end residue position
 106    */
 107   public Conservation(String name, int threshold,
 108           List<SequenceI> sequences, int start, int end)
 109   {
 110     this.name = name;
 111     this.threshold = threshold;
 112     this.start = start;
 113     this.end = end;
 114
 115     maxLength = end - start + 1; // default width includes bounds of
 116     // calculation
 117
 118     int s, sSize = sequences.size();
 119     SequenceI[] sarray = new SequenceI[sSize];
 120     this.sequences = sarray;
 121     try
 122     {
 123       for (s = 0; s < sSize; s++)
 124       {
 125         sarray[s] = sequences.get(s);
 126         if (sarray[s].getLength() > maxLength)
 127         {
 128           maxLength = sarray[s].getLength();
 129         }
 130       }
 131     } catch (ArrayIndexOutOfBoundsException ex)
 132     {
 133       // bail - another thread has modified the sequence array, so the current
 134       // calculation is probably invalid.
 135       this.sequences = new SequenceI[0];
 136       maxLength = 0;
 137     }
 138   }
 139
 140   /**
 141    * Translate sequence i into a numerical representation and store it in the
 142    * i'th position of the seqNums array.
 143    *
 144    * @param i
 145    */
 146   private void calcSeqNum(int i)
 147   {
 148     String sq = null; // for dumb jbuilder not-inited exception warning
 149     int[] sqnum = null;
 150
 151     int sSize = sequences.length;
 152
 153     if ((i > -1) && (i < sSize))
 154     {
 155       sq = sequences[i].getSequenceAsString();
 156
 157       if (seqNums.size() <= i)
 158       {
 159         seqNums.addElement(new int[sq.length() + 1]);
 160       }
 161
 162       if (sq.hashCode() != seqNums.elementAt(i)[0])
 163       {
 164         int j;
 165         int len;
 166         seqNumsChanged = true;
 167         len = sq.length();
 168
 169         if (maxLength < len)
 170         {
 171           maxLength = len;
 172         }
 173
 174         sqnum = new int[len + 1]; // better to always make a new array -
 175         // sequence can change its length
 176         sqnum[0] = sq.hashCode();
 177
 178         for (j = 1; j <= len; j++)
 179         {
 180           sqnum[j] = jalview.schemes.ResidueProperties.aaIndex[sq
 181                   .charAt(j - 1)];
 182         }
 183
 184         seqNums.setElementAt(sqnum, i);
 185       }
 186       else
 187       {
 188         System.out.println("SEQUENCE HAS BEEN DELETED!!!");
 189       }
 190     }
 191     else
 192     {
 193       // JBPNote INFO level debug
 194       System.err
 195               .println("ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
 196     }
 197   }
 198
 199   /**
 200    * Calculates the conservation values for given set of sequences
 201    */
 202   public void calculate()
 203   {
 204     int height = sequences.length;
 205
 206     total = new Map[maxLength];
 207
 208     for (int column = start; column <= end; column++)
 209     {
 210       ResidueCount values = countResidues(column);
 211
 212       // TODO is threshold a percentage or count value?
 213       int thresh = (threshold * height) / 100;
 214
 215       /*
 216        * check observed residues in column and record whether each
 217        * physico-chemical property is conserved (+1), negatively conserved (0),
 218        * or not conserved (-1)
 219        * Using TreeMap means properties are displayed in alphabetical order
 220        */
 221       Map<String, Integer> resultHash = new TreeMap<String, Integer>();
 222       SymbolCounts symbolCounts = values.getSymbolCounts();
 223       char[] symbols = symbolCounts.symbols;
 224       int[] counts = symbolCounts.values;
 225       for (int j = 0; j < symbols.length; j++)
 226       {
 227         char c = symbols[j];
 228         if (counts[j] > thresh)
 229         {
 230           recordConservation(resultHash, String.valueOf(c));
 231         }
 232       }
 233       if (values.getGapCount() > thresh)
 234       {
 235         recordConservation(resultHash, "-");
 236       }
 237
 238       if (total.length > 0)
 239       {
 240         total[column - start] = resultHash;
 241       }
 242     }
 243   }
 244
 245   /**
 246    * Updates the conservation results for an observed residue
 247    *
 248    * @param resultMap
 249    *          a map of {property, conservation} where conservation value is +1
 250    *          (all residues have the property), 0 (no residue has the property)
 251    *          or -1 (some do, some don't)
 252    * @param res
 253    */
 254   protected static void recordConservation(Map<String, Integer> resultMap,
 255           String res)
 256   {
 257     res = res.toUpperCase();
 258     for (Entry<String, Map<String, Integer>> property : ResidueProperties.propHash
 259             .entrySet())
 260     {
 261       String propertyName = property.getKey();
 262       Integer residuePropertyValue = property.getValue().get(res);
 263
 264       if (!resultMap.containsKey(propertyName))
 265       {
 266         /*
 267          * first time we've seen this residue - note whether it has this property
 268          */
 269         if (residuePropertyValue != null)
 270         {
 271           resultMap.put(propertyName, residuePropertyValue);
 272         }
 273         else
 274         {
 275           /*
 276            * unrecognised residue - use default value for property
 277            */
 278           resultMap.put(propertyName, property.getValue().get("-"));
 279         }
 280       }
 281       else
 282       {
 283         Integer currentResult = resultMap.get(propertyName);
 284         if (currentResult.intValue() != -1
 285                 && !currentResult.equals(residuePropertyValue))
 286         {
 287           /*
 288            * property is unconserved - residues seen both with and without it
 289            */
 290           resultMap.put(propertyName, Integer.valueOf(-1));
 291         }
 292       }
 293     }
 294   }
 295
 296   /**
 297    * Counts residues (upper-cased) and gaps in the given column
 298    *
 299    * @param column
 300    * @return
 301    */
 302   protected ResidueCount countResidues(int column)
 303   {
 304     ResidueCount values = new ResidueCount(false);
 305
 306     for (int row = 0; row < sequences.length; row++)
 307     {
 308       if (sequences[row].getLength() > column)
 309       {
 310         char c = sequences[row].getCharAt(column);
 311         if (canonicaliseAa)
 312         {
 313           int index = ResidueProperties.aaIndex[c];
 314           c = index > 20 ? '-' : ResidueProperties.aa[index].charAt(0);
 315         }
 316         else
 317         {
 318           c = toUpperCase(c);
 319         }
 320         if (Comparison.isGap(c))
 321         {
 322           values.addGap();
 323         }
 324         else
 325         {
 326           values.add(c);
 327         }
 328       }
 329       else
 330       {
 331         values.addGap();
 332       }
 333     }
 334     return values;
 335   }
 336
 337   /**
 338    * Counts conservation and gaps for a column of the alignment
 339    *
 340    * @return { 1 if fully conserved, else 0, gap count }
 341    */
 342   public int[] countConservationAndGaps(int column)
 343   {
 344     int gapCount = 0;
 345     boolean fullyConserved = true;
 346     int iSize = sequences.length;
 347
 348     if (iSize == 0)
 349     {
 350       return new int[] { 0, 0 };
 351     }
 352
 353     char lastRes = '0';
 354     for (int i = 0; i < iSize; i++)
 355     {
 356       if (column >= sequences[i].getLength())
 357       {
 358         gapCount++;
 359         continue;
 360       }
 361
 362       char c = sequences[i].getCharAt(column); // gaps do not have upper/lower case
 363
 364       if (Comparison.isGap((c)))
 365       {
 366         gapCount++;
 367       }
 368       else
 369       {
 370         c = toUpperCase(c);
 371         if (lastRes == '0')
 372         {
 373           lastRes = c;
 374         }
 375         if (c != lastRes)
 376         {
 377           fullyConserved = false;
 378         }
 379       }
 380     }
 381
 382     int[] r = new int[] { fullyConserved ? 1 : 0, gapCount };
 383     return r;
 384   }
 385
 386   /**
 387    * Returns the upper-cased character if between 'a' and 'z', else the
 388    * unchanged value
 389    *
 390    * @param c
 391    * @return
 392    */
 393   char toUpperCase(char c)
 394   {
 395     if ('a' <= c && c <= 'z')
 396     {
 397       c -= TOUPPERCASE;
 398     }
 399     return c;
 400   }
 401
 402   /**
 403    * Calculates the conservation sequence
 404    *
 405    * @param positiveOnly
 406    *          if true, calculate positive conservation; else calculate both
 407    *          positive and negative conservation
 408    * @param maxPercentageGaps
 409    *          the percentage of gaps in a column, at or above which no
 410    *          conservation is asserted
 411    */
 412   public void verdict(boolean positiveOnly, float maxPercentageGaps)
 413   {
 414     // TODO call this at the end of calculate(), should not be a public method
 415
 416     StringBuilder consString = new StringBuilder(end);
 417
 418     // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY
 419     // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE
 420     // DOES NOT EXIST IN JALVIEW 2.1.2
 421     for (int i = 0; i < start; i++)
 422     {
 423       consString.append('-');
 424     }
 425     consSymbs = new String[end - start + 1];
 426     for (int i = start; i <= end; i++)
 427     {
 428       int[] gapcons = countConservationAndGaps(i);
 429       boolean fullyConserved = gapcons[0] == 1;
 430       int totGaps = gapcons[1];
 431       float pgaps = (totGaps * 100f) / sequences.length;
 432
 433       if (maxPercentageGaps > pgaps)
 434       {
 435         Map<String, Integer> resultHash = total[i - start];
 436         int count = 0;
 437         StringBuilder positives = new StringBuilder(64);
 438         StringBuilder negatives = new StringBuilder(32);
 439         for (String type : resultHash.keySet())
 440         {
 441           int result = resultHash.get(type).intValue();
 442           if (result == -1)
 443           {
 444             /*
 445              * not conserved either positively or negatively
 446              */
 447             continue;
 448           }
 449           count++;
 450           if (result == 1)
 451           {
 452             /*
 453              * positively conserved property (all residues have it)
 454              */
 455             positives.append(positives.length() == 0 ? "" : " ");
 456             positives.append(type);
 457           }
 458           if (result == 0 && !positiveOnly)
 459           {
 460             /*
 461              * negatively conserved property (all residues lack it)
 462              */
 463             negatives.append(negatives.length() == 0 ? "" : " ");
 464             negatives.append("!").append(type);
 465           }
 466         }
 467         if (negatives.length() > 0)
 468         {
 469           positives.append(" ").append(negatives);
 470         }
 471         consSymbs[i - start] = positives.toString();
 472
 473         if (count < 10)
 474         {
 475           consString.append(count); // Conserved props!=Identity
 476         }
 477         else
 478         {
 479           consString.append(fullyConserved ? "*" : "+");
 480         }
 481       }
 482       else
 483       {
 484         consString.append('-');
 485       }
 486     }
 487
 488     consSequence = new Sequence(name, consString.toString(), start, end);
 489   }
 490
 491   /**
 492    *
 493    *
 494    * @return Conservation sequence
 495    */
 496   public Sequence getConsSequence()
 497   {
 498     return consSequence;
 499   }
 500
 501   // From Alignment.java in jalview118
 502   public void findQuality()
 503   {
 504     findQuality(0, maxLength - 1);
 505   }
 506
 507   /**
 508    * DOCUMENT ME!
 509    */
 510   private void percentIdentity2()
 511   {
 512     seqNums = new Vector<int[]>();
 513     // calcSeqNum(s);
 514     int i = 0, iSize = sequences.length;
 515     // Do we need to calculate this again?
 516     for (i = 0; i < iSize; i++)
 517     {
 518       calcSeqNum(i);
 519     }
 520
 521     if ((cons2 == null) || seqNumsChanged)
 522     {
 523       cons2 = new int[maxLength][24];
 524
 525       // Initialize the array
 526       for (int j = 0; j < 24; j++)
 527       {
 528         for (i = 0; i < maxLength; i++)
 529         {
 530           cons2[i][j] = 0;
 531         }
 532       }
 533
 534       int[] sqnum;
 535       int j = 0;
 536
 537       while (j < sequences.length)
 538       {
 539         sqnum = seqNums.elementAt(j);
 540
 541         for (i = 1; i < sqnum.length; i++)
 542         {
 543           cons2[i - 1][sqnum[i]]++;
 544         }
 545
 546         for (i = sqnum.length - 1; i < maxLength; i++)
 547         {
 548           cons2[i][23]++; // gap count
 549         }
 550
 551         j++;
 552       }
 553
 554       // unnecessary ?
 555
 556       /*
 557        * for (int i=start; i <= end; i++) { int max = -1000; int maxi = -1; int
 558        * maxj = -1;
 559        *
 560        * for (int j=0;j<24;j++) { if (cons2[i][j] > max) { max = cons2[i][j];
 561        * maxi = i; maxj = j; } } }
 562        */
 563     }
 564   }
 565
 566   /**
 567    * Calculates the quality of the set of sequences
 568    *
 569    * @param startRes
 570    *          Start residue
 571    * @param endRes
 572    *          End residue
 573    */
 574   public void findQuality(int startRes, int endRes)
 575   {
 576     quality = new Vector<Double>();
 577
 578     double max = -10000;
 579     int[][] BLOSUM62 = ResidueProperties.getBLOSUM62();
 580
 581     // Loop over columns // JBPNote Profiling info
 582     // long ts = System.currentTimeMillis();
 583     // long te = System.currentTimeMillis();
 584     percentIdentity2();
 585
 586     int size = seqNums.size();
 587     int[] lengths = new int[size];
 588     double tot, bigtot, sr, tmp;
 589     double[] x, xx;
 590     int l, j, i, ii, i2, k, seqNum;
 591
 592     for (l = 0; l < size; l++)
 593     {
 594       lengths[l] = seqNums.elementAt(l).length - 1;
 595     }
 596
 597     for (j = startRes; j <= endRes; j++)
 598     {
 599       bigtot = 0;
 600
 601       // First Xr = depends on column only
 602       x = new double[24];
 603
 604       for (ii = 0; ii < 24; ii++)
 605       {
 606         x[ii] = 0;
 607
 608         for (i2 = 0; i2 < 24; i2++)
 609         {
 610           x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4);
 611         }
 612
 613         x[ii] /= size;
 614       }
 615
 616       // Now calculate D for each position and sum
 617       for (k = 0; k < size; k++)
 618       {
 619         tot = 0;
 620         xx = new double[24];
 621         seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : 23; // Sequence,
 622                                                                       // or gap
 623                                                                       // at the
 624                                                                       // end
 625
 626         // This is a loop over r
 627         for (i = 0; i < 23; i++)
 628         {
 629           sr = 0;
 630
 631           sr = (double) BLOSUM62[i][seqNum] + 4;
 632
 633           // Calculate X with another loop over residues
 634           // System.out.println("Xi " + i + " " + x[i] + " " + sr);
 635           xx[i] = x[i] - sr;
 636
 637           tot += (xx[i] * xx[i]);
 638         }
 639
 640         bigtot += Math.sqrt(tot);
 641       }
 642
 643       // This is the quality for one column
 644       if (max < bigtot)
 645       {
 646         max = bigtot;
 647       }
 648
 649       // bigtot = bigtot * (size-cons2[j][23])/size;
 650       quality.addElement(new Double(bigtot));
 651
 652       // Need to normalize by gaps
 653     }
 654
 655     double newmax = -10000;
 656
 657     for (j = startRes; j <= endRes; j++)
 658     {
 659       tmp = quality.elementAt(j).doubleValue();
 660       tmp = ((max - tmp) * (size - cons2[j][23])) / size;
 661
 662       // System.out.println(tmp+ " " + j);
 663       quality.setElementAt(new Double(tmp), j);
 664
 665       if (tmp > newmax)
 666       {
 667         newmax = tmp;
 668       }
 669     }
 670
 671     // System.out.println("Quality " + s);
 672     qualityRange[0] = 0D;
 673     qualityRange[1] = newmax;
 674   }
 675
 676   /**
 677    * Complete the given consensus and quuality annotation rows. Note: currently
 678    * this method will enlarge the given annotation row if it is too small,
 679    * otherwise will leave its length unchanged.
 680    *
 681    * @param conservation
 682    *          conservation annotation row
 683    * @param quality2
 684    *          (optional - may be null)
 685    * @param istart
 686    *          first column for conservation
 687    * @param alWidth
 688    *          extent of conservation
 689    */
 690   public void completeAnnotations(AlignmentAnnotation conservation,
 691           AlignmentAnnotation quality2, int istart, int alWidth)
 692   {
 693     char[] sequence = getConsSequence().getSequence();
 694     float minR;
 695     float minG;
 696     float minB;
 697     float maxR;
 698     float maxG;
 699     float maxB;
 700     minR = 0.3f;
 701     minG = 0.0f;
 702     minB = 0f;
 703     maxR = 1.0f - minR;
 704     maxG = 0.9f - minG;
 705     maxB = 0f - minB; // scalable range for colouring both Conservation and
 706     // Quality
 707
 708     float min = 0f;
 709     float max = 11f;
 710     float qmin = 0f;
 711     float qmax = 0f;
 712
 713     char c;
 714
 715     if (conservation != null && conservation.annotations != null
 716             && conservation.annotations.length < alWidth)
 717     {
 718       conservation.annotations = new Annotation[alWidth];
 719     }
 720
 721     if (quality2 != null)
 722     {
 723       quality2.graphMax = (float) qualityRange[1];
 724       if (quality2.annotations != null
 725               && quality2.annotations.length < alWidth)
 726       {
 727         quality2.annotations = new Annotation[alWidth];
 728       }
 729       qmin = (float) qualityRange[0];
 730       qmax = (float) qualityRange[1];
 731     }
 732
 733     for (int i = istart; i < alWidth; i++)
 734     {
 735       float value = 0;
 736
 737       c = sequence[i];
 738
 739       if (Character.isDigit(c))
 740       {
 741         value = c - '0';
 742       }
 743       else if (c == '*')
 744       {
 745         value = 11;
 746       }
 747       else if (c == '+')
 748       {
 749         value = 10;
 750       }
 751
 752       if (conservation != null)
 753       {
 754         float vprop = value - min;
 755         vprop /= max;
 756         int consp = i - start;
 757         String conssym = (value > 0 && consp > -1 && consp < consSymbs.length) ? consSymbs[consp]
 758                 : "";
 759         conservation.annotations[i] = new Annotation(String.valueOf(c),
 760                 conssym, ' ', value, new Color(minR + (maxR * vprop), minG
 761                         + (maxG * vprop), minB + (maxB * vprop)));
 762       }
 763
 764       // Quality calc
 765       if (quality2 != null)
 766       {
 767         value = quality.elementAt(i).floatValue();
 768         float vprop = value - qmin;
 769         vprop /= qmax;
 770         quality2.annotations[i] = new Annotation(" ",
 771                 String.valueOf(value), ' ', value, new Color(minR
 772                         + (maxR * vprop), minG + (maxG * vprop), minB
 773                         + (maxB * vprop)));
 774       }
 775     }
 776   }
 777
 778   /**
 779    * construct and call the calculation methods on a new Conservation object
 780    *
 781    * @param name
 782    *          - name of conservation
 783    * @param threshold
 784    *          - minimum number of conserved residues needed to indicate
 785    *          conservation (typically 3)
 786    * @param seqs
 787    * @param start
 788    *          first column in calculation window
 789    * @param end
 790    *          last column in calculation window
 791    * @param positiveOnly
 792    *          calculate positive (true) or positive and negative (false)
 793    *          conservation
 794    * @param maxPercentGaps
 795    *          percentage of gaps tolerated in column
 796    * @param calcQuality
 797    *          flag indicating if alignment quality should be calculated
 798    * @return Conservation object ready for use in visualization
 799    */
 800   public static Conservation calculateConservation(String name,
 801           int threshold, List<SequenceI> seqs, int start, int end,
 802           boolean positiveOnly, int maxPercentGaps, boolean calcQuality)
 803   {
 804     Conservation cons = new Conservation(name, threshold, seqs, start, end);
 805     cons.calculate();
 806     cons.verdict(positiveOnly, maxPercentGaps);
 807
 808     if (calcQuality)
 809     {
 810       cons.findQuality();
 811     }
 812
 813     return cons;
 814   }
 815
 816   /**
 817    * Returns the computed tooltip (annotation description) for a given column.
 818    * The tip is empty if the conservation score is zero, otherwise holds the
 819    * positively (and, optionally, negatively) conserved properties.
 820    *
 821    * @param column
 822    * @return
 823    */
 824   String getTooltip(int column)
 825   {
 826     char[] sequence = getConsSequence().getSequence();
 827     char val = column < sequence.length ? sequence[column] : '-';
 828     boolean hasConservation = val != '-' && val != '0';
 829     int consp = column - start;
 830     String tip = (hasConservation && consp > -1 && consp < consSymbs.length) ? consSymbs[consp]
 831             : "";
 832     return tip;
 833   }
 834 }