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