830 - BestAlignment::MaqAlignment()

From Genome Analysis Wiki
Jump to navigationJump to search
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)));
}