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