Examples with Sequence beast.evolution.alignment.Sequence used on opensource projects

Example 6 with Sequence

use of beast.evolution.alignment.Sequence in project beast2 by CompEvol.

the class TreeLikelihoodTest method testMarginalisationOfLikelihoodNucleotide.

@Test
public void testMarginalisationOfLikelihoodNucleotide() throws Exception {
    // test summation over all patterns adds to 1 for nucleotide data
    Sequence German_ST = new Sequence("German_ST", "           AAAAAAAAAAAAAAAA CCCCCCCCCCCCCCCC GGGGGGGGGGGGGGGG TTTTTTTTTTTTTTTT");
    Sequence Dutch_List = new Sequence("Dutch_List", "          AAAACCCCGGGGTTTT AAAACCCCGGGGTTTT AAAACCCCGGGGTTTT AAAACCCCGGGGTTTT");
    Sequence English_ST = new Sequence("English_ST", "          ACGTACGTACGTACGT ACGTACGTACGTACGT ACGTACGTACGTACGT ACGTACGTACGTACGT");
    Alignment data = new Alignment();
    data.initByName("sequence", German_ST, "sequence", Dutch_List, "sequence", English_ST, "dataType", "nucleotide");
    Tree tree = BEASTTestCase.getTree(data, "(English_ST:0.22743347188019544,(German_ST:0.10557648379843088,Dutch_List:0.10557648379843088):0.12185698808176457):0.0;");
    RealParameter frequencies = new RealParameter("0.2 0.3 0.4 0.1");
    Frequencies freqs = new Frequencies();
    freqs.initByName("frequencies", frequencies);
    GeneralSubstitutionModel gsm = new GeneralSubstitutionModel();
    gsm.initByName("rates", "1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0", "frequencies", freqs);
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("mutationRate", "1.0", "gammaCategoryCount", 1, "substModel", gsm);
    TreeLikelihood likelihood = newTreeLikelihood();
    likelihood.initByName("data", data, "tree", tree, "siteModel", siteModel);
    likelihood.initByName("useAmbiguities", false, "data", data, "tree", tree, "siteModel", siteModel);
    likelihood.calculateLogP();
    double[] logPs = likelihood.getPatternLogLikelihoods();
    double P = 0;
    for (double d : logPs) {
        P += Math.exp(d);
    }
    assertEquals(P, 1.0, BEASTTestCase.PRECISION);
}

Example 7 with Sequence

use of beast.evolution.alignment.Sequence in project beast2 by CompEvol.

the class ExchangeOperatorTest method testNarrowExchange6Taxa.

@Test
public void testNarrowExchange6Taxa() throws Exception {
    int runs = 10000;
    Randomizer.setSeed(666);
    Sequence A = new Sequence("A", "A");
    Sequence B = new Sequence("B", "A");
    Sequence C = new Sequence("C", "A");
    Sequence D = new Sequence("D", "A");
    Sequence E = new Sequence("E", "A");
    Sequence F = new Sequence("F", "A");
    Alignment data = new Alignment();
    data.initByName("sequence", A, "sequence", B, "sequence", C, "sequence", D, "sequence", E, "sequence", F, "dataType", "nucleotide");
    // String sourceTree = "((((A:2.0,B:2.0):1.0,(C:1.0,D:1.0):2.0):1.0,E:4.0):1.0,F:5.0):0.0"; // ((((A,B),(C,D)),E),F)
    // String targetTree = "((((A:2.0,(C:1.0,D:1.0):1.0):1.0,B:3.0):1.0,E:4.0):1.0,F:5.0):0.0"; // ((((A,(C,D)),B),E),F)
    String sourceTree = "(((A:5.0,B:5.0):2.0,((C:5.0,D:5.0):1.0,E:6.0):1.0):1.0,F:8.0):0.0";
    String targetTree = "(((A:5.0,B:5.0):2.0,F:7.0):1.0,((C:5.0,D:5.0):1.0,E:6.0):2.0):0.0";
    testNarrowExchange(sourceTree, targetTree, runs, data);
}

