Search in sources :

Example 1 with GY94CodonModel

use of dr.evomodel.substmodel.codon.GY94CodonModel in project beast-mcmc by beast-dev.

the class RandomBranchModelParser method parseXMLObject.

@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Logger.getLogger("dr.evomodel").info("\nUsing random assignment branch model.");
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    GY94CodonModel baseSubstitutionModel = (GY94CodonModel) xo.getElementFirstChild(BASE_MODEL);
    long seed = -1;
    boolean hasSeed = false;
    if (xo.hasAttribute(SEED)) {
        seed = xo.getLongIntegerAttribute(SEED);
        hasSeed = true;
    }
    double rate = 1;
    if (xo.hasAttribute(RATE)) {
        rate = xo.getDoubleAttribute(RATE);
    }
    return new RandomBranchModel(treeModel, baseSubstitutionModel, rate, hasSeed, seed);
}
Also used : TreeModel(dr.evomodel.tree.TreeModel) RandomBranchModel(dr.evomodel.branchmodel.RandomBranchModel) GY94CodonModel(dr.evomodel.substmodel.codon.GY94CodonModel)

Example 2 with GY94CodonModel

use of dr.evomodel.substmodel.codon.GY94CodonModel in project beast-mcmc by beast-dev.

the class CompleteHistorySimulatorTest method testCodonSimulation.

public void testCodonSimulation() {
    Parameter kappa = new Parameter.Default(1, 2.0);
    // Expect many more non-syn changes
    Parameter omega = new Parameter.Default(1, 5.0);
    Codons codons = Codons.UNIVERSAL;
    int stateCount = codons.getStateCount();
    double[] p = new double[stateCount];
    for (int i = 0; i < stateCount; i++) {
        p[i] = 1.0 / (double) stateCount;
    }
    Parameter pi = new Parameter.Default(p);
    FrequencyModel f = new FrequencyModel(codons, pi);
    GY94CodonModel codonModel = new GY94CodonModel(codons, omega, kappa, f);
    Parameter mu = new Parameter.Default(1, 0.5);
    Parameter alpha = new Parameter.Default(1, 0.5);
    GammaSiteRateModel siteModel = new GammaSiteRateModel("gammaModel", mu, alpha, 4, null);
    siteModel.setSubstitutionModel(codonModel);
    BranchRateModel branchRateModel = new DefaultBranchRateModel();
    double analyticResult = TreeUtils.getTreeLength(tree, tree.getRoot()) * mu.getParameterValue(0);
    int nSites = 100;
    // use base 61
    double[] synRegMatrix = CodonLabeling.getRegisterMatrix(CodonLabeling.SYN, codons, false);
    // use base 61
    double[] nonSynRegMatrix = CodonLabeling.getRegisterMatrix(CodonLabeling.NON_SYN, codons, false);
    runSimulation(tree, siteModel, branchRateModel, nSites, new double[][] { synRegMatrix, nonSynRegMatrix }, analyticResult);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) Parameter(dr.inference.model.Parameter) GY94CodonModel(dr.evomodel.substmodel.codon.GY94CodonModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) Codons(dr.evolution.datatype.Codons) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel)

Example 3 with GY94CodonModel

use of dr.evomodel.substmodel.codon.GY94CodonModel in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateRandomBranchAssignment.

