File: scoring.cpp

package info (click to toggle)
centrifuge 1.0.3-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 11,864 kB
  • sloc: cpp: 51,936; perl: 1,919; python: 1,538; makefile: 618; sh: 352
file content (286 lines) | stat: -rw-r--r-- 9,495 bytes parent folder | download | duplicates (4)
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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
/*
 * Copyright 2011, Ben Langmead <langmea@cs.jhu.edu>
 *
 * This file is part of Bowtie 2.
 *
 * Bowtie 2 is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Bowtie 2 is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Bowtie 2.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <iostream>
#include "scoring.h"

using namespace std;

/**
 * Return true iff a read of length 'rdlen' passes the score filter, i.e.,
 * has enough characters to rise above the minimum score threshold.
 */
bool Scoring::scoreFilter(
	int64_t minsc,
	size_t rdlen) const
{
	int64_t sc = (int64_t)(rdlen * match(30));
	return sc >= minsc;
}

/**
 * Given the score floor for valid alignments and the length of the read,
 * calculate the maximum possible number of read gaps that could occur in a
 * valid alignment.
 */
int Scoring::maxReadGaps(
	int64_t minsc,
	size_t rdlen) const
{
	// Score if all characters match.  TODO: remove assumption that match bonus
	// is independent of quality value.
	int64_t sc = (int64_t)(rdlen * match(30));
	assert_geq(sc, minsc);
	// Now convert matches to read gaps until sc calls below minsc
	bool first = true;
	int num = 0;
	while(sc >= minsc) {
		if(first) {
			first = false;
			// Subtract both penalties
			sc -= readGapOpen();
		} else {
			// Subtract just the extension penalty
			sc -= readGapExtend();
		}
		num++;
	}
	assert_gt(num, 0);
	return num-1;
}

/**
 * Given the score floor for valid alignments and the length of the read,
 * calculate the maximum possible number of reference gaps that could occur
 * in a valid alignment.
 */
int Scoring::maxRefGaps(
	int64_t minsc,
	size_t rdlen) const
{
	// Score if all characters match.  TODO: remove assumption that match bonus
	// is independent of quality value.
	int64_t sc = (int64_t)(rdlen * match(30));
	assert_geq(sc, minsc);
	// Now convert matches to read gaps until sc calls below minsc
	bool first = true;
	int num = 0;
	while(sc >= minsc) {
		sc -= match(30);
		if(first) {
			first = false;
			// Subtract both penalties
			sc -= refGapOpen();
		} else {
			// Subtract just the extension penalty
			sc -= refGapExtend();
		}
		num++;
	}
	assert_gt(num, 0);
	return num-1;
}

/**
 * Given a read sequence, return true iff the read passes the N filter.
 * The N filter rejects reads with more than the number of Ns.
 */
bool Scoring::nFilter(const BTDnaString& rd, size_t& ns) const {
	size_t rdlen = rd.length();
	size_t maxns = nCeil.f<size_t>((double)rdlen);
	assert_geq(rd.length(), 0);
	for(size_t i = 0; i < rdlen; i++) {
		if(rd[i] == 4) {
			ns++;
			if(ns > maxns) {
				return false; // doesn't pass
			}
		}
	}
	return true; // passes
}

/**
 * Given a read sequence, return true iff the read passes the N filter.
 * The N filter rejects reads with more than the number of Ns.
 *
 * For paired-end reads, there is a	question of how to apply the filter.
 * The filter could be applied to both mates separately, which might then
 * prevent paired-end alignment.  Or the filter could be applied to the
 * reads as though they're concatenated together.  The latter approach has
 * pros and cons.  The pro is that we can use paired-end information to
 * recover alignments for mates that would not have passed the N filter on
 * their own.  The con is that we might not want to do that, since the
 * non-N portion of the bad mate might contain particularly unreliable
 * information.
 */
