Examples with JukesCantor beast.evolution.substitutionmodel.JukesCantor used on opensource projects

Example 1 with JukesCantor

use of beast.evolution.substitutionmodel.JukesCantor in project beast2 by CompEvol.

the class ClusterTreeTest method testUPGMA.

public void testUPGMA() throws Exception {
    Alignment alignment = BEASTTestCase.getAlignment();
    JukesCantor JC = new JukesCantor();
    JC.initAndValidate();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("substModel", JC);
    ClusterTree tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    String expectedNewick = "((((human:0.01903085702575253,(chimp:0.008560512208575313,bonobo:0.008560512208575313):0.010470344817177218):0.007962255880644985,gorilla:0.026993112906397516):0.019197419394211015,orangutan:0.04619053230060853):0.007214240662738673,siamang:0.053404772963347204):0.0";
    String actualNewick = tree.getRoot().toNewick();
    assertEquals(expectedNewick, actualNewick);
    // select 3 sequences
    List<Sequence> seqs = new ArrayList<Sequence>();
    seqs.addAll(alignment.sequenceInput.get());
    List<Sequence> newseqs = alignment.sequenceInput.get();
    newseqs.clear();
    newseqs.add(Sequence.getSequenceByTaxon("bonobo", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("chimp", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("human", seqs));
    alignment.initAndValidate();
    tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    expectedNewick = "((bonobo:0.008560512208575313,chimp:0.008560512208575313):0.010470344817177218,human:0.01903085702575253):0.0";
    System.err.println("Seqs:");
    for (Sequence s : seqs) {
        System.err.println(s.taxonInput.get());
    }
    System.err.println("Newseqs:");
    for (Sequence s : newseqs) {
        System.err.println(s.taxonInput.get());
    }
    actualNewick = tree.getRoot().toNewick();
    assertEquals(expectedNewick, actualNewick);
    // same sequences in different order
    newseqs.clear();
    newseqs.add(Sequence.getSequenceByTaxon("human", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("chimp", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("bonobo", seqs));
    alignment.initAndValidate();
    tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    actualNewick = tree.getRoot().toNewick();
    expectedNewick = "(human:0.01903085702575253,(chimp:0.008560512208575313,bonobo:0.008560512208575313):0.010470344817177218):0.0";
    assertEquals(expectedNewick, actualNewick);
}

Example 2 with JukesCantor

use of beast.evolution.substitutionmodel.JukesCantor in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodCaching.

@Test
public void testLikelihoodCaching() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    ConversionGraph acg = new SimulatedACG();
    acg.initByName("rho", 5.0 / locus.getSiteCount(), "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
    State state = new State();
    state.initByName("stateNode", acg);
    state.initialise();
    System.out.println(acg);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", new RealParameter("1"), "substModel", jc);
    // Simulate alignment:
    SimulatedAlignment alignment = new SimulatedAlignment();
    alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
    // Calculate likelihood 1:
    ACGLikelihood argLikelihood = new ACGLikelihood();
    argLikelihood.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    double logP1 = argLikelihood.calculateLogP();
    double logP1prime = argLikelihoodSlow.calculateLogP();
    double relError = 2.0 * Math.abs(logP1 - logP1prime) / Math.abs(logP1 + logP1prime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP1, logP1prime, relError);
    assertTrue(relError < 1e-13);
    // Add a single recombination event
    Node node1 = acg.getExternalNodes().get(0);
    Node node2 = node1.getParent();
    double height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
    double height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
    int startLocus = 500;
    int endLocus = 600;
    Conversion recomb1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(recomb1);
    double logP2 = argLikelihood.calculateLogP();
    double logP2prime = argLikelihoodSlow.calculateLogP();
    relError = 2.0 * Math.abs(logP2 - logP2prime) / Math.abs(logP2 + logP2prime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP2, logP2prime, relError);
    assertTrue(relError < 1e-13);
    state.restore();
    double logP3 = argLikelihood.calculateLogP();
    double logP3prime = argLikelihoodSlow.calculateLogP();
    relError = 2.0 * Math.abs(logP3 - logP3prime) / Math.abs(logP3 + logP3prime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP3, logP3prime, relError);
    assertTrue(relError < 1e-13);
}

Example 3 with JukesCantor

use of beast.evolution.substitutionmodel.JukesCantor in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodAmbiguities.

@Test
public void testLikelihoodAmbiguities() throws Exception {
    Locus locus = new Locus("locus", getAlignment());
    // ConversionGraph
    ConversionGraph acg = new ConversionGraph();
    ClusterTree tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", locus.getAlignment());
    acg.assignFrom(tree);
    acg.initByName("locus", locus);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("substModel", jc);
    // Likelihood
    ACGLikelihood argLikelihood = new ACGLikelihood();
    argLikelihood.initByName("locus", locus, "tree", acg, "siteModel", siteModel, "useAmbiguities", true);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel, "useAmbiguities", true);
    acg.setEverythingDirty(true);
    double logP = argLikelihood.calculateLogP();
    double logPtrue = argLikelihoodSlow.calculateLogP();
    double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    System.out.println("Test 1.  Truth: " + logPtrue + " Value: " + logP);
    assertTrue(relativeDiff < 1e-14);
    // Add a single recombination event
    Node node1 = acg.getExternalNodes().get(0);
    Node node2 = node1.getParent();
    double height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
    double height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
    int startLocus = 100;
    int endLocus = 200;
    Conversion recomb1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(recomb1);
    logP = argLikelihood.calculateLogP();
    logPtrue = argLikelihoodSlow.calculateLogP();
    relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    assertTrue(relativeDiff < 1e-14);
    // Add another recombination event
    node1 = acg.getExternalNodes().get(0);
    node2 = acg.getNode(20);
    height1 = 0.75 * (node1.getHeight() + node1.getParent().getHeight());
    height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
    startLocus = 250;
    endLocus = 300;
    Conversion recomb2 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(recomb2);
    logP = argLikelihood.calculateLogP();
    logPtrue = argLikelihoodSlow.calculateLogP();
    relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
    assertTrue(relativeDiff < 1e-14);
}

Example 4 with JukesCantor

use of beast.evolution.substitutionmodel.JukesCantor in project bacter by tgvaughan.

the class SimulatedAlignmentTest method test.

@Test
public void test() throws Exception {
    // Randomizer.setSeed(26);
    Randomizer.setSeed(7);
    Locus locus = new Locus("locus", 100000);
    TaxonSet taxonSet = getTaxonSet(10);
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    ConversionGraph acg = new SimulatedACG();
    acg.initByName("rho", 1.0 / 100000, "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
    System.out.println(acg);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", new RealParameter("1.0"), "substModel", jc);
    // Simulate alignment:
    SimulatedAlignment alignment = new SimulatedAlignment();
    alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
    for (Region region : acg.getRegions(locus)) System.out.println(new MarginalTree(acg, region));
    // (Should be enough info here for precise agreement)
    for (Region region : acg.getRegions(locus)) {
        Alignment margAlign = createMarginalAlignment(alignment, region);
        ClusterTree upgmaTree = new ClusterTree();
        upgmaTree.initByName("clusterType", "upgma", "taxa", margAlign);
        MarginalTree marginalTree = new MarginalTree(acg, region);
        // outfMarg.println(marginalTree.getRoot() + ";");
        // outfMarg.flush();
        // 
        // outfUPGMA.println(upgmaTree.getRoot() + ";");
        // outfUPGMA.flush();
        assertTrue(topologiesEquivalent(marginalTree.getRoot(), upgmaTree.getRoot()));
    }
// outfMarg.close();
// outfUPGMA.close();
}

Example 5 with JukesCantor

use of beast.evolution.substitutionmodel.JukesCantor in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodUsingSimulatedData.

@Test
public void testLikelihoodUsingSimulatedData() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    ConversionGraph acg = new SimulatedACG();
    acg.initByName("rho", 5.0 / locus.getSiteCount(), "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
    System.out.println(acg);
    // Site model:
    JukesCantor jc = new JukesCantor();
    jc.initByName();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", new RealParameter("1"), "substModel", jc);
    // Simulate alignment:
    SimulatedAlignment alignment = new SimulatedAlignment();
    alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
    // Calculate likelihood:
    ACGLikelihood argLikelihood = new ACGLikelihood();
    argLikelihood.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    double logP = argLikelihood.calculateLogP();
    // Compare product of likelihoods of "marginal alignments" with
    // likelihood computed using RGL.
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "data", alignment, "tree", acg, "siteModel", siteModel);
    double logPprime = argLikelihoodSlow.calculateLogP();
    double relError = 2.0 * Math.abs(logP - logPprime) / Math.abs(logP + logPprime);
    System.out.format("logP=%g\nlogPprime=%g\nrelError=%g\n", logP, logPprime, relError);
    assertTrue(relError < 1e-13);
}