diff --git a/alignment/substitution.cpp b/alignment/substitution.cpp index 953f41ae9..81f62681e 100644 --- a/alignment/substitution.cpp +++ b/alignment/substitution.cpp @@ -23,7 +23,7 @@ Substitution::Substitution(const std::string& sub_str, Alignment* const aln, con outError("Failed to parse the predefined mutation: '" + sub_str + "'. The new state is invalid."); // Parse the position - position = convert_int(sub_str.substr(num_sites_per_state, input_length - (num_sites_per_state + num_sites_per_state)).c_str()) - 1; + position = convert_int(sub_str.substr(num_sites_per_state, input_length - (num_sites_per_state + num_sites_per_state)).c_str()) - Params::getInstance().site_starting_index; if (aln->seq_type == SEQ_CODON) position = position * ONE_THIRD; diff --git a/main/alisim.cpp b/main/alisim.cpp index d5ebc313f..d9d37e549 100644 --- a/main/alisim.cpp +++ b/main/alisim.cpp @@ -364,7 +364,7 @@ std::vector> readMutations(const std::string& } // if not found -> all characters are spaces -> invalid else - outError("Invalid format '" + line + "'. Expected format should be "); + outError("Invalid format '" + line + "'. Expected format should be "); // remove spaces at the ending of the line pos = line.find_last_not_of(" "); @@ -377,11 +377,11 @@ std::vector> readMutations(const std::string& } // if not found -> all characters are spaces -> invalid else - outError("Invalid format '" + line + "'. Expected format should be "); + outError("Invalid format '" + line + "'. Expected format should be "); - // validate the input format + // validate the input format if (line.length() < 5) - outError("Invalid format '" + line + "'. Expected format should be "); + outError("Invalid format '" + line + "'. Expected format should be "); // replace \t by a space line = regex_replace(line, std::regex("\t"), " "); @@ -397,7 +397,7 @@ std::vector> readMutations(const std::string& transform(node_name.begin(), node_name.end(), node_name.begin(), ::toupper); } else - outError("Invalid format '" + line + "'. Expected format should be "); + outError("Invalid format '" + line + "'. Expected format should be "); // get mutation list at that node std::string mutation_str = ""; @@ -408,7 +408,7 @@ std::vector> readMutations(const std::string& if (line.length() - pos) mutation_str = line.substr(pos, line.length() - pos); else - outError("Invalid format '" + line + "'. Expected format should be "); + outError("Invalid format '" + line + "'. Expected format should be "); // record the mutations at that node node_mutations.push_back(std::pair(node_name, mutation_str)); @@ -646,7 +646,7 @@ void getLockedSites(Node* const node, Node* const dad, std::vector* const // vailidate position if (pos >= seq_length) { - outWarning("Ignore a predefined mutation " + aln->convertStateBackStr(mut_it->getOldState()) + convertIntToString((aln->seq_type == SEQ_CODON ? pos * 3 : pos) + 1) + aln->convertStateBackStr(mut_it->getNewState()) + ". Position exceeds the sequence length " + convertIntToString(aln->seq_type == SEQ_CODON ? seq_length * 3 : seq_length)); + outWarning("Ignore a predefined mutation " + aln->convertStateBackStr(mut_it->getOldState()) + convertIntToString((aln->seq_type == SEQ_CODON ? pos * 3 : pos) + Params::getInstance().site_starting_index) + aln->convertStateBackStr(mut_it->getNewState()) + ". Position exceeds the sequence length " + convertIntToString(aln->seq_type == SEQ_CODON ? seq_length * 3 : seq_length)); } // mark the site as locked else diff --git a/simulator/alisimulator.cpp b/simulator/alisimulator.cpp index f522d7afa..2d701472e 100644 --- a/simulator/alisimulator.cpp +++ b/simulator/alisimulator.cpp @@ -1591,12 +1591,12 @@ void AliSimulator::handlePreMutations(const NeighborVec::iterator& it, int& pred { // show infor if (verbose_mode >= VB_MED) - std::cout << "Apply a predefined mutation " << tree->aln->convertStateBackStr(mut_it->getOldState()) << convertIntToString(mut_it->getPosition() * num_sites_per_state + 1) << tree->aln->convertStateBackStr(mut_it->getNewState()) << std::endl; + std::cout << "Apply a predefined mutation " << tree->aln->convertStateBackStr(mut_it->getOldState()) << convertIntToString(mut_it->getPosition() * num_sites_per_state + params->site_starting_index) << tree->aln->convertStateBackStr(mut_it->getNewState()) << std::endl; // if the old state is different from the one observed at the specified position -> show a warning if ((*node_seq_chunk)[relative_pos] != mut_it->getOldState()) { - outWarning("State " + tree->aln->convertStateBackStr(mut_it->getOldState()) + " is not observed at position/site " + convertIntToString(mut_it->getPosition() * num_sites_per_state + 1) + ", we observe state " + tree->aln->convertStateBackStr((*node_seq_chunk)[relative_pos]) + " instead."); + outWarning("State " + tree->aln->convertStateBackStr(mut_it->getOldState()) + " is not observed at position/site " + convertIntToString(mut_it->getPosition() * num_sites_per_state + params->site_starting_index) + ", we observe state " + tree->aln->convertStateBackStr((*node_seq_chunk)[relative_pos]) + " instead."); } // apply the substitution @@ -1608,7 +1608,7 @@ void AliSimulator::handlePreMutations(const NeighborVec::iterator& it, int& pred // write log - for debugging only /*ofstream myfile; myfile.open ("thread_" + convertIntToString(omp_get_thread_num()) + ".txt", std::ostream::app); - myfile << tree->aln->convertStateBackStr(mut_it->getOldState()) << mut_it->getPosition() + 1 << tree->aln->convertStateBackStr(mut_it->getNewState())<< std::endl; + myfile << tree->aln->convertStateBackStr(mut_it->getOldState()) << mut_it->getPosition() + params->site_starting_index << tree->aln->convertStateBackStr(mut_it->getNewState())<< std::endl; myfile.close();*/ } } diff --git a/utils/tools.cpp b/utils/tools.cpp index c36f6af97..0ccbf1aa3 100644 --- a/utils/tools.cpp +++ b/utils/tools.cpp @@ -1487,6 +1487,7 @@ void parseArg(int argc, char *argv[], Params ¶ms) { params.alignment_id = 0; params.include_pre_mutations = false; params.mutation_file = ""; + params.site_starting_index = 0; // store original params for (cnt = 1; cnt < argc; cnt++) { @@ -5267,6 +5268,12 @@ void parseArg(int argc, char *argv[], Params ¶ms) { continue; } + if (strcmp(argv[cnt], "--index-from-one") == 0) { + params.site_starting_index = 1; + + continue; + } + if (strcmp(argv[cnt], "--distribution") == 0) { cnt++; if (cnt >= argc || argv[cnt][0] == '-') diff --git a/utils/tools.h b/utils/tools.h index e1d3fc716..4feccba8c 100644 --- a/utils/tools.h +++ b/utils/tools.h @@ -2601,6 +2601,11 @@ class Params { * Mutation file that specifies pre-defined mutations occurs at nodes */ std::string mutation_file; + + /** + * site starting index (for predefined mutations in AliSim) + */ + int site_starting_index; }; /**