// END: main
static void simulateRandomBranchAssignment() {
    try {
        System.out.println("Test case I dunno which: simulate random branch assignments");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        File treeFile = new File("/home/filip/Dropbox/BeagleSequenceSimulator/SimTree/SimTree.figtree");
        Tree tree = Utils.importTreeFromFile(treeFile);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(Utils.UNIFORM_CODON_FREQUENCIES);
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create base subst model
        Parameter omegaParameter = new Parameter.Default("omega", 1, 1.0);
        Parameter kappaParameter = new Parameter.Default("kappa", 1, 1.0);
        GY94CodonModel baseSubModel = new GY94CodonModel(Codons.UNIVERSAL, omegaParameter, kappaParameter, freqModel);
        RandomBranchModel substitutionModel = new RandomBranchModel(treeModel, baseSubModel, 0.25, false, -1);
        // 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);
        //			Sequence ancestralSequence = new Sequence();
        //			ancestralSequence.appendSequenceString("TCAAGTGAGG");
        //			partition1.setRootSequence(ancestralSequence);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        SimpleAlignment alignment = simulator.simulate(simulateInPar, true);
        // alignment.setOutputType(SimpleAlignment.OutputType.NEXUS);
        alignment.setOutputType(SimpleAlignment.OutputType.XML);
        System.out.println(alignment.toString());
    } catch (Exception e) {
        e.printStackTrace();
    }
// END: try-catch
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Partition(dr.app.beagle.tools.Partition) ArrayList(java.util.ArrayList) 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) RandomBranchModel(dr.evomodel.branchmodel.RandomBranchModel) SimpleAlignment(dr.evolution.alignment.SimpleAlignment) DefaultBranchRateModel(dr.evomodel.branchratemodel.DefaultBranchRateModel) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) Tree(dr.evolution.tree.Tree) Parameter(dr.inference.model.Parameter) GY94CodonModel(dr.evomodel.substmodel.codon.GY94CodonModel) File(java.io.File)

Example 4 with GY94CodonModel

use of dr.evomodel.substmodel.codon.GY94CodonModel in project beast-mcmc by beast-dev.

the class PartitionData method createBranchModel.

public BranchModel createBranchModel() {
    BranchModel branchModel = null;
    if (this.substitutionModelIndex == 0) {
        // HKY
        Parameter kappa = new Parameter.Default(1, substitutionParameterValues[0]);
        FrequencyModel frequencyModel = this.createFrequencyModel();
        HKY hky = new HKY(kappa, frequencyModel);
        branchModel = new HomogeneousBranchModel(hky);
    } else if (this.substitutionModelIndex == 1) {
        // GTR
        Parameter ac = new Parameter.Default(1, substitutionParameterValues[1]);
        Parameter ag = new Parameter.Default(1, substitutionParameterValues[2]);
        Parameter at = new Parameter.Default(1, substitutionParameterValues[3]);
        Parameter cg = new Parameter.Default(1, substitutionParameterValues[4]);
        Parameter ct = new Parameter.Default(1, substitutionParameterValues[5]);
        Parameter gt = new Parameter.Default(1, substitutionParameterValues[6]);
        FrequencyModel frequencyModel = this.createFrequencyModel();
        GTR gtr = new GTR(ac, ag, at, cg, ct, gt, frequencyModel);
        branchModel = new HomogeneousBranchModel(gtr);
    } else if (this.substitutionModelIndex == 2) {
        // TN93
        Parameter kappa1 = new Parameter.Default(1, substitutionParameterValues[7]);
        Parameter kappa2 = new Parameter.Default(1, substitutionParameterValues[8]);
        FrequencyModel frequencyModel = this.createFrequencyModel();
        TN93 tn93 = new TN93(kappa1, kappa2, frequencyModel);
        branchModel = new HomogeneousBranchModel(tn93);
    } else if (this.substitutionModelIndex == 3) {
        // Yang Codon Model
        FrequencyModel frequencyModel = this.createFrequencyModel();
        Parameter kappa = new Parameter.Default(1, substitutionParameterValues[9]);
        Parameter omega = new Parameter.Default(1, substitutionParameterValues[10]);
        GY94CodonModel yangCodonModel = new GY94CodonModel(Codons.UNIVERSAL, omega, kappa, frequencyModel);
        branchModel = new HomogeneousBranchModel(yangCodonModel);
    } else if (this.substitutionModelIndex == 4) {
        // MG94CodonModel
        FrequencyModel frequencyModel = this.createFrequencyModel();
        Parameter alpha = new Parameter.Default(1, substitutionParameterValues[11]);
        Parameter beta = new Parameter.Default(1, substitutionParameterValues[12]);
        Parameter kappa = new Parameter.Default(1, substitutionParameterValues[13]);
        MG94HKYCodonModel mg94 = new MG94HKYCodonModel(Codons.UNIVERSAL, alpha, beta, kappa, frequencyModel);
        branchModel = new HomogeneousBranchModel(mg94);
    } else if (this.substitutionModelIndex == 5) {
        // Blosum62
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = Blosum62.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 6) {
        // CPREV
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = CPREV.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 7) {
        // Dayhoff
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = Dayhoff.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 8) {
        // JTT
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = JTT.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 9) {
        // LG
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = LG.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 10) {
        // MTREV
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = MTREV.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 11) {
        // WAG
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = WAG.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else {
        System.out.println("Not yet implemented");
    }
    return branchModel;
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) EmpiricalRateMatrix(dr.evomodel.substmodel.EmpiricalRateMatrix) GTR(dr.evomodel.substmodel.nucleotide.GTR) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) BranchModel(dr.evomodel.branchmodel.BranchModel) TN93(dr.evomodel.substmodel.nucleotide.TN93) HKY(dr.evomodel.substmodel.nucleotide.HKY) Parameter(dr.inference.model.Parameter) GY94CodonModel(dr.evomodel.substmodel.codon.GY94CodonModel) MG94HKYCodonModel(dr.evomodel.substmodel.codon.MG94HKYCodonModel)

