diff --git a/modules/seq/alg/pymod/wrap_seq_alg.cc b/modules/seq/alg/pymod/wrap_seq_alg.cc index b9935d13d6f833cfb46caa8961ece096f676cf05..40a858f6058da4fcc5393af84950b8b72efed0ef 100644 --- a/modules/seq/alg/pymod/wrap_seq_alg.cc +++ b/modules/seq/alg/pymod/wrap_seq_alg.cc @@ -267,7 +267,10 @@ void export_hmm_algorithms() { arg("a")=0.9, arg("b")=4.0, arg("c")=1.0)); - def("AddTransitionPseudoCounts", &AddTransitionPseudoCounts, (arg("profile"))); + def("AddTransitionPseudoCounts", &AddTransitionPseudoCounts, (arg("profile"), + arg("gapb")=1.0, + arg("gapd")=0.15, + arg("gape")=1.0)); def("AddNullPseudoCounts", &AddNullPseudoCounts, (arg("profile"))); def("HMMScore", &HMMScore, (arg("profile_0"), arg("profile_1"), arg("alignment"), arg("s_0_idx"), arg("s_1_idx"), diff --git a/modules/seq/alg/src/hmm_pseudo_counts.cc b/modules/seq/alg/src/hmm_pseudo_counts.cc index 6c031e868fbc52d92610f9e6fb4528902c8d56a8..37d2d490b7a6c7bfdc4fea67dc84d289c77c2960 100644 --- a/modules/seq/alg/src/hmm_pseudo_counts.cc +++ b/modules/seq/alg/src/hmm_pseudo_counts.cc @@ -416,15 +416,16 @@ ContextProfileDBPtr ContextProfileDB::FromCRF(const String& filename) { } -void AddTransitionPseudoCounts(ost::seq::ProfileHandle& profile) { +void AddTransitionPseudoCounts(ost::seq::ProfileHandle& profile, + Real gapb, Real gapd, Real gape) { // a priori probabilities estimated with default values of HHblits - Real pM2D = 0.15 * 0.0286; + Real pM2D = gapd * 0.0286; Real pM2I = pM2D; Real pM2M = 1 - pM2D - pM2I; - Real pI2I = 0.75; + Real pI2I = 1.0 * gape / (gape - 1 + 1.0 / 0.75);; Real pI2M = 1 - pI2I; - Real pD2D = 0.75; + Real pD2D = pI2I; Real pD2M = 1 - pD2D; for (size_t col_idx = 0; col_idx < profile.size(); ++col_idx) { @@ -433,9 +434,9 @@ void AddTransitionPseudoCounts(ost::seq::ProfileHandle& profile) { // Transitions from M state Real neff = data->GetNeff(); - Real p0 = (neff - 1) * data->GetProb(ost::seq::HMM_M2M) + pM2M; - Real p1 = (neff - 1) * data->GetProb(ost::seq::HMM_M2D) + pM2D; - Real p2 = (neff - 1) * data->GetProb(ost::seq::HMM_M2I) + pM2I; + Real p0 = (neff - 1) * data->GetProb(ost::seq::HMM_M2M) + gapb*pM2M; + Real p1 = (neff - 1) * data->GetProb(ost::seq::HMM_M2D) + gapb*pM2D; + Real p2 = (neff - 1) * data->GetProb(ost::seq::HMM_M2I) + gapb*pM2I; Real sum = p0 + p1 + p2; data->SetProb(ost::seq::HMM_M2M, p0/sum); data->SetProb(ost::seq::HMM_M2D, p1/sum); @@ -443,16 +444,16 @@ void AddTransitionPseudoCounts(ost::seq::ProfileHandle& profile) { // Transitions from I state Real neff_i = data->GetNeff_I(); - p0 = neff_i * data->GetProb(ost::seq::HMM_I2M) + pI2M; - p1 = neff_i * data->GetProb(ost::seq::HMM_I2I) + pI2I; + p0 = neff_i * data->GetProb(ost::seq::HMM_I2M) + gapb*pI2M; + p1 = neff_i * data->GetProb(ost::seq::HMM_I2I) + gapb*pI2I; sum = p0 + p1; data->SetProb(ost::seq::HMM_I2M, p0/sum); data->SetProb(ost::seq::HMM_I2I, p1/sum); // Transitions from D state Real neff_d = data->GetNeff_D(); - p0 = neff_d * data->GetProb(ost::seq::HMM_D2M) + pD2M; - p1 = neff_d * data->GetProb(ost::seq::HMM_D2D) + pD2D; + p0 = neff_d * data->GetProb(ost::seq::HMM_D2M) + gapb*pD2M; + p1 = neff_d * data->GetProb(ost::seq::HMM_D2D) + gapb*pD2D; sum = p0 + p1; data->SetProb(ost::seq::HMM_D2M, p0/sum); data->SetProb(ost::seq::HMM_D2D, p1/sum); diff --git a/modules/seq/alg/src/hmm_pseudo_counts.hh b/modules/seq/alg/src/hmm_pseudo_counts.hh index 573f7638520f2790a362031253070b0595a73584..b967d9e6aeda7ddafdcfa1603a299b09ed4819ff 100644 --- a/modules/seq/alg/src/hmm_pseudo_counts.hh +++ b/modules/seq/alg/src/hmm_pseudo_counts.hh @@ -158,7 +158,9 @@ private: std::vector<ContextProfile> profiles_; }; -void AddTransitionPseudoCounts(ost::seq::ProfileHandle& profile); +void AddTransitionPseudoCounts(ost::seq::ProfileHandle& profile, + Real gapb = 1.0, Real gabd = 0.15, + Real gape = 1.0); void AddAAPseudoCounts(ost::seq::ProfileHandle& profile, Real a = 1.0, Real b = 1.5, Real c = 1.0);