Slice Tools
  • Home
  • SourceForge Page


  • libSlice
  • Home
  • Modules
  • Files
  • Examples
  • Additional Information


  • SourceForge.net Logo
     

    getConsQV.c

    Go to the documentation of this file.
    00001 // $Id: getConsQV.c,v 1.1 2005/07/29 02:55:17 mschatz Exp $
    00002 
    00008 #include <ctype.h>
    00009 #include <stdio.h>
    00010 #include <math.h>
    00011 #include <stdlib.h>
    00012 
    00013 #include "Slice.h"
    00014 
    00016 libSlice_BaseDistribution standardDistribution = {.20, .20, .20, .20, .20};
    00017 
    00018 
    00024 inline int prToQV(long double pr)
    00025 {
    00026   if (pr == 0.0L || pr < 0)
    00027   {
    00028     //fprintf (stderr, "floating point underflow in prToQV\n");
    00029     return MAX_QUALITY_VALUE;
    00030   }
    00031 
    00032   return floor(0.5 - 10.0 * log10(pr));
    00033 }
    00034 
    00039 inline void distributeQV(unsigned int qv, 
    00040                          long double * base,
    00041                          long double * error1, double frequency1,
    00042                          long double * error2, double frequency2,
    00043                          long double * error3, double frequency3,
    00044                          long double * error4, double frequency4)
    00045 {
    00046   if (qv)
    00047   {
    00048     long double ep = pow(0.1, (double) qv/10.0);
    00049     long double sum = frequency1 + frequency2 + frequency3 + frequency4;
    00050 
    00051     // error probability 1-Pr(x_ij=A | s_i=A)
    00052     *base = 1-ep;
    00053 
    00054     // Divide the remaining probability among the other possibilities
    00055     // based on their global frequency
    00056 
    00057     *error1 = ep * frequency1/sum;
    00058     *error2 = ep * frequency2/sum;
    00059     *error3 = ep * frequency3/sum;
    00060     *error4 = ep * frequency4/sum;
    00061   }
    00062 }
    00063 
    00069 
    00071 
    00091 char libSlice_calculateAmbiguityFlags(long double cpA,
    00092                                       long double cpC,
    00093                                       long double cpG,
    00094                                       long double cpT,
    00095                                       long double cpGap,
    00096                                       int highQualityThreshold,
    00097                                       int baseCount)
    00098 {
    00099   int codeSize = 0;
    00100   char ambiguityFlags = 0;
    00101   long double probSum = 0;
    00102   long double baseProb [5] = {cpA, cpC, cpG, cpT, cpGap};
    00103 
    00104   long double cpAmbiguity;
    00105 
    00106   if (highQualityThreshold == 0)
    00107   {
    00108     highQualityThreshold = HIGH_QUALITY_THRESHOLD;
    00109   }
    00110 
    00111   if (baseCount < 1 || baseCount > 5) baseCount = 5;
    00112 
    00113   /*
    00114   fprintf(stdout,
    00115           "cpA=%Le, cpC=%Le, cpG=%Le, cpT=%Le, cpGap=%Le, hq=%d\n",
    00116           cpA, cpC, cpG, cpT, cpGap, highQualityThreshold);
    00117   */
    00118 
    00119   // Find the minimum set of bases such that the sum of the probabilities
    00120   // exceeds the threshold and the size of the code is at most baseCount bases
    00121   do
    00122   {
    00123     int i;
    00124     int min = 0; // Assume position 0 is best
    00125 
    00126     // Find the base with the lowest error probability
    00127     for (i = 1; i < 5; i++)
    00128     {
    00129       if (baseProb[i] < baseProb[min]) min = i;
    00130     }
    00131 
    00132     // Add this base to the running sum of probability of success
    00133     probSum += (1.0 - baseProb[min]);
    00134     baseProb[min] = 2.0; // Guaranteed not to be included in next pass
    00135     codeSize++;
    00136 
    00137     switch (min)
    00138     {
    00139       case 0: ambiguityFlags |= AMBIGUITY_FLAGBIT_A;   break;
    00140       case 1: ambiguityFlags |= AMBIGUITY_FLAGBIT_C;   break;
    00141       case 2: ambiguityFlags |= AMBIGUITY_FLAGBIT_G;   break;
    00142       case 3: ambiguityFlags |= AMBIGUITY_FLAGBIT_T;   break;
    00143       case 4: ambiguityFlags |= AMBIGUITY_FLAGBIT_GAP; break;
    00144     };
    00145 
    00146     if (probSum == 1.0)
    00147     {
    00148       // A QV as low as 200 can be problematic because 1-baseProb[min] converts
    00149       // a probability close to 0 to one close to 1, and it takes more bits to
    00150       // represent a number close to 1.
    00151       //
    00152       // Fortunately, at that level, the consensus should not be ambigious
    00153       
    00154       //fprintf(stdout, "floating point underflow detected in calcAmbFlags\n");
    00155       break;
    00156     }
    00157 
    00158     cpAmbiguity = 1.0-probSum;
    00159   } 
    00160   while ((prToQV(cpAmbiguity) < highQualityThreshold) && 
    00161          (codeSize < baseCount));
    00162 
    00163   return ambiguityFlags;
    00164 }
    00165 
    00167 
    00173 char libSlice_convertAmbiguityFlags(char ambiguityFlags)
    00174 {
    00175   int addA   = (ambiguityFlags & AMBIGUITY_FLAGBIT_A);
    00176   int addC   = (ambiguityFlags & AMBIGUITY_FLAGBIT_C);
    00177   int addG   = (ambiguityFlags & AMBIGUITY_FLAGBIT_G);
    00178   int addT   = (ambiguityFlags & AMBIGUITY_FLAGBIT_T);
    00179   int addGap = (ambiguityFlags & AMBIGUITY_FLAGBIT_GAP);
    00180   
    00181   char consensus = '?';  
    00182 
    00183   // Pick IUPAC code 
    00184  
    00185   if      (addA && addC && addG && addT) consensus = 'N';
    00186   else if (addA && addC && addG)         consensus = 'V';
    00187   else if (addA && addC && addT)         consensus = 'H';
    00188   else if (addA && addG && addT)         consensus = 'D';
    00189   else if (addC && addG && addT)         consensus = 'B';
    00190   else if (addA && addT)                 consensus = 'W';
    00191   else if (addA && addC)                 consensus = 'M';
    00192   else if (addA && addG)                 consensus = 'R';
    00193   else if (addC && addG)                 consensus = 'S';
    00194   else if (addC && addT)                 consensus = 'Y';
    00195   else if (addG && addT)                 consensus = 'K';
    00196   else if (addA)                         consensus = 'A';
    00197   else if (addC)                         consensus = 'C';
    00198   else if (addG)                         consensus = 'G';
    00199   else if (addT)                         consensus = 'T';
    00200   else if (addGap)                       consensus = '-';
    00201 
    00202   // Encode extended ambiguity codes as well
    00203   if ((consensus != '-') && addGap)
    00204   {
    00205     consensus = tolower(consensus);
    00206   }
    00207 
    00208   return consensus;
    00209 }
    00210 
    00211 
    00213 
    00229 char libSlice_getAmbiguityCode(long double cpA,
    00230                                long double cpC,
    00231                                long double cpG,
    00232                                long double cpT,
    00233                                long double cpGap,
    00234                                int highQualityThreshold,
    00235                                int baseCount)
    00236 {
    00237   char ambiguityFlags = libSlice_calculateAmbiguityFlags(cpA, cpC, cpG, cpT, 
    00238                                                          cpGap, 
    00239                                                          highQualityThreshold, 
    00240                                                          baseCount);
    00241   return libSlice_convertAmbiguityFlags(ambiguityFlags);
    00242 }
    00243 
    00245 
    00250 
    00251 // Not thread safe, but protects the API
    00252 int m_recallEmpty = 0;
    00253 
    00255 
    00259 void libSlice_setRecallEmpty(int recallEmpty)
    00260 {
    00261   m_recallEmpty = recallEmpty;
    00262 }
    00263 
    00265 
    00280 int libSlice_getConsensusParam(const libSlice_Slice *s, 
    00281                                libSlice_Consensus * result, 
    00282                                const libSlice_BaseDistribution * dist,
    00283                                int highQualityThreshold,
    00284                                int doAmbiguity) 
    00285 {
    00286   int retval = 0;
    00287     
    00288   // quality value multiplicities
    00289   unsigned int qvMultA   = 0; 
    00290   unsigned int qvMultC   = 0;
    00291   unsigned int qvMultG   = 0;
    00292   unsigned int qvMultT   = 0;
    00293   unsigned int qvMultGap = 0;
    00294   
    00295   // consensus quality value parameters
    00296   // current consensus base: s_i, the column of base calls above it: x_ij
    00297   long double pr_A_A     = 1.0; // Pr(x_ij=A | s_i=A)
    00298   long double pr_C_A     = 1.0; // Pr(x_ij=C | s_i=A)
    00299   long double pr_G_A     = 1.0; // Pr(x_ij=G | s_i=A)
    00300   long double pr_T_A     = 1.0; // Pr(x_ij=T | s_i=A)
    00301   long double pr_GAP_A   = 1.0;
    00302   
    00303   long double pr_A_C     = 1.0;
    00304   long double pr_C_C     = 1.0;
    00305   long double pr_G_C     = 1.0;
    00306   long double pr_T_C     = 1.0;
    00307   long double pr_GAP_C   = 1.0;
    00308   
    00309   long double pr_A_G     = 1.0;
    00310   long double pr_C_G     = 1.0;
    00311   long double pr_G_G     = 1.0;
    00312   long double pr_T_G     = 1.0;
    00313   long double pr_GAP_G   = 1.0;
    00314 
    00315   long double pr_A_T     = 1.0;
    00316   long double pr_C_T     = 1.0;
    00317   long double pr_G_T     = 1.0;
    00318   long double pr_T_T     = 1.0;
    00319   long double pr_GAP_T   = 1.0;
    00320 
    00321   long double pr_A_GAP   = 1.0;
    00322   long double pr_C_GAP   = 1.0;
    00323   long double pr_G_GAP   = 1.0;
    00324   long double pr_T_GAP   = 1.0;
    00325   long double pr_GAP_GAP = 1.0;
    00326 
    00327   long double   pA,  pC,  pG,  pT,  pGap, sum;
    00328   long double  cpA, cpC, cpG, cpT, cpGap, cpConsensus;
    00329   unsigned int qvA, qvC, qvG, qvT, qvGap;
    00330 
    00331   char consensus;        // calculated consensus
    00332   int baseCount = 0;     // number of bases seen in read
    00333   unsigned int cqv = 0;  // calculated consensus quality value
    00334 
    00335   if (!result || !s)
    00336   {
    00337     // Result can't be NULL
    00338     return -1;
    00339   }
    00340 
    00341   if (highQualityThreshold < 0)
    00342   {
    00343     doAmbiguity = 0;
    00344   }
    00345 
    00346   if (!s->dcov && !m_recallEmpty)
    00347   {
    00348     result->qvA   = 0;
    00349     result->qvC   = 0;
    00350     result->qvG   = 0;
    00351     result->qvT   = 0;
    00352     result->qvGap = 0;
    00353 
    00354     result->cpA   = 0.2;
    00355     result->cpC   = 0.2;
    00356     result->cpG   = 0.2;
    00357     result->cpT   = 0.2;
    00358     result->cpGap = 0.2;
    00359 
    00360     result->qvConsensus = 0;
    00361     result->ambiguityFlags = 0;
    00362 
    00363     // Just set the consensensus to be the old consensus
    00364     result->consensus = s->c;
    00365   }
    00366   else
    00367   {
    00368     if (s->dcov)
    00369     {
    00370       // Loop through the bases and sum the qv for each component
    00371       
    00372       int qv;
    00373       int j;
    00374       char a;
    00375 
    00376       for(j = 0; j < s->dcov; j++)
    00377       {
    00378         a = toupper(s->bc[j]);
    00379         qv = (s->qv[j] < MIN_QUALITY_VALUE) ? MIN_QUALITY_VALUE : s->qv[j];
    00380 
    00381         switch(a)
    00382         {
    00383           case 'A': qvMultA   += qv; break;
    00384           case 'C': qvMultC   += qv; break;
    00385           case 'G': qvMultG   += qv; break;
    00386           case 'T': qvMultT   += qv; break;
    00387           case '-': qvMultGap += qv; break;
    00388 
    00389           default:
    00390               // Ambiguity codes cannot be processed
    00391               // fprintf(stderr, "Bad input=%c\n", a);
    00392               break;
    00393         };
    00394       }
    00395     }
    00396     else
    00397     {
    00398       // An empty slice is a gap (by definition)
    00399 
    00400       qvMultGap = GAP_QUALITY_VALUE_EMPTY_SLICE;
    00401     }
    00402 
    00403     if (qvMultA)   baseCount++;
    00404     if (qvMultC)   baseCount++;
    00405     if (qvMultG)   baseCount++;
    00406     if (qvMultT)   baseCount++;
    00407     if (qvMultGap) baseCount++;
    00408 
    00409     if (!dist) dist = &standardDistribution;
    00410 
    00411     distributeQV(qvMultA, 
    00412                  &pr_A_A,
    00413                  &pr_A_C,   dist->freqC, 
    00414                  &pr_A_G,   dist->freqG,
    00415                  &pr_A_T,   dist->freqT, 
    00416                  &pr_A_GAP, dist->freqGap);
    00417 
    00418     distributeQV(qvMultC, 
    00419                  &pr_C_C,
    00420                  &pr_C_A,   dist->freqA, 
    00421                  &pr_C_G,   dist->freqG,
    00422                  &pr_C_T,   dist->freqT, 
    00423                  &pr_C_GAP, dist->freqGap);
    00424 
    00425     distributeQV(qvMultG, 
    00426                  &pr_G_G,
    00427                  &pr_G_A,   dist->freqA, 
    00428                  &pr_G_C,   dist->freqC,
    00429                  &pr_G_T,   dist->freqT, 
    00430                  &pr_G_GAP, dist->freqGap);
    00431 
    00432     distributeQV(qvMultT, 
    00433                  &pr_T_T,
    00434                  &pr_T_A,   dist->freqA, 
    00435                  &pr_T_C,   dist->freqC,
    00436                  &pr_T_G,   dist->freqG, 
    00437                  &pr_T_GAP, dist->freqGap);
    00438 
    00439     distributeQV(qvMultGap, 
    00440                  &pr_GAP_GAP,
    00441                  &pr_GAP_A, dist->freqA, 
    00442                  &pr_GAP_C, dist->freqC,
    00443                  &pr_GAP_G, dist->freqG, 
    00444                  &pr_GAP_T, dist->freqT);
    00445 
    00446     // Calculate raw values
    00447     pA   = pr_A_A   * pr_C_A   * pr_G_A   * pr_T_A   * pr_GAP_A;
    00448     pC   = pr_A_C   * pr_C_C   * pr_G_C   * pr_T_C   * pr_GAP_C;
    00449     pG   = pr_A_G   * pr_C_G   * pr_G_G   * pr_T_G   * pr_GAP_G;
    00450     pT   = pr_A_T   * pr_C_T   * pr_G_T   * pr_T_T   * pr_GAP_T;
    00451     pGap = pr_A_GAP * pr_C_GAP * pr_G_GAP * pr_T_GAP * pr_GAP_GAP;
    00452 
    00453     sum = pA + pC + pG + pT + pGap;
    00454 
    00455     
    00456     /*
    00457     fprintf(stderr,
    00458             "sum=%Lg, pA=%Lg, pC=%Lg, pG=%Lg, pT=%Lg, pGap=%Lg\n",
    00459              sum,     pA,     pC,     pG,     pT,     pGap);
    00460     */
    00461     
    00462 
    00463     // Normalize values and computer error probabilities
    00464     // Note: Error probability is the sum of the probability of the other bases
    00465     
    00466     if (sum != 0.0)
    00467     {
    00468       cpA   = (0.0 + pC  + pG  + pT  + pGap) / sum;
    00469       cpC   = (pA  + 0.0 + pG  + pT  + pGap) / sum;
    00470       cpG   = (pA  + pC  + 0.0 + pT  + pGap) / sum;
    00471       cpT   = (pA  + pC  + pG  + 0.0 + pGap) / sum;
    00472       cpGap = (pA  + pC  + pG  + pT  + 0.0)  / sum;
    00473     }
    00474     else
    00475     {
    00476       // Sum can go to zero by floating point underflow
    00477       cpA = cpC = cpG = cpT = cpGap = .20;
    00478     }
    00479 
    00480     // Note: cpA + cpC + cpG + cpG + cpT + cpGap == 1
    00481 
    00482     
    00483     /*
    00484     fprintf(stderr,
    00485             "sum=%Lg, cpA=%Lg, cpC=%Lg, cpG=%Lg, cpL=%Lg, cpGap=%Lg\n",
    00486              sum,     cpA,     cpC,     cpG,     cpT,     cpGap);
    00487     */
    00488     
    00489 
    00490     // Convert error probabilities into quality values
    00491     qvA   = prToQV(cpA);
    00492     qvC   = prToQV(cpC);
    00493     qvG   = prToQV(cpG);
    00494     qvT   = prToQV(cpT);
    00495     qvGap = prToQV(cpGap);
    00496 
    00497     // The consensus quality value is the maximum of the quality 
    00498     // values for each base. The base that has the lowest probability of
    00499     // error is the unambiguous consensus of the slice.
    00500     
    00501     // Find the absolute winner
    00502     consensus = 'A'; cqv = qvA; cpConsensus = cpA;
    00503 
    00504     // Allow a small margin to handle floating point weirdness
    00505     if ((cpC - cpConsensus) < -0.000000001) 
    00506     {
    00507       consensus = 'C'; cqv = qvC; cpConsensus = cpC;
    00508     }
    00509     
    00510     if ((cpG - cpConsensus) < -0.000000001) 
    00511     {
    00512       consensus = 'G'; cqv = qvG; cpConsensus = cpG;
    00513     }
    00514     
    00515     if ((cpT - cpConsensus) < -0.000000001) 
    00516     {
    00517       consensus = 'T'; cqv = qvT; cpConsensus = cpT;
    00518     }
    00519     
    00520     if ((cpGap-cpConsensus) < -0.000000001) 
    00521     {
    00522       consensus = '-'; cqv=qvGap; cpConsensus = cpGap;
    00523     }
    00524 
    00525     /*
    00526     fprintf(stderr, "cp: A %s -\n", (cpA==cpGap)?"==":"!=");
    00527     fprintf(stderr, "qv: A %s -\n", (qvA==qvGap)?"==":"!=");
    00528 
    00529     fprintf(stderr, "%Lf\n", cpGap - cpA);
    00530     */
    00531 
    00532     result->qvA   = qvA;
    00533     result->qvC   = qvC;
    00534     result->qvG   = qvG;
    00535     result->qvT   = qvT;
    00536     result->qvGap = qvGap;
    00537 
    00538     result->cpA   = cpA;
    00539     result->cpC   = cpC;
    00540     result->cpG   = cpG;
    00541     result->cpT   = cpT;
    00542     result->cpGap = cpGap;
    00543 
    00544     result->qvConsensus = cqv;
    00545 
    00546     // This finds all ambiguity flags
    00547     result->ambiguityFlags = 
    00548             libSlice_calculateAmbiguityFlags(cpA, cpC, cpG, cpT, cpGap,
    00549                                              highQualityThreshold,
    00550                                              baseCount);
    00551     if (doAmbiguity)
    00552     {
    00553       consensus = libSlice_convertAmbiguityFlags(result->ambiguityFlags);
    00554     }
    00555 
    00556     result->consensus = consensus;
    00557   }
    00558 
    00559   return retval;
    00560 }
    00561 
    00563 int libSlice_updateAmbiguityConic(const libSlice_Slice * s, 
    00564                                   libSlice_Consensus * consensus,
    00565                                   int highQuality)
    00566 {
    00567   int retval = 0;
    00568 
    00569   if (highQuality == 0) { highQuality = HIGH_QUALITY_THRESHOLD; }
    00570 
    00571   if (!s || !consensus)
    00572   {
    00573     retval = -1;
    00574   }
    00575   else if (s->dcov == 0)
    00576   {
    00577     if (m_recallEmpty)
    00578     {
    00579       // Zero Coverage, assign gap
    00580       consensus->ambiguityFlags = AMBIGUITY_FLAGBIT_GAP;
    00581     }
    00582   }
    00583   else
    00584   {
    00585     int i;
    00586     char ambiguityFlags = consensus->ambiguityFlags;
    00587 
    00588     unsigned int qvArr[5]      = {0,0,0,0,0};
    00589     unsigned int countArr[5]   = {0,0,0,0,0};
    00590     unsigned int hqCountArr[5] = {0,0,0,0,0};
    00591 
    00592     unsigned int maxQV = 0;
    00593     unsigned int secondQV = -1;
    00594 
    00595     int base;
    00596     int qv;
    00597     char a;
    00598 
    00599     // Figure out base counts and sum quality values
    00600     for (i = 0; i < s->dcov; i++)
    00601     {
    00602       a = toupper(s->bc[i]);
    00603       switch(a)
    00604       {
    00605         case 'A': base = 0; break;
    00606         case 'C': base = 1; break;
    00607         case 'G': base = 2; break;
    00608         case 'T': base = 3; break;
    00609         case '-': base = 4; break;
    00610 
    00611         default:
    00612             // fprintf(stderr, "Bad input=%c\n", a);
    00613             // Ambiguity codes cannot be processed
    00614             base = -1;
    00615       };
    00616 
    00617       if (base != -1)
    00618       {
    00619         qv = (s->qv[i] < MIN_QUALITY_VALUE) ? MIN_QUALITY_VALUE : s->qv[i];
    00620         qvArr[base] += qv;
    00621 
    00622         countArr[base]++;
    00623         if (s->qv[i] >= highQuality) { hqCountArr[base]++; }
    00624       }
    00625     }
    00626 
    00627     // Find the two bases with the highest qv and top two highest counts
    00628     // Start at 1 by assuming 0 is max
    00629     for (i = 1; i < 5; i++)
    00630     {
    00631       if (qvArr[i] > qvArr[maxQV])
    00632       {
    00633         secondQV = maxQV;
    00634         maxQV = i;
    00635       }
    00636       else if (secondQV == -1 || qvArr[i] > qvArr[secondQV])
    00637       {
    00638         secondQV = i;
    00639       }
    00640     }
    00641 
    00642     if (1) // CONIC_SCORING
    00643     {
    00644       double ambiguityAngle = 36.8698977;
    00645       double effectiveAngle = ((ambiguityAngle / 2) <= 45) 
    00646                               ? ((double)ambiguityAngle)/2 : 45;
    00647 
    00648       double lowerRadians = (M_PI/(double) 180) * (45.0 - effectiveAngle);
    00649       double upperRadians = (M_PI/(double) 180) * (45.0 + effectiveAngle);
    00650 
    00651       double lowerlimit = tan(lowerRadians);
    00652       double upperlimit = tan(upperRadians);
    00653 
    00654       int component = 0;
    00655 
    00656       ambiguityFlags = 0;
    00657 
    00658       for (component = 0; component < 5; component++)
    00659       {
    00660         // Initialize this to be effectively infinite in case qvArr[maxQV] == 0
    00661         // but set to be zero if qvArr[component] == 0
    00662         double tangent = (qvArr[maxQV]) ? 
    00663                          ((double) qvArr[component]) / qvArr[maxQV] : 
    00664                          qvArr[component] * 10000000.0;
    00665         
    00666         // Everything is compared to the max, except itself
    00667         if (component == maxQV) { continue; }
    00668 
    00669         /*
    00670         fprintf(stderr, "CSS: ll = %0.5f, uu=%0.5f, tan=%0.5f\n",
    00671                 lowerlimit, upperlimit, tangent);
    00672 
    00673         */
    00674 
    00675         if (tangent < lowerlimit)
    00676         {
    00677           // Consensus in non-ambiguous for maxQV
    00678           // Normal case
    00679 
    00680           ambiguityFlags |= (1 << maxQV);
    00681         }
    00682         else if (tangent < upperlimit)
    00683         {
    00684           // Consensus is ambiguous between maxQV and component
    00685           ambiguityFlags |= (1 << maxQV);
    00686           ambiguityFlags |= (1 << component);
    00687         }
    00688         else
    00689         {
    00690           // Consensus is non-ambiguous for component
    00691           // This should never be the case
    00692           ambiguityFlags |= (1 << component);
    00693         }
    00694       }
    00695     }
    00696 
    00697     // Set the new ambiguity flags
    00698     consensus->ambiguityFlags = ambiguityFlags;
    00699   }
    00700 
    00701   return retval;
    00702 }
    00703 
    00704 
    00706 
    00716 int libSlice_updateAmbiguity(const libSlice_Slice * s, 
    00717                              libSlice_Consensus * consensus,
    00718                              int highQuality)
    00719 {
    00720   int retval = 0;
    00721 
    00722   if (!s || !consensus)
    00723   {
    00724     retval = -1;
    00725   }
    00726   else
    00727   {
    00728     char ambiguityFlags = 0;
    00729 
    00730     // Get the original flags
    00731     int addA   = ((consensus->ambiguityFlags & AMBIGUITY_FLAGBIT_A) != 0);
    00732     int addC   = ((consensus->ambiguityFlags & AMBIGUITY_FLAGBIT_C) != 0);
    00733     int addG   = ((consensus->ambiguityFlags & AMBIGUITY_FLAGBIT_G) != 0);
    00734     int addT   = ((consensus->ambiguityFlags & AMBIGUITY_FLAGBIT_T) != 0);
    00735     int addGap = ((consensus->ambiguityFlags & AMBIGUITY_FLAGBIT_GAP) != 0);
    00736     
    00737     int j;
    00738 
    00739     if (highQuality == 0) { highQuality = HIGH_QUALITY_THRESHOLD; }
    00740     
    00741     // Add high quality conflicts
    00742     for(j = 0; j < s->dcov; j++)
    00743     {
    00744       int qv = (s->qv[j] < MIN_QUALITY_VALUE) ? MIN_QUALITY_VALUE : s->qv[j];
    00745 
    00746       if (qv >= highQuality)
    00747       {
    00748         char a = toupper(s->bc[j]);
    00749         switch(a)
    00750         {
    00751           case 'A': addA = 1;   break;
    00752           case 'C': addC = 1;   break;
    00753           case 'G': addG = 1;   break;
    00754           case 'T': addT = 1;   break;
    00755           case '-': addGap = 1; break;
    00756 
    00757           default:
    00758               // fprintf(stderr, "Bad input=%c\n", a);
    00759               // Ambiguity codes cannot be processed
    00760               break;
    00761         };
    00762       }
    00763     }
    00764 
    00765     if (addA)   ambiguityFlags += AMBIGUITY_FLAGBIT_A;
    00766     if (addC)   ambiguityFlags += AMBIGUITY_FLAGBIT_C;
    00767     if (addG)   ambiguityFlags += AMBIGUITY_FLAGBIT_G;
    00768     if (addT)   ambiguityFlags += AMBIGUITY_FLAGBIT_T;
    00769     if (addGap) ambiguityFlags += AMBIGUITY_FLAGBIT_GAP;
    00770 
    00771     // Set the new ambiguity flags
    00772     consensus->ambiguityFlags = ambiguityFlags;
    00773   }
    00774 
    00775   return retval;
    00776 }
    00777 
    00779 
    00788 int libSlice_getConsensus(const libSlice_Slice *s, 
    00789                           libSlice_Consensus * result, 
    00790                           const libSlice_BaseDistribution * dist,
    00791                           int highQualityThreshold)
    00792 {
    00793   return libSlice_getConsensusParam(s, result, dist, highQualityThreshold, 0);
    00794 }
    00795 
    00797 
    00806 int libSlice_getConsensusWithAmbiguity(const libSlice_Slice *s, 
    00807                                        libSlice_Consensus * result, 
    00808                                        const libSlice_BaseDistribution * dist,
    00809                                        int highQualityThreshold)
    00810 {
    00811   return libSlice_getConsensusParam(s, result, dist, highQualityThreshold, 1);
    00812 }
    00813 
    00815 
    00816 
    00821 
    00823 
    00834 int libSlice_getConsensusRangeParam(const libSlice_Slice s [], 
    00835                                     libSlice_Consensus results [], 
    00836                                     int len,
    00837                                     const libSlice_BaseDistribution * dist,
    00838                                     int highQualityThreshold,
    00839                                     int doAmbiguity)
    00840 {
    00841   int retval = 0;
    00842   int i;
    00843 
    00844   for(i = 0; i < len; i++)
    00845   {
    00846     retval += libSlice_getConsensusParam(&s[i], &results[i], dist,
    00847                                          highQualityThreshold, doAmbiguity);
    00848   }
    00849   
    00850   return retval;
    00851 }
    00852 
    00854 
    00864 int libSlice_getConsensusRange(const libSlice_Slice s [], 
    00865                                libSlice_Consensus results [], 
    00866                                int len,
    00867                                const libSlice_BaseDistribution * dist,
    00868                                int highQualityThreshold)
    00869 {
    00870   return libSlice_getConsensusRangeParam(s, results, len, dist, 
    00871                                          highQualityThreshold, 0);
    00872 }
    00873 
    00874 
    00876 
    00886 int libSlice_getConsensusRangeWithAmbiguity(const libSlice_Slice s [], 
    00887                                          libSlice_Consensus results [], 
    00888                                          int len,
    00889                                          const libSlice_BaseDistribution * dist,
    00890                                          int highQualityThreshold)
    00891 {
    00892   return libSlice_getConsensusRangeParam(s, results, len, dist, 
    00893                                          highQualityThreshold, 1);
    00894 }
    00895 
    00897