Example 5 with GY94CodonModel

use of dr.evomodel.substmodel.codon.GY94CodonModel in project beast-mcmc by beast-dev.

the class GY94CodonModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Codons codons = Codons.UNIVERSAL;
    if (xo.hasAttribute(GeneticCode.GENETIC_CODE)) {
        String codeStr = xo.getStringAttribute(GeneticCode.GENETIC_CODE);
        codons = Codons.findByName(codeStr);
    }
    Parameter omegaParameter = (Parameter) xo.getElementFirstChild(OMEGA);
    int dim = omegaParameter.getDimension();
    double value = omegaParameter.getParameterValue(dim - 1);
    if (value < 0) {
        throw new RuntimeException("Negative Omega parameter value " + value);
    }
    //END: negative check
    Parameter kappaParameter = (Parameter) xo.getElementFirstChild(KAPPA);
    dim = kappaParameter.getDimension();
    value = kappaParameter.getParameterValue(dim - 1);
    if (value < 0) {
        throw new RuntimeException("Negative kappa parameter value value " + value);
    }
    //END: negative check
    FrequencyModel freqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    GY94CodonModel codonModel = new GY94CodonModel(codons, omegaParameter, kappaParameter, freqModel);
    return codonModel;
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Parameter(dr.inference.model.Parameter) GY94CodonModel(dr.evomodel.substmodel.codon.GY94CodonModel) Codons(dr.evolution.datatype.Codons)

Aggregations

GY94CodonModel (dr.evomodel.substmodel.codon.GY94CodonModel)6 FrequencyModel (dr.evomodel.substmodel.FrequencyModel)5 Parameter (dr.inference.model.Parameter)5 Codons (dr.evolution.datatype.Codons)2 RandomBranchModel (dr.evomodel.branchmodel.RandomBranchModel)2 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)2 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)2 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)2 TreeModel (dr.evomodel.tree.TreeModel)2 BeagleSequenceSimulator (dr.app.beagle.tools.BeagleSequenceSimulator)1 Partition (dr.app.beagle.tools.Partition)1 SimpleAlignment (dr.evolution.alignment.SimpleAlignment)1 DataType (dr.evolution.datatype.DataType)1 ImportException (dr.evolution.io.Importer.ImportException)1 NodeRef (dr.evolution.tree.NodeRef)1 Tree (dr.evolution.tree.Tree)1 BranchModel (dr.evomodel.branchmodel.BranchModel)1 HomogeneousBranchModel (dr.evomodel.branchmodel.HomogeneousBranchModel)1 EmpiricalRateMatrix (dr.evomodel.substmodel.EmpiricalRateMatrix)1 SubstitutionModel (dr.evomodel.substmodel.SubstitutionModel)1