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