Search in sources :

Example 1 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class CheckPointTreeModifier method incorporateAdditionalTaxa.

/**
     * Add the remaining taxa, which can be identified through the TreeDataLikelihood XML elements.
     */
public ArrayList<NodeRef> incorporateAdditionalTaxa(CheckPointUpdaterApp.UpdateChoice choice, BranchRates rateModel) {
    System.out.println("Tree before adding taxa:\n" + treeModel.toString() + "\n");
    ArrayList<NodeRef> newTaxaNodes = new ArrayList<NodeRef>();
    for (String str : newTaxaNames) {
        for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
            if (treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId().equals(str)) {
                newTaxaNodes.add(treeModel.getExternalNode(i));
                //always take into account Taxon dates vs. dates set through a TreeModel
                System.out.println(treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId() + " with height " + treeModel.getNodeHeight(treeModel.getExternalNode(i)) + " or " + treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getHeight());
            }
        }
    }
    System.out.println("newTaxaNodes length = " + newTaxaNodes.size());
    ArrayList<Taxon> currentTaxa = new ArrayList<Taxon>();
    for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
        boolean taxonFound = false;
        for (String str : newTaxaNames) {
            if (str.equals((treeModel.getNodeTaxon(treeModel.getExternalNode(i))).getId())) {
                taxonFound = true;
            }
        }
        if (!taxonFound) {
            System.out.println("Adding " + treeModel.getNodeTaxon(treeModel.getExternalNode(i)).getId() + " to list of current taxa");
            currentTaxa.add(treeModel.getNodeTaxon(treeModel.getExternalNode(i)));
        }
    }
    System.out.println("Current taxa count = " + currentTaxa.size());
    //iterate over both current taxa and to be added taxa
    boolean originTaxon = true;
    for (Taxon taxon : currentTaxa) {
        if (taxon.getHeight() == 0.0) {
            originTaxon = false;
            System.out.println("Current taxon " + taxon.getId() + " has node height 0.0");
        }
    }
    for (NodeRef newTaxon : newTaxaNodes) {
        if (treeModel.getNodeTaxon(newTaxon).getHeight() == 0.0) {
            originTaxon = false;
            System.out.println("New taxon " + treeModel.getNodeTaxon(newTaxon).getId() + " has node height 0.0");
        }
    }
    //check the Tree(Data)Likelihoods in the connected set of likelihoods
    //focus on TreeDataLikelihood, which has getTree() to get the tree for each likelihood
    //also get the DataLikelihoodDelegate from TreeDataLikelihood
    ArrayList<TreeDataLikelihood> likelihoods = new ArrayList<TreeDataLikelihood>();
    ArrayList<Tree> trees = new ArrayList<Tree>();
    ArrayList<DataLikelihoodDelegate> delegates = new ArrayList<DataLikelihoodDelegate>();
    for (Likelihood likelihood : Likelihood.CONNECTED_LIKELIHOOD_SET) {
        if (likelihood instanceof TreeDataLikelihood) {
            likelihoods.add((TreeDataLikelihood) likelihood);
            trees.add(((TreeDataLikelihood) likelihood).getTree());
            delegates.add(((TreeDataLikelihood) likelihood).getDataLikelihoodDelegate());
        }
    }
    //suggested to go through TreeDataLikelihoodParser and give it an extra option to create a HashMap
    //keyed by the tree; am currently not overly fond of this approach
    ArrayList<PatternList> patternLists = new ArrayList<PatternList>();
    for (DataLikelihoodDelegate del : delegates) {
        if (del instanceof BeagleDataLikelihoodDelegate) {
            patternLists.add(((BeagleDataLikelihoodDelegate) del).getPatternList());
        } else if (del instanceof MultiPartitionDataLikelihoodDelegate) {
            MultiPartitionDataLikelihoodDelegate mpdld = (MultiPartitionDataLikelihoodDelegate) del;
            List<PatternList> list = mpdld.getPatternLists();
            for (PatternList pList : list) {
                patternLists.add(pList);
            }
        }
    }
    if (patternLists.size() == 0) {
        throw new RuntimeException("No patterns detected. Please make sure the XML file is BEAST 1.9 compatible.");
    }
    //aggregate all patterns to create distance matrix
    //TODO What about different trees for different partitions?
    Patterns patterns = new Patterns(patternLists.get(0));
    if (patternLists.size() > 1) {
        for (int i = 1; i < patternLists.size(); i++) {
            patterns.addPatterns(patternLists.get(i));
        }
    }
    //set the patterns for the distance matrix computations
    choice.setPatterns(patterns);
    //add new taxa one at a time
    System.out.println("Adding " + newTaxaNodes.size() + " taxa ...");
    for (NodeRef newTaxon : newTaxaNodes) {
        treeModel.setNodeHeight(newTaxon, treeModel.getNodeTaxon(newTaxon).getHeight());
        System.out.println("\nadding Taxon: " + newTaxon + " (height = " + treeModel.getNodeHeight(newTaxon) + ")");
        //check if this taxon has a more recent sampling date than all other nodes in the current TreeModel
        double offset = checkCurrentTreeNodes(newTaxon, treeModel.getRoot());
        System.out.println("Sampling date offset when adding " + newTaxon + " = " + offset);
        //AND set its current node height to 0.0 IF no originTaxon has been found
        if (offset < 0.0) {
            if (!originTaxon) {
                System.out.println("Updating all node heights with offset " + Math.abs(offset));
                updateAllTreeNodes(Math.abs(offset), treeModel.getRoot());
                treeModel.setNodeHeight(newTaxon, 0.0);
            }
        } else if (offset == 0.0) {
            if (!originTaxon) {
                treeModel.setNodeHeight(newTaxon, 0.0);
            }
        }
        //get the closest Taxon to the Taxon that needs to be added
        //take into account which taxa can currently be chosen
        Taxon closest = choice.getClosestTaxon(treeModel.getNodeTaxon(newTaxon), currentTaxa);
        System.out.println("\nclosest Taxon: " + closest + " with original height: " + closest.getHeight());
        //get the distance between these two taxa
        double distance = choice.getDistance(treeModel.getNodeTaxon(newTaxon), closest);
        System.out.println("at distance: " + distance);
        //TODO what if distance == 0.0 ??? how to choose closest taxon then (in absence of geo info)?
        //find the NodeRef for the closest Taxon (do not rely on node numbering)
        NodeRef closestRef = null;
        //careful with node numbering and subtract number of new taxa
        for (int i = 0; i < treeModel.getExternalNodeCount(); i++) {
            if (treeModel.getNodeTaxon(treeModel.getExternalNode(i)) == closest) {
                closestRef = treeModel.getExternalNode(i);
            }
        }
        System.out.println(closestRef + " with height " + treeModel.getNodeHeight(closestRef));
        //System.out.println("trying to set node height: " + closestRef + " from " + treeModel.getNodeHeight(closestRef) + " to " + closest.getHeight());
        //treeModel.setNodeHeight(closestRef, closest.getHeight());
        double timeForDistance = distance / rateModel.getBranchRate(treeModel, closestRef);
        System.out.println("timeForDistance = " + timeForDistance);
        //get parent node of branch that will be split
        NodeRef parent = treeModel.getParent(closestRef);
        //determine height of new node
        double insertHeight;
        if (treeModel.getNodeHeight(closestRef) == treeModel.getNodeHeight(newTaxon)) {
            insertHeight = treeModel.getNodeHeight(closestRef) + timeForDistance / 2.0;
            System.out.println("treeModel.getNodeHeight(closestRef) == treeModel.getNodeHeight(newTaxon): " + insertHeight);
            if (insertHeight >= treeModel.getNodeHeight(parent)) {
                insertHeight = treeModel.getNodeHeight(closestRef) + EPSILON * (treeModel.getNodeHeight(parent) - treeModel.getNodeHeight(closestRef));
            }
        } else {
            double remainder = (timeForDistance - Math.abs(treeModel.getNodeHeight(closestRef) - treeModel.getNodeHeight(newTaxon))) / 2.0;
            if (remainder > 0) {
                insertHeight = Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon)) + remainder;
                System.out.println("remainder > 0: " + insertHeight);
                if (insertHeight >= treeModel.getNodeHeight(parent)) {
                    insertHeight = treeModel.getNodeHeight(closestRef) + EPSILON * (treeModel.getNodeHeight(parent) - treeModel.getNodeHeight(closestRef));
                }
            } else {
                insertHeight = EPSILON * (treeModel.getNodeHeight(parent) - Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon)));
                insertHeight += Math.max(treeModel.getNodeHeight(closestRef), treeModel.getNodeHeight(newTaxon));
                System.out.println("remainder <= 0: " + insertHeight);
            }
        }
        System.out.println("insert at height: " + insertHeight);
        //pass on all the necessary variables to a method that adds the new taxon to the tree
        addTaxonAlongBranch(newTaxon, parent, closestRef, insertHeight);
        //option to print tree after each taxon addition
        System.out.println("\nTree after adding taxon " + newTaxon + ":\n" + treeModel.toString());
        //add newly added Taxon to list of current taxa
        currentTaxa.add(treeModel.getNodeTaxon(newTaxon));
    }
    return newTaxaNodes;
}
Also used : Likelihood(dr.inference.model.Likelihood) TreeDataLikelihood(dr.evomodel.treedatalikelihood.TreeDataLikelihood) Taxon(dr.evolution.util.Taxon) ArrayList(java.util.ArrayList) PatternList(dr.evolution.alignment.PatternList) BeagleDataLikelihoodDelegate(dr.evomodel.treedatalikelihood.BeagleDataLikelihoodDelegate) NodeRef(dr.evolution.tree.NodeRef) TreeDataLikelihood(dr.evomodel.treedatalikelihood.TreeDataLikelihood) BeagleDataLikelihoodDelegate(dr.evomodel.treedatalikelihood.BeagleDataLikelihoodDelegate) MultiPartitionDataLikelihoodDelegate(dr.evomodel.treedatalikelihood.MultiPartitionDataLikelihoodDelegate) DataLikelihoodDelegate(dr.evomodel.treedatalikelihood.DataLikelihoodDelegate) Tree(dr.evolution.tree.Tree) PatternList(dr.evolution.alignment.PatternList) ArrayList(java.util.ArrayList) List(java.util.List) MultiPartitionDataLikelihoodDelegate(dr.evomodel.treedatalikelihood.MultiPartitionDataLikelihoodDelegate) Patterns(dr.evolution.alignment.Patterns)

