Difference between revisions of "830 - BestAlignment::MaqAlignment()"
From Genome Analysis Wiki
Jump to navigationJump to search (Created page with "void BestAlignment::MaqAlignment(MaqIndex & index1, MaqIndex & index2, MaqIndex & index3, FastQ & read) { position = -1; sumQ = mapQ = 0; mismatches = read.Length(); forwa...") |
|||
| Line 1: | Line 1: | ||
| + | <source lang="cpp"> | ||
| + | |||
void BestAlignment::MaqAlignment(MaqIndex & index1, MaqIndex & index2, MaqIndex & index3, FastQ & read) | void BestAlignment::MaqAlignment(MaqIndex & index1, MaqIndex & index2, MaqIndex & index3, FastQ & read) | ||
{ | { | ||
| Line 37: | Line 39: | ||
mapQ = quality < 0.5 ? 0 : (quality > 0.9999999998 ? 99 : (int) (0.5 + 10 * log(1.0 - quality) / log(0.1))); | mapQ = quality < 0.5 ? 0 : (quality > 0.9999999998 ? 99 : (int) (0.5 + 10 * log(1.0 - quality) / log(0.1))); | ||
} | } | ||
| + | |||
| + | </source> | ||
Latest revision as of 13:45, 28 October 2013
void BestAlignment::MaqAlignment(MaqIndex & index1, MaqIndex & index2, MaqIndex & index3, FastQ & read)
{
position = -1;
sumQ = mapQ = 0;
mismatches = read.Length();
forward = true;
for (unsigned int i = 0; i < read.quality.Length(); i++)
sumQ += read.quality[i];
mapQdenominator = 0;
ClearPositions();
EvaluateAlignment(index1, read, 0);
EvaluateAlignment(index1, read, index1.bases);
EvaluateAlignment(index1, read, index1.bases / 2);
EvaluateAlignment(index2, read, 0);
EvaluateAlignment(index2, read, index2.bases / 2);
EvaluateAlignment(index3, read, 0);
read.Reverse();
int offset = read.Length() - 24;
ClearPositions();
EvaluateAlignment(index1, read, offset + 0, false);
EvaluateAlignment(index1, read, offset + index1.bases, false);
EvaluateAlignment(index1, read, offset + index1.bases / 2, false);
EvaluateAlignment(index2, read, offset + 0, false);
EvaluateAlignment(index2, read, offset + index2.bases / 2, false);
EvaluateAlignment(index3, read, offset + 0, false);
if (forward) read.Reverse();
double quality = pow(10.0, -0.1 * sumQ) / (mapQdenominator + 1e-30);
mapQ = quality < 0.5 ? 0 : (quality > 0.9999999998 ? 99 : (int) (0.5 + 10 * log(1.0 - quality) / log(0.1)));
}
