src/jalview/analysis/Conservation.java

   1 /*
   2 * Jalview - A Sequence Alignment Editor and Viewer
   3 * Copyright (C) 2006 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 jalview.datamodel.*;
  22
  23 import java.util.*;
  24
  25
  26 /**
  27  * Calculates conservation values for a given set of sequences
  28  *
  29  * @author $author$
  30  * @version $Revision$
  31  */
  32 public class Conservation
  33 {
  34     SequenceI [] sequences;
  35     int start;
  36     int end;
  37     Vector seqNums; // vector of int vectors where first is sequence checksum
  38     int maxLength = 0; //  used by quality calcs
  39     boolean seqNumsChanged = false; // updated after any change via calcSeqNum;
  40     Hashtable [] total;
  41
  42     /** Stores calculated quality values */
  43     public Vector quality;
  44
  45     /** Stores maximum and minimum values of quality values  */
  46     public Double[] qualityRange = new Double[2];
  47     String consString = "";
  48     Sequence consSequence;
  49     Hashtable propHash;
  50     int threshold;
  51     String name = "";
  52     int[][] cons2;
  53
  54     /**
  55      * Creates a new Conservation object.
  56      *
  57      * @param name Name of conservation
  58      * @param propHash DOCUMENT ME!
  59      * @param threshold to count the residues in residueHash(). commonly used value is 3
  60      * @param sequences sequences to be used in calculation
  61      * @param start start residue position
  62      * @param end end residue position
  63      */
  64     public Conservation(String name, Hashtable propHash, int threshold,
  65         Vector sequences, int start, int end)
  66     {
  67
  68         this.name = name;
  69         this.propHash = propHash;
  70         this.threshold = threshold;
  71         this.start = start;
  72         this.end = end;
  73
  74         maxLength=end-start+1; // default width includes bounds of calculation
  75
  76         int s, sSize = sequences.size();
  77         SequenceI[] sarray = new SequenceI[sSize];
  78         this.sequences = sarray;
  79
  80         for (s = 0; s < sSize; s++)
  81         {
  82           sarray[s] = (SequenceI) sequences.elementAt(s);
  83           if(sarray[s].getLength()>maxLength)
  84             maxLength = sarray[s].getLength();
  85         }
  86     }
  87
  88
  89     /**
  90      * DOCUMENT ME!
  91      *
  92      * @param i DOCUMENT ME!
  93      */
  94     private void calcSeqNum(int i)
  95     {
  96         String sq = null; // for dumb jbuilder not-inited exception warning
  97         int[] sqnum = null;
  98
  99         int sSize = sequences.length;
 100
 101         if ((i > -1) && (i < sSize))
 102         {
 103             sq = sequences[i].getSequence();
 104
 105             if (seqNums.size() <= i)
 106             {
 107                 seqNums.addElement(new int[sq.length() + 1]);
 108             }
 109
 110             if (sq.hashCode() != ((int[]) seqNums.elementAt(i))[0])
 111             {
 112                 int j;
 113                 int len;
 114                 seqNumsChanged = true;
 115                 len = sq.length();
 116
 117                 if (maxLength < len)
 118                 {
 119                     maxLength = len;
 120                 }
 121
 122                 sqnum = new int[len + 1]; // better to always make a new array - sequence can change its length
 123                 sqnum[0] = sq.hashCode();
 124
 125                 for (j = 1; j <= len; j++)
 126                 {
 127                     sqnum[j] = jalview.schemes.ResidueProperties.aaIndex[sq.charAt(j-1)];
 128                 }
 129
 130
 131                 seqNums.setElementAt(sqnum, i);
 132             }
 133             else
 134               System.out.println("SEQUENCE HAS BEEN DELETED!!!");
 135         }
 136         else
 137         {
 138             // JBPNote INFO level debug
 139             System.err.println(
 140                 "ERROR: calcSeqNum called with out of range sequence index for Alignment\n");
 141         }
 142     }
 143
 144     /**
 145      * Calculates the conservation values for given set of sequences
 146      */
 147     public void calculate()
 148     {
 149       Hashtable resultHash, ht;
 150       int thresh, j, jSize = sequences.length;
 151       int[] values; // Replaces residueHash
 152       String type, res=null;
 153       char c;
 154       Enumeration enumeration2;
 155
 156       total = new Hashtable[maxLength];
 157
 158       for (int i = start; i <= end; i++)
 159         {
 160             values = new int[132];
 161
 162             for (j = 0; j < jSize; j++)
 163             {
 164               if (sequences[j].getLength() > i)
 165               {
 166                 c = sequences[j].getCharAt(i);
 167
 168                 // No need to check if its a '-'
 169                 if (c == '.' || c == ' ')
 170                   c = '-';
 171
 172                 if ('a' <= c && c <= 'z')
 173                 {
 174                   c -= (32);// 32 = 'a' - 'A'
 175                 }
 176
 177                 values[c]++;
 178               }
 179               else
 180               {
 181                 values['-']++;
 182               }
 183             }
 184
 185             //What is the count threshold to count the residues in residueHash()
 186             thresh = (threshold * (jSize)) / 100;
 187
 188             //loop over all the found residues
 189             resultHash = new Hashtable();
 190             for (int v = '-'; v < 'Z'; v++)
 191             {
 192
 193                 if (values[v] > thresh)
 194                 {
 195                   res =  String.valueOf( (char) v);
 196
 197                     //Now loop over the properties
 198                     enumeration2 = propHash.keys();
 199
 200                     while (enumeration2.hasMoreElements())
 201                     {
 202                         type = (String) enumeration2.nextElement();
 203                         ht = (Hashtable) propHash.get(type);
 204
 205                         //Have we ticked this before?
 206                         if (!resultHash.containsKey(type))
 207                         {
 208                             if (ht.containsKey(res))
 209                             {
 210                                 resultHash.put(type, ht.get(res));
 211                             }
 212                             else
 213                             {
 214                                 resultHash.put(type, ht.get("-"));
 215                             }
 216                         }
 217                         else if (((Integer) resultHash.get(type)).equals(
 218                                     (Integer) ht.get(res)) == false)
 219                         {
 220                             resultHash.put(type, new Integer(-1));
 221                         }
 222                     }
 223                 }
 224             }
 225
 226             total[i-start] = resultHash;
 227         }
 228     }
 229
 230
 231     /***
 232      * countConsNGaps
 233      * returns gap count in int[0], and conserved residue count in int[1]
 234      */
 235     public int[] countConsNGaps(int j)
 236     {
 237         int count = 0;
 238         int cons = 0;
 239         int nres = 0;
 240         int[] r = new int[2];
 241         char f = '$';
 242         int i, iSize = sequences.length;
 243         char c;
 244
 245         for (i = 0; i < iSize; i++)
 246         {
 247             if (j >= sequences[i].getLength())
 248             {
 249                 count++;
 250                 continue;
 251             }
 252
 253             c = sequences[i].getCharAt(j); // gaps do not have upper/lower case
 254
 255             if (jalview.util.Comparison.isGap((c)))
 256             {
 257                 count++;
 258             }
 259             else
 260             {
 261                 nres++;
 262
 263                 if (nres == 1)
 264                 {
 265                     f = c;
 266                     cons++;
 267                 }
 268                 else if (f == c)
 269                 {
 270                     cons++;
 271                 }
 272             }
 273         }
 274
 275         r[0] = (nres == cons) ? 1 : 0;
 276         r[1] = count;
 277
 278         return r;
 279     }
 280
 281     /**
 282      * Calculates the conservation sequence
 283      *
 284      * @param consflag if true, poitiveve conservation; false calculates negative conservation
 285      * @param percentageGaps commonly used value is 25
 286      */
 287     public void verdict(boolean consflag, float percentageGaps)
 288     {
 289         StringBuffer consString = new StringBuffer();
 290         String type;
 291         Integer result;
 292         int[] gapcons;
 293         int totGaps, count;
 294         float pgaps;
 295         Hashtable resultHash ;
 296         Enumeration enumeration;
 297
 298
 299         for (int i = start; i <= end; i++)
 300         {
 301             gapcons = countConsNGaps(i);
 302             totGaps = gapcons[1];
 303             pgaps = ((float) totGaps * 100) / (float) sequences.length;
 304
 305             if (percentageGaps > pgaps)
 306             {
 307                 resultHash =  total[i - start];
 308
 309                 //Now find the verdict
 310                 count = 0;
 311                 enumeration = resultHash.keys();
 312
 313                 while (enumeration.hasMoreElements())
 314                 {
 315                    type = (String) enumeration.nextElement();
 316                    result = (Integer) resultHash.get(type);
 317
 318                     //Do we want to count +ve conservation or +ve and -ve cons.?
 319                     if (consflag)
 320                     {
 321                         if (result.intValue() == 1)
 322                         {
 323                             count++;
 324                         }
 325                     }
 326                     else
 327                     {
 328                         if (result.intValue() != -1)
 329                         {
 330                             count++;
 331                         }
 332                     }
 333                 }
 334
 335                 if (count < 10)
 336                 {
 337                     consString.append(count); // Conserved props!=Identity
 338                 }
 339                 else
 340                 {
 341                     consString.append((gapcons[0] == 1) ? "*" : "+");
 342                 }
 343             }
 344             else
 345             {
 346                 consString.append("-");
 347             }
 348         }
 349
 350         consSequence = new Sequence(name, consString.toString(), start, end);
 351     }
 352
 353     /**
 354      *
 355      *
 356      * @return Conservation sequence
 357      */
 358     public Sequence getConsSequence()
 359     {
 360         return consSequence;
 361     }
 362
 363     // From Alignment.java in jalview118
 364     public void findQuality()
 365     {
 366         findQuality(0, maxLength - 1);
 367     }
 368
 369     /**
 370      * DOCUMENT ME!
 371      */
 372     private void percentIdentity2()
 373     {
 374       seqNums = new Vector();
 375      // calcSeqNum(s);
 376       int i = 0, iSize = sequences.length;
 377     //Do we need to calculate this again?
 378       for (i = 0; i < iSize; i++)
 379       {
 380        calcSeqNum(i);
 381       }
 382
 383
 384         if ((cons2 == null) || seqNumsChanged)
 385         {
 386             cons2 = new int[maxLength][24];
 387
 388             // Initialize the array
 389             for (int j = 0; j < 24; j++)
 390             {
 391                 for (i = 0; i < maxLength; i++)
 392                 {
 393                     cons2[i][j] = 0;
 394                 }
 395             }
 396
 397             int[] sqnum;
 398             int j = 0;
 399
 400             while (j < sequences.length)
 401             {
 402                 sqnum = (int[]) seqNums.elementAt(j);
 403
 404                 for (i = 1; i < sqnum.length; i++)
 405                 {
 406                     cons2[i - 1][sqnum[i]]++;
 407                 }
 408
 409                 for (i = sqnum.length - 1; i < maxLength; i++)
 410                 {
 411                     cons2[i][23]++; // gap count
 412                 }
 413
 414                 j++;
 415             }
 416
 417             // unnecessary ?
 418
 419             /* for (int i=start; i <= end; i++) {
 420                  int max = -1000;
 421             int maxi = -1;
 422             int maxj = -1;
 423
 424             for (int j=0;j<24;j++) {
 425               if (cons2[i][j] > max) {
 426               max = cons2[i][j];
 427               maxi = i;
 428               maxj = j;
 429             }
 430
 431             }
 432             } */
 433         }
 434     }
 435
 436     /**
 437      * Calculates the quality of the set of sequences
 438      *
 439      * @param start Start residue
 440      * @param end End residue
 441      */
 442     public void findQuality(int start, int end)
 443     {
 444         quality = new Vector();
 445
 446         double max = -10000;
 447         int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
 448
 449         //Loop over columns // JBPNote Profiling info
 450         //long ts = System.currentTimeMillis();
 451         //long te = System.currentTimeMillis();
 452         percentIdentity2();
 453
 454         int size = seqNums.size();
 455         int[] lengths = new int[size];
 456         double tot, bigtot, sr, tmp;
 457         double [] x, xx;
 458         int l, j, i, ii, i2, k, seqNum;
 459
 460         for (l = 0; l < size; l++)
 461             lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;
 462
 463         for (j = start; j <= end; j++)
 464         {
 465             bigtot = 0;
 466
 467             // First Xr = depends on column only
 468             x = new double[24];
 469
 470             for (ii = 0; ii < 24; ii++)
 471             {
 472                 x[ii] = 0;
 473
 474                 for (i2 = 0; i2 < 24; i2++)
 475                 {
 476                   x[ii] += ( ( (double) cons2[j][i2] * BLOSUM62[ii][i2]) +
 477                             4);
 478                 }
 479
 480                 x[ii] /= size;
 481             }
 482
 483             // Now calculate D for each position and sum
 484             for (k = 0; k < size; k++)
 485             {
 486                 tot = 0;
 487                 xx = new double[24];
 488                 seqNum = (j < lengths[k])
 489                     ? ((int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end
 490
 491                 // This is a loop over r
 492                 for (i = 0; i < 23; i++)
 493                 {
 494                     sr = 0;
 495
 496                     sr = (double) BLOSUM62[i][seqNum] + 4;
 497
 498                     //Calculate X with another loop over residues
 499                     //  System.out.println("Xi " + i + " " + x[i] + " " + sr);
 500                     xx[i] = x[i] - sr;
 501
 502                     tot += (xx[i] * xx[i]);
 503                 }
 504
 505                 bigtot += Math.sqrt(tot);
 506             }
 507
 508             // This is the quality for one column
 509             if (max < bigtot)
 510             {
 511                 max = bigtot;
 512             }
 513
 514             //      bigtot  = bigtot * (size-cons2[j][23])/size;
 515             quality.addElement(new Double(bigtot));
 516
 517             // Need to normalize by gaps
 518         }
 519
 520         double newmax = -10000;
 521
 522         for (j = start; j <= end; j++)
 523         {
 524             tmp = ((Double) quality.elementAt(j)).doubleValue();
 525             tmp = ((max - tmp) * (size - cons2[j][23])) / size;
 526
 527             //     System.out.println(tmp+ " " + j);
 528             quality.setElementAt(new Double(tmp), j);
 529
 530             if (tmp > newmax)
 531             {
 532                 newmax = tmp;
 533             }
 534         }
 535
 536         //    System.out.println("Quality " + s);
 537         qualityRange[0] = new Double(0);
 538         qualityRange[1] = new Double(newmax);
 539     }
 540 }