Example 2 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class BeagleTreeLikelihood method main.

public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 1: simulateOnePartition");
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        Parameter kappa2 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        HKY hky2 = new HKY(kappa2, freqModel);
        HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        substitutionModels.add(hky2);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        Parameter epochTimes = new Parameter.Default(1, 20);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 0.001);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogenousBranchSubstitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        Alignment alignment = simulator.simulate(false, false);
        BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
        nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Partition(dr.app.beagle.tools.Partition) EpochBranchModel(dr.evomodel.branchmodel.EpochBranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) BranchModel(dr.evomodel.branchmodel.BranchModel) EpochBranchModel(dr.evomodel.branchmodel.EpochBranchModel) StrictClockBranchRates(dr.evomodel.branchratemodel.StrictClockBranchRates) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) MarkovModulatedSubstitutionModel(dr.evomodel.substmodel.MarkovModulatedSubstitutionModel) BeagleSequenceSimulator(dr.app.beagle.tools.BeagleSequenceSimulator) TreeModel(dr.evomodel.tree.TreeModel) Alignment(dr.evolution.alignment.Alignment) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) NewickImporter(dr.evolution.io.NewickImporter) HKY(dr.evomodel.substmodel.nucleotide.HKY) Tree(dr.evolution.tree.Tree) Parameter(dr.inference.model.Parameter)

