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.List;
  32 import java.util.Map;
  33 import java.util.TreeMap;
  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 jSize = sequences.length;
 193     // int[] values; // Replaces residueHash
 194     SparseIntArray values = new SparseIntArray();
 195
 196     total = new Map[maxLength];
 197
 198     for (int i = start; i <= end; i++)
 199     {
 200       // values = new int[255];
 201       values.clear();
 202
 203       for (int j = 0; j < jSize; j++)
 204       {
 205         if (sequences[j].getLength() > i)
 206         {
 207           char 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.add(c, 1);
 235         }
 236         else
 237         {
 238           // values['-']++;
 239           values.add('-', 1);
 240         }
 241       }
 242
 243       // What is the count threshold to count the residues in residueHash()
 244       int thresh = (threshold * jSize) / 100;
 245
 246       // loop over all the found residues
 247       // Hashtable<String, Integer> resultHash = new Hashtable<String,
 248       // Integer>();
 249       Map<String, Integer> resultHash = new TreeMap<String, Integer>();
 250       // for (char v = '-'; v < 'Z'; v++)
 251       for (int key = 0; key < values.size(); key++)
 252       {
 253         char v = (char) values.keyAt(key);
 254         // if (values[v] > thresh)
 255         if (values.valueAt(key) > thresh)
 256         {
 257           String res = String.valueOf(v);
 258
 259           // Now loop over the properties
 260           for (String type : ResidueProperties.propHash.keySet())
 261           {
 262             Map<String, Integer> ht = ResidueProperties.propHash.get(type);
 263
 264             // Have we ticked this before?
 265             if (!resultHash.containsKey(type))
 266             {
 267               if (ht.containsKey(res))
 268               {
 269                 resultHash.put(type, ht.get(res));
 270               }
 271               else
 272               {
 273                 resultHash.put(type, ht.get("-"));
 274               }
 275             }
 276             else if (!resultHash.get(type).equals(ht.get(res)))
 277             {
 278               resultHash.put(type, new Integer(-1));
 279             }
 280           }
 281         }
 282       }
 283
 284       if (total.length > 0)
 285       {
 286         total[i - start] = resultHash;
 287       }
 288     }
 289   }
 290
 291   /*****************************************************************************
 292    * count conservation for the j'th column of the alignment
 293    *
 294    * @return { gap count, conserved residue count}
 295    */
 296   public int[] countConsNGaps(int j)
 297   {
 298     int count = 0;
 299     int cons = 0;
 300     int nres = 0;
 301     int[] r = new int[2];
 302     char f = '$';
 303     int i, iSize = sequences.length;
 304     char c;
 305
 306     for (i = 0; i < iSize; i++)
 307     {
 308       if (j >= sequences[i].getLength())
 309       {
 310         count++;
 311         continue;
 312       }
 313
 314       c = sequences[i].getCharAt(j); // gaps do not have upper/lower case
 315
 316       if (jalview.util.Comparison.isGap((c)))
 317       {
 318         count++;
 319       }
 320       else
 321       {
 322         c = toUpperCase(c);
 323         nres++;
 324
 325         if (nres == 1)
 326         {
 327           f = c;
 328           cons++;
 329         }
 330         else if (f == c)
 331         {
 332           cons++;
 333         }
 334       }
 335     }
 336
 337     r[0] = (nres == cons) ? 1 : 0;
 338     r[1] = count;
 339
 340     return r;
 341   }
 342
 343   /**
 344    * Returns the upper-cased character if between 'a' and 'z', else the
 345    * unchanged value
 346    *
 347    * @param c
 348    * @return
 349    */
 350   char toUpperCase(char c)
 351   {
 352     if ('a' <= c && c <= 'z')
 353     {
 354       c -= (32); // 32 = 'a' - 'A'
 355     }
 356     return c;
 357   }
 358
 359   /**
 360    * Calculates the conservation sequence
 361    *
 362    * @param consflag
 363    *          if true, positive conservation; false calculates negative
 364    *          conservation
 365    * @param percentageGaps
 366    *          commonly used value is 25
 367    */
 368   public void verdict(boolean consflag, float percentageGaps)
 369   {
 370     StringBuilder consString = new StringBuilder(end);
 371
 372     // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY
 373     // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE
 374     // DOES NOT EXIST IN JALVIEW 2.1.2
 375     for (int i = 0; i < start; i++)
 376     {
 377       consString.append('-');
 378     }
 379     consSymbs = new String[end - start + 1];
 380     for (int i = start; i <= end; i++)
 381     {
 382       int[] gapcons = countConsNGaps(i);
 383       int totGaps = gapcons[1];
 384       float pgaps = ((float) totGaps * 100) / sequences.length;
 385       StringBuilder positives = new StringBuilder(64);
 386       StringBuilder negatives = new StringBuilder(32);
 387       // consSymbs[i - start] = "";
 388
 389       if (percentageGaps > pgaps)
 390       {
 391         Map<String, Integer> resultHash = total[i - start];
 392         // Now find the verdict
 393         int count = 0;
 394         for (String type : resultHash.keySet())
 395         {
 396           int result = resultHash.get(type).intValue();
 397           // Do we want to count +ve conservation or +ve and -ve cons.?
 398           if (consflag)
 399           {
 400             if (result == 1)
 401             {
 402               // consSymbs[i - start] = type + " " + consSymbs[i - start];
 403               positives.append(positives.length() == 0 ? "" : " ");
 404               positives.append(type);
 405               count++;
 406             }
 407           }
 408           else
 409           {
 410             if (result != -1)
 411             {
 412               if (result == 0)
 413               {
 414                 /*
 415                  * add negatively conserved properties on the end
 416                  */
 417                 // consSymbs[i - start] = consSymbs[i - start] + " !" + type;
 418                 negatives.append(negatives.length() == 0 ? "" : " ");
 419                 negatives.append("!").append(type);
 420               }
 421               else
 422               {
 423                 /*
 424                  * put positively conserved properties on the front
 425                  */
 426                 // consSymbs[i - start] = type + " " + consSymbs[i - start];
 427                 positives.append(positives.length() == 0 ? "" : " ");
 428                 positives.append(type);
 429               }
 430               count++;
 431             }
 432           }
 433         }
 434         if (negatives.length() > 0)
 435         {
 436           positives.append(" ").append(negatives);
 437         }
 438         consSymbs[i - start] = positives.toString();
 439
 440         if (count < 10)
 441         {
 442           consString.append(count); // Conserved props!=Identity
 443         }
 444         else
 445         {
 446           consString.append((gapcons[0] == 1) ? "*" : "+");
 447         }
 448       }
 449       else
 450       {
 451         consString.append('-');
 452       }
 453     }
 454
 455     consSequence = new Sequence(name, consString.toString(), start, end);
 456   }
 457
 458   /**
 459    *
 460    *
 461    * @return Conservation sequence
 462    */
 463   public Sequence getConsSequence()
 464   {
 465     return consSequence;
 466   }
 467
 468   // From Alignment.java in jalview118
 469   public void findQuality()
 470   {
 471     findQuality(0, maxLength - 1);
 472   }
 473
 474   /**
 475    * DOCUMENT ME!
 476    */
 477   private void percentIdentity2()
 478   {
 479     seqNums = new Vector<int[]>();
 480     // calcSeqNum(s);
 481     int i = 0, iSize = sequences.length;
 482     // Do we need to calculate this again?
 483     for (i = 0; i < iSize; i++)
 484     {
 485       calcSeqNum(i);
 486     }
 487
 488     if ((cons2 == null) || seqNumsChanged)
 489     {
 490       cons2 = new int[maxLength][24];
 491
 492       // Initialize the array
 493       for (int j = 0; j < 24; j++)
 494       {
 495         for (i = 0; i < maxLength; i++)
 496         {
 497           cons2[i][j] = 0;
 498         }
 499       }
 500
 501       int[] sqnum;
 502       int j = 0;
 503
 504       while (j < sequences.length)
 505       {
 506         sqnum = seqNums.elementAt(j);
 507
 508         for (i = 1; i < sqnum.length; i++)
 509         {
 510           cons2[i - 1][sqnum[i]]++;
 511         }
 512
 513         for (i = sqnum.length - 1; i < maxLength; i++)
 514         {
 515           cons2[i][23]++; // gap count
 516         }
 517
 518         j++;
 519       }
 520
 521       // unnecessary ?
 522
 523       /*
 524        * for (int i=start; i <= end; i++) { int max = -1000; int maxi = -1; int
 525        * maxj = -1;
 526        *
 527        * for (int j=0;j<24;j++) { if (cons2[i][j] > max) { max = cons2[i][j];
 528        * maxi = i; maxj = j; } } }
 529        */
 530     }
 531   }
 532
 533   /**
 534    * Calculates the quality of the set of sequences
 535    *
 536    * @param startRes
 537    *          Start residue
 538    * @param endRes
 539    *          End residue
 540    */
 541   public void findQuality(int startRes, int endRes)
 542   {
 543     quality = new Vector<Double>();
 544
 545     double max = -10000;
 546     int[][] BLOSUM62 = ResidueProperties.getBLOSUM62();
 547
 548     // Loop over columns // JBPNote Profiling info
 549     // long ts = System.currentTimeMillis();
 550     // long te = System.currentTimeMillis();
 551     percentIdentity2();
 552
 553     int size = seqNums.size();
 554     int[] lengths = new int[size];
 555     double tot, bigtot, sr, tmp;
 556     double[] x, xx;
 557     int l, j, i, ii, i2, k, seqNum;
 558
 559     for (l = 0; l < size; l++)
 560     {
 561       lengths[l] = seqNums.elementAt(l).length - 1;
 562     }
 563
 564     for (j = startRes; j <= endRes; j++)
 565     {
 566       bigtot = 0;
 567
 568       // First Xr = depends on column only
 569       x = new double[24];
 570
 571       for (ii = 0; ii < 24; ii++)
 572       {
 573         x[ii] = 0;
 574
 575         for (i2 = 0; i2 < 24; i2++)
 576         {
 577           x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4);
 578         }
 579
 580         x[ii] /= size;
 581       }
 582
 583       // Now calculate D for each position and sum
 584       for (k = 0; k < size; k++)
 585       {
 586         tot = 0;
 587         xx = new double[24];
 588         seqNum = (j < lengths[k]) ? seqNums.elementAt(k)[j + 1] : 23; // Sequence,
 589                                                                       // or gap
 590                                                                       // at the
 591                                                                       // end
 592
 593         // This is a loop over r
 594         for (i = 0; i < 23; i++)
 595         {
 596           sr = 0;
 597
 598           sr = (double) BLOSUM62[i][seqNum] + 4;
 599
 600           // Calculate X with another loop over residues
 601           // System.out.println("Xi " + i + " " + x[i] + " " + sr);
 602           xx[i] = x[i] - sr;
 603
 604           tot += (xx[i] * xx[i]);
 605         }
 606
 607         bigtot += Math.sqrt(tot);
 608       }
 609
 610       // This is the quality for one column
 611       if (max < bigtot)
 612       {
 613         max = bigtot;
 614       }
 615
 616       // bigtot = bigtot * (size-cons2[j][23])/size;
 617       quality.addElement(new Double(bigtot));
 618
 619       // Need to normalize by gaps
 620     }
 621
 622     double newmax = -10000;
 623
 624     for (j = startRes; j <= endRes; j++)
 625     {
 626       tmp = quality.elementAt(j).doubleValue();
 627       tmp = ((max - tmp) * (size - cons2[j][23])) / size;
 628
 629       // System.out.println(tmp+ " " + j);
 630       quality.setElementAt(new Double(tmp), j);
 631
 632       if (tmp > newmax)
 633       {
 634         newmax = tmp;
 635       }
 636     }
 637
 638     // System.out.println("Quality " + s);
 639     qualityRange[0] = 0D;
 640     qualityRange[1] = newmax;
 641   }
 642
 643   /**
 644    * Complete the given consensus and quuality annotation rows. Note: currently
 645    * this method will enlarge the given annotation row if it is too small,
 646    * otherwise will leave its length unchanged.
 647    *
 648    * @param conservation
 649    *          conservation annotation row
 650    * @param quality2
 651    *          (optional - may be null)
 652    * @param istart
 653    *          first column for conservation
 654    * @param alWidth
 655    *          extent of conservation
 656    */
 657   public void completeAnnotations(AlignmentAnnotation conservation,
 658           AlignmentAnnotation quality2, int istart, int alWidth)
 659   {
 660     char[] sequence = getConsSequence().getSequence();
 661     float minR;
 662     float minG;
 663     float minB;
 664     float maxR;
 665     float maxG;
 666     float maxB;
 667     minR = 0.3f;
 668     minG = 0.0f;
 669     minB = 0f;
 670     maxR = 1.0f - minR;
 671     maxG = 0.9f - minG;
 672     maxB = 0f - minB; // scalable range for colouring both Conservation and
 673     // Quality
 674
 675     float min = 0f;
 676     float max = 11f;
 677     float qmin = 0f;
 678     float qmax = 0f;
 679
 680     char c;
 681
 682     if (conservation != null && conservation.annotations != null
 683             && conservation.annotations.length < alWidth)
 684     {
 685       conservation.annotations = new Annotation[alWidth];
 686     }
 687
 688     if (quality2 != null)
 689     {
 690       quality2.graphMax = (float) qualityRange[1];
 691       if (quality2.annotations != null
 692               && quality2.annotations.length < alWidth)
 693       {
 694         quality2.annotations = new Annotation[alWidth];
 695       }
 696       qmin = (float) qualityRange[0];
 697       qmax = (float) qualityRange[1];
 698     }
 699
 700     for (int i = istart; i < alWidth; i++)
 701     {
 702       float value = 0;
 703
 704       c = sequence[i];
 705
 706       if (Character.isDigit(c))
 707       {
 708         value = c - '0';
 709       }
 710       else if (c == '*')
 711       {
 712         value = 11;
 713       }
 714       else if (c == '+')
 715       {
 716         value = 10;
 717       }
 718
 719       if (conservation != null)
 720       {
 721         float vprop = value - min;
 722         vprop /= max;
 723         int consp = i - start;
 724         String conssym = (value > 0 && consp > -1 && consp < consSymbs.length) ? consSymbs[consp]
 725                 : "";
 726         conservation.annotations[i] = new Annotation(String.valueOf(c),
 727                 conssym, ' ', value, new Color(minR + (maxR * vprop), minG
 728                         + (maxG * vprop), minB + (maxB * vprop)));
 729       }
 730
 731       // Quality calc
 732       if (quality2 != null)
 733       {
 734         value = quality.elementAt(i).floatValue();
 735         float vprop = value - qmin;
 736         vprop /= qmax;
 737         quality2.annotations[i] = new Annotation(" ",
 738                 String.valueOf(value), ' ', value, new Color(minR
 739                         + (maxR * vprop), minG + (maxG * vprop), minB
 740                         + (maxB * vprop)));
 741       }
 742     }
 743   }
 744
 745   /**
 746    * construct and call the calculation methods on a new Conservation object
 747    *
 748    * @param name
 749    *          - name of conservation
 750    * @param threshold
 751    *          - minimum number of conserved residues needed to indicate
 752    *          conservation (typically 3)
 753    * @param seqs
 754    * @param start
 755    *          first column in calculation window
 756    * @param end
 757    *          last column in calculation window
 758    * @param posOrNeg
 759    *          positive (true) or negative (false) conservation
 760    * @param consPercGaps
 761    *          percentage of gaps tolerated in column
 762    * @param calcQuality
 763    *          flag indicating if alignment quality should be calculated
 764    * @return Conservation object ready for use in visualization
 765    */
 766   public static Conservation calculateConservation(String name,
 767           int threshold, List<SequenceI> seqs, int start, int end,
 768           boolean posOrNeg, int consPercGaps, boolean calcQuality)
 769   {
 770     Conservation cons = new Conservation(name, threshold, seqs, start, end);
 771     cons.calculate();
 772     cons.verdict(posOrNeg, consPercGaps);
 773
 774     if (calcQuality)
 775     {
 776       cons.findQuality();
 777     }
 778
 779     return cons;
 780   }
 781 }