Enable stored E-value params for any genetic code
authorMartin C. Frith
Mon Oct 07 12:02:11 2019 +0900 (6 weeks ago)
changeset 987f757e6ca4edb
parent 986 88adecc284af
child 988 ca6613dd1f06
Enable stored E-value params for any genetic code
src/LastEvaluer.cc
src/LastEvaluer.hh
src/lastal.cc
     1.1 --- a/src/LastEvaluer.cc	Mon Oct 07 10:57:27 2019 +0900
     1.2 +++ b/src/LastEvaluer.cc	Mon Oct 07 12:02:11 2019 +0900
     1.3 @@ -31,6 +31,7 @@
     1.4  };
     1.5  
     1.6  struct FrameshiftEvalueParameters {
     1.7 +  const char *geneticCodeName;
     1.8    const char *matrixName;
     1.9    int gapOpen;
    1.10    int gapEpen;
    1.11 @@ -129,41 +130,41 @@
    1.12  };
    1.13  
    1.14  const FrameshiftEvalueParameters frameshiftEvalueParameters[] = {
    1.15 -  {"BL62", 11, 1, 15, {0.31457181182385774, 0.077024909125411836,
    1.16 -		       3.5355057419386005, -39.014329998056937,
    1.17 -		       1.1739847579695837, -12.896364921780187,
    1.18 -		       160.29749789587885, -3200.2716722761552,
    1.19 -		       17.539096792506459, -349.06406999516196,
    1.20 -		       51.79266617536797, -1027.401229133852}},
    1.21 +  {"1", "BL62", 11, 1, 15, {0.31457181182385774, 0.077024909125411836,
    1.22 +			    3.5355057419386005, -39.014329998056937,
    1.23 +			    1.1739847579695837, -12.896364921780187,
    1.24 +			    160.29749789587885, -3200.2716722761552,
    1.25 +			    17.539096792506459, -349.06406999516196,
    1.26 +			    51.79266617536797, -1027.401229133852}},
    1.27  
    1.28    // this one is included only to make lest-test.sh faster:
    1.29 -  {"BL62", 11, 2, 12, {0.32704828292493776, 0.10543687626821494,
    1.30 -		       2.7166240938670798, -19.361170444340445,
    1.31 -		       0.89899140520122844, -6.29652445533966,
    1.32 -		       74.744763449386667, -1147.0060455603425,
    1.33 -		       8.3059555754127565, -127.46868078491309,
    1.34 -		       24.78179014970371, -379.14020451790986}},
    1.35 +  {"1", "BL62", 11, 2, 12, {0.32704828292493776, 0.10543687626821494,
    1.36 +			    2.7166240938670798, -19.361170444340445,
    1.37 +			    0.89899140520122844, -6.29652445533966,
    1.38 +			    74.744763449386667, -1147.0060455603425,
    1.39 +			    8.3059555754127565, -127.46868078491309,
    1.40 +			    24.78179014970371, -379.14020451790986}},
    1.41  
    1.42 -  {"BL62", 11, 2, 15, {0.33388770870821022, 0.11532516007803961,
    1.43 -		       2.4678058049483518, -14.50532580281522,
    1.44 -		       0.82160372991753583, -4.8091552692419572,
    1.45 -		       55.103072639059761, -731.90592162187147,
    1.46 -		       6.1010321043683131, -80.763060603166991,
    1.47 -		       18.237750047538203, -240.59017890476582}},
    1.48 +  {"1", "BL62", 11, 2, 15, {0.33388770870821022, 0.11532516007803961,
    1.49 +			    2.4678058049483518, -14.50532580281522,
    1.50 +			    0.82160372991753583, -4.8091552692419572,
    1.51 +			    55.103072639059761, -731.90592162187147,
    1.52 +			    6.1010321043683131, -80.763060603166991,
    1.53 +			    18.237750047538203, -240.59017890476582}},
    1.54  
    1.55 -  {"BL80", 11, 1, 15, {0.35400649542314511, 0.13270256108942211,
    1.56 -		       1.8960749679829285, -14.061923223904673,
    1.57 -		       0.62940827123451903, -4.6245065070665863,
    1.58 -		       35.18772909081801, -583.19242886423649,
    1.59 -		       3.8260214679558033, -62.789729751450864,
    1.60 -		       11.072568656496113, -178.63729131744145}},
    1.61 +  {"1", "BL80", 11, 1, 15, {0.35400649542314511, 0.13270256108942211,
    1.62 +			    1.8960749679829285, -14.061923223904673,
    1.63 +			    0.62940827123451903, -4.6245065070665863,
    1.64 +			    35.18772909081801, -583.19242886423649,
    1.65 +			    3.8260214679558033, -62.789729751450864,
    1.66 +			    11.072568656496113, -178.63729131744145}},
    1.67  
    1.68 -  {"BL80", 11, 2, 15, {0.3652492341706855, 0.16850422398182774,
    1.69 -		       1.5316138005575799, -5.7577598061709985,
    1.70 -		       0.5101720323233776, -1.9097398376324572,
    1.71 -		       17.427364219333899, -170.0223112776693,
    1.72 -		       1.9259860816444827, -18.621287186644096,
    1.73 -		       5.7329546520583801, -54.693768145180513}},
    1.74 +  {"1", "BL80", 11, 2, 15, {0.3652492341706855, 0.16850422398182774,
    1.75 +			    1.5316138005575799, -5.7577598061709985,
    1.76 +			    0.5101720323233776, -1.9097398376324572,
    1.77 +			    17.427364219333899, -170.0223112776693,
    1.78 +			    1.9259860816444827, -18.621287186644096,
    1.79 +			    5.7329546520583801, -54.693768145180513}},
    1.80  };
    1.81  
    1.82  static bool isEqual(const char *x, const char *y) {
    1.83 @@ -182,8 +183,8 @@
    1.84  }
    1.85  
    1.86  static bool isHit(const FrameshiftEvalueParameters &p,
    1.87 -		  const char *n, int a, int b, int f) {
    1.88 -  return isEqual(p.matrixName, n) &&
    1.89 +		  const char *g, const char *n, int a, int b, int f) {
    1.90 +  return isEqual(p.geneticCodeName, g) && isEqual(p.matrixName, n) &&
    1.91      p.gapOpen == a && p.gapEpen == b && p.frameshiftCost == f;
    1.92  }
    1.93  
    1.94 @@ -264,7 +265,7 @@
    1.95  		       int insEpen,
    1.96  		       int frameshiftCost,
    1.97  		       const GeneticCode &geneticCode,
    1.98 -		       bool isStandardGeneticCode,
    1.99 +		       const char *geneticCodeName,
   1.100  		       int verbosity) {
   1.101    const double lambdaTolerance = 0.01;
   1.102    const double kTolerance = 0.05;
   1.103 @@ -276,10 +277,11 @@
   1.104  
   1.105    if (frameshiftCost >= 0) {  // DNA-versus-protein alignment:
   1.106      if (isGapped && insOpen == delOpen && insEpen == delEpen) {
   1.107 -      if (isProtein(alphabet) && isStandardGeneticCode) {
   1.108 +      if (isProtein(alphabet)) {
   1.109  	for (size_t i = 0; i < COUNTOF(frameshiftEvalueParameters); ++i) {
   1.110  	  const FrameshiftEvalueParameters &p = frameshiftEvalueParameters[i];
   1.111 -	  if (isHit(p, matrixName, delOpen, delEpen, frameshiftCost))
   1.112 +	  if (isHit(p, geneticCodeName, matrixName, delOpen, delEpen,
   1.113 +		    frameshiftCost))
   1.114  	    return frameshiftEvaluer.initParameters(p.parameters);
   1.115  	}
   1.116        }
     2.1 --- a/src/LastEvaluer.hh	Mon Oct 07 10:57:27 2019 +0900
     2.2 +++ b/src/LastEvaluer.hh	Mon Oct 07 12:02:11 2019 +0900
     2.3 @@ -32,7 +32,7 @@
     2.4    // alignment parameters.  It may fail, i.e. set the object to the
     2.5    // "bad" state and throw an Sls::error.
     2.6    // These arguments are only used to lookup pre-calculated cases:
     2.7 -  // matrixName, matchScore, mismatchCost, isStandardGeneticCode.
     2.8 +  // matrixName, matchScore, mismatchCost, geneticCodeName.
     2.9    // DNA-versus-protein alignment is indicated by: frameshiftCost >= 0.
    2.10    // As a special case, frameshiftCost==0 means no frameshifts.
    2.11    // For DNA-versus-protein alignment, letterFreqs2 is not used.
    2.12 @@ -50,7 +50,7 @@
    2.13  	    int insEpen,
    2.14  	    int frameshiftCost,
    2.15  	    const GeneticCode &geneticCode,
    2.16 -	    bool isStandardGeneticCode,
    2.17 +	    const char *geneticCodeName,
    2.18  	    int verbosity);
    2.19  
    2.20    void setSearchSpace(double databaseLength,  // number of database letters
     3.1 --- a/src/lastal.cc	Mon Oct 07 10:57:27 2019 +0900
     3.2 +++ b/src/lastal.cc	Mon Oct 07 12:02:11 2019 +0900
     3.3 @@ -259,7 +259,6 @@
     3.4    const char *canonicalMatrixName = ScoreMatrix::canonicalName( matrixName );
     3.5    if (args.temperature > 0 && !matrixName.empty()) canonicalMatrixName = " ";
     3.6    bool isGapped = (args.outputType > 1);
     3.7 -  bool isStandardGeneticCode = (args.geneticCodeFile == "1");
     3.8    LOG( "getting E-value parameters..." );
     3.9    try{
    3.10      const mcf::GapCosts::Piece &del = gapCosts.delPieces[0];
    3.11 @@ -268,8 +267,8 @@
    3.12                    alph.letters.c_str(), fwdMatrices.scoresMasked,
    3.13  		  stats.letterProbs1(), stats.letterProbs2(), isGapped,
    3.14  		  del.openCost, del.growCost, ins.openCost, ins.growCost,
    3.15 -                  args.frameshiftCost, geneticCode, isStandardGeneticCode,
    3.16 -		  args.verbosity );
    3.17 +		  args.frameshiftCost, geneticCode,
    3.18 +		  args.geneticCodeFile.c_str(), args.verbosity );
    3.19      countT m = std::min(refMaxSeqLen, refLetters);
    3.20      evaluer.setSearchSpace(refLetters, m, args.numOfStrands());
    3.21      if( args.verbosity > 0 ) evaluer.writeParameters( std::cerr );