Example 3 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class TiterImporter method importTrees.

/**
     * importTrees.
     */
public Tree[] importTrees(TaxonList taxonList) throws IOException, ImportException {
    boolean done = false;
    ArrayList<FlexibleTree> array = new ArrayList<FlexibleTree>();
    do {
        try {
            skipUntil("(");
            unreadCharacter('(');
            FlexibleNode root = readInternalNode(taxonList);
            FlexibleTree tree = new FlexibleTree(root, false, true);
            array.add(tree);
            if (taxonList == null) {
                taxonList = tree;
            }
            if (readCharacter() != ';') {
                throw new BadFormatException("Expecting ';' after tree");
            }
        } catch (EOFException e) {
            done = true;
        }
    } while (!done);
    Tree[] trees = new Tree[array.size()];
    array.toArray(trees);
    return trees;
}
Also used : FlexibleTree(dr.evolution.tree.FlexibleTree) FlexibleNode(dr.evolution.tree.FlexibleNode) ArrayList(java.util.ArrayList) EOFException(java.io.EOFException) FlexibleTree(dr.evolution.tree.FlexibleTree) Tree(dr.evolution.tree.Tree)

Example 4 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class OldTreesPanel method createTree.

private void createTree() {
    if (generateTreeDialog == null) {
        generateTreeDialog = new GenerateTreeDialog(frame);
    }
    int result = generateTreeDialog.showDialog(options);
    if (result != JOptionPane.CANCEL_OPTION) {
        GenerateTreeDialog.MethodTypes methodType = generateTreeDialog.getMethodType();
        PartitionData partition = generateTreeDialog.getDataPartition();
        Patterns patterns = new Patterns(partition.getAlignment());
        DistanceMatrix distances = new F84DistanceMatrix(patterns);
        Tree tree;
        TemporalRooting temporalRooting;
        switch(methodType) {
            case NJ:
                tree = new NeighborJoiningTree(distances);
                temporalRooting = new TemporalRooting(tree);
                tree = temporalRooting.findRoot(tree, TemporalRooting.RootingFunction.CORRELATION);
                break;
            case UPGMA:
                tree = new UPGMATree(distances);
                temporalRooting = new TemporalRooting(tree);
                break;
            default:
                throw new IllegalArgumentException("unknown method type");
        }
        tree.setId(generateTreeDialog.getName());
        options.userTrees.add(tree);
        treesTableModel.fireTableDataChanged();
        int row = options.userTrees.size() - 1;
        treesTable.getSelectionModel().setSelectionInterval(row, row);
    }
    fireTreePriorsChanged();
}
Also used : F84DistanceMatrix(dr.evolution.distance.F84DistanceMatrix) NeighborJoiningTree(dr.evolution.tree.NeighborJoiningTree) PartitionData(dr.app.beauti.options.PartitionData) UPGMATree(dr.evolution.tree.UPGMATree) NeighborJoiningTree(dr.evolution.tree.NeighborJoiningTree) Tree(dr.evolution.tree.Tree) TemporalRooting(dr.app.pathogen.TemporalRooting) F84DistanceMatrix(dr.evolution.distance.F84DistanceMatrix) DistanceMatrix(dr.evolution.distance.DistanceMatrix) UPGMATree(dr.evolution.tree.UPGMATree) Patterns(dr.evolution.alignment.Patterns)

