Skip to content

Rate variation and MAT updates #50

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
SR-Martin wants to merge 8 commits into
base: main_clone_rate_variation_MAT
Choose a base branch
from
Open

Rate variation and MAT updates #50

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
a1add9a
Fix bug
SR-Martin Feb 11, 2026
70a0eca
MATs
SR-Martin Feb 24, 2026
ba03b9b
Updates to rate-variation models
SR-Martin Feb 25, 2026
dc6ddc3
Updates to scalar rate variation model
SR-Martin Mar 16, 2026
21319c7
Update model_dna_rate_variation.cpp
SR-Martin Mar 18, 2026
fd95a17
Updates to MAT estimates
SR-Martin Mar 18, 2026
f8fc7c2
Update usage_cmaple
SR-Martin Mar 23, 2026
54a0a29
Update usage_cmaple 2
SR-Martin Mar 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 16 additions & 4 deletions maple/cmaple.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,8 +122,8 @@ void cmaple::runCMAPLE(cmaple::Params &params)
throw std::invalid_argument("Unknown Model " + params.sub_model_str);
}
assert(sub_model != cmaple::ModelBase::UNKNOWN);
bool useRateVariationModel = params.rate_variation || params.site_specific_rates;
Model model(aln.ref_seq.size(), useRateVariationModel, params.rate_variation, params.wt_pseudocount, params.rates_filename, sub_model, aln.getSeqType());
bool useRateVariationModel = params.rate_variation || params.site_specific_rate_matrix;
Model model(aln.ref_seq.size(), useRateVariationModel, params.rate_variation, params.wt_pseudocount, params.rates_filename, params.rate_variation_max_num_EM_steps, sub_model, aln.getSeqType());

// If users only want to convert the alignment to another format -> convert it and terminate
if (params.output_aln.length())
Expand Down Expand Up @@ -197,6 +197,16 @@ void cmaple::runCMAPLE(cmaple::Params &params)
out << tree.exportTSV();
out.close();
}

// export MAT if selected
if(params.output_MAT)
{
std::string filename = params.output_prefix + "_MAT.nex";
std::cout << "Writing MAT to file " << filename << std::endl;
ofstream out = ofstream(filename);
out << tree.exportNexus(tree_format, false, true);
out.close();
}

// output log-likelihood of the tree
if (cmaple::verbose_mode > cmaple::VB_QUIET) {
Expand All @@ -218,10 +228,12 @@ void cmaple::runCMAPLE(cmaple::Params &params)
// Show information about output files
std::cout << "Analysis results written to:" << std::endl;
std::cout << "Maximum-likelihood tree: " << output_treefile << std::endl;
if (params.output_MAT)
std::cout << "Estimated mutation-annotated tree (MAT): " << output_treefile + "_MAT.nwk" << std::endl;
if (params.output_NEXUS || params.compute_SPRTA)
std::cout << "Tree in NEXUS format: " << output_treefile + ".nex" << std::endl;
std::cout << "Tree in NEXUS format: " << output_treefile + ".nex" << std::endl;
if (params.compute_SPRTA && params.output_alternative_spr)
std::cout << "Meta data in TSV format: " << output_treefile + ".tsv" << std::endl;
std::cout << "Meta data in TSV format: " << output_treefile + ".tsv" << std::endl;
/*if (params.compute_aLRT_SH) {
std::cout << "Tree with aLRT-SH values: "
<< prefix + ".aLRT_SH.treefile" << std::endl;
Expand Down
3 changes: 2 additions & 1 deletion model/model.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ cmaple::Model::Model( cmaple::PositionType ref_genome_size,
bool _siteRates,
cmaple::RealNumType wt_pseudocount,
const std::string _rates_filename,
int _max_num_EM_steps,
const cmaple::ModelBase::SubModel sub_model,
const cmaple::SeqRegion::SeqType seqtype)
: model_base(nullptr) {
Expand Down Expand Up @@ -75,7 +76,7 @@ cmaple::Model::Model( cmaple::PositionType ref_genome_size,
}
case cmaple::SeqRegion::SEQ_DNA: {
if(rate_variation){
model_base = new ModelDNARateVariation(n_sub_model, ref_genome_size, _siteRates, wt_pseudocount, _rates_filename);
model_base = new ModelDNARateVariation(n_sub_model, ref_genome_size, _siteRates, wt_pseudocount, _rates_filename, _max_num_EM_steps);
} else {
model_base = new ModelDNA(n_sub_model);
}
Expand Down
1 change: 1 addition & 0 deletions model/model.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ class Model {
bool _siteRates,
cmaple::RealNumType wt_pseudocount,
const std::string _rates_filename,
int _max_num_EM_steps,
const cmaple::ModelBase::SubModel sub_model = cmaple::ModelBase::DEFAULT,
const cmaple::SeqRegion::SeqType seqtype = cmaple::SeqRegion::SEQ_AUTO);

Expand Down
482 changes: 335 additions & 147 deletions model/model_dna_rate_variation.cpp
View file
Open in desktop

Large diffs are not rendered by default.

28 changes: 24 additions & 4 deletions model/model_dna_rate_variation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,14 @@ class Tree;
/** Class of DNA evolutionary models with rate variation */
class ModelDNARateVariation : public ModelDNA {
public:
ModelDNARateVariation( const cmaple::ModelBase::SubModel sub_model, PositionType _genome_size,
bool _use_site_rates, cmaple::RealNumType _wt_pseudocount, std::string _rates_filename);
ModelDNARateVariation(
const cmaple::ModelBase::SubModel sub_model,
PositionType _genome_size,
bool _use_site_rates,
cmaple::RealNumType _wt_pseudocount,
std::string _rates_filename,
int _max_num_EM_steps);

virtual ~ModelDNARateVariation();

void estimateRates(cmaple::Tree* tree);
Expand Down Expand Up @@ -79,14 +85,27 @@ class ModelDNARateVariation : public ModelDNA {

private:

void updateCountsAndWaitingTimesAcrossRoot( PositionType start, PositionType end,
void updateCountsAndWaitingTimesAcrossRoot( PositionType genome_pos,
StateType parent_state, StateType child_state,
RealNumType dist_to_root, RealNumType dist_to_observed,
RealNumType* waiting_times, RealNumType* counts,
RealNumType weight = 1.);

void readRatesFile();

std::vector<RealNumType> getRelativeProbabilityOfParentOStatesForRegion( const cmaple::SeqRegion* seqP_region,
StateType child_state,
RealNumType branch_length_to_obs,
PositionType genome_pos);
std::vector<RealNumType> getRelativeProbabilityOfChildOStatesForRegion( const cmaple::SeqRegion* seqC_region,
StateType parent_state,
RealNumType branch_length_to_obs,
PositionType genome_pos);
std::vector<RealNumType> getRelativeProbabilityOfParentOChildOStatesForRegion(const cmaple::SeqRegion* seqP_region,
const cmaple::SeqRegion* seqC_region,
RealNumType branch_length_to_obs,
PositionType genome_pos);

cmaple::PositionType genome_size;

cmaple::RealNumType* mutation_matrices = nullptr;
Expand All @@ -96,12 +115,13 @@ class ModelDNARateVariation : public ModelDNA {
cmaple::RealNumType* freqj_transposedijs = nullptr;
cmaple::RealNumType* rates = nullptr;
uint16_t mat_size;
bool use_site_rates = false;
bool scalar_rate_model = false;
bool rates_estimated = false;

cmaple::RealNumType waiting_time_pseudocount;

std::string rates_filename;
int max_num_EM_steps;

};
}
Loading