Search in sources :

Example 6 with Partition

use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.

the class BeagleSequenceSimulatorParser method parseXMLObject.

// END: getSyntaxRules
@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    String msg = "";
    boolean parallel = false;
    boolean outputAncestralSequences = false;
    if (xo.hasAttribute(PARALLEL)) {
        parallel = xo.getBooleanAttribute(PARALLEL);
    }
    if (xo.hasAttribute(OUTPUT_ANCESTRAL_SEQUENCES)) {
        outputAncestralSequences = xo.getBooleanAttribute(OUTPUT_ANCESTRAL_SEQUENCES);
    }
    SimpleAlignment.OutputType output = SimpleAlignment.OutputType.FASTA;
    if (xo.hasAttribute(OUTPUT)) {
        output = SimpleAlignment.OutputType.parseFromString(xo.getStringAttribute(OUTPUT));
    }
    int siteCount = 0;
    int to = 0;
    for (int i = 0; i < xo.getChildCount(); i++) {
        Partition partition = (Partition) xo.getChild(i);
        to = partition.to + 1;
        if (to > siteCount) {
            siteCount = to;
        }
    }
    // END: partitions loop
    ArrayList<Partition> partitionsList = new ArrayList<Partition>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        Partition partition = (Partition) xo.getChild(i);
        if (partition.from > siteCount) {
            throw new XMLParseException("Illegal 'from' attribute in " + PartitionParser.PARTITION + " element");
        }
        if (partition.to > siteCount) {
            throw new XMLParseException("Illegal 'to' attribute in " + PartitionParser.PARTITION + " element");
        }
        if (partition.to == -1) {
            partition.to = siteCount - 1;
        }
        if (partition.getRootSequence() != null) {
            //            	TODO: what about 'every'?
            int partitionSiteCount = (partition.to - partition.from) + 1;
            if (partition.getRootSequence().getLength() != 3 * partitionSiteCount && partition.getFreqModel().getDataType() instanceof Codons) {
                throw new RuntimeException("Root codon sequence " + "for partition " + (i + 1) + " has " + partition.getRootSequence().getLength() + " characters " + "expecting " + 3 * partitionSiteCount + " characters");
            } else if (partition.getRootSequence().getLength() != partitionSiteCount && partition.getFreqModel().getDataType() instanceof Nucleotides) {
                throw new RuntimeException("Root nuleotide sequence " + "for partition " + (i + 1) + " has " + partition.getRootSequence().getLength() + " characters " + "expecting " + partitionSiteCount + " characters");
            }
        // END: dataType check
        //                System.exit(-1);
        }
        // END: ancestralSequence check
        partitionsList.add(partition);
    }
    // END: partitions loop
    msg += "\n\t" + siteCount + ((siteCount > 1) ? " replications " : " replication");
    if (msg.length() > 0) {
        Logger.getLogger("dr.app.beagle.tools").info("Using Beagle Sequence Simulator: " + msg);
    }
    BeagleSequenceSimulator s = new BeagleSequenceSimulator(partitionsList);
    SimpleAlignment alignment = s.simulate(parallel, outputAncestralSequences);
    alignment.setOutputType(output);
    return alignment;
}
Also used : Partition(dr.app.beagle.tools.Partition) ArrayList(java.util.ArrayList) Nucleotides(dr.evolution.datatype.Nucleotides) BeagleSequenceSimulator(dr.app.beagle.tools.BeagleSequenceSimulator) SimpleAlignment(dr.evolution.alignment.SimpleAlignment) Codons(dr.evolution.datatype.Codons)

Example 7 with Partition

use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.

the class PartitionParser method parseXMLObject.

