Function find

Map/map.cc:346–475

Return the mapping of the specified sequence. */

/** Return the mapping of the specified sequence. */
static void find(const FastaIndex& faIndex, const FMIndex& fmIndex,
		const FastqRecord& rec)
{
	if (rec.seq.empty()) {
		cerr << PROGRAM ": error: "
			"the sequence `" << rec.id << "' is empty\n";
		exit(EXIT_FAILURE);
	}

	Match m, rcm;
	tie(m, rcm) = findMatch(fmIndex, rec.seq);

	if (opt::dup) {
		printDuplicates(m, rcm, faIndex, fmIndex, rec);
		return;
	}

	bool rc;
	if (opt::ss) {
		rc = rec.id.size() > 2
			&& rec.id.substr(rec.id.size()-2) == "/1";
		bool prc = rcm.qspan() > m.qspan();
		if (prc != rc && ((rc && rcm.size() > 0)
					|| (!rc && m.size() > 0)))
#pragma omp atomic
			g_count.suboptimal++;
		if (prc != rc && ((rc && rcm.size() == 0 && m.size() > 0)
				|| (!rc && m.size() == 0 && rcm.size() > 0)))
#pragma omp atomic
			g_count.subunmapped++;
	} else {
		rc = rcm.qspan() > m.qspan();

		// if both matches are the same length, sum up the number of times
		// each were seen.
		if (rcm.qspan() == m.qspan())
			rc ? rcm.num += m.num : m.num += rcm.num;
	}

	vector<string> alts;
	Match mm = rc ? rcm : m;
	string mseq = rc ? reverseComplement(rec.seq) : rec.seq;

#if SAM_SEQ_QUAL
	if (opt::multi) {
		if (mm.qstart > 0) {
			string seq = mseq.substr(0, mm.qstart);
			Match m1, rcm1;
			tie(m1, rcm1) = findMatch(fmIndex, seq);
			bool rc1 = rcm1.qspan() > m1.qspan();
			string xa = toXA(faIndex, fmIndex, rc1 ? rcm1 : m1,
					rc ^ rc1, mseq.size(), 0);
			if (xa != "")
				alts.push_back(xa);
		}

		if (mm.qend < mseq.size()) {
			string seq =