website/archive/binaries/mac/src/clustalw/src/multipleAlign/ProfileBase.cpp

   1 /**
   2  * Author: Mark Larkin
   3  *
   4  * Copyright (c) 2007 Des Higgins, Julie Thompson and Toby Gibson.
   5  */
   6 #ifdef HAVE_CONFIG_H
   7     #include "config.h"
   8 #endif
   9 #include "ProfileBase.h"
  10
  11 namespace clustalw
  12 {
  13
  14 /**
  15  *
  16  * @param prfLen
  17  * @param firstS
  18  * @param lastS
  19  */
  20 ProfileBase::ProfileBase(int prfLen, int firstS, int lastS)
  21  : vwindow(5),
  22    vll(50),
  23    reducedGap(1.0),
  24    prfLength(prfLen),
  25    firstSeq(firstS),
  26    lastSeq(lastS)
  27 {
  28     vlut.resize(numLetters, vector<int>(numLetters));
  29     for(int i = 0; i < numLetters; i++)
  30     {
  31         vlut[i][i] = 1;
  32     }
  33
  34     pascarellaRes = "ACDEFGHKILMNPQRSTVYW";
  35     int pasprob[] = { 87, 87, 104, 69, 80, 139, 100, 104, 68, 79,
  36                      71, 137, 126, 93, 128, 124, 111, 75, 100, 77};
  37     pascarellaProb = vector<int>(pasprob, pasprob + 20);
  38
  39     profile.resize(prfLength + 2, vector<int>(LENCOL + 2));
  40 }
  41
  42
  43 /**
  44  *
  45  * @param seqArray
  46  * @param gaps
  47  * @param useStructPenalties
  48  * @param gapPenaltyMask
  49  * @param gapCoef
  50  * @param lenCoef
  51  */
  52 void ProfileBase::calcGapCoeff(SeqArray* seqArray, vector<int>* gaps,
  53          bool useStructPenalties, vector<char>* gapPenaltyMask, int gapCoef, int lenCoef)
  54 {
  55     int c;
  56     int i, j;
  57     int is, ie;
  58     int _numSeq, val, pcid;
  59     vector<int> gapPos;
  60     vector<int> vWeight, resWeight, hydWeight;
  61     float scale;
  62     int _maxAA = userParameters->getMaxAA();
  63
  64     _numSeq = lastSeq - firstSeq;
  65     if(_numSeq == 2)
  66     {
  67         pcid = static_cast<int>(percentId(&(*seqArray)[firstSeq], &(*seqArray)[firstSeq + 1]));
  68     }
  69     else
  70         pcid = 0;
  71
  72     for (j = 0; j < prfLength; j++)
  73     {
  74         (*gaps)[j] = 0;
  75     }
  76
  77     // Check for a gap penalty mask
  78
  79     if (useStructPenalties != false)
  80     {
  81         nVarPen = nHydPen = nPrefPen = true;
  82         gdist = 0;
  83     }
  84     else if (userParameters->getNoVarPenalties() == false && pcid > 60)
  85     {
  86         nHydPen = nPrefPen = true;
  87         nVarPen = false;
  88     }
  89     else
  90     {
  91         nVarPen = true;
  92         nHydPen = userParameters->getNoHydPenalties();
  93         nPrefPen = userParameters->getNoPrefPenalties();
  94         gdist = userParameters->getGapDist();
  95     }
  96
  97     for (i = firstSeq; i < lastSeq; i++)
  98     {
  99         // Include end gaps as gaps ?
 100         is = 0;
 101         ie = prfLength;
 102         if (userParameters->getUseEndGaps() == false &&
 103             userParameters->getEndGapPenalties() == false)
 104         {
 105             for (j = 0; j < prfLength; j++)
 106             {
 107                 c = (*seqArray)[i][j];
 108                 if ((c < 0) || (c > _maxAA))
 109                     is++;
 110                 else
 111                     break;
 112             }
 113             for (j = prfLength - 1; j >= 0; j--)
 114             {
 115                 c = (*seqArray)[i][j];
 116                 if ((c < 0) || (c > _maxAA))
 117                     ie--;
 118                 else
 119                     break;
 120             }
 121         }
 122
 123         for (j = is; j < ie; j++)
 124         {
 125             if (((*seqArray)[i][j] < 0) || ((*seqArray)[i][j] > _maxAA))
 126             {
 127                 (*gaps)[j]++;
 128             }
 129         }
 130     }
 131
 132     int _DNAFlag = userParameters->getDNAFlag();
 133
 134     if ((!_DNAFlag) && (nVarPen == false))
 135     {
 136         vWeight.resize(prfLength + 2);
 137         calcVPenalties(seqArray, &vWeight);
 138     }
 139
 140
 141     if ((!_DNAFlag) && (nPrefPen == false))
 142     {
 143         resWeight.resize(prfLength + 2);
 144         calcResidueSpecificPen(seqArray, &resWeight);
 145     }
 146
 147     if ((!_DNAFlag) && (nHydPen == false))
 148     {
 149         hydWeight.resize(prfLength + 2);
 150         calcHydrophilicPen(seqArray, &hydWeight);
 151     }
 152
 153     gapPos.resize(prfLength + 2);
 154
 155     // mark the residues close to an existing gap (set gaps[i] = -ve)
 156
 157     if (_DNAFlag || (gdist <= 0))
 158     {
 159         for (i = 0; i < prfLength; i++)
 160         {
 161             gapPos[i] = (*gaps)[i];
 162         }
 163     }
 164     else
 165     {
 166         i = 0;
 167         while (i < prfLength)
 168         {
 169             if ((*gaps)[i] <= 0)
 170             {
 171                 gapPos[i] = (*gaps)[i];
 172                 i++;
 173             }
 174             else
 175             {
 176                 for (j = -gdist + 1; j < 0; j++)
 177                 {
 178                     if ((i + j >= 0) && (i + j < prfLength) &&
 179                        (((*gaps)[i + j] == 0) || ((*gaps)[i + j] < j)))
 180                     {
 181                         gapPos[i + j] = j;
 182                     }
 183                 }
 184                 while ((*gaps)[i] > 0)
 185                 {
 186                     if (i >= prfLength)
 187                     {
 188                         break;
 189                     }
 190                     gapPos[i] = (*gaps)[i];
 191                     i++;
 192                 }
 193                 for (j = 0; j < gdist; j++)
 194                 {
 195                     if ((*gaps)[i + j] > 0)
 196                     {
 197                         break;
 198                     }
 199                     if ((i + j >= 0) && (i + j < prfLength) &&
 200                        (((*gaps)[i + j] == 0) || ((*gaps)[i + j] < -j)))
 201                     {
 202                          gapPos[i + j] = -j-1;
 203                     }
 204                 }
 205                 i += j;
 206             }
 207         }
 208     }
 209
 210     for (j = 0; j < prfLength; j++)
 211     {
 212         if (gapPos[j] <= 0)
 213         {
 214
 215          // apply residue-specific and hydrophilic gap penalties.
 216
 217             if (!_DNAFlag)
 218             {
 219                 profile[j + 1][GAPCOL] = localPenalty(gapCoef, j, &resWeight, &hydWeight,
 220                                                        &vWeight);
 221                 profile[j+1][LENCOL] = lenCoef;
 222             }
 223             else
 224             {
 225                 profile[j + 1][GAPCOL] = gapCoef;
 226                 profile[j + 1][LENCOL] = lenCoef;
 227             }
 228
 229          // increase gap penalty near to existing gaps.
 230
 231             if (gapPos[j] < 0)
 232             {
 233                 profile[j + 1][GAPCOL] = static_cast<int>((profile[j + 1][GAPCOL] * (2.0 + 2.0 * (gdist + gapPos[j]) / gdist)));
 234             }
 235         }
 236         else
 237         {
 238             scale = ((float)(_numSeq - (*gaps)[j]) / (float)_numSeq) * reducedGap;
 239             profile[j + 1][GAPCOL] = static_cast<int>(scale * gapCoef);
 240             profile[j + 1][LENCOL] = static_cast<int>(0.5 * lenCoef);
 241         }
 242
 243         // apply the gap penalty mask
 244
 245         if (useStructPenalties != NONE)
 246         {
 247             val = (*gapPenaltyMask)[j] - '0';
 248             if (val > 0 && val < 10)
 249             {
 250                 profile[j + 1][GAPCOL] *= val;
 251                 profile[j + 1][LENCOL] *= val;
 252             }
 253         }
 254
 255         // make sure no penalty is zero - even for all-gap positions
 256
 257         if (profile[j + 1][GAPCOL] <= 0)
 258         {
 259             profile[j + 1][GAPCOL] = 1;
 260         }
 261         if (profile[j + 1][LENCOL] <= 0)
 262         {
 263             profile[j + 1][LENCOL] = 1;
 264         }
 265     }
 266
 267     // set the penalties at the beginning and end of the profile
 268     if(userParameters->getEndGapPenalties() == true)
 269     {
 270         profile[0][GAPCOL] = gapCoef;
 271         profile[0][LENCOL] = lenCoef;
 272     }
 273     else
 274     {
 275         profile[0][GAPCOL] = 0;
 276         profile[0][LENCOL] = 0;
 277         profile[prfLength][GAPCOL] = 0;
 278         profile[prfLength][LENCOL] = 0;
 279     }
 280     if (userParameters->getDebug() > 0)
 281     {
 282         cout << "Opening penalties:\n";
 283         for(i = 0; i <= prfLength; i++)
 284         {
 285             cout <<" " << i << ":" << profile[i][GAPCOL]<< " ";
 286         }
 287         cout << "\n";
 288     }
 289     if (userParameters->getDebug() > 0)
 290     {
 291         cout << "Extension penalties:\n";
 292         for(i = 0; i <= prfLength; i++)
 293         {
 294             cout << i << ":" << profile[i][LENCOL] << " ";
 295         }
 296         cout << "\n";
 297     }
 298 }
 299
 300 /** **************************************************************************************
 301  *                               Protected functions                                     *
 302  *****************************************************************************************/
 303
 304
 305 /**
 306  *
 307  * @param aln
 308  * @param weight
 309  */
 310 void ProfileBase::calcVPenalties(SeqArray* aln, vector<int>* weight)
 311 {
 312     int ix1, ix2;
 313     int i, j, t;
 314     int _maxAA = userParameters->getMaxAA();
 315     int aminoCodeix1, aminoCodeix2;
 316
 317     for (i = 0; i < prfLength; i++)
 318     {
 319         (*weight)[i] = 0;
 320         t = 0;
 321         for(j = i - vwindow;j < i + vwindow; j++)
 322         {
 323             if(j >= 0 && j < prfLength)
 324             {
 325                 ix1 = (*aln)[firstSeq][j];
 326                 ix2 = (*aln)[firstSeq + 1][j];
 327                 if ((ix1 < 0) || (ix1 > _maxAA) || (ix2 < 0) || (ix2 > _maxAA))
 328                 {
 329                     continue;
 330                 }
 331                 aminoCodeix1 = userParameters->getAminoAcidCode(ix1);
 332                 aminoCodeix2 = userParameters->getAminoAcidCode(ix2);
 333                 (*weight)[i] += vlut[aminoCodeix1 - 'A'][aminoCodeix2 - 'A'];
 334                 t++;
 335             }
 336         }
 337         /* now we have a weight -t < w < t */
 338         (*weight)[i] +=t;
 339         if(t > 0)
 340             (*weight)[i] = ((*weight)[i] * 100)/(2 * t);
 341         else
 342             (*weight)[i] = 100;
 343         /* now we have a weight vll < w < 100 */
 344         if ((*weight)[i] < vll)
 345             (*weight)[i] = vll;
 346     }
 347 }
 348
 349
 350 /**
 351  *
 352  * @param aln
 353  * @param weight
 354  */
 355 void ProfileBase::calcResidueSpecificPen(SeqArray* aln, vector<int>* weight)
 356 {
 357     char ix;
 358     int j, k, _numSeq;
 359     int i;
 360     int _maxAA = userParameters->getMaxAA();
 361     int _pascarellaNumRes = pascarellaRes.size();
 362
 363     _numSeq = lastSeq - firstSeq;
 364     for (i = 0; i < prfLength; i++)
 365     {
 366         (*weight)[i] = 0;
 367         for (k = firstSeq; k < lastSeq; k++)
 368         {
 369             for (j = 0; j < _pascarellaNumRes; j++)
 370             {
 371                 ix = (*aln)[k][i];
 372                 if ((ix < 0) || (ix > _maxAA))
 373                     continue;
 374                 if (userParameters->getAminoAcidCode(ix) == pascarellaRes[j])
 375                 {
 376                     (*weight)[i] += (180 - pascarellaProb[j]);
 377                     break;
 378                 }
 379             }
 380         }
 381         (*weight)[i] /= _numSeq;
 382     }
 383 }
 384
 385 /**
 386  *
 387  * @param aln
 388  * @param weight
 389  */
 390 void ProfileBase::calcHydrophilicPen(SeqArray* aln, vector<int>* weight)
 391 {
 392     int res;
 393     int numHydResidues, j, k;
 394     int i, e, s;
 395     vector<int> hyd;
 396     float scale;
 397     int _maxAA = userParameters->getMaxAA();
 398
 399     hyd.resize(prfLength + 2);
 400     string _hydResidues(userParameters->getHydResidues());
 401     numHydResidues = _hydResidues.size();
 402
 403     for (i = 0; i < prfLength; i++)
 404     {
 405         (*weight)[i] = 0;
 406     }
 407
 408     for (k = firstSeq; k < lastSeq; k++)
 409     {
 410         for (i = 0; i < prfLength; i++)
 411         {
 412             hyd[i] = 0;
 413             for (j = 0; j < numHydResidues; j++)
 414             {
 415                 res = (*aln)[k][i];
 416                 if ((res < 0) || (res > _maxAA))
 417                     continue;
 418                 if (userParameters->getAminoAcidCode(res) == _hydResidues[j])
 419                 {
 420                     hyd[i] = 1;
 421                     break;
 422                 }
 423             }
 424         }
 425         i = 0;
 426         while (i < prfLength)
 427         {
 428             if (hyd[i] == 0)
 429                 i++;
 430             else
 431             {
 432                 s = i;
 433                 while ((hyd[i] != 0) && (i < prfLength))
 434                 {
 435                     i++;
 436                 }
 437                 e = i;
 438                 if (e - s > 3)
 439                 {
 440                     for (j = s; j < e; j++)
 441                     {
 442                         (*weight)[j] += 100;
 443                     }
 444                 }
 445             }
 446         }
 447     }
 448
 449     scale = lastSeq - firstSeq;
 450     for (i = 0; i < prfLength; i++)
 451     {
 452         (*weight)[i] = static_cast<int>(((*weight)[i] / scale)); // Mark change 17-5-07
 453     }
 454 }
 455
 456
 457 /**
 458  *
 459  * @param penalty
 460  * @param n
 461  * @param resWeight
 462  * @param hydWeight
 463  * @param vWeight
 464  * @return
 465  */
 466 int ProfileBase::localPenalty(int penalty, int n, vector<int>* resWeight,
 467                               vector<int>* hydWeight, vector<int>* vWeight)
 468 {
 469     bool h = false;
 470     float gw;
 471
 472     if (userParameters->getDNAFlag())
 473     {
 474         return(1);
 475     }
 476
 477     gw = 1.0;
 478     if (nVarPen == false)
 479     {
 480         gw *= (*vWeight)[n] / 100.0;
 481     }
 482
 483     if (nHydPen == false)
 484     {
 485         if ((*hydWeight)[n] > 0)
 486         {
 487             gw *= 0.5;
 488             h = true;
 489         }
 490     }
 491     if ((nPrefPen == false) && (h == false))
 492     {
 493          gw *= ((*resWeight)[n] / 100.0);
 494     }
 495
 496     gw *= penalty;
 497     return (int)gw;
 498
 499 }
 500
 501 /** *********************************************************************
 502  * Note dont think this will work. Dont have -3. Need to use lengths!   *
 503  ************************************************************************/
 504
 505 /**
 506  *
 507  * @param s1
 508  * @param s2
 509  * @return
 510  */
 511 float ProfileBase::percentId(vector<int>* s1, vector<int>* s2)
 512 {
 513     int i;
 514     int count, total;
 515     float score;
 516
 517     count = total = 0;
 518     for (i = 0; i < prfLength; i++)
 519     {
 520         if (((*s1)[i] >= 0) && ((*s1)[i] < userParameters->getMaxAA()))
 521         {
 522             total++;
 523             if ((*s1)[i] == (*s2)[i])
 524             {
 525                 count++;
 526             }
 527         }
 528         if ((*s1)[i]==(-3) || (*s2)[i]==(-3))
 529         {
 530             break; // I dont have -3 at the end!
 531         }
 532
 533     }
 534
 535     if(total == 0)
 536     {
 537         score = 0;
 538     }
 539     else
 540     {
 541         score = 100.0 * (float)count / (float)total;
 542     }
 543     return (score);
 544 }
 545
 546
 547 }