Examples with Codons dr.evolution.datatype.Codons used on opensource projects

Example 1 with Codons

use of dr.evolution.datatype.Codons in project beast-mcmc by beast-dev.

the class CodonPartitionedRobustCountingParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    AncestralStateBeagleTreeLikelihood[] partition = new AncestralStateBeagleTreeLikelihood[3];
    String[] labels = new String[] { FIRST, SECOND, THIRD };
    int patternCount = -1;
    BranchRateModel testBranchRateModel = null;
    for (int i = 0; i < 3; i++) {
        partition[i] = (AncestralStateBeagleTreeLikelihood) xo.getChild(labels[i]).getChild(AncestralStateBeagleTreeLikelihood.class);
        if (i == 0) {
            patternCount = partition[i].getPatternCount();
        } else {
            if (partition[i].getPatternCount() != patternCount) {
                throw new XMLParseException("Codon-partitioned robust counting requires all partitions to have the same length." + " Make sure that partitions include all unique sites and do not strip gaps.");
            }
        }
        // Ensure that siteRateModel has one category
        if (partition[i].getSiteRateModel().getCategoryCount() > 1) {
            throw new XMLParseException("Robust counting currently only implemented for single category models");
        }
        // Ensure that branchRateModel is the same across all partitions
        if (testBranchRateModel == null) {
            testBranchRateModel = partition[i].getBranchRateModel();
        } else if (testBranchRateModel != partition[i].getBranchRateModel()) {
            throw new XMLParseException("Robust counting currently requires the same branch rate model for all partitions");
        }
    }
    TreeModel tree = (TreeModel) xo.getChild(TreeModel.class);
    Codons codons = Codons.UNIVERSAL;
    if (xo.hasAttribute(GeneticCode.GENETIC_CODE)) {
        String codeStr = xo.getStringAttribute(GeneticCode.GENETIC_CODE);
        codons = Codons.findByName(codeStr);
    }
    String labelingString = (String) xo.getAttribute(LABELING);
    CodonLabeling codonLabeling = CodonLabeling.parseFromString(labelingString);
    if (codonLabeling == null) {
        throw new XMLParseException("Unrecognized codon labeling '" + labelingString + "'");
    }
    String branchFormatString = xo.getAttribute(BRANCH_FORMAT, StratifiedTraitOutputFormat.SUM_OVER_SITES.getText());
    StratifiedTraitOutputFormat branchFormat = StratifiedTraitOutputFormat.parseFromString(branchFormatString);
    if (branchFormat == null) {
        throw new XMLParseException("Unrecognized branch output format '" + branchFormat + "'");
    }
    String logFormatString = xo.getAttribute(LOG_FORMAT, StratifiedTraitOutputFormat.SUM_OVER_SITES.getText());
    StratifiedTraitOutputFormat logFormat = StratifiedTraitOutputFormat.parseFromString(logFormatString);
    if (logFormat == null) {
        throw new XMLParseException("Unrecognized log output format '" + branchFormat + "'");
    }
    boolean useUniformization = xo.getAttribute(USE_UNIFORMIZATION, false);
    boolean includeExternalBranches = xo.getAttribute(INCLUDE_EXTERNAL, true);
    boolean includeInternalBranches = xo.getAttribute(INCLUDE_INTERNAL, true);
    boolean doUnconditionedPerBranch = xo.getAttribute(DO_UNCONDITIONED_PER_BRANCH, false);
    boolean averageRates = xo.getAttribute(AVERAGE_RATES, true);
    boolean saveCompleteHistory = xo.getAttribute(SAVE_HISTORY, false);
    boolean useNewNeutralModel = xo.getAttribute(USE_NEW_NEUTRAL_MODEL, false);
    String prefix = xo.hasAttribute(PREFIX) ? xo.getStringAttribute(PREFIX) : null;
    return new CodonPartitionedRobustCounting(xo.getId(), tree, partition, codons, codonLabeling, useUniformization, includeExternalBranches, includeInternalBranches, doUnconditionedPerBranch, saveCompleteHistory, averageRates, useNewNeutralModel, branchFormat, logFormat, prefix);
}

Example 2 with Codons

use of dr.evolution.datatype.Codons 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;
}

Example 3 with Codons

use of dr.evolution.datatype.Codons 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 4 with Codons

use of dr.evolution.datatype.Codons in project beast-mcmc by beast-dev.

the class FrequencyModelParser method getEmpirical3x4Freqs.