void Scoring::nFilterPair(
	const BTDnaString* rd1, // mate 1
	const BTDnaString* rd2, // mate 2
	size_t& ns1,            // # Ns in mate 1
	size_t& ns2,            // # Ns in mate 2
	bool& filt1,            // true -> mate 1 rejected by filter
	bool& filt2)            // true -> mate 2 rejected by filter
	const
{
	// Both fail to pass by default
	filt1 = filt2 = false;
	if(rd1 != NULL && rd2 != NULL && ncatpair) {
		size_t rdlen1 = rd1->length();
		size_t rdlen2 = rd2->length();
		size_t maxns = nCeil.f<size_t>((double)(rdlen1 + rdlen2));
		for(size_t i = 0; i < rdlen1; i++) {
			if((*rd1)[i] == 4) ns1++;
			if(ns1 > maxns) {
				// doesn't pass
				return;
			}
		}
		for(size_t i = 0; i < rdlen2; i++) {
			if((*rd2)[i] == 4) ns2++;
			if(ns2 > maxns) {
				// doesn't pass
				return;
			}
		}
		// Both pass
		filt1 = filt2 = true;
	} else {
		if(rd1 != NULL) filt1 = nFilter(*rd1, ns1);
		if(rd2 != NULL) filt2 = nFilter(*rd2, ns2);
	}
}

#ifdef SCORING_MAIN

