Giter VIP home page Giter VIP logo

Comments (10)

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024 1

@pedrotaucce will work on this tomorrow :-)

from beautier.

pedrotaucce avatar pedrotaucce commented on August 21, 2024 1

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.

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024

Hi @pedrotaucce, thanks for this excellent bug report (and sorry for replying so late). I will take a look at it now :-)

from beautier.

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024

Plan:

  • Recreate XML in BEAUti: maybe it cannot be done at all
  • If it can be done: compare XML to babette XML

from beautier.

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024

Recreate XML:

  • Install BEAST2 by doing beastierinstall::install_beast2()
  • Loaded the alignments

Screenshot from 2023-02-15 11-57-29

  • Change the clock rate

Screenshot from 2023-02-15 12-04-45

Aha, silly me: a strict clock is always normalized to 1.0 🤦

from beautier.

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024

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.

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024

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.

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024

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:

  1. No mrca prior and a single value at clock rate (0.00277)
  2. No mrca prior, estimating clock rate from a uniform prior (starting value: 0.0035, lower: 0.0027, upper: 0.00542)
  3. Mrca prior with monophyletic constrain and a single value at clock rate (0.00277)
  4. 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:

  1. 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:

image

  1. 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:

image

  1. 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

image

image

  1. 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:

image

issue_135.zip
image

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.

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024

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.

richelbilderbeek avatar richelbilderbeek commented on August 21, 2024

Here I see a first problem:

Screenshot from 2023-03-06 13-03-49

from beautier.

Related Issues (20)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo 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.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.