Example 5 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class TimeSlicer method readAndAnalyzeTrees.

// private List<Tree> importTrees(String treeFileName, int burnin) throws IOException, Importer.ImportException {
// 
// int totalTrees = 10000;
// 
// progressStream.println("Reading trees (bar assumes 10,000 trees)...");
// progressStream.println("0              25             50             75            100");
// progressStream.println("|--------------|--------------|--------------|--------------|");
// 
// int stepSize = totalTrees / 60;
// if (stepSize < 1) stepSize = 1;
// 
// List<Tree> treeList = new ArrayList<Tree>();
// BufferedReader reader1 = new BufferedReader(new FileReader(treeFileName));
// 
// String line1 = reader1.readLine();
// TreeImporter importer1;
// if (line1.toUpperCase().startsWith("#NEXUS")) {
// importer1 = new NexusImporter(new FileReader(treeFileName));
// } else {
// importer1 = new NewickImporter(new FileReader(treeFileName));
// }
// totalTrees = 0;
// while (importer1.hasTree()) {
// Tree treeTime = importer1.importNextTree();
// 
// if (totalTrees > burnin)
// treeList.add(treeTime);
// 
// if (totalTrees > 0 && totalTrees % stepSize == 0) {
// progressStream.print("*");
// progressStream.flush();
// }
// totalTrees++;
// }
// return treeList;
// }
private void readAndAnalyzeTrees(String treeFileName, int burnin, int skipEvery, String[] traits, double[] slices, boolean impute, boolean trueNoise, Normalization normalize, boolean divideByBranchLength, BranchSet branchset, Set taxaSet) throws IOException, Importer.ImportException {
    int totalTrees = 10000;
    int totalStars = 0;
    progressStream.println("Reading and analyzing trees (bar assumes 10,000 trees)...");
    progressStream.println("0              25             50             75            100");
    progressStream.println("|--------------|--------------|--------------|--------------|");
    int stepSize = totalTrees / 60;
    if (stepSize < 1)
        stepSize = 1;
    BufferedReader reader1 = new BufferedReader(new FileReader(treeFileName));
    String line1 = reader1.readLine();
    TreeImporter importer1;
    if (line1.toUpperCase().startsWith("#NEXUS")) {
        importer1 = new NexusImporter(new FileReader(treeFileName));
    } else {
        importer1 = new NewickImporter(new FileReader(treeFileName));
    }
    totalTrees = 0;
    while (importer1.hasTree()) {
        Tree treeTime = importer1.importNextTree();
        if (totalTrees % skipEvery == 0) {
            treesRead++;
            if (totalTrees >= burnin) {
                analyzeTree(treeTime, traits, slices, impute, trueNoise, normalize, divideByBranchLength, branchset, taxaSet);
            }
        }
        if (totalTrees > 0 && totalTrees % stepSize == 0) {
            progressStream.print("*");
            totalStars++;
            if (totalStars % 61 == 0)
                progressStream.print("\n");
            progressStream.flush();
        }
        totalTrees++;
    }
    progressStream.print("\n");
}
Also used : NexusImporter(dr.evolution.io.NexusImporter) NewickImporter(dr.evolution.io.NewickImporter) TreeImporter(dr.evolution.io.TreeImporter) Tree(dr.evolution.tree.Tree)

Aggregations

Tree (dr.evolution.tree.Tree)167 NewickImporter (dr.evolution.io.NewickImporter)43 ArrayList (java.util.ArrayList)35 Parameter (dr.inference.model.Parameter)27 DefaultTreeModel (dr.evomodel.tree.DefaultTreeModel)22 NexusImporter (dr.evolution.io.NexusImporter)20 TaxonList (dr.evolution.util.TaxonList)20 Taxa (dr.evolution.util.Taxa)19 Taxon (dr.evolution.util.Taxon)19 FlexibleTree (dr.evolution.tree.FlexibleTree)18 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)18 NodeRef (dr.evolution.tree.NodeRef)17 Importer (dr.evolution.io.Importer)15 TreeModel (dr.evomodel.tree.TreeModel)13 ImportException (dr.evolution.io.Importer.ImportException)12 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)12 SimpleTree (dr.evolution.tree.SimpleTree)11 TreeUtils (dr.evolution.tree.TreeUtils)11 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)11 IOException (java.io.IOException)11