Examples with MG94K80CodonModel dr.evomodel.substmodel.codon.MG94K80CodonModel used on opensource projects

Example 1 with MG94K80CodonModel

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

the class MG94CodonModelParser 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 alphaParam = (Parameter) xo.getElementFirstChild(ALPHA);
    Parameter betaParam = (Parameter) xo.getElementFirstChild(BETA);
    FrequencyModel freqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    boolean isParameterTotalRate = xo.getAttribute(TOTAL_RATES, true);
    MG94HKYCodonModel codonModel;
    if (xo.hasChildNamed(GTR_MODEL)) {
        // TODO: change this into constructing a MG94HKYCodonModel (needs to be written), which is started underneath
        codonModel = new MG94K80CodonModel(codons, alphaParam, betaParam, freqModel, new CodonOptions(isParameterTotalRate));
        Parameter rateACValue = null;
        if (xo.hasChildNamed(A_TO_C)) {
            rateACValue = (Parameter) xo.getElementFirstChild(A_TO_C);
        }
        Parameter rateAGValue = null;
        if (xo.hasChildNamed(A_TO_G)) {
            rateAGValue = (Parameter) xo.getElementFirstChild(A_TO_G);
        }
        Parameter rateATValue = null;
        if (xo.hasChildNamed(A_TO_T)) {
            rateATValue = (Parameter) xo.getElementFirstChild(A_TO_T);
        }
        Parameter rateCGValue = null;
        if (xo.hasChildNamed(C_TO_G)) {
            rateCGValue = (Parameter) xo.getElementFirstChild(C_TO_G);
        }
        Parameter rateCTValue = null;
        if (xo.hasChildNamed(C_TO_T)) {
            rateCTValue = (Parameter) xo.getElementFirstChild(C_TO_T);
        }
        Parameter rateGTValue = null;
        if (xo.hasChildNamed(G_TO_T)) {
            rateGTValue = (Parameter) xo.getElementFirstChild(G_TO_T);
        }
        int countNull = 0;
        if (rateACValue == null)
            countNull++;
        if (rateAGValue == null)
            countNull++;
        if (rateATValue == null)
            countNull++;
        if (rateCGValue == null)
            countNull++;
        if (rateCTValue == null)
            countNull++;
        if (rateGTValue == null)
            countNull++;
        if (countNull != 1)
            throw new XMLParseException("Only five parameters may be specified in GTR, leave exactly one out, the others will be specifed relative to the one left out.");
        if (xo.hasAttribute(KAPPA)) {
            System.err.print("using GTR rates -- overrides KAPPA");
        }
    } else if (xo.hasChildNamed(KAPPA)) {
        Parameter kappaParam = (Parameter) xo.getElementFirstChild(KAPPA);
        codonModel = new MG94HKYCodonModel(codons, alphaParam, betaParam, kappaParam, freqModel, new CodonOptions(isParameterTotalRate));
    // System.err.println("setting up MG94K80CodonModel");
    } else {
        // resort to standard MG94 without nucleotide rate bias
        codonModel = new MG94K80CodonModel(codons, alphaParam, betaParam, freqModel, new CodonOptions(isParameterTotalRate));
    }
    if (!xo.getAttribute(NORMALIZED, true)) {
    // codonModel.setNormalization(false);
    // Logger.getLogger("dr.app.beagle.evomodel").info("MG94HKYCodonModel normalization: false");
    }
    return codonModel;
}

Example 2 with MG94K80CodonModel

use of dr.evomodel.substmodel.codon.MG94K80CodonModel 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);
        DefaultTreeModel treeModel = new DefaultTreeModel(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);
        MG94HKYCodonModel mg94 = new MG94K80CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel, new CodonOptions());
        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
}

Example 3 with MG94K80CodonModel

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

the class LineageSpecificBranchModel method main.

// END: acceptState
public static void main(String[] args) {
    try {
        // the seed of the BEAST
        MathUtils.setSeed(666);
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        TreeModel tree = new DefaultTreeModel(importer.importTree(null));
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { // 
        0.0163936, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344, // 
        0.01639344 });
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create substitution model
        Parameter alpha = new Parameter.Default(1, 10);
        Parameter beta = new Parameter.Default(1, 5);
        MG94HKYCodonModel mg94 = new MG94K80CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel, new CodonOptions());
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(mg94);
        // create partition
        Partition partition1 = new // 
        Partition(// 
        tree, // 
        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(false, false);
        ConvertAlignment convert = new ConvertAlignment(Nucleotides.INSTANCE, GeneticCode.UNIVERSAL, alignment);
        List<SubstitutionModel> substModels = new ArrayList<SubstitutionModel>();
        for (int i = 0; i < 2; i++) {
            // alpha = new Parameter.Default(1, 10 );
            // beta = new Parameter.Default(1, 5 );
            // mg94 = new MG94HKYCodonModel(Codons.UNIVERSAL, alpha, beta,
            // freqModel);
            substModels.add(mg94);
        }
        Parameter uCategories = new Parameter.Default(2, 0);
        // CountableBranchCategoryProvider provider = new CountableBranchCategoryProvider.IndependentBranchCategoryModel(tree, uCategories);
        LineageSpecificBranchModel branchSpecific = new // provider,
        LineageSpecificBranchModel(// provider,
        tree, // provider,
        freqModel, // provider,
        substModels, uCategories);
        BeagleTreeLikelihood like = new // 
        BeagleTreeLikelihood(// 
        convert, // 
        tree, // 
        branchSpecific, // 
        siteRateModel, // 
        branchRateModel, // 
        null, // 
        false, PartialsRescalingScheme.DEFAULT, true);
        BeagleTreeLikelihood gold = new // 
        BeagleTreeLikelihood(// 
        convert, // 
        tree, // 
        substitutionModel, // 
        siteRateModel, // 
        branchRateModel, // 
        null, // 
        false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("likelihood (gold) = " + gold.getLogLikelihood());
        System.out.println("likelihood = " + like.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
    }
}