src/jalview/analysis/Conservation.java

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