/*
* 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 <string>
#include <iostream>
#include <sstream>
#include <limits>
#include "ds.h"
#include "aligner_seed_policy.h"
#include "mem_ids.h"
using namespace std;
static int parseFuncType(const std::string& otype) {
string type = otype;
if(type == "C" || type == "Constant") {
return SIMPLE_FUNC_CONST;
} else if(type == "L" || type == "Linear") {
return SIMPLE_FUNC_LINEAR;
} else if(type == "S" || type == "Sqrt") {
return SIMPLE_FUNC_SQRT;
} else if(type == "G" || type == "Log") {
return SIMPLE_FUNC_LOG;
}
std::cerr << "Error: Bad function type '" << otype.c_str()
<< "'. Should be C (constant), L (linear), "
<< "S (square root) or G (natural log)." << std::endl;
throw 1;
}
#define PARSE_FUNC(fv) { \
if(ctoks.size() >= 1) { \
fv.setType(parseFuncType(ctoks[0])); \
} \
if(ctoks.size() >= 2) { \
double co; \
istringstream tmpss(ctoks[1]); \
tmpss >> co; \
fv.setConst(co); \
} \
if(ctoks.size() >= 3) { \
double ce; \
istringstream tmpss(ctoks[2]); \
tmpss >> ce; \
fv.setCoeff(ce); \
} \
if(ctoks.size() >= 4) { \
double mn; \
istringstream tmpss(ctoks[3]); \
tmpss >> mn; \
fv.setMin(mn); \
} \
if(ctoks.size() >= 5) { \
double mx; \
istringstream tmpss(ctoks[4]); \
tmpss >> mx; \
fv.setMin(mx); \
} \
}
/**
* Parse alignment policy when provided in this format:
* <lab>=<val>;<lab>=<val>;<lab>=<val>...
*
* And label=value possibilities are:
*
* Bonus for a match
* -----------------
*
* MA=xx (default: MA=0, or MA=2 if --local is set)
*
* xx = Each position where equal read and reference characters match up
* in the alignment contriubtes this amount to the total score.
*
* Penalty for a mismatch
* ----------------------
*
* MMP={Cxx|Q|RQ} (default: MMP=C6)
*
* Cxx = Each mismatch costs xx. If MMP=Cxx is specified, quality
* values are ignored when assessing penalities for mismatches.
* Q = Each mismatch incurs a penalty equal to the mismatched base's
* value.
* R = Each mismatch incurs a penalty equal to the mismatched base's
* rounded quality value. Qualities are rounded off to the
* nearest 10, and qualities greater than 30 are rounded to 30.
*
* Penalty for position with N (in either read or reference)
* ---------------------------------------------------------
*
* NP={Cxx|Q|RQ} (default: NP=C1)
*
* Cxx = Each alignment position with an N in either the read or the
* reference costs xx. If NP=Cxx is specified, quality values are
* ignored when assessing penalities for Ns.
* Q = Each alignment position with an N in either the read or the
* reference incurs a penalty equal to the read base's quality
* value.
* R = Each alignment position with an N in either the read or the
* reference incurs a penalty equal to the read base's rounded
* quality value. Qualities are rounded off to the nearest 10,
* and qualities greater than 30 are rounded to 30.
*
* Penalty for a read gap
* ----------------------
*
* RDG=xx,yy (default: RDG=5,3)
*
* xx = Read gap open penalty.
* yy = Read gap extension penalty.
*
* Total cost incurred by a read gap = xx + (yy * gap length)
*
* Penalty for a reference gap
* ---------------------------
*
* RFG=xx,yy (default: RFG=5,3)
*
* xx = Reference gap open penalty.
* yy = Reference gap extension penalty.
*
* Total cost incurred by a reference gap = xx + (yy * gap length)
*
* Minimum score for valid alignment
* ---------------------------------
*
* MIN=xx,yy (defaults: MIN=-0.6,-0.6, or MIN=0.0,0.66 if --local is set)
*
* xx,yy = For a read of length N, the total score must be at least
* xx + (read length * yy) for the alignment to be valid. The
* total score is the sum of all negative penalties (from
* mismatches and gaps) and all positive bonuses. The minimum
* can be negative (and is by default in global alignment mode).
*
* Score floor for local alignment
* -------------------------------
*
* FL=xx,yy (defaults: FL=-Infinity,0.0, or FL=0.0,0.0 if --local is set)
*
* xx,yy = If a cell in the dynamic programming table has a score less
* than xx + (read length * yy), then no valid alignment can go
* through it. Defaults are highly recommended.
*
* N ceiling
* ---------
*
* NCEIL=xx,yy (default: NCEIL=0.0,0.15)
*
* xx,yy = For a read of length N, the number of alignment
* positions with an N in either the read or the
* reference cannot exceed
* ceiling = xx + (read length * yy). If the ceiling is
* exceeded, the alignment is considered invalid.
*
* Seeds
* -----
*
* SEED=mm,len,ival (default: SEED=0,22)
*
* mm = Maximum number of mismatches allowed within a seed.
* Must be >= 0 and <= 2. Note that 2-mismatch mode is
* not fully sensitive; i.e. some 2-mismatch seed
* alignments may be missed.
* len = Length of seed.
* ival = Interval between seeds. If not specified, seed
* interval is determined by IVAL.
*
* Seed interval
* -------------
*
* IVAL={L|S|C},xx,yy (default: IVAL=S,1.0,0.0)
*
* L = let interval between seeds be a linear function of the
* read length. xx and yy are the constant and linear
* coefficients respectively. In other words, the interval
* equals a * len + b, where len is the read length.
* Intervals less than 1 are rounded up to 1.
* S = let interval between seeds be a function of the sqaure
* root of the read length. xx and yy are the
* coefficients. In other words, the interval equals
* a * sqrt(len) + b, where len is the read length.
* Intervals less than 1 are rounded up to 1.
* C = Like S but uses cube root of length instead of square
* root.
*
* Example 1:
*
* SEED=1,10,5 and read sequence is TGCTATCGTACGATCGTAC:
*
* The following seeds are extracted from the forward
* representation of the read and aligned to the reference
* allowing up to 1 mismatch:
*
* Read: TGCTATCGTACGATCGTACA
*
* Seed 1+: TGCTATCGTA
* Seed 2+: TCGTACGATC
* Seed 3+: CGATCGTACA
*
* ...and the following are extracted from the reverse-complement
* representation of the read and align to the reference allowing
* up to 1 mismatch:
*
* Seed 1-: TACGATAGCA
* Seed 2-: GATCGTACGA
* Seed 3-: TGTACGATCG
*
* Example 2:
*
* SEED=1,20,20 and read sequence is TGCTATCGTACGATC. The seed
* length is 20 but the read is only 15 characters long. In this
* case, Bowtie2 automatically shrinks the seed length to be equal
* to the read length.
*
* Read: TGCTATCGTACGATC
*
* Seed 1+: TGCTATCGTACGATC
* Seed 1-: GATCGTACGATAGCA
*
* Example 3:
*
* SEED=1,10,10 and read sequence is TGCTATCGTACGATC. Only one seed
* fits on the read; a second seed would overhang the end of the read
* by 5 positions. In this case, Bowtie2 extracts one seed.
*
* Read: TGCTATCGTACGATC
*
* Seed 1+: TGCTATCGTA
* Seed 1-: TACGATAGCA
*/
void SeedAlignmentPolicy::parseString(
const std::string& s,
bool local,
bool noisyHpolymer,
bool ignoreQuals,
int& bonusMatchType,
int& bonusMatch,
int& penMmcType,
int& penMmcMax,
int& penMmcMin,
int& penNType,
int& penN,
int& penRdExConst,
int& penRfExConst,
int& penRdExLinear,
int& penRfExLinear,
SimpleFunc& costMin,
SimpleFunc& nCeil,
bool& nCatPair,
int& multiseedMms,
int& multiseedLen,
SimpleFunc& multiseedIval,
size_t& failStreak,
size_t& seedRounds,
SimpleFunc* penCanSplice,
SimpleFunc* penNoncanSplice,
SimpleFunc* penIntronLen)
{
bonusMatchType = local ? DEFAULT_MATCH_BONUS_TYPE_LOCAL : DEFAULT_MATCH_BONUS_TYPE;
bonusMatch = local ? DEFAULT_MATCH_BONUS_LOCAL : DEFAULT_MATCH_BONUS;
penMmcType = ignoreQuals ? DEFAULT_MM_PENALTY_TYPE_IGNORE_QUALS :
DEFAULT_MM_PENALTY_TYPE;
penMmcMax = DEFAULT_MM_PENALTY_MAX;
penMmcMin = DEFAULT_MM_PENALTY_MIN;
penNType = DEFAULT_N_PENALTY_TYPE;
penN = DEFAULT_N_PENALTY;
const double DMAX = std::numeric_limits<double>::max();
/*
costMin.init(
local ? SIMPLE_FUNC_LOG : SIMPLE_FUNC_LINEAR,
local ? DEFAULT_MIN_CONST_LOCAL : DEFAULT_MIN_CONST,
local ? DEFAULT_MIN_LINEAR_LOCAL : DEFAULT_MIN_LINEAR);
*/
costMin.init(
local ? SIMPLE_FUNC_LOG : SIMPLE_FUNC_CONST,
local ? DEFAULT_MIN_CONST_LOCAL : -18,
local ? DEFAULT_MIN_LINEAR_LOCAL : 0);
nCeil.init(
SIMPLE_FUNC_LINEAR, 0.0f, DMAX,
DEFAULT_N_CEIL_CONST, DEFAULT_N_CEIL_LINEAR);
multiseedIval.init(
DEFAULT_IVAL, 1.0f, DMAX,
DEFAULT_IVAL_B, DEFAULT_IVAL_A);
nCatPair = DEFAULT_N_CAT_PAIR;
if(!noisyHpolymer) {
penRdExConst = DEFAULT_READ_GAP_CONST;
penRdExLinear = DEFAULT_READ_GAP_LINEAR;
penRfExConst = DEFAULT_REF_GAP_CONST;
penRfExLinear = DEFAULT_REF_GAP_LINEAR;
} else {
penRdExConst = DEFAULT_READ_GAP_CONST_BADHPOLY;
penRdExLinear = DEFAULT_READ_GAP_LINEAR_BADHPOLY;
penRfExConst = DEFAULT_REF_GAP_CONST_BADHPOLY;
penRfExLinear = DEFAULT_REF_GAP_LINEAR_BADHPOLY;
}
multiseedMms = DEFAULT_SEEDMMS;
multiseedLen = DEFAULT_SEEDLEN;
EList<string> toks(MISC_CAT);
string tok;
istringstream ss(s);
int setting = 0;
// Get each ;-separated token
while(getline(ss, tok, ';')) {
setting++;
EList<string> etoks(MISC_CAT);
string etok;
// Divide into tokens on either side of =
istringstream ess(tok);
while(getline(ess, etok, '=')) {
etoks.push_back(etok);
}
// Must be exactly 1 =
if(etoks.size() != 2) {
cerr << "Error parsing alignment policy setting " << setting
<< "; must be bisected by = sign" << endl
<< "Policy: " << s.c_str() << endl;
assert(false); throw 1;
}
// LHS is tag, RHS value
string tag = etoks[0], val = etoks[1];
// Separate value into comma-separated tokens
EList<string> ctoks(MISC_CAT);
string ctok;
istringstream css(val);
while(getline(css, ctok, ',')) {
ctoks.push_back(ctok);
}
if(ctoks.size() == 0) {
cerr << "Error parsing alignment policy setting " << setting
<< "; RHS must have at least 1 token" << endl
<< "Policy: " << s.c_str() << endl;
assert(false); throw 1;
}
for(size_t i = 0; i < ctoks.size(); i++) {
if(ctoks[i].length() == 0) {
cerr << "Error parsing alignment policy setting " << setting
<< "; token " << i+1 << " on RHS had length=0" << endl
<< "Policy: " << s.c_str() << endl;
assert(false); throw 1;
}
}
// Bonus for a match
// MA=xx (default: MA=0, or MA=10 if --local is set)
if(tag == "MA") {
if(ctoks.size() != 1) {
cerr << "Error parsing alignment policy setting " << setting
<< "; RHS must have 1 token" << endl
<< "Policy: " << s.c_str() << endl;
assert(false); throw 1;
}
string tmp = ctoks[0];
istringstream tmpss(tmp);
tmpss >> bonusMatch;
}
// Scoring for mismatches
// MMP={Cxx|Q|RQ}
// Cxx = constant, where constant is integer xx
// Qxx = equal to quality, scaled
// R = equal to maq-rounded quality value (rounded to nearest
// 10, can't be greater than 30)
else if(tag == "MMP") {
if(ctoks.size() > 3) {
cerr << "Error parsing alignment policy setting "
<< "'" << tag.c_str() << "'"
<< "; RHS must have at most 3 tokens" << endl
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
if(ctoks[0][0] == 'C') {
string tmp = ctoks[0].substr(1);
// Parse constant penalty
istringstream tmpss(tmp);
tmpss >> penMmcMax;
penMmcMin = penMmcMax;
// Parse constant penalty
penMmcType = COST_MODEL_CONSTANT;
} else if(ctoks[0][0] == 'Q') {
if(ctoks.size() >= 2) {
string tmp = ctoks[1];
istringstream tmpss(tmp);
tmpss >> penMmcMax;
} else {
penMmcMax = DEFAULT_MM_PENALTY_MAX;
}
if(ctoks.size() >= 3) {
string tmp = ctoks[2];
istringstream tmpss(tmp);
tmpss >> penMmcMin;
} else {
penMmcMin = DEFAULT_MM_PENALTY_MIN;
}
if(penMmcMin > penMmcMax) {
cerr << "Error: Maximum mismatch penalty (" << penMmcMax
<< ") is less than minimum penalty (" << penMmcMin
<< endl;
throw 1;
}
// Set type to =quality
penMmcType = COST_MODEL_QUAL;
} else if(ctoks[0][0] == 'R') {
// Set type to=Maq-quality
penMmcType = COST_MODEL_ROUNDED_QUAL;
} else {
cerr << "Error parsing alignment policy setting "
<< "'" << tag.c_str() << "'"
<< "; RHS must start with C, Q or R" << endl
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
}
// Scoring for mismatches where read char=N
// NP={Cxx|Q|RQ}
// Cxx = constant, where constant is integer xx
// Q = equal to quality
// R = equal to maq-rounded quality value (rounded to nearest
// 10, can't be greater than 30)
else if(tag == "NP") {
if(ctoks.size() != 1) {
cerr << "Error parsing alignment policy setting "
<< "'" << tag.c_str() << "'"
<< "; RHS must have 1 token" << endl
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
if(ctoks[0][0] == 'C') {
string tmp = ctoks[0].substr(1);
// Parse constant penalty
istringstream tmpss(tmp);
tmpss >> penN;
// Parse constant penalty
penNType = COST_MODEL_CONSTANT;
} else if(ctoks[0][0] == 'Q') {
// Set type to =quality
penNType = COST_MODEL_QUAL;
} else if(ctoks[0][0] == 'R') {
// Set type to=Maq-quality
penNType = COST_MODEL_ROUNDED_QUAL;
} else {
cerr << "Error parsing alignment policy setting "
<< "'" << tag.c_str() << "'"
<< "; RHS must start with C, Q or R" << endl
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
}
// Scoring for read gaps
// RDG=xx,yy,zz
// xx = read gap open penalty
// yy = read gap extension penalty constant coefficient
// (defaults to open penalty)
// zz = read gap extension penalty linear coefficient
// (defaults to 0)
else if(tag == "RDG") {
if(ctoks.size() >= 1) {
istringstream tmpss(ctoks[0]);
tmpss >> penRdExConst;
} else {
penRdExConst = noisyHpolymer ?
DEFAULT_READ_GAP_CONST_BADHPOLY :
DEFAULT_READ_GAP_CONST;
}
if(ctoks.size() >= 2) {
istringstream tmpss(ctoks[1]);
tmpss >> penRdExLinear;
} else {
penRdExLinear = noisyHpolymer ?
DEFAULT_READ_GAP_LINEAR_BADHPOLY :
DEFAULT_READ_GAP_LINEAR;
}
}
// Scoring for reference gaps
// RFG=xx,yy,zz
// xx = ref gap open penalty
// yy = ref gap extension penalty constant coefficient
// (defaults to open penalty)
// zz = ref gap extension penalty linear coefficient
// (defaults to 0)
else if(tag == "RFG") {
if(ctoks.size() >= 1) {
istringstream tmpss(ctoks[0]);
tmpss >> penRfExConst;
} else {
penRfExConst = noisyHpolymer ?
DEFAULT_REF_GAP_CONST_BADHPOLY :
DEFAULT_REF_GAP_CONST;
}
if(ctoks.size() >= 2) {
istringstream tmpss(ctoks[1]);
tmpss >> penRfExLinear;
} else {
penRfExLinear = noisyHpolymer ?
DEFAULT_REF_GAP_LINEAR_BADHPOLY :
DEFAULT_REF_GAP_LINEAR;
}
}
// Minimum score as a function of read length
// MIN=xx,yy
// xx = constant coefficient
// yy = linear coefficient
else if(tag == "MIN") {
PARSE_FUNC(costMin);
}
// Per-read N ceiling as a function of read length
// NCEIL=xx,yy
// xx = N ceiling constant coefficient
// yy = N ceiling linear coefficient (set to 0 if unspecified)
else if(tag == "NCEIL") {
PARSE_FUNC(nCeil);
}
/*
* Seeds
* -----
*
* SEED=mm,len,ival (default: SEED=0,22)
*
* mm = Maximum number of mismatches allowed within a seed.
* Must be >= 0 and <= 2. Note that 2-mismatch mode is
* not fully sensitive; i.e. some 2-mismatch seed
* alignments may be missed.
* len = Length of seed.
* ival = Interval between seeds. If not specified, seed
* interval is determined by IVAL.
*/
else if(tag == "SEED") {
if(ctoks.size() > 2) {
cerr << "Error parsing alignment policy setting "
<< "'" << tag.c_str() << "'; RHS must have 1 or 2 tokens, "
<< "had " << ctoks.size() << ". "
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
if(ctoks.size() >= 1) {
istringstream tmpss(ctoks[0]);
tmpss >> multiseedMms;
if(multiseedMms > 1) {
cerr << "Error: -N was set to " << multiseedMms << ", but cannot be set greater than 1" << endl;
throw 1;
}
if(multiseedMms < 0) {
cerr << "Error: -N was set to a number less than 0 (" << multiseedMms << ")" << endl;
throw 1;
}
}
if(ctoks.size() >= 2) {
istringstream tmpss(ctoks[1]);
tmpss >> multiseedLen;
} else {
multiseedLen = DEFAULT_SEEDLEN;
}
}
else if(tag == "SEEDLEN") {
if(ctoks.size() > 1) {
cerr << "Error parsing alignment policy setting "
<< "'" << tag.c_str() << "'; RHS must have 1 token, "
<< "had " << ctoks.size() << ". "
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
if(ctoks.size() >= 1) {
istringstream tmpss(ctoks[0]);
tmpss >> multiseedLen;
}
}
else if(tag == "DPS") {
if(ctoks.size() > 1) {
cerr << "Error parsing alignment policy setting "
<< "'" << tag.c_str() << "'; RHS must have 1 token, "
<< "had " << ctoks.size() << ". "
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
if(ctoks.size() >= 1) {
istringstream tmpss(ctoks[0]);
tmpss >> failStreak;
}
}
else if(tag == "ROUNDS") {
if(ctoks.size() > 1) {
cerr << "Error parsing alignment policy setting "
<< "'" << tag.c_str() << "'; RHS must have 1 token, "
<< "had " << ctoks.size() << ". "
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
if(ctoks.size() >= 1) {
istringstream tmpss(ctoks[0]);
tmpss >> seedRounds;
}
}
/*
* Seed interval
* -------------
*
* IVAL={L|S|C},a,b (default: IVAL=S,1.0,0.0)
*
* L = let interval between seeds be a linear function of the
* read length. xx and yy are the constant and linear
* coefficients respectively. In other words, the interval
* equals a * len + b, where len is the read length.
* Intervals less than 1 are rounded up to 1.
* S = let interval between seeds be a function of the sqaure
* root of the read length. xx and yy are the
* coefficients. In other words, the interval equals
* a * sqrt(len) + b, where len is the read length.
* Intervals less than 1 are rounded up to 1.
* C = Like S but uses cube root of length instead of square
* root.
*/
else if(tag == "IVAL") {
PARSE_FUNC(multiseedIval);
}
else if(tag == "INTRONLEN") {
assert(penIntronLen != NULL);
PARSE_FUNC((*penIntronLen));
}
else {
// Unknown tag
cerr << "Unexpected alignment policy setting "
<< "'" << tag.c_str() << "'" << endl
<< "Policy: '" << s.c_str() << "'" << endl;
assert(false); throw 1;
}
}
}
#ifdef ALIGNER_SEED_POLICY_MAIN
int main() {
int bonusMatchType;
int bonusMatch;
int penMmcType;
int penMmc;
int penNType;
int penN;
int penRdExConst;
int penRfExConst;
int penRdExLinear;
int penRfExLinear;
SimpleFunc costMin;
SimpleFunc costFloor;
SimpleFunc nCeil;
bool nCatPair;
int multiseedMms;
int multiseedLen;
SimpleFunc msIval;
SimpleFunc posfrac;
SimpleFunc rowmult;
uint32_t mhits;
{
cout << "Case 1: Defaults 1 ... ";
const char *pol = "";
SeedAlignmentPolicy::parseString(
string(pol),
false, // --local?
false, // noisy homopolymers a la 454?
false, // ignore qualities?
bonusMatchType,
bonusMatch,
penMmcType,
penMmc,
penNType,
penN,
penRdExConst,
penRfExConst,
penRdExLinear,
penRfExLinear,
costMin,
costFloor,
nCeil,
nCatPair,
multiseedMms,
multiseedLen,
msIval,
mhits);
assert_eq(DEFAULT_MATCH_BONUS_TYPE, bonusMatchType);
assert_eq(DEFAULT_MATCH_BONUS, bonusMatch);
assert_eq(DEFAULT_MM_PENALTY_TYPE, penMmcType);
assert_eq(DEFAULT_MM_PENALTY_MAX, penMmcMax);
assert_eq(DEFAULT_MM_PENALTY_MIN, penMmcMin);
assert_eq(DEFAULT_N_PENALTY_TYPE, penNType);
assert_eq(DEFAULT_N_PENALTY, penN);
assert_eq(DEFAULT_MIN_CONST, costMin.getConst());
assert_eq(DEFAULT_MIN_LINEAR, costMin.getCoeff());
assert_eq(DEFAULT_FLOOR_CONST, costFloor.getConst());
assert_eq(DEFAULT_FLOOR_LINEAR, costFloor.getCoeff());
assert_eq(DEFAULT_N_CEIL_CONST, nCeil.getConst());
assert_eq(DEFAULT_N_CAT_PAIR, nCatPair);
assert_eq(DEFAULT_READ_GAP_CONST, penRdExConst);
assert_eq(DEFAULT_READ_GAP_LINEAR, penRdExLinear);
assert_eq(DEFAULT_REF_GAP_CONST, penRfExConst);
assert_eq(DEFAULT_REF_GAP_LINEAR, penRfExLinear);
assert_eq(DEFAULT_SEEDMMS, multiseedMms);
assert_eq(DEFAULT_SEEDLEN, multiseedLen);
assert_eq(DEFAULT_IVAL, msIval.getType());
assert_eq(DEFAULT_IVAL_A, msIval.getCoeff());
assert_eq(DEFAULT_IVAL_B, msIval.getConst());
cout << "PASSED" << endl;
}
{
cout << "Case 2: Defaults 2 ... ";
const char *pol = "";
SeedAlignmentPolicy::parseString(
string(pol),
false, // --local?
true, // noisy homopolymers a la 454?
false, // ignore qualities?
bonusMatchType,
bonusMatch,
penMmcType,
penMmc,
penNType,
penN,
penRdExConst,
penRfExConst,
penRdExLinear,
penRfExLinear,
costMin,
costFloor,
nCeil,
nCatPair,
multiseedMms,
multiseedLen,
msIval,
mhits);
assert_eq(DEFAULT_MATCH_BONUS_TYPE, bonusMatchType);
assert_eq(DEFAULT_MATCH_BONUS, bonusMatch);
assert_eq(DEFAULT_MM_PENALTY_TYPE, penMmcType);
assert_eq(DEFAULT_MM_PENALTY_MAX, penMmc);
assert_eq(DEFAULT_MM_PENALTY_MIN, penMmc);
assert_eq(DEFAULT_N_PENALTY_TYPE, penNType);
assert_eq(DEFAULT_N_PENALTY, penN);
assert_eq(DEFAULT_MIN_CONST, costMin.getConst());
assert_eq(DEFAULT_MIN_LINEAR, costMin.getCoeff());
assert_eq(DEFAULT_FLOOR_CONST, costFloor.getConst());
assert_eq(DEFAULT_FLOOR_LINEAR, costFloor.getCoeff());
assert_eq(DEFAULT_N_CEIL_CONST, nCeil.getConst());
assert_eq(DEFAULT_N_CAT_PAIR, nCatPair);
assert_eq(DEFAULT_READ_GAP_CONST_BADHPOLY, penRdExConst);
assert_eq(DEFAULT_READ_GAP_LINEAR_BADHPOLY, penRdExLinear);
assert_eq(DEFAULT_REF_GAP_CONST_BADHPOLY, penRfExConst);
assert_eq(DEFAULT_REF_GAP_LINEAR_BADHPOLY, penRfExLinear);
assert_eq(DEFAULT_SEEDMMS, multiseedMms);
assert_eq(DEFAULT_SEEDLEN, multiseedLen);
assert_eq(DEFAULT_IVAL, msIval.getType());
assert_eq(DEFAULT_IVAL_A, msIval.getCoeff());
assert_eq(DEFAULT_IVAL_B, msIval.getConst());
cout << "PASSED" << endl;
}
{
cout << "Case 3: Defaults 3 ... ";
const char *pol = "";
SeedAlignmentPolicy::parseString(
string(pol),
true, // --local?
false, // noisy homopolymers a la 454?
false, // ignore qualities?
bonusMatchType,
bonusMatch,
penMmcType,
penMmc,
penNType,
penN,
penRdExConst,
penRfExConst,
penRdExLinear,
penRfExLinear,
costMin,
costFloor,
nCeil,
nCatPair,
multiseedMms,
multiseedLen,
msIval,
mhits);
assert_eq(DEFAULT_MATCH_BONUS_TYPE_LOCAL, bonusMatchType);
assert_eq(DEFAULT_MATCH_BONUS_LOCAL, bonusMatch);
assert_eq(DEFAULT_MM_PENALTY_TYPE, penMmcType);
assert_eq(DEFAULT_MM_PENALTY_MAX, penMmcMax);
assert_eq(DEFAULT_MM_PENALTY_MIN, penMmcMin);
assert_eq(DEFAULT_N_PENALTY_TYPE, penNType);
assert_eq(DEFAULT_N_PENALTY, penN);
assert_eq(DEFAULT_MIN_CONST_LOCAL, costMin.getConst());
assert_eq(DEFAULT_MIN_LINEAR_LOCAL, costMin.getCoeff());
assert_eq(DEFAULT_FLOOR_CONST_LOCAL, costFloor.getConst());
assert_eq(DEFAULT_FLOOR_LINEAR_LOCAL, costFloor.getCoeff());
assert_eq(DEFAULT_N_CEIL_CONST, nCeil.getConst());
assert_eq(DEFAULT_N_CEIL_LINEAR, nCeil.getCoeff());
assert_eq(DEFAULT_N_CAT_PAIR, nCatPair);
assert_eq(DEFAULT_READ_GAP_CONST, penRdExConst);
assert_eq(DEFAULT_READ_GAP_LINEAR, penRdExLinear);
assert_eq(DEFAULT_REF_GAP_CONST, penRfExConst);
assert_eq(DEFAULT_REF_GAP_LINEAR, penRfExLinear);
assert_eq(DEFAULT_SEEDMMS, multiseedMms);
assert_eq(DEFAULT_SEEDLEN, multiseedLen);
assert_eq(DEFAULT_IVAL, msIval.getType());
assert_eq(DEFAULT_IVAL_A, msIval.getCoeff());
assert_eq(DEFAULT_IVAL_B, msIval.getConst());
cout << "PASSED" << endl;
}
{
cout << "Case 4: Simple string 1 ... ";
const char *pol = "MMP=C44;MA=4;RFG=24,12;FL=C,8;RDG=2;NP=C4;MIN=C,7";
SeedAlignmentPolicy::parseString(
string(pol),
true, // --local?
false, // noisy homopolymers a la 454?
false, // ignore qualities?
bonusMatchType,
bonusMatch,
penMmcType,
penMmc,
penNType,
penN,
penRdExConst,
penRfExConst,
penRdExLinear,
penRfExLinear,
costMin,
costFloor,
nCeil,
nCatPair,
multiseedMms,
multiseedLen,
msIval,
mhits);
assert_eq(COST_MODEL_CONSTANT, bonusMatchType);
assert_eq(4, bonusMatch);
assert_eq(COST_MODEL_CONSTANT, penMmcType);
assert_eq(44, penMmc);
assert_eq(COST_MODEL_CONSTANT, penNType);
assert_eq(4.0f, penN);
assert_eq(7, costMin.getConst());
assert_eq(DEFAULT_MIN_LINEAR_LOCAL, costMin.getCoeff());
assert_eq(8, costFloor.getConst());
assert_eq(DEFAULT_FLOOR_LINEAR_LOCAL, costFloor.getCoeff());
assert_eq(DEFAULT_N_CEIL_CONST, nCeil.getConst());
assert_eq(DEFAULT_N_CEIL_LINEAR, nCeil.getCoeff());
assert_eq(DEFAULT_N_CAT_PAIR, nCatPair);
assert_eq(2.0f, penRdExConst);
assert_eq(DEFAULT_READ_GAP_LINEAR, penRdExLinear);
assert_eq(24.0f, penRfExConst);
assert_eq(12.0f, penRfExLinear);
assert_eq(DEFAULT_SEEDMMS, multiseedMms);
assert_eq(DEFAULT_SEEDLEN, multiseedLen);
assert_eq(DEFAULT_IVAL, msIval.getType());
assert_eq(DEFAULT_IVAL_A, msIval.getCoeff());
assert_eq(DEFAULT_IVAL_B, msIval.getConst());
cout << "PASSED" << endl;
}
}
#endif /*def ALIGNER_SEED_POLICY_MAIN*/