// END: parseXMLObject
private double[] getEmpirical3x4Freqs(PatternList patternList) {
    DataType nucleotideDataType = Nucleotides.INSTANCE;
    Codons codonDataType = Codons.UNIVERSAL;
    List<String> stopCodonsList = Arrays.asList(STOP_CODONS);
    int cStateCount = codonDataType.getStateCount();
    int nStateCount = nucleotideDataType.getStateCount();
    int nPosition = 3;
    double[] stopCodonFreqs = new double[STOP_CODONS.length];
    double[][] counts = new double[nStateCount][nPosition];
    int[] countsPos = new int[nPosition];
    int patternCount = patternList.getPatternCount();
    for (int i = 0; i < patternCount; i++) {
        int[] sitePatterns = patternList.getPattern(i);
        for (int codonState : sitePatterns) {
            int[] nucleotideStates = codonDataType.getTripletStates(codonState);
            String triplet = codonDataType.getTriplet(codonState);
            if (stopCodonsList.contains(triplet)) {
                int stopCodonIndex = stopCodonsList.indexOf(triplet);
                stopCodonFreqs[stopCodonIndex]++;
            }
            for (int pos = 0; pos < nPosition; pos++) {
                int nucleotideState = nucleotideStates[pos];
                counts[nucleotideState][pos]++;
                countsPos[pos]++;
            }
        // END: nucleotide positions loop
        }
    // END: sitePatterns loop
    }
    // sites loop
    int total = 0;
    for (int pos = 0; pos < nPosition; pos++) {
        int totalPos = countsPos[pos];
        for (int s = 0; s < nStateCount; s++) {
            counts[s][pos] = counts[s][pos] / totalPos;
        }
        // END: nucleotide states loop
        total += totalPos;
    }
    // END: nucleotide positions loop
    // Utils.printArray(stopCodonFreqs);
    // System.out.println(stopCodonFreqs.length);
    // add stop codon frequencies
    double pi = 0.0;
    for (double stopCodonFreq : stopCodonFreqs) {
        double freq = stopCodonFreq / total;
        pi += freq;
    }
    // System.out.println(pi);
    double[] freqs = new double[cStateCount];
    Arrays.fill(freqs, 1.0);
    for (int codonState = 0; codonState < cStateCount; codonState++) {
        int[] nucleotideStates = codonDataType.getTripletStates(codonState);
        for (int pos = 0; pos < nPosition; pos++) {
            int nucleotide = nucleotideStates[pos];
            freqs[codonState] *= counts[nucleotide][pos];
        }
        // END: nucleotide positions loop
        // TODO: stop codons freqs
        freqs[codonState] = freqs[codonState] / (1 - pi);
    }
    return freqs;
}

Example 5 with Codons

use of dr.evolution.datatype.Codons in project beast-mcmc by beast-dev.

the class Partition method simulatePartition.

// END: loadBeagleInstance
public void simulatePartition() {
    try {
        NodeRef root = treeModel.getRoot();
        // gamma category rates
        double[] categoryRates = siteRateModel.getCategoryRates();
        beagle.setCategoryRates(categoryRates);
        // probabilities for gamma category rates
        double[] categoryProbs = siteRateModel.getCategoryProportions();
        // beagle.setCategoryWeights(0, categoryProbs);
        // Utils.printArray(categoryRates);
        // Utils.printArray(categoryProbs);
        int[] category = new int[partitionSiteCount];
        for (int i = 0; i < partitionSiteCount; i++) {
            category[i] = randomChoicePDF(categoryProbs, partitionNumber, "categories");
        }
        if (DEBUG) {
            System.out.println("category for each site:");
            Utils.printArray(category);
        }
        // END: DEBUG
        int[] parentSequence = new int[partitionSiteCount];
        // set ancestral sequence for partition if it exists
        if (hasRootSequence) {
            if (rootSequence.getLength() == partitionSiteCount) {
                parentSequence = sequence2intArray(rootSequence);
            } else if (dataType instanceof Codons && rootSequence.getLength() == 3 * partitionSiteCount) {
                parentSequence = sequence2intArray(rootSequence);
            } else {
                throw new RuntimeException("Ancestral sequence length of " + rootSequence.getLength() + " does not match partition site count of " + partitionSiteCount + ".");
            }
        } else {
            double[] frequencies = freqModel.getFrequencies();
            for (int i = 0; i < partitionSiteCount; i++) {
                parentSequence[i] = randomChoicePDF(frequencies, partitionNumber, "root");
            }
        }
        if (DEBUG) {
            synchronized (this) {
                System.out.println();
                System.out.println("root Sequence:");
                Utils.printArray(parentSequence);
            }
        }
        // END: DEBUG
        substitutionModelDelegate.updateSubstitutionModels(beagle);
        traverse(root, parentSequence, category);
        if (DEBUG) {
            synchronized (this) {
                System.out.println("Simulated alignment:");
                printSequences();
            }
        }
        // END: DEBUG
        beagle.finalize();
    } catch (Exception e) {
        e.printStackTrace();
    } catch (Throwable e) {
        System.err.println("BeagleException: " + e.getMessage());
        System.exit(-1);
    }
}