Examples with Locus bacter.Locus used on opensource projects

Example 1 with Locus

use of bacter.Locus 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 2 with Locus

use of bacter.Locus 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 3 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class AddRemoveConversionTest method testHR.

/**
 * Tests that probability density of forward move calculated
 * by drawNewConversion() matches probability density of backward
 * move calculated by getConversionProb().
 *
 * @throws Exception
 */
@Test
public void testHR() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    SimulatedACG acg = new SimulatedACG();
    acg.initByName("rho", 1.0 / locus.getSiteCount(), "delta", 50.0, "locus", locus, "taxonset", taxonSet, "populationModel", popFunc);
    AddRemoveConversion operator = new AddRemoveConversion();
    // Loop until a valid proposal is made
    double logP1;
    List<Conversion> oldConversions;
    do {
        operator.initByName("weight", 1.0, "acg", acg, "delta", new RealParameter("50.0"), "populationModel", popFunc);
        oldConversions = Lists.newArrayList(acg.getConversions(locus));
        logP1 = operator.drawNewConversion();
    } while (Double.isInfinite(logP1));
    System.out.println("logP1 = " + logP1);
    // Identify new recomination
    Conversion newRecomb = null;
    for (Conversion recomb : acg.getConversions(locus)) {
        if (!oldConversions.contains(recomb))
            newRecomb = recomb;
    }
    assertNotNull(newRecomb);
    double logP2 = operator.getConversionProb(newRecomb);
    System.out.println("logP2 = " + logP2);
    assertTrue(Math.abs(logP1 - logP2) < 1e-10);
}

Example 4 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class SimulatedACG method initAndValidate.

@Override
public void initAndValidate() {
    rho = rhoInput.get();
    delta = deltaInput.get();
    popFunc = popFuncInput.get();
    // Need to do this here as Tree.processTraits(), which is called
    // by hasDateTrait() and hence simulateClonalFrame(), expects a
    // tree with nodes.
    super.initAndValidate();
    if (clonalFrameInput.get() == null)
        simulateClonalFrame();
    else
        assignFromWithoutID(clonalFrameInput.get());
    // Need to do this here as this sets the tree object that the nodes
    // point to, so without it they point to the dummy tree created by
    // super.initAndValidate().
    initArrays();
    // Generate recombinations
    generateConversions();
    // Write output file
    if (outputFileNameInput.get() != null) {
        NexusBuilder nexusBuilder = new NexusBuilder();
        nexusBuilder.append(new TaxaBlock(m_taxonset.get()));
        nexusBuilder.append((new TreesBlock() {

            @Override
            public String getTreeString(Tree tree) {
                return ((ConversionGraph) tree).getExtendedNewick();
            }
        }).addTree(this, "simulatedARG"));
        nexusBuilder.append(new NexusBlock() {

            @Override
            public String getBlockName() {
                return "bacter";
            }

            @Override
            public List<String> getBlockLines() {
                List<String> lines = new ArrayList<>();
                String lociLine = "loci";
                for (Locus locus : getLoci()) lociLine += " " + locus.getID() + ":" + locus.getSiteCount();
                lines.add(lociLine);
                lines.add("clonalframe_labeled " + root.toNewick());
                lines.add("clonalframe_numbered " + root.toShortNewick(true));
                for (Locus locus : getLoci()) {
                    for (Conversion conv : getConversions(locus)) {
                        lines.add("conversion node1=" + conv.getNode1().getNr() + " node2=" + conv.getNode2().getNr() + " site1=" + conv.getStartSite() + " site2=" + conv.getEndSite() + " locus=" + locus.getID());
                    }
                }
                return lines;
            }
        });
        try (PrintStream pstream = new PrintStream(outputFileNameInput.get())) {
            nexusBuilder.write(pstream);
        } catch (FileNotFoundException e) {
            e.printStackTrace();
        }
    }
}

Example 5 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ACGOperator method connectEdge.

/**
 * Connect edge <node, node.parent> above destEdgeBase, forming new
 * edge <destEdgeBase, node.parent> and <node.parent, destEdgeBase.parent>.
 * All conversions above destEdgeBase that are older than destTime
 * are transferred to the new edge above node.parent.
 *
 * Conversions on edge above node are not modified.
 *
 * @param node base of edge to attach
 * @param destEdgeBase base of edge to be bisected
 * @param destTime height at which bisection will occur
 */
protected void connectEdge(Node node, Node destEdgeBase, double destTime) {
    if (node.isRoot())
        throw new IllegalArgumentException("Programmer error: " + "root argument passed to connectEdge().");
    Node parent = node.getParent();
    if (destEdgeBase.isRoot()) {
        parent.addChild(destEdgeBase);
    } else {
        Node grandParent = destEdgeBase.getParent();
        grandParent.removeChild(destEdgeBase);
        grandParent.addChild(parent);
        parent.addChild(destEdgeBase);
    }
    parent.setHeight(destTime);
    for (Locus locus : acg.getLoci()) {
        for (Conversion conv : acg.getConversions(locus)) {
            if (conv.getNode1() == destEdgeBase && conv.getHeight1() > destTime)
                conv.setNode1(parent);
            if (conv.getNode2() == destEdgeBase && conv.getHeight2() > destTime)
                conv.setNode2(parent);
        }
    }
}