src/jalview/analysis/Conservation.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer
   3  * Copyright (C) 2007 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                                   // map all symbols to canonical aa numbering
  49                                   // rather than consider conservation of that
  50                                   // symbol
  51
  52   /** Stores calculated quality values */
  53   public Vector quality;
  54
  55   /** Stores maximum and minimum values of quality values */
  56   public Double[] qualityRange = new Double[2];
  57
  58   String consString = "";
  59
  60   Sequence consSequence;
  61
  62   Hashtable propHash;
  63
  64   int threshold;
  65
  66   String name = "";
  67
  68   int[][] cons2;
  69
  70   /**
  71    * Creates a new Conservation object.
  72    *
  73    * @param name
  74    *                Name of conservation
  75    * @param propHash
  76    *                hash of properties for each symbol
  77    * @param threshold
  78    *                to count the residues in residueHash(). commonly used value
  79    *                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
 116    * representation and store it in the 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    * @return { gap count, conserved residue count}
 281    */
 282   public int[] countConsNGaps(int j)
 283   {
 284     int count = 0;
 285     int cons = 0;
 286     int nres = 0;
 287     int[] r = new int[2];
 288     char f = '$';
 289     int i, iSize = sequences.length;
 290     char c;
 291
 292     for (i = 0; i < iSize; i++)
 293     {
 294       if (j >= sequences[i].getLength())
 295       {
 296         count++;
 297         continue;
 298       }
 299
 300       c = sequences[i].getCharAt(j); // gaps do not have upper/lower case
 301
 302       if (jalview.util.Comparison.isGap((c)))
 303       {
 304         count++;
 305       }
 306       else
 307       {
 308         nres++;
 309
 310         if (nres == 1)
 311         {
 312           f = c;
 313           cons++;
 314         }
 315         else if (f == c)
 316         {
 317           cons++;
 318         }
 319       }
 320     }
 321
 322     r[0] = (nres == cons) ? 1 : 0;
 323     r[1] = count;
 324
 325     return r;
 326   }
 327
 328   /**
 329    * Calculates the conservation sequence
 330    *
 331    * @param consflag
 332    *                if true, poitiveve conservation; false calculates negative
 333    *                conservation
 334    * @param percentageGaps
 335    *                commonly used value is 25
 336    */
 337   public void verdict(boolean consflag, float percentageGaps)
 338   {
 339     StringBuffer consString = new StringBuffer();
 340     String type;
 341     Integer result;
 342     int[] gapcons;
 343     int totGaps, count;
 344     float pgaps;
 345     Hashtable resultHash;
 346     Enumeration enumeration;
 347
 348     // NOTE THIS SHOULD CHECK IF THE CONSEQUENCE ALREADY
 349     // EXISTS AND NOT OVERWRITE WITH '-', BUT THIS CASE
 350     // DOES NOT EXIST IN JALVIEW 2.1.2
 351     for (int i = 0; i < start; i++)
 352     {
 353       consString.append('-');
 354     }
 355
 356     for (int i = start; i <= end; i++)
 357     {
 358       gapcons = countConsNGaps(i);
 359       totGaps = gapcons[1];
 360       pgaps = ((float) totGaps * 100) / (float) sequences.length;
 361
 362       if (percentageGaps > pgaps)
 363       {
 364         resultHash = total[i - start];
 365
 366         // Now find the verdict
 367         count = 0;
 368         enumeration = resultHash.keys();
 369
 370         while (enumeration.hasMoreElements())
 371         {
 372           type = (String) enumeration.nextElement();
 373           result = (Integer) resultHash.get(type);
 374
 375           // Do we want to count +ve conservation or +ve and -ve cons.?
 376           if (consflag)
 377           {
 378             if (result.intValue() == 1)
 379             {
 380               count++;
 381             }
 382           }
 383           else
 384           {
 385             if (result.intValue() != -1)
 386             {
 387               count++;
 388             }
 389           }
 390         }
 391
 392         if (count < 10)
 393         {
 394           consString.append(count); // Conserved props!=Identity
 395         }
 396         else
 397         {
 398           consString.append((gapcons[0] == 1) ? "*" : "+");
 399         }
 400       }
 401       else
 402       {
 403         consString.append('-');
 404       }
 405     }
 406
 407     consSequence = new Sequence(name, consString.toString(), start, end);
 408   }
 409
 410   /**
 411    *
 412    *
 413    * @return Conservation sequence
 414    */
 415   public Sequence getConsSequence()
 416   {
 417     return consSequence;
 418   }
 419
 420   // From Alignment.java in jalview118
 421   public void findQuality()
 422   {
 423     findQuality(0, maxLength - 1);
 424   }
 425
 426   /**
 427    * DOCUMENT ME!
 428    */
 429   private void percentIdentity2()
 430   {
 431     seqNums = new Vector();
 432     // calcSeqNum(s);
 433     int i = 0, iSize = sequences.length;
 434     // Do we need to calculate this again?
 435     for (i = 0; i < iSize; i++)
 436     {
 437       calcSeqNum(i);
 438     }
 439
 440     if ((cons2 == null) || seqNumsChanged)
 441     {
 442       cons2 = new int[maxLength][24];
 443
 444       // Initialize the array
 445       for (int j = 0; j < 24; j++)
 446       {
 447         for (i = 0; i < maxLength; i++)
 448         {
 449           cons2[i][j] = 0;
 450         }
 451       }
 452
 453       int[] sqnum;
 454       int j = 0;
 455
 456       while (j < sequences.length)
 457       {
 458         sqnum = (int[]) seqNums.elementAt(j);
 459
 460         for (i = 1; i < sqnum.length; i++)
 461         {
 462           cons2[i - 1][sqnum[i]]++;
 463         }
 464
 465         for (i = sqnum.length - 1; i < maxLength; i++)
 466         {
 467           cons2[i][23]++; // gap count
 468         }
 469
 470         j++;
 471       }
 472
 473       // unnecessary ?
 474
 475       /*
 476        * for (int i=start; i <= end; i++) { int max = -1000; int maxi = -1; int
 477        * maxj = -1;
 478        *
 479        * for (int j=0;j<24;j++) { if (cons2[i][j] > max) { max = cons2[i][j];
 480        * maxi = i; maxj = j; }
 481        *  } }
 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 }