File: search_seeded_phase2.c

package info (click to toggle)
bowtie 1.2.3+dfsg-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 16,888 kB
  • sloc: cpp: 35,784; perl: 5,903; ansic: 1,247; sh: 1,180; python: 487; makefile: 426
file content (99 lines) | stat: -rw-r--r-- 2,976 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
/*
 * This is a fragment, included from multiple places in ebwt_search.cpp.
 * It implements the logic of the second phase of the seeded, quality-
 * aware search routine.  It is implemented as a code fragment so that
 * it can be reused in both the half-index-in-memory and full-index-in-
 * memory situations.
 */
{
	if(!nofw) {
		// If we reach here, then cases 1R, 2R, and 3R have been
		// eliminated.  The next most likely cases are 1F, 2F and
		// 3F...
		params.setFw(true);  // looking at forward strand
		btf2.setReportExacts(false);
		btr2.setReportExacts(false);
		btf2.setQuery(patsrc->bufa());
		// Set up seed bounds
		if(qs < s) {
			btf2.setOffs(0, 0,
						(seedMms > 0)? qs5 : qs,
						(seedMms > 1)? qs5 : qs,
						(seedMms > 2)? qs5 : qs,
						(seedMms > 3)? qs5 : qs);
		} else {
			btf2.setOffs(0, 0,
						(seedMms > 0)? s5 : s,
						(seedMms > 1)? s5 : s,
						(seedMms > 2)? s5 : s,
						(seedMms > 3)? s5 : s);
		}
		// Do a 12/24 backtrack on the forward-strand read using
		// the mirror index.  This will find all case 1F, 2F
		// and 3F hits.
		if(btf2.backtrack()) {
			// The reverse complement hit, so we're done with this
			// read
			DONEMASK_SET(patid);
			continue;
		}

		if(sink->finishedWithStratum(0)) { // no more exact hits are possible
			DONEMASK_SET(patid);
			continue;
		}
	}

	// No need to collect partial alignments if we're not
	// allowing mismatches in the 5' seed half
	if(seedMms == 0) continue;

	if(!norc) {
		// If we reach here, then cases 1F, 2F, 3F, 1R, 2R, and 3R
		// have been eliminated, leaving us with cases 4F and 4R
		// (the cases with 1 mismatch in the 5' half of the seed)
		params.setFw(false);  // looking at reverse-comp strand
		// Set up seed bounds
		if(qs < s) {
			btr2.setOffs(0, 0,
						qs3,
						(seedMms > 1)? qs3 : qs,
						(seedMms > 2)? qs3 : qs,
						(seedMms > 3)? qs3 : qs);
		} else {
			btr2.setOffs(0, 0,
						s3,
						(seedMms > 1)? s3 : s,
						(seedMms > 2)? s3 : s,
						(seedMms > 3)? s3 : s);
		}
		btr2.setQuery(patsrc->bufa());
		btr2.setQlen(s); // just look at the seed
		// Find partial alignments for case 4R
		ASSERT_ONLY(bool done =) btr2.backtrack();
#ifndef NDEBUG
		vector<PartialAlignment> partials;
		assert(pamRc != NULL);
		pamRc->getPartials(patid, partials);
		if(done) assert_gt(partials.size(), 0);
		for(size_t i = 0; i < partials.size(); i++) {
			uint32_t pos0 = partials[i].entry.pos0;
			assert_lt(pos0, s5);
			uint8_t oldChar = (uint8_t)patRcRev[pos0];
			assert_neq(oldChar, partials[i].entry.char0);
			if(partials[i].entry.pos1 != 0xffff) {
				uint32_t pos1 = partials[i].entry.pos1;
				assert_lt(pos1, s5);
				oldChar = (uint8_t)patRcRev[pos1];
				assert_neq(oldChar, partials[i].entry.char1);
				if(partials[i].entry.pos2 != 0xffff) {
					uint32_t pos2 = partials[i].entry.pos2;
					assert_lt(pos2, s5);
					oldChar = (uint8_t)patRcRev[pos2];
					assert_neq(oldChar, partials[i].entry.char2);
				}
			}
		}
#endif
	}
}