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...")
(No difference)

Revision as of 13:44, 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))); }