Search in sources :

Example 16 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateOnePartition.

// END: simulate topology
static void simulateOnePartition() {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 2: simulateOnePartition");
        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 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
        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, false);
        // alignment.setOutputType(SimpleAlignment.OutputType.NEXUS);
        alignment.setOutputType(SimpleAlignment.OutputType.XML);
        System.out.println(alignment.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) ArrayList(java.util.ArrayList) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) Sequence(dr.evolution.sequence.Sequence) 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) SimpleAlignment(dr.evolution.alignment.SimpleAlignment) 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 17 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class IstvanOperator method doOperation.

public double doOperation() {
    Tree tree = likelihood.getTreeModel();
    alignment = likelihood.getAlignment();
    //System.out.println("Incoming alignment");
    //System.out.println(alignment);
    //System.out.println();
    SubstitutionModel substModel = likelihood.getSiteModel().getSubstitutionModel();
    // initialize the iParent and iTau arrays based on the given tree.
    initTree(tree, likelihood.getSiteModel().getMu());
    int[] treeIndex = new int[tree.getTaxonCount()];
    for (int i = 0; i < treeIndex.length; i++) {
        treeIndex[i] = tree.getTaxonIndex(alignment.getTaxonId(i));
    }
    // initialize the iAlignment array from the given alignment.
    initAlignment(alignment, treeIndex);
    // initialize the iProbs array from the substitution model -- must be called after populating tree!
    initSubstitutionModel(substModel);
    DataType dataType = substModel.getDataType();
    proposal.setGapSymbol(dataType.getGapState());
    int[][] returnedAlignment = new int[iAlignment.length][];
    //System.out.println("Initialization done, starting proposal proper...");
    double logq = proposal.propose(iAlignment, iProbs, iBaseFreqs, iParent, iTau, returnedAlignment, tuning, exponent, gapPenalty);
    //System.out.println("Proposal finished, logq=" + logq);
    //create new alignment object
    SimpleAlignment newAlignment = new SimpleAlignment();
    for (int i = 0; i < alignment.getTaxonCount(); i++) {
        StringBuffer seqBuffer = new StringBuffer();
        for (int j = 0; j < returnedAlignment[i].length; j++) {
            seqBuffer.append(dataType.getChar(returnedAlignment[treeIndex[i]][j]));
        }
        // add sequences in order of tree
        String seqString = seqBuffer.toString();
        Sequence sequence = new Sequence(alignment.getTaxon(i), seqString);
        newAlignment.addSequence(sequence);
        String oldunaligned = alignment.getUnalignedSequenceString(i);
        String unaligned = newAlignment.getUnalignedSequenceString(i);
        if (!unaligned.equals(oldunaligned)) {
            System.err.println("Sequence changed from:");
            System.err.println("old:'" + oldunaligned + "'");
            System.err.println("new:'" + unaligned + "'");
            throw new RuntimeException();
        }
    }
    //System.out.println("Outgoing alignment");
    //System.out.println(newAlignment);
    //System.out.println();
    likelihood.setAlignment(newAlignment);
    return logq;
}
Also used : SimpleAlignment(dr.evolution.alignment.SimpleAlignment) Tree(dr.evolution.tree.Tree) DataType(dr.evolution.datatype.DataType) Sequence(dr.evolution.sequence.Sequence) SubstitutionModel(dr.oldevomodel.substmodel.SubstitutionModel)

Example 18 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class AncestralStateBeagleTreeLikelihoodTest method testJointLikelihood.

