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

Example 71 with Alignment

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

the class FilteredAlignmentTest method testWeightedSitesReordered.

public void testWeightedSitesReordered() throws Exception {
    // reorder taxa
    Alignment data = getAlignmentNoTInHuman();
    data.setID("data");
    List<Taxon> taxa = new ArrayList<>();
    taxa.add(new Taxon("1chimp"));
    taxa.add(new Taxon("0human"));
    TaxonSet set = new TaxonSet(taxa);
    data.taxonSetInput.setValue(set, data);
    data.siteWeightsInput.setValue("11232, 2, 3, 4 ,1123,2,3,4,112,2,3,4,11,2,3,	4 ", data);
    data.initAndValidate();
    String weights = Arrays.toString(data.getWeights());
    System.out.println(weights + "\n0human\t" + alignmentToString(data, data.getTaxonIndex("0human")) + "\n1chimp\t" + alignmentToString(data, data.getTaxonIndex("1chimp")));
    assertEquals("[11243, 1123, 112, 4, 2, 2, 6, 3, 3, 8, 4, 4]", weights);
}

Example 72 with Alignment

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

the class FilteredAlignmentTest method getAlignment3.

public static Alignment getAlignment3() throws Exception {
    // reordered from getAlignment() & with duplicates
    Sequence human = new Sequence("0human", "GGGAAA");
    Sequence chimp = new Sequence("1chimp", "AGGACA");
    Alignment data = new Alignment();
    data.initByName("sequence", human, "sequence", chimp, "dataType", "nucleotide");
    return data;
}

Example 73 with Alignment

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

the class UnorderedAlignmentsTest method testUnorderedAlignment.

@Test
public void testUnorderedAlignment() throws Exception {
    TaxonSet taxa = getTaxa();
    Tree tree = getTree(taxa);
    SiteModel siteModel = getSiteModel();
    double logP = 0.0;
    double shuffledLogP = 0.0;
    for (int i = 0; i < 3; ++i) {
        Alignment data = getAlignment(taxa, tree, siteModel);
        // First calculate in order
        TreeLikelihood likelihood = new TreeLikelihood();
        likelihood.initByName("data", data, "tree", tree, "siteModel", siteModel);
        logP += likelihood.calculateLogP();
        // Now calculate again, with shuffled taxon order
        Collections.shuffle(data.sequenceInput.get());
        likelihood = new TreeLikelihood();
        likelihood.initByName("data", data, "tree", tree, "siteModel", siteModel);
        shuffledLogP += likelihood.calculateLogP();
    }
    assertEquals(logP, shuffledLogP, 1E-9);
}

Example 74 with Alignment

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

the class UnorderedAlignmentsTest method getAlignment.

public static Alignment getAlignment(TaxonSet taxa, Tree tree, SiteModel siteModel) throws Exception {
    Alignment dummy = new Alignment();
    Object[] args = new String[2 * taxa.getTaxonCount() + 2];
    args[args.length - 2] = "dataType";
    args[args.length - 1] = "nucleotide";
    for (int i = 0; i < taxa.getTaxonCount(); ++i) {
        args[2 * i] = taxa.getTaxonId(i);
        args[2 * i + 1] = "N";
    }
    dummy.initByName(args);
    SimulatedAlignment data = new SimulatedAlignment();
    data.initByName("data", dummy, "tree", tree, "siteModel", siteModel);
    return data;
}

Example 75 with Alignment

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

the class StarBeastStartState method fullInit.

