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