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