public void testJointLikelihood() {
    TreeModel treeModel = new TreeModel("treeModel", tree);
    Sequence[] sequence = new Sequence[3];
    sequence[0] = new Sequence(new Taxon("0"), "A");
    sequence[1] = new Sequence(new Taxon("1"), "C");
    sequence[2] = new Sequence(new Taxon("2"), "C");
    Taxa taxa = new Taxa();
    for (Sequence s : sequence) {
        taxa.addTaxon(s.getTaxon());
    }
    SimpleAlignment alignment = new SimpleAlignment();
    for (Sequence s : sequence) {
        alignment.addSequence(s);
    }
    Parameter mu = new Parameter.Default(1, 1.0);
    Parameter kappa = new Parameter.Default(1, 1.0);
    double[] pi = { 0.25, 0.25, 0.25, 0.25 };
    Parameter freqs = new Parameter.Default(pi);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    GammaSiteRateModel siteRateModel = new GammaSiteRateModel("gammaModel", mu, null, -1, null);
    siteRateModel.setSubstitutionModel(hky);
    BranchModel branchModel = new HomogeneousBranchModel(siteRateModel.getSubstitutionModel());
    BranchRateModel branchRateModel = null;
    AncestralStateBeagleTreeLikelihood treeLikelihood = new AncestralStateBeagleTreeLikelihood(alignment, treeModel, branchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true, null, hky.getDataType(), "stateTag", // useMap = true
    true, false);
    double logLike = treeLikelihood.getLogLikelihood();
    StringBuffer buffer = new StringBuffer();
    //        Tree.Utils.newick(treeModel, treeModel.getRoot(), false, Tree.BranchLengthType.LENGTHS_AS_TIME,
    //                null, null, new NodeAttributeProvider[]{treeLikelihood}, null, null, buffer);
    TreeUtils.newick(treeModel, treeModel.getRoot(), false, TreeUtils.BranchLengthType.LENGTHS_AS_TIME, null, null, new TreeTraitProvider[] { treeLikelihood }, null, buffer);
    System.out.println(buffer);
    System.out.println("t_CA(2) = " + t(false, 2.0));
    System.out.println("t_CC(1) = " + t(true, 1.0));
    double trueValue = 0.25 * t(false, 2.0) * Math.pow(t(true, 1.0), 3.0);
    assertEquals(logLike, Math.log(trueValue), 1e-6);
}
Also used : FrequencyModel(dr.evomodel.substmodel.FrequencyModel) Taxon(dr.evolution.util.Taxon) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) GammaSiteRateModel(dr.evomodel.siteratemodel.GammaSiteRateModel) Sequence(dr.evolution.sequence.Sequence) BranchModel(dr.evomodel.branchmodel.BranchModel) HomogeneousBranchModel(dr.evomodel.branchmodel.HomogeneousBranchModel) TreeModel(dr.evomodel.tree.TreeModel) Taxa(dr.evolution.util.Taxa) SimpleAlignment(dr.evolution.alignment.SimpleAlignment) BranchRateModel(dr.evomodel.branchratemodel.BranchRateModel) HKY(dr.evomodel.substmodel.nucleotide.HKY) Parameter(dr.inference.model.Parameter) AncestralStateBeagleTreeLikelihood(dr.evomodel.treelikelihood.AncestralStateBeagleTreeLikelihood)

Example 19 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class NexusImporter method readCharactersBlock.

/**
     * Reads a 'CHARACTERS' block.
     */
private Alignment readCharactersBlock(TaxonList taxonList) throws ImportException, IOException {
    siteCount = 0;
    dataType = null;
    readDataBlockHeader("MATRIX", CHARACTERS_BLOCK);
    SimpleAlignment alignment = new SimpleAlignment();
    readSequenceData(alignment, taxonList);
    alignment.updateSiteCount();
    findEndBlock();
    return alignment;
}
Also used : SimpleAlignment(dr.evolution.alignment.SimpleAlignment)

Example 20 with SimpleAlignment

use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.

the class PhylipSequentialImporter method importAlignment.

/**
	 * importAlignment. 
	 */
public Alignment importAlignment() throws IOException, Importer.ImportException {
    SimpleAlignment alignment = null;
    try {
        int taxonCount = readInteger();
        int siteCount = readInteger();
        String firstSeq = null;
        for (int i = 0; i < taxonCount; i++) {
            StringBuffer name = new StringBuffer();
            char ch = read();
            int n = 0;
            while (!Character.isWhitespace(ch) && (maxNameLength < 1 || n < maxNameLength)) {
                name.append(ch);
                ch = read();
                n++;
            }
            StringBuffer seq = new StringBuffer(siteCount);
            readSequence(seq, dataType, "", siteCount, "-", "?", ".", firstSeq);
            if (firstSeq == null) {
                firstSeq = seq.toString();
            }
            if (alignment == null) {
                alignment = new SimpleAlignment();
            }
            alignment.addSequence(new Sequence(new Taxon(name.toString()), seq.toString()));
        }
    } catch (EOFException e) {
    }
    return alignment;
}
Also used : SimpleAlignment(dr.evolution.alignment.SimpleAlignment) Taxon(dr.evolution.util.Taxon) EOFException(java.io.EOFException) Sequence(dr.evolution.sequence.Sequence)

Aggregations

SimpleAlignment (dr.evolution.alignment.SimpleAlignment)28 Sequence (dr.evolution.sequence.Sequence)15 Taxon (dr.evolution.util.Taxon)10 ArrayList (java.util.ArrayList)9 TreeModel (dr.evomodel.tree.TreeModel)7 Parameter (dr.inference.model.Parameter)7 Tree (dr.evolution.tree.Tree)6 BranchRateModel (dr.evomodel.branchratemodel.BranchRateModel)6 GammaSiteRateModel (dr.evomodel.siteratemodel.GammaSiteRateModel)6 BeagleSequenceSimulator (dr.app.beagle.tools.BeagleSequenceSimulator)5 Partition (dr.app.beagle.tools.Partition)5 ImportException (dr.evolution.io.Importer.ImportException)5 FrequencyModel (dr.evomodel.substmodel.FrequencyModel)5 Date (dr.evolution.util.Date)4 HomogeneousBranchModel (dr.evomodel.branchmodel.HomogeneousBranchModel)4 DefaultBranchRateModel (dr.evomodel.branchratemodel.DefaultBranchRateModel)4 HKY (dr.evomodel.substmodel.nucleotide.HKY)4 IOException (java.io.IOException)4 Taxa (dr.evolution.util.Taxa)3 BranchModel (dr.evomodel.branchmodel.BranchModel)3