Example 8 with Sequence

use of beast.evolution.alignment.Sequence in project beast2 by CompEvol.

the class RandomTreeTest method testCoalescentTimes.

@Test
public void testCoalescentTimes() throws Exception {
    Randomizer.setSeed(53);
    int Nleaves = 10;
    int Niter = 5000;
    // (Serially sampled) coalescent time means and variances
    // estimated from 50000 trees simulated using MASTER
    double[] coalTimeMeansTruth = { 1.754662, 2.833337, 3.843532, 4.850805, 5.849542, 6.847016, 7.8482, 8.855137, 10.15442 };
    double[] coalTimeVarsTruth = { 0.2751625, 0.2727121, 0.2685172, 0.2705117, 0.2678611, 0.2671793, 0.2686952, 0.2828477, 1.076874 };
    // Assemble BEASTObjects needed by RandomTree
    StringBuilder traitSB = new StringBuilder();
    List<Sequence> seqList = new ArrayList<Sequence>();
    for (int i = 0; i < Nleaves; i++) {
        String taxonID = "t " + i;
        seqList.add(new Sequence(taxonID, "?"));
        if (i > 0)
            traitSB.append(",");
        traitSB.append(taxonID).append("=").append(i);
    }
    Alignment alignment = new Alignment(seqList, "nucleotide");
    TaxonSet taxonSet = new TaxonSet(alignment);
    TraitSet timeTrait = new TraitSet();
    timeTrait.initByName("traitname", "date-backward", "taxa", taxonSet, "value", traitSB.toString());
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    // Create RandomTree and TreeInterval instances
    RandomTree tree = new RandomTree();
    TreeIntervals intervals = new TreeIntervals();
    // Estimate coalescence time moments
    double[] coalTimeMeans = new double[Nleaves - 1];
    double[] coalTimeVars = new double[Nleaves - 1];
    double[] coalTimes = new double[Nleaves - 1];
    for (int i = 0; i < Niter; i++) {
        tree.initByName("taxa", alignment, "populationModel", popFunc, "trait", timeTrait);
        intervals.initByName("tree", tree);
        intervals.getCoalescentTimes(coalTimes);
        for (int j = 0; j < Nleaves - 1; j++) {
            coalTimeMeans[j] += coalTimes[j];
            coalTimeVars[j] += coalTimes[j] * coalTimes[j];
        }
    }
    // Normalise means and variances
    for (int j = 0; j < Nleaves - 1; j++) {
        coalTimeMeans[j] /= Niter;
        coalTimeVars[j] /= Niter;
        coalTimeVars[j] -= coalTimeMeans[j] * coalTimeMeans[j];
    }
    // Test means and variances against independently estimated values
    for (int j = 0; j < Nleaves - 1; j++) {
        assert (relError(coalTimeMeans[j], coalTimeMeansTruth[j]) < 5e-3);
        assert (relError(coalTimeVars[j], coalTimeVarsTruth[j]) < 5e-2);
    }
}

Example 9 with Sequence

use of beast.evolution.alignment.Sequence in project beast2 by CompEvol.

the class ClusterTreeTest method testUPGMA.

