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