@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    int from = 0;
    int to = -1;
    int every = xo.getAttribute(EVERY, 1);
    if (xo.hasAttribute(FROM)) {
        from = xo.getIntegerAttribute(FROM) - 1;
        if (from < 0) {
            throw new XMLParseException("Illegal 'from' attribute in patterns element");
        }
    }
    if (xo.hasAttribute(TO)) {
        to = xo.getIntegerAttribute(TO) - 1;
        if (to < 0 || to < from) {
            throw new XMLParseException("Illegal 'to' attribute in patterns element");
        }
    }
    if (every <= 0) {
        throw new XMLParseException("Illegal 'every' attribute in patterns element");
    }
    // END: every check
    TreeModel tree = (TreeModel) xo.getChild(TreeModel.class);
    GammaSiteRateModel siteModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    FrequencyModel freqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    Sequence rootSequence = (Sequence) xo.getChild(Sequence.class);
    BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (rateModel == null) {
        rateModel = new DefaultBranchRateModel();
    }
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        branchModel = new HomogeneousBranchModel(substitutionModel);
    }
    Partition partition = new Partition(tree, branchModel, siteModel, rateModel, freqModel, from, to, every);
    if (rootSequence != null) {
        partition.setRootSequence(rootSequence);
    }
    return partition;
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) TreeModel(dr.evomodel.tree.TreeModel) Partition(dr.app.beagle.tools.Partition) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) Sequence(dr.evolution.sequence.Sequence) XMLParseException(dr.xml.XMLParseException) BranchModel(dr.evomodel.branchmodel.BranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) SubstitutionModel(dr.evomodel.substmodel.SubstitutionModel) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel)

Example 8 with Partition

use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateCodon.

// END: simulateAminoAcid
static void simulateCodon() {
    try {
        boolean calculateLikelihood = true;
        System.out.println("Test case 6: simulate codons");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        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 site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(Utils.UNIFORM_CODON_FREQUENCIES);
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create substitution model
        Parameter alpha = new Parameter.Default(1, 10);
        Parameter beta = new Parameter.Default(1, 5);
        //			Parameter kappa = new Parameter.Default(1, 1);
        MG94CodonModel mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(mg94);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        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(simulateInPar, false);
        System.out.println(alignment.toString());
        if (calculateLikelihood) {
            // NewBeagleSequenceLikelihood nbtl = new
            // NewBeagleSequenceLikelihood(alignment, treeModel,
            // substitutionModel, (SiteModel) siteRateModel,
            // branchRateModel, null, false,
            // PartialsRescalingScheme.DEFAULT);
            ConvertAlignment convert = new ConvertAlignment(Nucleotides.INSTANCE, GeneticCode.UNIVERSAL, alignment);
            BeagleTreeLikelihood nbtl = new //
            BeagleTreeLikelihood(//
            convert, //
            treeModel, //
            substitutionModel, //
            siteRateModel, //
            branchRateModel, //
            null, //
            false, PartialsRescalingScheme.DEFAULT, true);
            System.out.println("likelihood = " + nbtl.getLogLikelihood());
        }
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Partition(dr.app.beagle.tools.Partition) BeagleTreeLikelihood(dr.evomodel.treelikelihood.BeagleTreeLikelihood) MG94CodonModel(dr.evomodel.substmodel.codon.MG94CodonModel) ArrayList(java.util.ArrayList) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) BeagleSequenceSimulator(dr.app.beagle.tools.BeagleSequenceSimulator) ImportException(dr.evolution.io.Importer.ImportException) IOException(java.io.IOException) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) TreeModel(dr.evomodel.tree.TreeModel) ConvertAlignment(dr.evolution.alignment.ConvertAlignment) Alignment(dr.evolution.alignment.Alignment) SimpleAlignment(dr.evolution.alignment.SimpleAlignment) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) NewickImporter(dr.evolution.io.NewickImporter) ConvertAlignment(dr.evolution.alignment.ConvertAlignment) Tree(dr.evolution.tree.Tree) Parameter(dr.inference.model.Parameter)

Example 9 with Partition

use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateThreePartitions.

