src/jalview/analysis/Conservation.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.4)
   3  * Copyright (C) 2008 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle
   4  *
   5  * This program is free software; you can redistribute it and/or
   6  * modify it under the terms of the GNU General Public License
   7  * as published by the Free Software Foundation; either version 2
   8  * of the License, or (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  18  */
  19 package jalview.analysis;
  20
  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
  80    *                is 3
  81    * @param sequences
  82    *                sequences to be used in calculation
  83    * @param start
  84    *                start residue position
  85    * @param end
  86    *                end residue position
  87    */
  88   public Conservation(String name, Hashtable propHash, int threshold,
  89           Vector sequences, int start, int end)
  90   {
  91
  92     this.name = name;
  93     this.propHash = propHash;
  94     this.threshold = threshold;
  95     this.start = start;
  96     this.end = end;
  97
  98     maxLength = end - start + 1; // default width includes bounds of
  99     // calculation
 100
 101     int s, sSize = sequences.size();
 102     SequenceI[] sarray = new SequenceI[sSize];
 103     this.sequences = sarray;
 104
 105     for (s = 0; s < sSize; s++)
 106     {
 107       sarray[s] = (SequenceI) sequences.elementAt(s);
 108       if (sarray[s].getLength() > maxLength)
 109       {
 110         maxLength = sarray[s].getLength();
 111       }
 112     }
 113   }
 114
 115   /**
 116    * Translate sequence i into a numerical representation and store it in the
 117    * i'th position of the seqNums array.
 118    *
 119    * @param i
 120    */
 121   private void calcSeqNum(int i)
 122   {
 123     String sq = null; // for dumb jbuilder not-inited exception warning
 124     int[] sqnum = null;
 125
 126     int sSize = sequences.length;
 127
 128     if ((i > -1) && (i < sSize))
 129     {
 130       sq = sequences[i].getSequenceAsString();
 131
 132       if (seqNums.size() <= i)
 133       {
 134         seqNums.addElement(new int[sq.length() + 1]);
 135       }
 136
 137       if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0])
 138       {
 139         int j;
 140         int len;
 141         seqNumsChanged = true;
 142         len = sq.length();
 143
 144         if (maxLength < len)
 145         {
 146           maxLength = len;
 147         }
 148
 149         sqnum = new int[len + 1]; // better to always make a new array -
 150         // sequence can change its length
 151         sqnum[0] = sq.hashCode();
 152
 153         for (j = 1; j <= len; j++)
 154         {
 155           sqnum[j] = jalview.schemes.ResidueProperties.aaIndex[sq
 156                   .charAt(j - 1)];
 157         }
 158
 159         seqNums.setElementAt(sqnum, i);
 160       }
 161       else
 162       {
 163         System.out.println("SEQUENCE HAS BEEN DELETED!!!");
 164       }
 165     }
 166     else
 167     {
 168       // JBPNote INFO level debug
 169       System.err
 170               .println("ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
 171     }
 172   }
 173
 174   /**
 175    * Calculates the conservation values for given set of sequences
 176    */
 177   public void calculate()
 178   {
 179     Hashtable resultHash, ht;
 180     int thresh, j, jSize = sequences.length;
 181     int[] values; // Replaces residueHash
 182     String type, res = null;
 183     char c;
 184     Enumeration enumeration2;
 185
 186     total = new Hashtable[maxLength];
 187
 188     for (int i = start; i <= end; i++)
 189     {
 190       values = new int[255];
 191
 192       for (j = 0; j < jSize; j++)
 193       {
 194         if (sequences[j].getLength() > i)
 195         {
 196           c = sequences[j].getCharAt(i);
 197
 198           if (canonicaliseAa)
 199           { // lookup the base aa code symbol
 200             c = (char) jalview.schemes.ResidueProperties.aaIndex[sequences[j]
 201                     .getCharAt(i)];
 202             if (c > 20)
 203             {
 204               c = '-';
 205             }
 206             else
 207             {
 208               // recover canonical aa symbol
 209               c = jalview.schemes.ResidueProperties.aa[c].charAt(0);
 210             }
 211           }
 212           else
 213           {
 214             // original behaviour - operate on ascii symbols directly
 215             // No need to check if its a '-'
 216             if (c == '.' || c == ' ')
 217             {
 218               c = '-';
 219             }
 220
 221             if (!canonicaliseAa && 'a' <= c && c <= 'z')
 222             {
 223               c -= (32); // 32 = 'a' - 'A'
 224             }
 225           }
 226           values[c]++;
 227         }
 228         else
 229         {
 230           values['-']++;
 231         }
 232       }
 233
 234       // What is the count threshold to count the residues in residueHash()
 235       thresh = (threshold * (jSize)) / 100;
 236
 237       // loop over all the found residues
 238       resultHash = new Hashtable();
 239       for (char v = '-'; v < 'Z'; v++)
 240       {
 241
 242         if (values[v] > thresh)
 243         {
 244           res = String.valueOf(v);
 245
 246           // Now loop over the properties
 247           enumeration2 = propHash.keys();
 248
 249           while (enumeration2.hasMoreElements())
 250           {
 251             type = (String) enumeration2.nextElement();
 252             ht = (Hashtable) propHash.get(type);
 253
 254             // Have we ticked this before?
 255             if (!resultHash.containsKey(type))
 256             {
 257               if (ht.containsKey(res))
 258               {
 259                 resultHash.put(type, ht.get(res));
 260               }
 261               else
 262               {
 263                 resultHash.put(type, ht.get("-"));
 264               }
 265             }
 266             else if (((Integer) resultHash.get(type)).equals((Integer) ht
 267                     .get(res)) == false)
 268             {
 269               resultHash.put(type, new Integer(-1));
 270             }
 271           }
 272         }
 273       }
 274
 275       total[i - start] = resultHash;
 276     }
 277   }
 278
 279   /*****************************************************************************
 280    * count conservation for the j'th column of the alignment
 281    *
 282    * @return { gap count, conserved residue count}
 283    */
 284   public int[] countConsNGaps(int j)
 285   {
 286     int count = 0;
 287     int cons = 0;
 288     int nres = 0;
 289     int[] r = new int[2];
 290     char f = '$';
 291     int i, iSize = sequences.length;
 292     char c;
 293
 294     for (i = 0; i < iSize; i++)
 295     {
 296       if (j >= sequences[i].getLength())
 297       {
 298         count++;
 299         continue;
 300       }
 301
 302       c = sequences[i].getCharAt(j); // gaps do not have upper/lower case
 303
 304       if (jalview.util.Comparison.isGap((c)))
 305       {
 306         count++;
 307       }
 308       else
 309       {
 310         nres++;
 311
 312         if (nres == 1)
 313         {
 314           f = c;
 315           cons++;
 316         }
 317         else if (f == c)
 318         {
 319           cons++;
 320         }
 321       }
 322     }
 323
 324     r[0] = (nres == cons) ? 1 : 0;
 325     r[1] = count;
 326
 327     return r;
 328   }
 329
 330   /**
 331    * Calculates the conservation sequence
 332    *
 333    * @param consflag
 334    *                if true, poitiveve conservation; false calculates negative
 335    *                conservation
 336    * @param percentageGaps
 337    *                commonly used value is 25
 338    */
 339   public void verdict(boolean consflag, float percentageGaps)
 340   {
 341     StringBuffer consString = new StringBuffer();
 342     String type;
 343     Integer result;
 344     int[] gapcons;
 345     int totGaps, count;
 346     float pgaps;
 347     Hashtable resultHash;
 348     Enumeration enumeration;
 349
 350     // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY
 351     // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE
 352     // DOES NOT EXIST IN JALVIEW 2.1.2
 353     for (int i = 0; i < start; i++)
 354     {
 355       consString.append('-');
 356     }
 357
 358     for (int i = start; i <= end; i++)
 359     {
 360       gapcons = countConsNGaps(i);
 361       totGaps = gapcons[1];
 362       pgaps = ((float) totGaps * 100) / (float) sequences.length;
 363
 364       if (percentageGaps > pgaps)
 365       {
 366         resultHash = total[i - start];
 367
 368         // Now find the verdict
 369         count = 0;
 370         enumeration = resultHash.keys();
 371
 372         while (enumeration.hasMoreElements())
 373         {
 374           type = (String) enumeration.nextElement();
 375           result = (Integer) resultHash.get(type);
 376
 377           // Do we want to count +ve conservation or +ve and -ve cons.?
 378           if (consflag)
 379           {
 380             if (result.intValue() == 1)
 381             {
 382               count++;
 383             }
 384           }
 385           else
 386           {
 387             if (result.intValue() != -1)
 388             {
 389               count++;
 390             }
 391           }
 392         }
 393
 394         if (count < 10)
 395         {
 396           consString.append(count); // Conserved props!=Identity
 397         }
 398         else
 399         {
 400           consString.append((gapcons[0] == 1) ? "*" : "+");
 401         }
 402       }
 403       else
 404       {
 405         consString.append('-');
 406       }
 407     }
 408
 409     consSequence = new Sequence(name, consString.toString(), start, end);
 410   }
 411
 412   /**
 413    *
 414    *
 415    * @return Conservation sequence
 416    */
 417   public Sequence getConsSequence()
 418   {
 419     return consSequence;
 420   }
 421
 422   // From Alignment.java in jalview118
 423   public void findQuality()
 424   {
 425     findQuality(0, maxLength - 1);
 426   }
 427
 428   /**
 429    * DOCUMENT ME!
 430    */
 431   private void percentIdentity2()
 432   {
 433     seqNums = new Vector();
 434     // calcSeqNum(s);
 435     int i = 0, iSize = sequences.length;
 436     // Do we need to calculate this again?
 437     for (i = 0; i < iSize; i++)
 438     {
 439       calcSeqNum(i);
 440     }
 441
 442     if ((cons2 == null) || seqNumsChanged)
 443     {
 444       cons2 = new int[maxLength][24];
 445
 446       // Initialize the array
 447       for (int j = 0; j < 24; j++)
 448       {
 449         for (i = 0; i < maxLength; i++)
 450         {
 451           cons2[i][j] = 0;
 452         }
 453       }
 454
 455       int[] sqnum;
 456       int j = 0;
 457
 458       while (j < sequences.length)
 459       {
 460         sqnum = (int[]) seqNums.elementAt(j);
 461
 462         for (i = 1; i < sqnum.length; i++)
 463         {
 464           cons2[i - 1][sqnum[i]]++;
 465         }
 466
 467         for (i = sqnum.length - 1; i < maxLength; i++)
 468         {
 469           cons2[i][23]++; // gap count
 470         }
 471
 472         j++;
 473       }
 474
 475       // unnecessary ?
 476
 477       /*
 478        * for (int i=start; i <= end; i++) { int max = -1000; int maxi = -1; int
 479        * maxj = -1;
 480        *
 481        * for (int j=0;j<24;j++) { if (cons2[i][j] > max) { max = cons2[i][j];
 482        * maxi = i; maxj = j; } } }
 483        */
 484     }
 485   }
 486
 487   /**
 488    * Calculates the quality of the set of sequences
 489    *
 490    * @param start
 491    *                Start residue
 492    * @param end
 493    *                End residue
 494    */
 495   public void findQuality(int start, int end)
 496   {
 497     quality = new Vector();
 498
 499     double max = -10000;
 500     int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
 501
 502     // Loop over columns // JBPNote Profiling info
 503     // long ts = System.currentTimeMillis();
 504     // long te = System.currentTimeMillis();
 505     percentIdentity2();
 506
 507     int size = seqNums.size();
 508     int[] lengths = new int[size];
 509     double tot, bigtot, sr, tmp;
 510     double[] x, xx;
 511     int l, j, i, ii, i2, k, seqNum;
 512
 513     for (l = 0; l < size; l++)
 514     {
 515       lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;
 516     }
 517
 518     for (j = start; j <= end; j++)
 519     {
 520       bigtot = 0;
 521
 522       // First Xr = depends on column only
 523       x = new double[24];
 524
 525       for (ii = 0; ii < 24; ii++)
 526       {
 527         x[ii] = 0;
 528
 529         for (i2 = 0; i2 < 24; i2++)
 530         {
 531           x[ii] += (((double) cons2[j][i2] * BLOSUM62[ii][i2]) + 4);
 532         }
 533
 534         x[ii] /= size;
 535       }
 536
 537       // Now calculate D for each position and sum
 538       for (k = 0; k < size; k++)
 539       {
 540         tot = 0;
 541         xx = new double[24];
 542         seqNum = (j < lengths[k]) ? ((int[]) seqNums.elementAt(k))[j + 1]
 543                 : 23; // Sequence, or gap at the end
 544
 545         // This is a loop over r
 546         for (i = 0; i < 23; i++)
 547         {
 548           sr = 0;
 549
 550           sr = (double) BLOSUM62[i][seqNum] + 4;
 551
 552           // Calculate X with another loop over residues
 553           // System.out.println("Xi " + i + " " + x[i] + " " + sr);
 554           xx[i] = x[i] - sr;
 555
 556           tot += (xx[i] * xx[i]);
 557         }
 558
 559         bigtot += Math.sqrt(tot);
 560       }
 561
 562       // This is the quality for one column
 563       if (max < bigtot)
 564       {
 565         max = bigtot;
 566       }
 567
 568       // bigtot = bigtot * (size-cons2[j][23])/size;
 569       quality.addElement(new Double(bigtot));
 570
 571       // Need to normalize by gaps
 572     }
 573
 574     double newmax = -10000;
 575
 576     for (j = start; j <= end; j++)
 577     {
 578       tmp = ((Double) quality.elementAt(j)).doubleValue();
 579       tmp = ((max - tmp) * (size - cons2[j][23])) / size;
 580
 581       // System.out.println(tmp+ " " + j);
 582       quality.setElementAt(new Double(tmp), j);
 583
 584       if (tmp > newmax)
 585       {
 586         newmax = tmp;
 587       }
 588     }
 589
 590     // System.out.println("Quality " + s);
 591     qualityRange[0] = new Double(0);
 592     qualityRange[1] = new Double(newmax);
 593   }
 594 }