Examples with Node beast.evolution.tree.Node used on opensource projects

Example 66 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ConversionGraph method fromExtendedNewick.

/**
 * Read in an ACG from a string in extended newick format.  Assumes
 * that the network is stored with exactly the same metadata as written
 * by the getExtendedNewick() method.
 *
 * @param string extended newick representation of ACG
 * @param numbered true indicates that the ACG is numbered.
 */
public void fromExtendedNewick(String string, boolean numbered, int nodeNumberoffset) {
    // Spin up ANTLR
    CharStream input = CharStreams.fromString(string);
    ExtendedNewickLexer lexer = new ExtendedNewickLexer(input);
    CommonTokenStream tokens = new CommonTokenStream(lexer);
    ExtendedNewickParser parser = new ExtendedNewickParser(tokens);
    ParseTree parseTree = parser.tree();
    Map<String, Conversion> convIDMap = new HashMap<>();
    Node root = new ExtendedNewickBaseVisitor<Node>() {

        /**
         * Convert branch lengths to node heights for all nodes in clade.
         *
         * @param node clade parent
         * @return minimum height assigned in clade.
         */
        private double branchLengthsToHeights(Node node) {
            if (node.isRoot())
                node.setHeight(0.0);
            else
                node.setHeight(node.getParent().getHeight() - node.getHeight());
            double minHeight = node.getHeight();
            for (Node child : node.getChildren()) {
                minHeight = Math.min(minHeight, branchLengthsToHeights(child));
            }
            return minHeight;
        }

        /**
         * Remove height offset from all nodes in clade
         * @param node parent of clade
         * @param offset offset to remove
         */
        private void removeOffset(Node node, double offset) {
            node.setHeight(node.getHeight() - offset);
            for (Node child : node.getChildren()) removeOffset(child, offset);
        }

        private Node getTrueNode(Node node) {
            if (node.isLeaf()) {
                assert !convIDMap.containsKey(node.getID());
                return node;
            }
            if (convIDMap.containsKey(node.getID()))
                return getTrueNode(node.getChild(0));
            int hybridIdx = -1;
            int nonHybridIdx = -1;
            for (int i = 0; i < node.getChildCount(); i++) {
                if (node.getChild(i).isLeaf() && convIDMap.containsKey(node.getChild(i).getID()))
                    hybridIdx = i;
                else
                    nonHybridIdx = i;
            }
            if (hybridIdx > 0)
                return getTrueNode(node.getChild(nonHybridIdx));
            return node;
        }

        /**
         * Traverse the newly constructed tree looking for
         * hybrid nodes and using these to set the heights of
         * Conversion objects.
         *
         * @param node parent of clade
         */
        private void findConversionAttachments(Node node) {
            if (convIDMap.containsKey(node.getID())) {
                Conversion conv = convIDMap.get(node.getID());
                if (node.isLeaf()) {
                    conv.setHeight1(node.getHeight());
                    conv.setHeight2(node.getParent().getHeight());
                    conv.setNode2(getTrueNode(node.getParent()));
                } else
                    conv.setNode1(getTrueNode(node));
            }
            for (Node child : node.getChildren()) findConversionAttachments(child);
        }

        /**
         * Remove all conversion-associated nodes, leaving only
         * the clonal frame.
         *
         * @param node parent of clade
         * @return new parent of same clade
         */
        private Node stripHybridNodes(Node node) {
            Node trueNode = getTrueNode(node);
            List<Node> trueChildren = new ArrayList<>();
            for (Node child : trueNode.getChildren()) {
                trueChildren.add(stripHybridNodes(child));
            }
            trueNode.removeAllChildren(false);
            for (Node trueChild : trueChildren) trueNode.addChild(trueChild);
            return trueNode;
        }

        private int numberInternalNodes(Node node, int nextNr) {
            if (node.isLeaf())
                return nextNr;
            for (Node child : node.getChildren()) nextNr = numberInternalNodes(child, nextNr);
            node.setNr(nextNr);
            return nextNr + 1;
        }

        @Override
        public Node visitTree(ExtendedNewickParser.TreeContext ctx) {
            Node root = visitNode(ctx.node());
            double minHeight = branchLengthsToHeights(root);
            removeOffset(root, minHeight);
            findConversionAttachments(root);
            root = stripHybridNodes(root);
            root.setParent(null);
            if (!numbered)
                numberInternalNodes(root, root.getAllLeafNodes().size());
            return root;
        }

        @Override
        public Node visitNode(ExtendedNewickParser.NodeContext ctx) {
            Node node = new Node();
            if (ctx.post().hybrid() != null) {
                String convID = ctx.post().hybrid().getText();
                node.setID(convID);
                Conversion conv;
                if (convIDMap.containsKey(convID))
                    conv = convIDMap.get(convID);
                else {
                    conv = new Conversion();
                    convIDMap.put(convID, conv);
                }
                if (ctx.node().isEmpty()) {
                    String locusID;
                    for (ExtendedNewickParser.AttribContext attribCtx : ctx.post().meta().attrib()) {
                        switch(attribCtx.attribKey.getText()) {
                            case "region":
                                conv.setStartSite(Integer.parseInt(attribCtx.attribValue().vector().attribValue(0).getText()));
                                conv.setEndSite(Integer.parseInt(attribCtx.attribValue().vector().attribValue(1).getText()));
                                break;
                            case "locus":
                                locusID = attribCtx.attribValue().getText();
                                if (locusID.startsWith("\""))
                                    locusID = locusID.substring(1, locusID.length() - 1);
                                Locus locus = null;
                                for (Locus thisLocus : getLoci()) {
                                    if (thisLocus.getID().equals(locusID))
                                        locus = thisLocus;
                                }
                                if (locus == null)
                                    throw new IllegalArgumentException("Locus with ID " + locusID + " not found.");
                                conv.setLocus(locus);
                                break;
                            default:
                                break;
                        }
                    }
                }
            }
            for (ExtendedNewickParser.NodeContext childCtx : ctx.node()) node.addChild(visitNode(childCtx));
            if (ctx.post().label() != null) {
                node.setID(ctx.post().label().getText());
                node.setNr(Integer.parseInt(ctx.post().label().getText()) - nodeNumberoffset);
            }
            node.setHeight(Double.parseDouble(ctx.post().length.getText()));
            return node;
        }
    }.visit(parseTree);
    m_nodes = root.getAllChildNodesAndSelf().toArray(m_nodes);
    nodeCount = m_nodes.length;
    leafNodeCount = root.getAllLeafNodes().size();
    setRoot(root);
    initArrays();
    for (Locus locus : getLoci()) convs.get(locus).clear();
    for (Conversion conv : convIDMap.values()) addConversion(conv);
}

