Comments (10)
@pedrotaucce will work on this tomorrow :-)
from beautier.
Thanks! I really appreciate that. I have more than 40 alignments and will have to different beast runs. Your package is being a lifesaver!
from beautier.
Hi @pedrotaucce, thanks for this excellent bug report (and sorry for replying so late). I will take a look at it now :-)
from beautier.
Plan:
- Recreate XML in BEAUti: maybe it cannot be done at all
- If it can be done: compare XML to babette XML
from beautier.
Recreate XML:
- Install BEAST2 by doing
beastierinstall::install_beast2()
- Loaded the alignments
- Change the clock rate
Aha, silly me: a strict clock is always normalized to 1.0 🤦
from beautier.
For completeness, this is the XML generated by BEAUti:
<?xml version="1.0" encoding="UTF-8" standalone="no"?><beast beautitemplate='Standard' beautistatus='' namespace="beast.core:beast.evolution.alignment:beast.evolution.tree.coalescent:beast.core.util:beast.evolution.nuc:beast.evolution.operators:beast.evolution.sitemodel:beast.evolution.substitutionmodel:beast.evolution.likelihood" required="" version="2.6">
<data
id="anthus_aco_sub"
spec="Alignment"
name="alignment">
<sequence id="seq_61430_aco" spec="Sequence" taxon="61430_aco" totalcount="4" value="ACAGGTTAGAAACTACTCTGTTTTCTGGCTCCTTGTTTAATGCCCTGTCCTATTTTATTGCGAAAATTGTCTGTTTTT"/>
<sequence id="seq_626029_aco" spec="Sequence" taxon="626029_aco" totalcount="4" value="ACAGGTTAGAAACTACTCTGTTTTCTGGCTGCTTGTTTAATGCCCTCTCCTATTTTATTGTGACGATTGTCTGTTTTT"/>
<sequence id="seq_630116_aco" spec="Sequence" taxon="630116_aco" totalcount="4" value="ACAGGTTAGAAACTACTCTGTTTTCTGGCTCCTTGTTTAATGCCCTGTCCTATTTTATTGCGAAAATTGTCTGTTTTT"/>
<sequence id="seq_630210_aco" spec="Sequence" taxon="630210_aco" totalcount="4" value="ACAGGTTAGAAACTACTCTGTTTTCTGGCTCCTTGTTTAATGCCCTGTCCTATTTTATTGTGACAATTGTCTGTTTTT"/>
<sequence id="seq_B25702_aco" spec="Sequence" taxon="B25702_aco" totalcount="4" value="ACAGGTTAGAAACTACTCTGTTTTCTGGCTCCTTGTTTAATGCCCTGTCCTATTTTATTGTGACAATTGTCTGTTTTT"/>
</data>
<map name="Uniform" >beast.math.distributions.Uniform</map>
<map name="Exponential" >beast.math.distributions.Exponential</map>
<map name="LogNormal" >beast.math.distributions.LogNormalDistributionModel</map>
<map name="Normal" >beast.math.distributions.Normal</map>
<map name="Beta" >beast.math.distributions.Beta</map>
<map name="Gamma" >beast.math.distributions.Gamma</map>
<map name="LaplaceDistribution" >beast.math.distributions.LaplaceDistribution</map>
<map name="prior" >beast.math.distributions.Prior</map>
<map name="InverseGamma" >beast.math.distributions.InverseGamma</map>
<map name="OneOnX" >beast.math.distributions.OneOnX</map>
<run id="mcmc" spec="MCMC" chainLength="10000000">
<state id="state" spec="State" storeEvery="5000">
<tree id="Tree.t:anthus_aco_sub" spec="beast.evolution.tree.Tree" name="stateNode">
<taxonset id="TaxonSet.anthus_aco_sub" spec="TaxonSet">
<alignment idref="anthus_aco_sub"/>
</taxonset>
</tree>
<parameter id="birthRate.t:anthus_aco_sub" spec="parameter.RealParameter" name="stateNode">1.0</parameter>
</state>
<init id="RandomTree.t:anthus_aco_sub" spec="beast.evolution.tree.RandomTree" estimate="false" initial="@Tree.t:anthus_aco_sub" taxa="@anthus_aco_sub">
<populationModel id="ConstantPopulation0.t:anthus_aco_sub" spec="ConstantPopulation">
<parameter id="randomPopSize.t:anthus_aco_sub" spec="parameter.RealParameter" name="popSize">1.0</parameter>
</populationModel>
</init>
<distribution id="posterior" spec="util.CompoundDistribution">
<distribution id="prior" spec="util.CompoundDistribution">
<distribution id="YuleModel.t:anthus_aco_sub" spec="beast.evolution.speciation.YuleModel" birthDiffRate="@birthRate.t:anthus_aco_sub" tree="@Tree.t:anthus_aco_sub"/>
<prior id="YuleBirthRatePrior.t:anthus_aco_sub" name="distribution" x="@birthRate.t:anthus_aco_sub">
<Uniform id="Uniform.1" name="distr" upper="Infinity"/>
</prior>
</distribution>
<distribution id="likelihood" spec="util.CompoundDistribution" useThreads="true">
<distribution id="treeLikelihood.anthus_aco_sub" spec="ThreadedTreeLikelihood" data="@anthus_aco_sub" tree="@Tree.t:anthus_aco_sub">
<siteModel id="SiteModel.s:anthus_aco_sub" spec="SiteModel">
<parameter id="mutationRate.s:anthus_aco_sub" spec="parameter.RealParameter" estimate="false" name="mutationRate">1.0</parameter>
<parameter id="gammaShape.s:anthus_aco_sub" spec="parameter.RealParameter" estimate="false" name="shape">1.0</parameter>
<parameter id="proportionInvariant.s:anthus_aco_sub" spec="parameter.RealParameter" estimate="false" lower="0.0" name="proportionInvariant" upper="1.0">0.0</parameter>
<substModel id="JC69.s:anthus_aco_sub" spec="JukesCantor"/>
</siteModel>
<branchRateModel id="StrictClock.c:anthus_aco_sub" spec="beast.evolution.branchratemodel.StrictClockModel">
<parameter id="clockRate.c:anthus_aco_sub" spec="parameter.RealParameter" estimate="false" lower="0.1" name="clock.rate" upper="0.3">0.2</parameter>
</branchRateModel>
</distribution>
</distribution>
</distribution>
<operator id="YuleBirthRateScaler.t:anthus_aco_sub" spec="ScaleOperator" parameter="@birthRate.t:anthus_aco_sub" weight="3.0"/>
<operator id="YuleModelTreeScaler.t:anthus_aco_sub" spec="ScaleOperator" scaleFactor="0.5" tree="@Tree.t:anthus_aco_sub" weight="3.0"/>
<operator id="YuleModelTreeRootScaler.t:anthus_aco_sub" spec="ScaleOperator" rootOnly="true" scaleFactor="0.5" tree="@Tree.t:anthus_aco_sub" weight="3.0"/>
<operator id="YuleModelUniformOperator.t:anthus_aco_sub" spec="Uniform" tree="@Tree.t:anthus_aco_sub" weight="30.0"/>
<operator id="YuleModelSubtreeSlide.t:anthus_aco_sub" spec="SubtreeSlide" tree="@Tree.t:anthus_aco_sub" weight="15.0"/>
<operator id="YuleModelNarrow.t:anthus_aco_sub" spec="Exchange" tree="@Tree.t:anthus_aco_sub" weight="15.0"/>
<operator id="YuleModelWide.t:anthus_aco_sub" spec="Exchange" isNarrow="false" tree="@Tree.t:anthus_aco_sub" weight="3.0"/>
<operator id="YuleModelWilsonBalding.t:anthus_aco_sub" spec="WilsonBalding" tree="@Tree.t:anthus_aco_sub" weight="3.0"/>
<logger id="tracelog" spec="Logger" fileName="anthus_aco_sub.log" logEvery="1000" model="@posterior" sanitiseHeaders="true" sort="smart">
<log idref="posterior"/>
<log idref="likelihood"/>
<log idref="prior"/>
<log idref="treeLikelihood.anthus_aco_sub"/>
<log id="TreeHeight.t:anthus_aco_sub" spec="beast.evolution.tree.TreeHeightLogger" tree="@Tree.t:anthus_aco_sub"/>
<log idref="YuleModel.t:anthus_aco_sub"/>
<log idref="birthRate.t:anthus_aco_sub"/>
</logger>
<logger id="screenlog" spec="Logger" logEvery="1000">
<log idref="posterior"/>
<log idref="likelihood"/>
<log idref="prior"/>
</logger>
<logger id="treelog.t:anthus_aco_sub" spec="Logger" fileName="$(tree).trees" logEvery="1000" mode="tree">
<log id="TreeWithMetaDataLogger.t:anthus_aco_sub" spec="beast.evolution.tree.TreeWithMetaDataLogger" tree="@Tree.t:anthus_aco_sub"/>
</logger>
<operatorschedule id="OperatorSchedule" spec="OperatorSchedule"/>
</run>
</beast>
Take special note at ...
<branchRateModel id="StrictClock.c:anthus_aco_sub" spec="beast.evolution.branchratemodel.StrictClockModel">
<parameter id="clockRate.c:anthus_aco_sub" spec="parameter.RealParameter" estimate="false" lower="0.1" name="clock.rate" upper="0.3">0.2</parameter>
</branchRateModel>
... especially the estimate="false"
.
This proves -in a way- that time is normalized to 1.0
from beautier.
Here I write down the solution to this Issue:
Dear user,
Thanks for this excellent bug report. It is clear the problem is (1) you want to estimate a clock rate, (2) the output of BEAST2 does not show an estimated clock rate.
This is not a bug in nor babette
/beautier
, nor BEAST2: it is a feature: BEAST2 normalizes time from time 0 (at the root) to time 1 (the present). This causes -I forgot the details!- that estimating a clock rate is nonsense; it simple is. To find out a more detailed answer, read the BEAST2 book or ask the BEAST2 people; you can use the pictures and XML file above.
I predict you want to, however, get an actual clock rate. The way to do this is simple: if you know the phylogeny is -say- 15 million, then one can multiply/divide things by 15.
Another way is to take a known timespan within the phylogeny and rescale based on that:
+--- A
+---+
---+ +--- B
+------- C
Total phylogeny time: 1.0 (as normalized by BEAST2)
Known time between A and B: 10 million years
Concluded time: Total phylogeny is 20 million years
I hope this answer helps you and future users!
I'll close this Issue if you agree or after a month :-)
from beautier.
Some emails later we've zoomed in on the Issue, with again an awesome bug by @pedrotaucce:
[...] I wrote to you at first because I believe beautier and BEAUTi are resulting in different .xml files, if I am not mistaken.
Earlier I misunderstood this. Thanks Pedro for correcting me
I have built 4 xml files using BEAUTi 2.6.7 and beautier (in R 4.2.2) with different priors:
- No mrca prior and a single value at clock rate (0.00277)
- No mrca prior, estimating clock rate from a uniform prior (starting value: 0.0035, lower: 0.0027, upper: 0.00542)
- Mrca prior with monophyletic constrain and a single value at clock rate (0.00277)
- Mrca prior with monophyletic constrain and a uniform prior (starting value: 0.0035, lower: 0.0027, upper: 0.00542)
BEAUTi and beautier yielded the same xml file in number 1, but differed in the other three scenarios.
In 2, beautier xml file estimated the clock rate with a starting value of 0.0035, but without the specified uniform prior.
In 3 and 4, the clock rate was 1.0 and estimate = false in beautier xml file.
Here the code and BEAUTi printscreens for each scenario:
- beautier:
#without mrca prior, single value at clock rate
fasta_filename <- get_babette_path("anthus_aco_sub.fas")
clock.rate <- beautier::create_clock_rate_param(value = 0.00277,estimate=FALSE)
inference_model <- create_inference_model(
site_model = beautier::create_hky_site_model(),
clock_model = beautier::create_strict_clock_model(id = NA,clock.rate),
tree_prior = create_yule_tree_prior(),
beauti_options = beautier::create_beauti_options_v2_6()
)
create_beast2_input_file_from_model(
input_filename = fasta_filename,
output_filename = "no_mrca_no_estimate_beautier.xml",
inference_model= inference_model
)
BEAUTi:
- beautier:
#134 without mrca prior, estimating clock rate from a uniform prior
fasta_filename <- get_babette_path("anthus_aco_sub.fas")
clock.rate <- beautier::create_clock_rate_param(value = "0.0035",estimate=TRUE)
clock.uniform<-beautier::create_uniform_distr(value = 0.0035,lower = 0.00277, upper = 0.00542)
inference_model <- create_inference_model(
site_model = beautier::create_hky_site_model(),
clock_model = beautier::create_strict_clock_model(id = NA,clock.rate, clock.uniform),
tree_prior = create_yule_tree_prior(),
beauti_options = beautier::create_beauti_options_v2_6()
)
create_beast2_input_file_from_model(
input_filename = fasta_filename,
output_filename = "no_mrca_estimate_beautier.xml",
inference_model= inference_model
)
BEAUTi:
- beautier:
#With mrca prior, single value at clock rate
fasta_filename <- get_babette_path("anthus_aco_sub.fas")
mrca.taxa <- get_taxa_names(fasta_filename)
mrca.taxa <- mrca.taxa[2:length(mrca.taxa)]
mrca.prior <- create_mrca_prior(taxa_names=mrca.taxa,is_monophyletic = T)
clock.rate <- beautier::create_clock_rate_param(value = 0.00277,estimate=FALSE)
inference_model <- create_inference_model(
site_model = beautier::create_hky_site_model(),
clock_model = beautier::create_strict_clock_model(id = NA,clock.rate),
tree_prior = create_yule_tree_prior(),
mrca_prior = mrca.prior,
beauti_options = beautier::create_beauti_options_v2_6()
)
create_beast2_input_file_from_model(
input_filename = fasta_filename,
output_filename = "mrca_no_estimate_beautier.xml",
inference_model= inference_model
)
BEAUTi: same as 1, but with MRCA prior
- beautier:
#With mrca prior, estimating clock rate from a uniform prior
fasta_filename <- get_babette_path("anthus_aco_sub.fas")
mrca.taxa <- get_taxa_names(fasta_filename)
mrca.taxa <- mrca.taxa[2:length(mrca.taxa)]
mrca.prior <- create_mrca_prior(taxa_names=mrca.taxa,is_monophyletic = T)
clock.rate <- beautier::create_clock_rate_param(value = "0.0035",estimate=TRUE)
clock.uniform<-beautier::create_uniform_distr(value = 0.0035,lower = 0.00277, upper = 0.00542)
inference_model <- create_inference_model(
site_model = beautier::create_hky_site_model(),
clock_model = beautier::create_strict_clock_model(id = NA,clock.rate,clock.uniform),
tree_prior = create_yule_tree_prior(),
mrca_prior = mrca.prior,
beauti_options = beautier::create_beauti_options_v2_6()
)
create_beast2_input_file_from_model(
input_filename = fasta_filename,
output_filename = "mrca_estimate_beautier.xml",
inference_model= inference_model
)
BEAUTi: the same as number 2, but with MRCA prior just like number 3.
I ran all the xml files, and numbers 2 and 4 did estimate the clock rate according to the prior:
Also, on page 107 of the BEAST2 book, the authors talk about clock rate priors and apparently there is nothing wrong about it and BEAST2 can deal with that, although I am not using the best prior and a lognormal distribution would be probably better.
I am attaching all the xml files and their names are what I expected them to be. I hope I am not being too annoying and I am sorry if I am wrong or making mistakes regarding BEAUTi and/or beautier. [...]
Nope, please keep sending these awesome bug reports!
from beautier.
I'll first concentrate on 2.
The error is quite clear: this is what beautier
gives:
<branchRateModel id="StrictClock.c:anthus_aco_sub" spec="beast.evolution.branchratemodel.StrictClockModel">
<parameter id="clockRate.c:anthus_aco_sub" spec="parameter.RealParameter" estimate="true" name="clock.rate">0.0035</parameter>
</branchRateModel>
This is what BEAUti gives:
<branchRateModel id="StrictClock.c:anthus_aco_sub" spec="beast.evolution.branchratemodel.StrictClockModel" clock.rate="@clockRate.c:anthus_aco_sub"/>
Note the lack of estimate="true"
in the latter, which is the BEAUti (hence: the truth) text :-/ ?
I will reproduce the BEAUti file with the beautier
-friendly BEAUti, i.e. the one installed by beastierinstall::install_beast2
.
Problem is that setting 1 already fails, so fix that one first :-)
from beautier.
Here I see a first problem:
from beautier.
Related Issues (20)
- Add value, upper and lower HOT 1
- Add missing data doc HOT 2
- Broken link found (/join_next?ref_cta=Sign+up&ref_loc=header+logged+out&ref_page=%2F%3Cuser-name%3E%2F%3Crepo-name%3E&source=header-repo) HOT 1
- Broken link found (/join_next?ref_cta=Sign+up&ref_loc=header+logged+out&ref_page=%2F%3Cuser-name%3E%2F%3Crepo-name%3E&source=header-repo&source_repo=ropensci%2Fbeautier) HOT 1
- Broken link found (/join_next?ref_cta=Sign+up&ref_loc=header+logged+out&ref_page=%2F%3Cuser-name%3E%2F%3Crepo-name%3E&source=header-repo) HOT 1
- Broken link found (/join_next?ref_cta=Sign+up&ref_loc=header+logged+out&ref_page=%2F%3Cuser-name%3E%2F%3Crepo-name%3E&source=header-repo&source_repo=ropensci%2Fbeautier) HOT 1
- Clock rate added twice HOT 2
- Doubling of clockPrior ID HOT 1
- beautier writes and leaves behind in .cache files HOT 2
- Feature request: add offset HOT 14
- Feature request: allow multiple MRCA priors HOT 7
- Fix CRAN issues HOT 3
- Implementing options to run Birth Death Skyline Serial HOT 1
- Incompatibility between BEAST 2.7 and babette HOT 1
- kappa in HKY site model must be a param
- Resulting xml cannot estimate clock rate in a lognormal distribution prior HOT 3
- Bug report: clock rate cannot be set to follow a normal distribution HOT 5
- Add `rate_scaler_factor` to have `scaleFactor` in XML HOT 1
- Fix https://github.com/ropensci/babette/issues/106 HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from beautier.