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         for(int i=0; i<start; i++)
 299           consString.append('-');
 300
 301
 302         for (int i = start; i <= end; i++)
 303         {
 304             gapcons = countConsNGaps(i);
 305             totGaps = gapcons[1];
 306             pgaps = ((float) totGaps * 100) / (float) sequences.length;
 307
 308             if (percentageGaps > pgaps)
 309             {
 310                 resultHash =  total[i - start];
 311
 312                 //Now find the verdict
 313                 count = 0;
 314                 enumeration = resultHash.keys();
 315
 316                 while (enumeration.hasMoreElements())
 317                 {
 318                    type = (String) enumeration.nextElement();
 319                    result = (Integer) resultHash.get(type);
 320
 321                     //Do we want to count +ve conservation or +ve and -ve cons.?
 322                     if (consflag)
 323                     {
 324                         if (result.intValue() == 1)
 325                         {
 326                             count++;
 327                         }
 328                     }
 329                     else
 330                     {
 331                         if (result.intValue() != -1)
 332                         {
 333                             count++;
 334                         }
 335                     }
 336                 }
 337
 338                 if (count < 10)
 339                 {
 340                     consString.append(count); // Conserved props!=Identity
 341                 }
 342                 else
 343                 {
 344                     consString.append((gapcons[0] == 1) ? "*" : "+");
 345                 }
 346             }
 347             else
 348             {
 349                 consString.append('-');
 350             }
 351         }
 352
 353         consSequence = new Sequence(name, consString.toString(), start, end);
 354     }
 355
 356     /**
 357      *
 358      *
 359      * @return Conservation sequence
 360      */
 361     public Sequence getConsSequence()
 362     {
 363         return consSequence;
 364     }
 365
 366     // From Alignment.java in jalview118
 367     public void findQuality()
 368     {
 369         findQuality(0, maxLength - 1);
 370     }
 371
 372     /**
 373      * DOCUMENT ME!
 374      */
 375     private void percentIdentity2()
 376     {
 377       seqNums = new Vector();
 378      // calcSeqNum(s);
 379       int i = 0, iSize = sequences.length;
 380     //Do we need to calculate this again?
 381       for (i = 0; i < iSize; i++)
 382       {
 383        calcSeqNum(i);
 384       }
 385
 386
 387         if ((cons2 == null) || seqNumsChanged)
 388         {
 389             cons2 = new int[maxLength][24];
 390
 391             // Initialize the array
 392             for (int j = 0; j < 24; j++)
 393             {
 394                 for (i = 0; i < maxLength; i++)
 395                 {
 396                     cons2[i][j] = 0;
 397                 }
 398             }
 399
 400             int[] sqnum;
 401             int j = 0;
 402
 403             while (j < sequences.length)
 404             {
 405                 sqnum = (int[]) seqNums.elementAt(j);
 406
 407                 for (i = 1; i < sqnum.length; i++)
 408                 {
 409                     cons2[i - 1][sqnum[i]]++;
 410                 }
 411
 412                 for (i = sqnum.length - 1; i < maxLength; i++)
 413                 {
 414                     cons2[i][23]++; // gap count
 415                 }
 416
 417                 j++;
 418             }
 419
 420             // unnecessary ?
 421
 422             /* for (int i=start; i <= end; i++) {
 423                  int max = -1000;
 424             int maxi = -1;
 425             int maxj = -1;
 426
 427             for (int j=0;j<24;j++) {
 428               if (cons2[i][j] > max) {
 429               max = cons2[i][j];
 430               maxi = i;
 431               maxj = j;
 432             }
 433
 434             }
 435             } */
 436         }
 437     }
 438
 439     /**
 440      * Calculates the quality of the set of sequences
 441      *
 442      * @param start Start residue
 443      * @param end End residue
 444      */
 445     public void findQuality(int start, int end)
 446     {
 447         quality = new Vector();
 448
 449         double max = -10000;
 450         int[][] BLOSUM62 = jalview.schemes.ResidueProperties.getBLOSUM62();
 451
 452         //Loop over columns // JBPNote Profiling info
 453         //long ts = System.currentTimeMillis();
 454         //long te = System.currentTimeMillis();
 455         percentIdentity2();
 456
 457         int size = seqNums.size();
 458         int[] lengths = new int[size];
 459         double tot, bigtot, sr, tmp;
 460         double [] x, xx;
 461         int l, j, i, ii, i2, k, seqNum;
 462
 463         for (l = 0; l < size; l++)
 464             lengths[l] = ((int[]) seqNums.elementAt(l)).length - 1;
 465
 466         for (j = start; j <= end; j++)
 467         {
 468             bigtot = 0;
 469
 470             // First Xr = depends on column only
 471             x = new double[24];
 472
 473             for (ii = 0; ii < 24; ii++)
 474             {
 475                 x[ii] = 0;
 476
 477                 for (i2 = 0; i2 < 24; i2++)
 478                 {
 479                   x[ii] += ( ( (double) cons2[j][i2] * BLOSUM62[ii][i2]) +
 480                             4);
 481                 }
 482
 483                 x[ii] /= size;
 484             }
 485
 486             // Now calculate D for each position and sum
 487             for (k = 0; k < size; k++)
 488             {
 489                 tot = 0;
 490                 xx = new double[24];
 491                 seqNum = (j < lengths[k])
 492                     ? ((int[]) seqNums.elementAt(k))[j + 1] : 23; // Sequence, or gap at the end
 493
 494                 // This is a loop over r
 495                 for (i = 0; i < 23; i++)
 496                 {
 497                     sr = 0;
 498
 499                     sr = (double) BLOSUM62[i][seqNum] + 4;
 500
 501                     //Calculate X with another loop over residues
 502                     //  System.out.println("Xi " + i + " " + x[i] + " " + sr);
 503                     xx[i] = x[i] - sr;
 504
 505                     tot += (xx[i] * xx[i]);
 506                 }
 507
 508                 bigtot += Math.sqrt(tot);
 509             }
 510
 511             // This is the quality for one column
 512             if (max < bigtot)
 513             {
 514                 max = bigtot;
 515             }
 516
 517             //      bigtot  = bigtot * (size-cons2[j][23])/size;
 518             quality.addElement(new Double(bigtot));
 519
 520             // Need to normalize by gaps
 521         }
 522
 523         double newmax = -10000;
 524
 525         for (j = start; j <= end; j++)
 526         {
 527             tmp = ((Double) quality.elementAt(j)).doubleValue();
 528             tmp = ((max - tmp) * (size - cons2[j][23])) / size;
 529
 530             //     System.out.println(tmp+ " " + j);
 531             quality.setElementAt(new Double(tmp), j);
 532
 533             if (tmp > newmax)
 534             {
 535                 newmax = tmp;
 536             }
 537         }
 538
 539         //    System.out.println("Quality " + s);
 540         qualityRange[0] = new Double(0);
 541         qualityRange[1] = new Double(newmax);
 542     }
 543 }