Example 67 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class SerializationDeserializationTest method testExtendedNewick.

@Test
public void testExtendedNewick() throws Exception {
    Alignment alignment = getAlignment();
    alignment.setID("alignment");
    Locus locus = new Locus("locus", alignment.getSiteCount());
    // ConversionGraph
    ConversionGraph acg = new ConversionGraph();
    ClusterTree tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    acg.assignFrom(tree);
    acg.initByName("locus", locus);
    // Add recombination event 1
    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 conv1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(conv1);
    // Add recombination event 2
    node1 = acg.getExternalNodes().get(8);
    node2 = acg.getRoot();
    height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
    height2 = node2.getHeight() + 1.0;
    startLocus = 300;
    endLocus = 400;
    Conversion conv2 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
    acg.addConversion(conv2);
    String newickString = acg.getExtendedNewick();
    System.out.println(newickString);
    ConversionGraph acgNew = new ConversionGraph();
    acgNew.initByName("locus", locus);
    acgNew.fromExtendedNewick(newickString);
    // Check that new ACG matches old
    Conversion newConv1 = acgNew.getConversions(locus).get(0);
    assertEquals(newConv1.getNode1().getNr(), conv1.getNode1().getNr());
    assertEquals(newConv1.getNode2().getNr(), conv1.getNode2().getNr());
    assertEquals(newConv1.getHeight1(), conv1.getHeight1(), 1e-15);
    assertEquals(newConv1.getHeight2(), conv1.getHeight2(), 1e-15);
    assertEquals(newConv1.getStartSite(), conv1.getStartSite());
    assertEquals(newConv1.getEndSite(), conv1.getEndSite());
    Conversion newConv2 = acgNew.getConversions(locus).get(1);
    assertEquals(newConv2.getNode1().getNr(), conv2.getNode1().getNr());
    assertEquals(newConv2.getNode2().getNr(), conv2.getNode2().getNr());
    assertEquals(newConv2.getHeight1(), conv2.getHeight1(), 1e-15);
    assertEquals(newConv2.getHeight2(), conv2.getHeight2(), 1e-15);
    assertEquals(newConv2.getStartSite(), conv2.getStartSite());
    assertEquals(newConv2.getEndSite(), conv2.getEndSite());
}

Example 68 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ACGLikelihoodTest method testBeagleLikelihood.

@Test
public void testBeagleLikelihood() 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
    ACGLikelihoodBeagle argLikelihood = new ACGLikelihoodBeagle();
    argLikelihood.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
    acg.setEverythingDirty(true);
    try {
        double logP = argLikelihood.calculateLogP();
        double logPtrue = argLikelihoodSlow.calculateLogP();
        double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + 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);
    } catch (RuntimeException ex) {
        System.err.println("Beagle library not found: skipping beagle likelihood test.");
    }
}

Example 69 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class ACGLikelihoodTest method testLikelihoodFixedData.

@Test
public void testLikelihoodFixedData() 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);
    ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
    argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
    acg.setEverythingDirty(true);
    double logP = argLikelihood.calculateLogP();
    double logPtrue = argLikelihoodSlow.calculateLogP();
    double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + 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 70 with Node

use of beast.evolution.tree.Node in project bacter by tgvaughan.

the class DifferenceFromTrueACG method getClades.

public static Clade getClades(Clade[] clades, Node node) {
    Clade clade = new Clade();
    clade.age = node.getHeight();
    if (node.isLeaf()) {
        clade.set(node.getNr());
    } else {
        for (Node child : node.getChildren()) clade.or(getClades(clades, child));
    }
    clades[node.getNr()] = clade;
    return clade;
}