Changes

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..."
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)));
}

Navigation menu