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