public void testUPGMA() throws Exception {
    Alignment alignment = BEASTTestCase.getAlignment();
    JukesCantor JC = new JukesCantor();
    JC.initAndValidate();
    SiteModel siteModel = new SiteModel();
    siteModel.initByName("substModel", JC);
    ClusterTree tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    String expectedNewick = "((((human:0.01903085702575253,(chimp:0.008560512208575313,bonobo:0.008560512208575313):0.010470344817177218):0.007962255880644985,gorilla:0.026993112906397516):0.019197419394211015,orangutan:0.04619053230060853):0.007214240662738673,siamang:0.053404772963347204):0.0";
    String actualNewick = tree.getRoot().toNewick();
    assertEquals(expectedNewick, actualNewick);
    // select 3 sequences
    List<Sequence> seqs = new ArrayList<Sequence>();
    seqs.addAll(alignment.sequenceInput.get());
    List<Sequence> newseqs = alignment.sequenceInput.get();
    newseqs.clear();
    newseqs.add(Sequence.getSequenceByTaxon("bonobo", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("chimp", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("human", seqs));
    alignment.initAndValidate();
    tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    expectedNewick = "((bonobo:0.008560512208575313,chimp:0.008560512208575313):0.010470344817177218,human:0.01903085702575253):0.0";
    System.err.println("Seqs:");
    for (Sequence s : seqs) {
        System.err.println(s.taxonInput.get());
    }
    System.err.println("Newseqs:");
    for (Sequence s : newseqs) {
        System.err.println(s.taxonInput.get());
    }
    actualNewick = tree.getRoot().toNewick();
    assertEquals(expectedNewick, actualNewick);
    // same sequences in different order
    newseqs.clear();
    newseqs.add(Sequence.getSequenceByTaxon("human", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("chimp", seqs));
    newseqs.add(Sequence.getSequenceByTaxon("bonobo", seqs));
    alignment.initAndValidate();
    tree = new ClusterTree();
    tree.initByName("clusterType", "upgma", "taxa", alignment);
    actualNewick = tree.getRoot().toNewick();
    expectedNewick = "(human:0.01903085702575253,(chimp:0.008560512208575313,bonobo:0.008560512208575313):0.010470344817177218):0.0";
    assertEquals(expectedNewick, actualNewick);
}

Example 10 with Sequence

use of beast.evolution.alignment.Sequence in project beast2 by CompEvol.

the class TaxonSetInputEditor method guessTaxonSets.

/**
 * guesses taxon sets based on pattern in regExp based on the taxa in
 * m_rawData
 */
public int guessTaxonSets(String regexp, int minSize) {
    m_taxonset.clear();
    HashMap<String, TaxonSet> map = new HashMap<>();
    Pattern m_pattern = Pattern.compile(regexp);
    Set<Taxon> taxa = new HashSet<>();
    Set<String> taxonIDs = new HashSet<>();
    for (Alignment alignment : getDoc().alignments) {
        for (String id : alignment.getTaxaNames()) {
            if (!taxonIDs.contains(id)) {
                Taxon taxon = getDoc().getTaxon(id);
                taxa.add(taxon);
                taxonIDs.add(id);
            }
        }
        for (Sequence sequence : alignment.sequenceInput.get()) {
            String id = sequence.taxonInput.get();
            if (!taxonIDs.contains(id)) {
                Taxon taxon = getDoc().getTaxon(sequence.taxonInput.get());
                // ensure sequence and taxon do not get same ID
                if (sequence.getID().equals(sequence.taxonInput.get())) {
                    sequence.setID("_" + sequence.getID());
                }
                taxa.add(taxon);
                taxonIDs.add(id);
            }
        }
    }
    List<String> unknowns = new ArrayList<>();
    for (Taxon taxon : taxa) {
        if (!(taxon instanceof TaxonSet)) {
            Matcher matcher = m_pattern.matcher(taxon.getID());
            String match;
            if (matcher.find()) {
                match = matcher.group(1);
            } else {
                match = "UNKNOWN";
                unknowns.add(taxon.getID());
            }
            try {
                if (map.containsKey(match)) {
                    TaxonSet set = map.get(match);
                    set.taxonsetInput.setValue(taxon, set);
                } else {
                    TaxonSet set = newTaxonSet(match);
                    set.taxonsetInput.setValue(taxon, set);
                    map.put(match, set);
                }
            } catch (Exception ex) {
                ex.printStackTrace();
            }
        }
    }
    if (unknowns.size() > 0) {
        showMisMatchMessage(unknowns);
    }
    // add taxon sets
    int ignored = 0;
    for (TaxonSet set : map.values()) {
        if (set.taxonsetInput.get().size() > minSize) {
            m_taxonset.add(set);
        } else {
            ignored += set.taxonsetInput.get().size();
        }
    }
    return ignored;
}