// END: simulateTwoPartitions
static void simulateThreePartitions(int i, int N) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 3: simulateThreePartitions");
        int sequenceLength = 100000;
        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 substitution model
        Parameter kappa = new Parameter.Default(1, 10);
        HKY hky = new HKY(kappa, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(hky);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        3);
        // create partition
        Partition Partition = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        1, // to
        sequenceLength - 1, // every
        3);
        // create partition
        Partition partition3 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        2, // to
        sequenceLength - 1, // every
        3);
        partitionsList.add(partition1);
        partitionsList.add(Partition);
        partitionsList.add(partition3);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        if (i == (N - 1)) {
            System.out.println(simulator.simulate(simulateInPar, false).toString());
        } else {
            simulator.simulate(simulateInPar, false);
        }
    } 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) ArrayList(java.util.ArrayList) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) BeagleSequenceSimulator(dr.app.beagle.tools.BeagleSequenceSimulator) ImportException(dr.evolution.io.Importer.ImportException) IOException(java.io.IOException) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) TreeModel(dr.evomodel.tree.TreeModel) 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 10 with Partition

use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateTopology.

// END: annotateTree
static void simulateTopology() {
    try {
        System.out.println("Test case 1: simulateTopology");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        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);
        // set demographic function
        ExponentialGrowth exponentialGrowth = new ExponentialGrowth(Units.Type.YEARS);
        exponentialGrowth.setN0(10);
        exponentialGrowth.setGrowthRate(0.5);
        Taxa taxa = new Taxa();
        for (Taxon taxon : tree.asList()) {
            double absoluteHeight = Utils.getAbsoluteTaxonHeight(taxon, tree);
            taxon.setAttribute(Utils.ABSOLUTE_HEIGHT, absoluteHeight);
            // taxon.setAttribute("date", new Date(absoluteHeight,
            // Units.Type.YEARS, true));
            taxa.addTaxon(taxon);
        }
        // END: taxon loop
        CoalescentSimulator topologySimulator = new CoalescentSimulator();
        TreeModel treeModel = new TreeModel(topologySimulator.simulateTree(taxa, exponentialGrowth));
        System.out.println(treeModel.toString());
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create substitution model
        Parameter kappa = new Parameter.Default(1, 10);
        HKY hky = new HKY(kappa, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(hky);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        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);
        System.out.println(simulator.simulate(simulateInPar, false).toString());
    } 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) ExponentialGrowth(dr.evolution.coalescent.ExponentialGrowth) Taxon(dr.evolution.util.Taxon) ArrayList(java.util.ArrayList) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) CoalescentSimulator(dr.evolution.coalescent.CoalescentSimulator) BeagleSequenceSimulator(dr.app.beagle.tools.BeagleSequenceSimulator) ImportException(dr.evolution.io.Importer.ImportException) IOException(java.io.IOException) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) Taxa(dr.evolution.util.Taxa) TreeModel(dr.evomodel.tree.TreeModel) 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)

Aggregations

Partition (dr.app.beagle.tools.Partition)14 BeagleSequenceSimulator (dr.app.beagle.tools.BeagleSequenceSimulator)13 TreeModel (dr.evomodel.tree.TreeModel)12 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)11 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)11 FrequencyModel (dr.evomodel.substmodel.FrequencyModel)11 Tree (dr.evolution.tree.Tree)10 HomogeneousBranchModel (dr.evomodel.branchmodel.HomogeneousBranchModel)10 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)10 Parameter (dr.inference.model.Parameter)10 ArrayList (java.util.ArrayList)10 NewickImporter (dr.evolution.io.NewickImporter)9 IOException (java.io.IOException)9 ImportException (dr.evolution.io.Importer.ImportException)8 SimpleAlignment (dr.evolution.alignment.SimpleAlignment)6 HKY (dr.evomodel.substmodel.nucleotide.HKY)6 Alignment (dr.evolution.alignment.Alignment)4 SubstitutionModel (dr.evomodel.substmodel.SubstitutionModel)4 Sequence (dr.evolution.sequence.Sequence)3 BranchModel (dr.evomodel.branchmodel.BranchModel)3