1 /* $Id: scoreblk_unit_test.cpp 462040 2015-03-16 13:14:29Z vasilche $
2 * ===========================================================================
3 *
4 * PUBLIC DOMAIN NOTICE
5 * National Center for Biotechnology Information
6 *
7 * This software/database is a "United States Government Work" under the
8 * terms of the United States Copyright Act. It was written as part of
9 * the author's official duties as a United States Government employee and
10 * thus cannot be copyrighted. This software/database is freely available
11 * to the public for use. The National Library of Medicine and the U.S.
12 * Government have not placed any restriction on its use or reproduction.
13 *
14 * Although all reasonable efforts have been taken to ensure the accuracy
15 * and reliability of the software and data, the NLM and the U.S.
16 * Government do not and cannot warrant the performance or results that
17 * may be obtained by using this software or data. The NLM and the U.S.
18 * Government disclaim all warranties, express or implied, including
19 * warranties of performance, merchantability or fitness for any particular
20 * purpose.
21 *
22 * Please cite the author in any work or product based on this material.
23 *
24 * ===========================================================================
25 *
26 * Author: Tom Madden
27 *
28 * File Description:
29 * Unit test module for ScoreBlk related functions.
30 *
31 * ===========================================================================
32 */
33 #include <ncbi_pch.hpp>
34 #include <corelib/test_boost.hpp>
35
36 #include <corelib/metareg.hpp>
37 #include <objmgr/util/sequence.hpp>
38 #include <objmgr/bioseq_handle.hpp>
39 #include <objmgr/seq_vector.hpp>
40
41 #include <algo/blast/api/blast_types.hpp>
42 #include <algo/blast/api/blast_aux.hpp>
43 #include <algo/blast/api/blast_exception.hpp>
44 #include <algo/blast/api/blast_options_handle.hpp>
45 #include <blast_objmgr_priv.hpp>
46
47 #include <algo/blast/core/ncbi_math.h>
48 #include <algo/blast/core/blast_setup.h>
49 #include <algo/blast/core/blast_stat.h>
50 #include <algo/blast/core/blast_encoding.h>
51
52 #include "test_objmgr.hpp"
53
54 #include <string>
55 #include <vector>
56
57 using namespace std;
58 using namespace ncbi;
59 using namespace ncbi::objects;
60 using namespace ncbi::blast;
61
62 BOOST_AUTO_TEST_SUITE(scoreblk)
63
BOOST_AUTO_TEST_CASE(GetScoreBlockNucl)64 BOOST_AUTO_TEST_CASE(GetScoreBlockNucl) {
65 const EBlastProgramType kProgram = eBlastTypeBlastn;
66 CSeq_id id("gi|1945388");
67 auto_ptr<SSeqLoc> qsl(CTestObjMgr::Instance().CreateSSeqLoc(id, eNa_strand_both));
68 TSeqLocVector query_v;
69 query_v.push_back(*qsl);
70 CBlastQueryInfo query_info;
71 CBLAST_SequenceBlk query_blk;
72 TSearchMessages blast_msg;
73 CRef<CBlastOptionsHandle> opts(CBlastOptionsFactory::Create(eBlastn));
74
75 const CBlastOptions& kOpts = opts->GetOptions();
76 EBlastProgramType prog = kOpts.GetProgramType();
77 ENa_strand strand_opt = kOpts.GetStrandOption();
78
79 SetupQueryInfo(query_v, prog, strand_opt, &query_info);
80 SetupQueries(query_v, query_info, &query_blk,
81 prog, strand_opt, blast_msg);
82 ITERATE(TSearchMessages, m, blast_msg) {
83 BOOST_REQUIRE(m->empty());
84 }
85
86 CBlastScoringOptions scoring_opts;
87 Int2 rv = BlastScoringOptionsNew(kProgram, &scoring_opts);
88 BOOST_REQUIRE(rv == 0);
89 BlastScoreBlk* sbp;
90 Blast_Message* blast_message = NULL;
91 Int2 status = BlastSetup_ScoreBlkInit(query_blk, query_info, scoring_opts,
92 kProgram, &sbp, 1.0, &blast_message, &BlastFindMatrixPath);
93 BOOST_REQUIRE(status == 0);
94
95 BOOST_REQUIRE_EQUAL(0, (int) sbp->protein_alphabet);
96 BOOST_REQUIRE_EQUAL(99, (int) sbp->alphabet_code);
97 BOOST_REQUIRE_EQUAL(16, (int) sbp->alphabet_size);
98 BOOST_REQUIRE_EQUAL(0, (int) sbp->alphabet_start);
99 BOOST_REQUIRE_EQUAL(-3, (int) sbp->loscore);
100 BOOST_REQUIRE_EQUAL(1, (int) sbp->hiscore);
101 BOOST_REQUIRE_EQUAL(-3, (int) sbp->penalty);
102 BOOST_REQUIRE_EQUAL(1, (int) sbp->reward);
103
104 sbp = BlastScoreBlkFree(sbp);
105 BOOST_REQUIRE(sbp == NULL);
106
107 }
108
BOOST_AUTO_TEST_CASE(GetScoreBlockProtein)109 BOOST_AUTO_TEST_CASE(GetScoreBlockProtein) {
110 const EBlastProgramType kProgram = eBlastTypeBlastp;
111 CSeq_id id("gi|3091");
112 auto_ptr<SSeqLoc> qsl(CTestObjMgr::Instance().CreateSSeqLoc(id));
113 TSeqLocVector query_v;
114 query_v.push_back(*qsl);
115 CBlastQueryInfo query_info;
116 CBLAST_SequenceBlk query_blk;
117 TSearchMessages blast_msg;
118 CRef<CBlastOptionsHandle> opts(CBlastOptionsFactory::Create(eBlastp));
119
120 const CBlastOptions& kOpts = opts->GetOptions();
121 EBlastProgramType prog = kOpts.GetProgramType();
122 ENa_strand strand_opt = kOpts.GetStrandOption();
123
124 SetupQueryInfo(query_v, prog, strand_opt, &query_info);
125 SetupQueries(query_v, query_info, &query_blk,
126 prog, strand_opt, blast_msg);
127 ITERATE(TSearchMessages, m, blast_msg) {
128 BOOST_REQUIRE(m->empty());
129 }
130
131 CBlastScoringOptions scoring_opts;
132 Int2 rv = BlastScoringOptionsNew(kProgram, &scoring_opts);
133 BOOST_REQUIRE(rv == 0);
134 BlastScoreBlk* sbp;
135 CBlast_Message error_msg;
136 Int2 status = BlastSetup_ScoreBlkInit(query_blk, query_info, scoring_opts,
137 kProgram, &sbp, 1.0, &error_msg, &BlastFindMatrixPath);
138 BOOST_REQUIRE(status == 0);
139
140 BOOST_REQUIRE_EQUAL(1, (int) sbp->protein_alphabet);
141 BOOST_REQUIRE_EQUAL(11, (int) sbp->alphabet_code);
142 BOOST_REQUIRE_EQUAL(BLASTAA_SIZE, (int) sbp->alphabet_size);
143 BOOST_REQUIRE_EQUAL(0, (int) sbp->alphabet_start);
144 BOOST_REQUIRE_EQUAL(-4, (int) sbp->loscore);
145 BOOST_REQUIRE_EQUAL(11, (int) sbp->hiscore);
146
147 sbp = BlastScoreBlkFree(sbp);
148 BOOST_REQUIRE(sbp == NULL);
149 }
150
BOOST_AUTO_TEST_CASE(GetScoreBlockPHI)151 BOOST_AUTO_TEST_CASE(GetScoreBlockPHI) {
152 const EBlastProgramType kProgram = eBlastTypePhiBlastp;
153 const string kPhiPattern("EVALNAEGWQSSG");
154 const string kMatrix("BLOSUM45");
155 const int kGapOpenBad = 16; // Unsupported value
156 const int kGapOpenGood = 14;// Supported value
157 const int kGapExtend = 2;
158 const string kErrorMsg("The combination 16 for gap opening cost");
159 const double kPhiLambda = 0.199;
160 const double kPhiK = 0.040;
161
162 CSeq_id id("gi|3091");
163 auto_ptr<SSeqLoc> qsl(CTestObjMgr::Instance().CreateSSeqLoc(id));
164 TSeqLocVector query_v;
165 query_v.push_back(*qsl);
166 CBlastQueryInfo query_info;
167 CBLAST_SequenceBlk query_blk;
168 TSearchMessages blast_msg;
169
170 CRef<CBlastOptionsHandle> opts(CBlastOptionsFactory::Create(eBlastp));
171 opts->SetOptions().SetPHIPattern(kPhiPattern.c_str(), false);
172
173 const CBlastOptions& kOpts = opts->GetOptions();
174 EBlastProgramType prog = kOpts.GetProgramType();
175 ENa_strand strand_opt = kOpts.GetStrandOption();
176
177 SetupQueryInfo(query_v, prog, strand_opt, &query_info);
178 SetupQueries(query_v, query_info, &query_blk,
179 prog, strand_opt, blast_msg);
180 ITERATE(TSearchMessages, m, blast_msg) {
181 BOOST_REQUIRE(m->empty());
182 }
183
184 CBlastScoringOptions scoring_opts;
185 Int2 rv = BlastScoringOptionsNew(kProgram, &scoring_opts);
186 BOOST_REQUIRE(rv == 0);
187 sfree(scoring_opts->matrix);
188 scoring_opts->matrix = strdup(kMatrix.c_str());
189 scoring_opts->gap_open = kGapOpenBad;
190 scoring_opts->gap_extend = kGapExtend;
191
192 BlastScoreBlk* sbp;
193 Blast_Message* blast_message=NULL;
194
195 Int2 status =
196 BlastSetup_ScoreBlkInit(query_blk, query_info,
197 scoring_opts,
198 kProgram, &sbp, 1.0,
199 &blast_message,
200 &BlastFindMatrixPath);
201
202 BOOST_REQUIRE_EQUAL(-1, (int)status);
203 sbp = BlastScoreBlkFree(sbp);
204 BOOST_REQUIRE(sbp == NULL);
205
206 BOOST_REQUIRE(!strncmp(kErrorMsg.c_str(), blast_message->message,
207 kErrorMsg.size()));
208 // blast_message will be reused in the next call, so this message
209 // must be freed.
210 Blast_MessageFree(blast_message);
211
212 scoring_opts->gap_open = kGapOpenGood;
213 status =
214 BlastSetup_ScoreBlkInit(query_blk, query_info,
215 scoring_opts,
216 kProgram, &sbp, 1.0,
217 &blast_message,
218 &BlastFindMatrixPath);
219 BOOST_REQUIRE_EQUAL(0, (int) status);
220 BOOST_REQUIRE(sbp->kbp_std != sbp->kbp_gap_std);
221 BOOST_REQUIRE(sbp->kbp == sbp->kbp_std);
222 BOOST_REQUIRE(sbp->kbp_gap == sbp->kbp_gap_std);
223 BOOST_REQUIRE(sbp->kbp_gap[0]->Lambda == kPhiLambda);
224 BOOST_REQUIRE(sbp->kbp_gap[0]->K == kPhiK);
225 BOOST_REQUIRE(sbp->kbp_gap[0]->H > 0);
226 BOOST_REQUIRE(sbp->kbp[0]->Lambda == kPhiLambda);
227 BOOST_REQUIRE(sbp->kbp[0]->K == kPhiK);
228 BOOST_REQUIRE(sbp->kbp[0]->H > 0);
229
230 sbp = BlastScoreBlkFree(sbp);
231 BOOST_REQUIRE(sbp == NULL);
232 }
233
BOOST_AUTO_TEST_CASE(GetScoreBlockForFullyMaskedProtein)234 BOOST_AUTO_TEST_CASE(GetScoreBlockForFullyMaskedProtein) {
235 const EBlastProgramType kProgram = eBlastTypeBlastp;
236 CSeq_id id("gi|3091");
237 const int start = 0;
238 const int stop = 27;
239 pair<TSeqPos, TSeqPos> range(start, stop);
240 auto_ptr<SSeqLoc> qsl(CTestObjMgr::Instance().CreateSSeqLoc(id, range));
241 TSeqLocVector query_v;
242 query_v.push_back(*qsl);
243 CBlastQueryInfo query_info;
244 CBLAST_SequenceBlk query_blk;
245 TSearchMessages blast_msg;
246 CRef<CBlastOptionsHandle> opts(CBlastOptionsFactory::Create(eBlastp));
247
248 const CBlastOptions& kOpts = opts->GetOptions();
249 EBlastProgramType prog = kOpts.GetProgramType();
250 ENa_strand strand_opt = kOpts.GetStrandOption();
251
252 SetupQueryInfo(query_v, prog, strand_opt, &query_info);
253 SetupQueries(query_v, query_info, &query_blk,
254 prog, strand_opt, blast_msg);
255 ITERATE(TSearchMessages, m, blast_msg) {
256 BOOST_REQUIRE(m->empty());
257 }
258
259 CBlastScoringOptions scoring_opts;
260 Int2 rv = BlastScoringOptionsNew(kProgram, &scoring_opts);
261 BOOST_REQUIRE(rv == 0);
262
263 BlastSeqLoc *loc = BlastSeqLocNew(NULL, start, stop);
264 BlastMaskLoc* filter_maskloc = BlastMaskLocNew(1);
265 filter_maskloc->seqloc_array[0] = loc;
266
267 BlastSetUp_MaskQuery(query_blk, query_info, filter_maskloc, kProgram);
268 filter_maskloc = BlastMaskLocFree(filter_maskloc);
269 BOOST_REQUIRE(filter_maskloc == NULL);
270
271
272 BlastScoreBlk* sbp;
273 Blast_Message* blast_message = NULL;
274 Int2 status = BlastSetup_ScoreBlkInit(query_blk, query_info,
275 scoring_opts, kProgram, &sbp,
276 1.0, &blast_message,
277 &BlastFindMatrixPath);
278
279 // Note that errors will come in the Blast_Message structure
280 BOOST_REQUIRE_EQUAL(0, (int) status);
281 BOOST_REQUIRE(blast_message != NULL);
282
283 BOOST_REQUIRE_EQUAL(1, (int) sbp->protein_alphabet);
284 BOOST_REQUIRE_EQUAL(11, (int) sbp->alphabet_code);
285 BOOST_REQUIRE_EQUAL(BLASTAA_SIZE, (int) sbp->alphabet_size);
286 BOOST_REQUIRE_EQUAL(0, (int) sbp->alphabet_start);
287 BOOST_REQUIRE_EQUAL(-4, (int) sbp->loscore);
288 BOOST_REQUIRE_EQUAL(11, (int) sbp->hiscore);
289
290 blast_message = Blast_MessageFree(blast_message);
291 BOOST_REQUIRE(blast_message == NULL);
292 sbp = BlastScoreBlkFree(sbp);
293 BOOST_REQUIRE(sbp == NULL);
294 }
295
BOOST_AUTO_TEST_CASE(BlastResFreqStdCompProteinTest)296 BOOST_AUTO_TEST_CASE(BlastResFreqStdCompProteinTest) {
297
298 BlastScoreBlk sbp;
299 sbp.alphabet_code = 11;
300 sbp.protein_alphabet = TRUE;
301 sbp.alphabet_start = 0;
302 sbp.alphabet_size = 26;
303
304 Blast_ResFreq* stdrfp = Blast_ResFreqNew(&sbp);
305 Blast_ResFreqStdComp(&sbp, stdrfp);
306
307 BOOST_REQUIRE_EQUAL(11, (int) stdrfp->alphabet_code);
308 // All frequencies multiplied by 100000 and truncated to nearest integer
309 // to avoid rounding error of floating points.
310 // some ambiguity codes, should be zero.
311 BOOST_REQUIRE_EQUAL(0, (int)BLAST_Nint(100000*stdrfp->prob[2])); // Asp or Asn
312 BOOST_REQUIRE_EQUAL(0, (int)BLAST_Nint(100000*stdrfp->prob[21])); // X
313 // some "real" residues we use.
314 BOOST_REQUIRE_EQUAL(3856, (int)BLAST_Nint(100000*stdrfp->prob[6]));
315 BOOST_REQUIRE_EQUAL(2243, (int)BLAST_Nint(100000*stdrfp->prob[12]));
316 BOOST_REQUIRE_EQUAL(3216, (int)BLAST_Nint(100000*stdrfp->prob[22]));
317
318 stdrfp = Blast_ResFreqFree(stdrfp);
319 BOOST_REQUIRE(stdrfp == NULL);
320 }
321
BOOST_AUTO_TEST_CASE(BlastResFreqStdCompNucleotideTest)322 BOOST_AUTO_TEST_CASE(BlastResFreqStdCompNucleotideTest) {
323
324 BlastScoreBlk sbp;
325 sbp.alphabet_code = 99;
326 sbp.protein_alphabet = FALSE;
327 sbp.alphabet_start = 0;
328 sbp.alphabet_size = 16;
329
330 Blast_ResFreq* stdrfp = Blast_ResFreqNew(&sbp);
331 Blast_ResFreqStdComp(&sbp, stdrfp);
332
333 BOOST_REQUIRE_EQUAL(99, (int) stdrfp->alphabet_code);
334 const int num_real_bases = 4;
335 // Multiplied by 100 and truncated to avoid rounding error of floating point.
336 for (int index=0; index<num_real_bases; index++) // A,C,G,T
337 BOOST_REQUIRE_EQUAL(25, (int)(100*stdrfp->prob[index]));
338
339 for (int index=num_real_bases; index<sbp.alphabet_size; index++) // ambiguity codes, are all zero.
340 BOOST_REQUIRE_EQUAL(0, (int)(100*stdrfp->prob[index]));
341
342 stdrfp = Blast_ResFreqFree(stdrfp);
343 BOOST_REQUIRE(stdrfp == NULL);
344 }
345
BOOST_AUTO_TEST_CASE(EqualRewardPenaltyLHtoK)346 BOOST_AUTO_TEST_CASE(EqualRewardPenaltyLHtoK)
347 {
348 const EBlastProgramType kProgram = eBlastTypeBlastn;
349 BlastScoringOptions* score_opts = NULL;
350 BlastScoringOptionsNew(kProgram, &score_opts);
351 score_opts->reward = 2;
352 score_opts->penalty = -2;
353
354 BlastScoreBlk* sbp = BlastScoreBlkNew(BLASTNA_SEQ_CODE, 1);
355 Blast_ScoreBlkMatrixInit(kProgram, score_opts, sbp, NULL);
356 Blast_ScoreBlkKbpIdealCalc(sbp);
357
358 BOOST_REQUIRE(fabs(sbp->kbp_ideal->K - 1.0/3) < 1e-6);
359 BlastScoreBlkFree(sbp);
360 BlastScoringOptionsFree(score_opts);
361 }
362
BOOST_AUTO_TEST_CASE(NuclGappedCalc)363 BOOST_AUTO_TEST_CASE(NuclGappedCalc)
364 {
365 const EBlastProgramType kProgram = eBlastTypeBlastn;
366 CBlastScoringOptions score_opts;
367 BlastScoringOptionsNew(kProgram, &score_opts);
368 score_opts->reward = 1;
369 score_opts->penalty = -2;
370 score_opts->gap_open = 3;
371 score_opts->gap_extend = 1;
372
373 CBlastScoreBlk sbp(BlastScoreBlkNew(BLASTNA_SEQ_CODE, 1));
374 Blast_ScoreBlkMatrixInit(kProgram, score_opts, sbp, NULL);
375 Blast_ScoreBlkKbpIdealCalc(sbp);
376
377 Blast_KarlinBlk* kbp = Blast_KarlinBlkNew();
378 CBlast_Message error_msg;
379 Int2 status = 0;
380 status =
381 Blast_KarlinBlkNuclGappedCalc(kbp, score_opts->gap_open,
382 score_opts->gap_extend, score_opts->reward,
383 score_opts->penalty, sbp->kbp_ideal,
384 &(sbp->round_down), &error_msg);
385
386 BOOST_REQUIRE_EQUAL(0, (int) status);
387 BOOST_REQUIRE_EQUAL(false, (bool) sbp->round_down);
388 BOOST_REQUIRE_CLOSE(1.32, kbp->Lambda, 0.001);
389 BOOST_REQUIRE_CLOSE(0.57, kbp->K, 0.001);
390 BOOST_REQUIRE_CLOSE(-0.562, kbp->logK, 0.1);
391 BOOST_REQUIRE(error_msg.Get() == NULL);
392 error_msg.Reset();
393
394 // Check values of alpha and beta parameters
395 double alpha, beta;
396 Blast_GetNuclAlphaBeta(score_opts->reward, score_opts->penalty,
397 score_opts->gap_open, score_opts->gap_extend,
398 sbp->kbp_ideal, TRUE, &alpha, &beta);
399 BOOST_REQUIRE_CLOSE(1.3, alpha, 0.001);
400 BOOST_REQUIRE_CLOSE(-1.0, beta, 0.001);
401
402 // Check gap costs for which Karlin-Altschul paramters are copied
403 // from the ungapped block.
404 score_opts->gap_open = 4;
405 score_opts->gap_extend = 2;
406
407 status =
408 Blast_KarlinBlkNuclGappedCalc(kbp, score_opts->gap_open,
409 score_opts->gap_extend, score_opts->reward,
410 score_opts->penalty, sbp->kbp_ideal,
411 &(sbp->round_down), &error_msg);
412 BOOST_REQUIRE_EQUAL(0, (int) status);
413 BOOST_REQUIRE_EQUAL(false, (bool) sbp->round_down);
414 BOOST_REQUIRE_EQUAL(sbp->kbp_ideal->Lambda, kbp->Lambda);
415 BOOST_REQUIRE_EQUAL(sbp->kbp_ideal->K, kbp->K);
416 BOOST_REQUIRE_EQUAL(sbp->kbp_ideal->logK, kbp->logK);
417 BOOST_REQUIRE(error_msg.Get() == NULL);
418 error_msg.Reset();
419
420 Blast_GetNuclAlphaBeta(score_opts->reward, score_opts->penalty,
421 score_opts->gap_open, score_opts->gap_extend,
422 sbp->kbp_ideal, TRUE, &alpha, &beta);
423 BOOST_REQUIRE_CLOSE(sbp->kbp_ideal->Lambda/sbp->kbp_ideal->H,
424 alpha, 1e-10);
425 BOOST_REQUIRE_EQUAL(0.0, beta);
426
427 // Check for scaled up values.
428 score_opts->reward = 10;
429 score_opts->penalty = -20;
430 score_opts->gap_open = 30;
431 score_opts->gap_extend = 10;
432
433 status =
434 Blast_KarlinBlkNuclGappedCalc(kbp, score_opts->gap_open,
435 score_opts->gap_extend, score_opts->reward,
436 score_opts->penalty, sbp->kbp_ideal,
437 &(sbp->round_down), &error_msg);
438
439 BOOST_REQUIRE_EQUAL(0, (int) status);
440 BOOST_REQUIRE_EQUAL(false, (bool) sbp->round_down);
441 BOOST_REQUIRE_CLOSE(0.132, kbp->Lambda, 0.001);
442 BOOST_REQUIRE_CLOSE(0.57, kbp->K, 0.001);
443 BOOST_REQUIRE_CLOSE(-0.562, kbp->logK, 0.1);
444
445 // For this set of values the score needs to be rounded down, mostly checking
446 // here that the round_down bool is set properly
447 score_opts->reward = 2;
448 score_opts->penalty = -7;
449 score_opts->gap_open = 4;
450 score_opts->gap_extend = 2;
451
452 status =
453 Blast_KarlinBlkNuclGappedCalc(kbp, score_opts->gap_open,
454 score_opts->gap_extend, score_opts->reward,
455 score_opts->penalty, sbp->kbp_ideal,
456 &(sbp->round_down), &error_msg);
457 BOOST_REQUIRE_EQUAL(0, (int) status);
458 BOOST_REQUIRE_EQUAL(true, (bool) sbp->round_down);
459 BOOST_REQUIRE_CLOSE(0.675, kbp->Lambda, 0.001);
460 BOOST_REQUIRE_CLOSE(0.62, kbp->K, 0.001);
461 BOOST_REQUIRE_CLOSE(-0.478036, kbp->logK, 0.001);
462 BOOST_REQUIRE(error_msg.Get() == NULL);
463 error_msg.Reset();
464
465 // Check invalid gap costs that were permitted owing to a bug.
466 score_opts->reward = 4;
467 score_opts->penalty = -5;
468 score_opts->gap_open = 3;
469 score_opts->gap_extend = 2;
470
471 status =
472 Blast_KarlinBlkNuclGappedCalc(kbp, score_opts->gap_open,
473 score_opts->gap_extend, score_opts->reward,
474 score_opts->penalty, sbp->kbp_ideal,
475 &(sbp->round_down), &error_msg);
476 BOOST_REQUIRE_EQUAL(1, (int) status);
477 BOOST_REQUIRE(!strncmp("Gap existence and extension values 3 and 2 are not supported for substitution scores 4 and -5",
478 error_msg->message, 90));
479 error_msg.Reset();
480
481 // Check invalid gap costs
482 score_opts->reward = 1;
483 score_opts->penalty = -2;
484 score_opts->gap_open = 1;
485 score_opts->gap_extend = 3;
486
487 status =
488 Blast_KarlinBlkNuclGappedCalc(kbp, score_opts->gap_open,
489 score_opts->gap_extend, score_opts->reward,
490 score_opts->penalty, sbp->kbp_ideal,
491 &(sbp->round_down), &error_msg);
492 BOOST_REQUIRE_EQUAL(1, (int) status);
493 BOOST_REQUIRE(!strncmp("Gap existence and extension values 1 and 3 are not supported for substitution scores 1 and -2",
494 error_msg->message, 90));
495 error_msg.Reset();
496
497 // Alpha and beta can be returned even for unsupported gap costs,
498 // because this function should work for an ungapped search.
499 Blast_GetNuclAlphaBeta(score_opts->reward, score_opts->penalty,
500 score_opts->gap_open, score_opts->gap_extend,
501 sbp->kbp_ideal, TRUE, &alpha, &beta);
502 BOOST_REQUIRE_CLOSE(sbp->kbp_ideal->Lambda/sbp->kbp_ideal->H,
503 alpha, 1e-10);
504 BOOST_REQUIRE_EQUAL(0.0, beta);
505
506 // Check invalid substitution scores
507 score_opts->reward = 2;
508 score_opts->penalty = -1;
509 status =
510 Blast_KarlinBlkNuclGappedCalc(kbp, score_opts->gap_open,
511 score_opts->gap_extend, score_opts->reward,
512 score_opts->penalty, sbp->kbp_ideal,
513 &(sbp->round_down), &error_msg);
514 BOOST_REQUIRE_EQUAL(-1, (int) status);
515 BOOST_REQUIRE(!strcmp("Substitution scores 2 and -1 are not supported",
516 error_msg->message));
517 error_msg.Reset();
518 // Alpha and beta would still be returned in this case as for an
519 // ungapped search.
520 Blast_GetNuclAlphaBeta(score_opts->reward, score_opts->penalty,
521 score_opts->gap_open, score_opts->gap_extend,
522 sbp->kbp_ideal, TRUE, &alpha, &beta);
523 BOOST_REQUIRE_CLOSE(sbp->kbp_ideal->Lambda/sbp->kbp_ideal->H,
524 alpha, 1e-10);
525 BOOST_REQUIRE_EQUAL(0.0, beta);
526
527 // Check alpha and beta for an ungapped search, with a different
528 // pair of substitution scores.
529 score_opts->penalty = -3;
530 Blast_GetNuclAlphaBeta(score_opts->reward, score_opts->penalty,
531 0, 0, sbp->kbp_ideal, FALSE, &alpha, &beta);
532
533 BOOST_REQUIRE_CLOSE(sbp->kbp_ideal->Lambda/sbp->kbp_ideal->H,
534 alpha, 1e-10);
535 BOOST_REQUIRE_EQUAL(-2.0, beta);
536
537 kbp = Blast_KarlinBlkFree(kbp);
538 BOOST_REQUIRE(kbp == NULL);
539 }
540
541 BOOST_AUTO_TEST_SUITE_END()
542