int main() {
	{
		cout << "Case 1: Simple 1 ... ";
		Scoring sc = Scoring::base1();
		assert_eq(COST_MODEL_CONSTANT, sc.matchType);
		
		assert_eq(0, sc.maxRefGaps(0, 10));  // 10 - 1 - 15 = -6
		assert_eq(0, sc.maxRefGaps(0, 11));  // 11 - 1 - 15 = -5
		assert_eq(0, sc.maxRefGaps(0, 12));  // 12 - 1 - 15 = -4
		assert_eq(0, sc.maxRefGaps(0, 13));  // 13 - 1 - 15 = -3
		assert_eq(0, sc.maxRefGaps(0, 14));  // 14 - 1 - 15 = -2
		assert_eq(0, sc.maxRefGaps(0, 15));  // 15 - 1 - 15 = -1
		assert_eq(1, sc.maxRefGaps(0, 16));  // 16 - 1 - 15 =  0
		assert_eq(1, sc.maxRefGaps(0, 17));  // 17 - 2 - 19 = -4
		assert_eq(1, sc.maxRefGaps(0, 18));  // 18 - 2 - 19 = -3
		assert_eq(1, sc.maxRefGaps(0, 19));  // 19 - 2 - 19 = -2
		assert_eq(1, sc.maxRefGaps(0, 20));  // 20 - 2 - 19 = -1
		assert_eq(2, sc.maxRefGaps(0, 21));  // 21 - 2 - 19 =  0
		
		assert_eq(0, sc.maxReadGaps(0, 10));   // 10 - 0 - 15 = -5
		assert_eq(0, sc.maxReadGaps(0, 11));   // 11 - 0 - 15 = -4
		assert_eq(0, sc.maxReadGaps(0, 12));   // 12 - 0 - 15 = -3
		assert_eq(0, sc.maxReadGaps(0, 13));   // 13 - 0 - 15 = -2
		assert_eq(0, sc.maxReadGaps(0, 14));   // 14 - 0 - 15 = -1
		assert_eq(1, sc.maxReadGaps(0, 15));   // 15 - 0 - 15 =  0
		assert_eq(1, sc.maxReadGaps(0, 16));   // 16 - 0 - 19 = -3
		assert_eq(1, sc.maxReadGaps(0, 17));   // 17 - 0 - 19 = -2
		assert_eq(1, sc.maxReadGaps(0, 18));   // 18 - 0 - 19 = -1
		assert_eq(2, sc.maxReadGaps(0, 19));   // 19 - 0 - 19 =  0
		assert_eq(2, sc.maxReadGaps(0, 20));   // 20 - 0 - 23 = -3
		assert_eq(2, sc.maxReadGaps(0, 21));   // 21 - 0 - 23 = -2
		
		// N ceiling: const=2, linear=0.1
		assert_eq(1, sc.nCeil(1));
		assert_eq(2, sc.nCeil(3));
		assert_eq(2, sc.nCeil(5));
		assert_eq(2, sc.nCeil(7));
		assert_eq(2, sc.nCeil(9));
		assert_eq(3, sc.nCeil(10));
		for(int i = 0; i < 30; i++) {
			assert_eq(3, sc.n(i));
			assert_eq(3, sc.mm(i));
		}
		assert_eq(5, sc.gapbar);
		cout << "PASSED" << endl;
	}
	{
		cout << "Case 2: Simple 2 ... ";
		Scoring sc(
			4,               // reward for a match
			COST_MODEL_QUAL, // how to penalize mismatches
			0,               // constant if mm pelanty is a constant
			30,              // penalty for nuc mm in decoded colorspace als
			-3.0f,           // constant coeff for minimum score
			-3.0f,           // linear coeff for minimum score
			DEFAULT_FLOOR_CONST,  // constant coeff for score floor
			DEFAULT_FLOOR_LINEAR, // linear coeff for score floor
			3.0f,            // max # ref Ns allowed in alignment; const coeff
			0.4f,            // max # ref Ns allowed in alignment; linear coeff
			COST_MODEL_QUAL, // how to penalize Ns in the read
			0,               // constant if N pelanty is a constant
			true,            // whether to concatenate mates before N filtering
			25,              // constant coeff for cost of gap in the read
			25,              // constant coeff for cost of gap in the ref
			10,              // coeff of linear term for cost of gap in read
			10,              // coeff of linear term for cost of gap in ref
			5,               // 5 rows @ top/bot diagonal-entrance-only
			-1,              // no restriction on row
			false            // score prioritized over row
		);

		assert_eq(COST_MODEL_CONSTANT, sc.matchType);
		assert_eq(4, sc.matchConst);
		assert_eq(COST_MODEL_QUAL, sc.mmcostType);
		assert_eq(COST_MODEL_QUAL, sc.npenType);
		
		assert_eq(0, sc.maxRefGaps(0, 8));  // 32 - 4 - 35 = -7
		assert_eq(0, sc.maxRefGaps(0, 9));  // 36 - 4 - 35 = -3
		assert_eq(1, sc.maxRefGaps(0, 10)); // 40 - 4 - 35 =  1
		assert_eq(1, sc.maxRefGaps(0, 11)); // 44 - 8 - 45 = -9
		assert_eq(1, sc.maxRefGaps(0, 12)); // 48 - 8 - 45 = -5
		assert_eq(1, sc.maxRefGaps(0, 13)); // 52 - 8 - 45 = -1
		assert_eq(2, sc.maxRefGaps(0, 14)); // 56 - 8 - 45 =  3
		
		assert_eq(0, sc.maxReadGaps(0, 8));   // 32 - 0 - 35 = -3
		assert_eq(1, sc.maxReadGaps(0, 9));   // 36 - 0 - 35 =  1
		assert_eq(1, sc.maxReadGaps(0, 10));  // 40 - 0 - 45 = -5
		assert_eq(1, sc.maxReadGaps(0, 11));  // 44 - 0 - 45 = -1
		assert_eq(2, sc.maxReadGaps(0, 12));  // 48 - 0 - 45 =  3
		assert_eq(2, sc.maxReadGaps(0, 13));  // 52 - 0 - 55 = -3
		assert_eq(3, sc.maxReadGaps(0, 14));  // 56 - 0 - 55 =  1

		// N ceiling: const=3, linear=0.4
		assert_eq(1, sc.nCeil(1));
		assert_eq(2, sc.nCeil(2));
		assert_eq(3, sc.nCeil(3));
		assert_eq(4, sc.nCeil(4));
		assert_eq(5, sc.nCeil(5));
		assert_eq(5, sc.nCeil(6));
		assert_eq(5, sc.nCeil(7));
		assert_eq(6, sc.nCeil(8));
		assert_eq(6, sc.nCeil(9));

		for(int i = 0; i < 256; i++) {
			assert_eq(i, sc.n(i));
			assert_eq(i, sc.mm(i));
		}

		assert_eq(5, sc.gapbar);

		cout << "PASSED" << endl;
	}
}

#endif /*def SCORING_MAIN*/