private void fullInit() {
    // Build gene trees from  alignments
    final Function muInput = this.muInput.get();
    final double mu = (muInput != null) ? muInput.getArrayValue() : 1;
    final Tree stree = speciesTreeInput.get();
    final TaxonSet species = stree.m_taxonset.get();
    final List<String> speciesNames = species.asStringList();
    final int speciesCount = speciesNames.size();
    final List<Tree> geneTrees = genes.get();
    // final List<Alignment> alignments = genes.get();
    // final List<Tree> geneTrees = new ArrayList<>(alignments.size());
    double maxNsites = 0;
    // for( final Alignment alignment : alignments)  {
    for (final Tree gtree : geneTrees) {
        // final Tree gtree = new Tree();
        final Alignment alignment = gtree.m_taxonset.get().alignmentInput.get();
        final ClusterTree ctree = new ClusterTree();
        ctree.initByName("initial", gtree, "clusterType", "upgma", "taxa", alignment);
        gtree.scale(1 / mu);
        maxNsites = max(maxNsites, alignment.getSiteCount());
    }
    final Map<String, Integer> geneTips2Species = new HashMap<>();
    final List<Taxon> taxonSets = species.taxonsetInput.get();
    for (int k = 0; k < speciesNames.size(); ++k) {
        final Taxon nx = taxonSets.get(k);
        final List<Taxon> taxa = ((TaxonSet) nx).taxonsetInput.get();
        for (final Taxon n : taxa) {
            geneTips2Species.put(n.getID(), k);
        }
    }
    final double[] dg = new double[(speciesCount * (speciesCount - 1)) / 2];
    final double[][] genesDmins = new double[geneTrees.size()][];
    for (int ng = 0; ng < geneTrees.size(); ++ng) {
        final Tree g = geneTrees.get(ng);
        final double[] dmin = firstMeetings(g, geneTips2Species, speciesCount);
        genesDmins[ng] = dmin;
        for (int i = 0; i < dmin.length; ++i) {
            dg[i] += dmin[i];
            if (dmin[i] == Double.MAX_VALUE) {
                // this happens when a gene tree has no taxa for some species-tree taxon.
                // TODO: ensure that if this happens, there will always be an "infinite"
                // distance between species-taxon 0 and the species-taxon with missing lineages,
                // so i < speciesCount - 1.
                // What if lineages for species-taxon 0 are missing? Then all entries will be 'infinite'.
                String id = (i < speciesCount - 1 ? stree.getExternalNodes().get(i + 1).getID() : "unknown taxon");
                if (i == 0) {
                    // test that all entries are 'infinite', which implies taxon 0 has lineages missing
                    boolean b = true;
                    for (int k = 1; b && k < speciesCount - 1; k++) {
                        b = (dmin[k] == Double.MAX_VALUE);
                    }
                    if (b) {
                        // if all entries have 'infinite' distances, it is probably the first taxon that is at fault
                        id = stree.getExternalNodes().get(0).getID();
                    }
                }
                throw new RuntimeException("Gene tree " + g.getID() + " has no lineages for species taxon " + id + " ");
            }
        }
    }
    for (int i = 0; i < dg.length; ++i) {
        double d = dg[i] / geneTrees.size();
        if (d == 0) {
            d = (0.5 / maxNsites) * (1 / mu);
        } else {
            // heights to distances
            d *= 2;
        }
        dg[i] = d;
    }
    final ClusterTree ctree = new ClusterTree();
    final Distance distance = new Distance() {

        @Override
        public double pairwiseDistance(final int s1, final int s2) {
            final int i = getDMindex(speciesCount, s1, s2);
            return dg[i];
        }
    };
    ctree.initByName("initial", stree, "taxonset", species, "clusterType", "upgma", "distance", distance);
    final Map<String, Integer> sptips2SpeciesIndex = new HashMap<>();
    for (int i = 0; i < speciesNames.size(); ++i) {
        sptips2SpeciesIndex.put(speciesNames.get(i), i);
    }
    final double[] spmin = firstMeetings(stree, sptips2SpeciesIndex, speciesCount);
    for (int ng = 0; ng < geneTrees.size(); ++ng) {
        final double[] dmin = genesDmins[ng];
        boolean compatible = true;
        for (int i = 0; i < spmin.length; ++i) {
            if (dmin[i] <= spmin[i]) {
                compatible = false;
                break;
            }
        }
        if (!compatible) {
            final Tree gtree = geneTrees.get(ng);
            final TaxonSet gtreeTaxa = gtree.m_taxonset.get();
            final Alignment alignment = gtreeTaxa.alignmentInput.get();
            final List<String> taxaNames = alignment.getTaxaNames();
            final int taxonCount = taxaNames.size();
            // speedup
            final Map<Integer, Integer> g2s = new HashMap<>();
            for (int i = 0; i < taxonCount; ++i) {
                g2s.put(i, geneTips2Species.get(taxaNames.get(i)));
            }
            final JukesCantorDistance jc = new JukesCantorDistance();
            jc.setPatterns(alignment);
            final Distance gdistance = new Distance() {

                @Override
                public double pairwiseDistance(final int t1, final int t2) {
                    final int s1 = g2s.get(t1);
                    final int s2 = g2s.get(t2);
                    double d = jc.pairwiseDistance(t1, t2) / mu;
                    if (s1 != s2) {
                        final int i = getDMindex(speciesCount, s1, s2);
                        final double minDist = 2 * spmin[i];
                        if (d <= minDist) {
                            d = minDist * 1.001;
                        }
                    }
                    return d;
                }
            };
            final ClusterTree gtreec = new ClusterTree();
            gtreec.initByName("initial", gtree, "taxonset", gtreeTaxa, "clusterType", "upgma", "distance", gdistance);
        }
    }
    {
        final RealParameter lambda = birthRate.get();
        if (lambda != null) {
            final double rh = stree.getRoot().getHeight();
            double l = 0;
            for (int i = 2; i < speciesCount + 1; ++i) {
                l += 1. / i;
            }
            setParameterValue(lambda, (1 / rh) * l);
        }
        double totBranches = 0;
        final Node[] streeNodeas = stree.getNodesAsArray();
        for (final Node n : streeNodeas) {
            if (!n.isRoot()) {
                totBranches += n.getLength();
            }
        }
        totBranches /= 2 * (streeNodeas.length - 1);
        final RealParameter popm = popMean.get();
        if (popm != null) {
            setParameterValue(popm, totBranches);
        }
        final SpeciesTreePrior speciesTreePrior = speciesTreePriorInput.get();
        if (speciesTreePrior != null) {
            final RealParameter popb = speciesTreePrior.popSizesBottomInput.get();
            if (popb != null) {
                for (int i = 0; i < popb.getDimension(); ++i) {
                    setParameterValue(popb, i, 2 * totBranches);
                }
            }
            final RealParameter popt = speciesTreePrior.popSizesTopInput.get();
            if (popt != null) {
                for (int i = 0; i < popt.getDimension(); ++i) {
                    setParameterValue(popt, i, totBranches);
                }
            }
        }
    }
}