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

Example 26 with Alignment

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

the class HiddenLinkageModel method computeTipPartials.

private void computeTipPartials(int nodeIndex) {
    if (!dirtyTipPartials[nodeIndex]) {
        swapTipPartials(nodeIndex);
        dirtyTipPartials[nodeIndex] = true;
    }
    double[] tipPartials = this.tipPartials[nodeIndex];
    // if this is one of the reference organisms, then return the resolved partials
    // if this is a linkage group, return partials that correspond to probabilities of each nucleotide.
    Alignment aln = data.getAlignment();
    int sc = aln.getStateCount();
    for (int i = 0; i < tipPartials.length; i++) {
        tipPartials[i] = 0.0;
    }
    if (nodeIndex < data.getReferenceTaxa().getTaxonCount()) {
        int j = 0;
        for (int i = 0; i < aln.getSiteCount(); i++) {
            int s = aln.getState(nodeIndex, i);
            if (s >= sc) {
                for (int k = 0; k < sc; k++) tipPartials[j + k] = 1.0;
            } else
                tipPartials[j + s] = 1.0;
            j += sc;
        }
    } else {
        int gI = nodeIndex - data.getReferenceTaxa().getTaxonCount();
        HashSet<Taxon> group = groups.get(gI);
        int internalNum = data.getReadsTaxa().getTaxonCount();
        Taxon firstTax = null;
        boolean peeled = false;
        for (Taxon tax : group) {
            if (firstTax == null) {
                firstTax = tax;
                continue;
            }
            int c2 = data.getReadsTaxa().getTaxonIndex(tax);
            if (!peeled) {
                int c1 = data.getReadsTaxa().getTaxonIndex(firstTax);
                core.setNodePartialsForUpdate(internalNum);
                core.calculatePartials(c1, c2, internalNum);
            } else {
                core.setNodePartialsForUpdate(internalNum);
                core.calculatePartials(internalNum - 1, c2, internalNum);
            }
            internalNum++;
            peeled = true;
        }
        if (group.size() == 0) {
            for (int i = 0; i < tipPartials.length; i++) tipPartials[i] = 1.0;
        } else if (!peeled)
            getPartialsForGroupSizeOne(firstTax, tipPartials);
        else
            core.getPartials(internalNum - 1, tipPartials);
    }
}

Example 27 with Alignment

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

the class ConvertAlignmentParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Alignment alignment = (Alignment) xo.getChild(Alignment.class);
    // Old parser always returned UNIVERSAL type for codon conversion
    DataType dataType = DataTypeUtils.getDataType(xo);
    GeneticCode geneticCode = GeneticCode.UNIVERSAL;
    if (dataType instanceof Codons) {
        geneticCode = ((Codons) dataType).getGeneticCode();
    }
    ConvertAlignment convert = new ConvertAlignment(dataType, geneticCode, alignment);
    Logger.getLogger("dr.evoxml").info("Converted alignment, '" + xo.getId() + "', from " + alignment.getDataType().getDescription() + " to " + dataType.getDescription());
    return convert;
}

Example 28 with Alignment

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

the class DataModelImporter method importBEASTFile.

// xml
private void importBEASTFile(File file, DateGuesser guesser, Map dataModel) throws IOException, ImportException {
    try {
        FileReader reader = new FileReader(file);
        BeastImporter importer = new BeastImporter(reader);
        List<TaxonList> taxonLists = new ArrayList<TaxonList>();
        List<Alignment> alignments = new ArrayList<Alignment>();
        importer.importBEAST(taxonLists, alignments);
        TaxonList taxa = taxonLists.get(0);
        int count = 1;
        for (Alignment alignment : alignments) {
            String name = file.getName();
            if (alignment.getId() != null && alignment.getId().length() > 0) {
                name = alignment.getId();
            } else {
                if (alignments.size() > 1) {
                    name += count;
                }
            }
            setData(dataModel, guesser, name, taxa, taxonLists, alignment, null, null, null);
            count++;
        }
        reader.close();
    } catch (JDOMException e) {
        throw new ImportException(e.getMessage());
    } catch (ImportException e) {
        throw new ImportException(e.getMessage());
    } catch (IOException e) {
        throw new IOException(e.getMessage());
    }
}

Example 29 with Alignment

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

the class RecomboGen method generate.

public Alignment generate() {
    List<Node> tips = new ArrayList<Node>();
    // add all the tips
    for (int i = 0; i < taxa.getTaxonCount(); i++) {
        tips.add(new Node(taxa.getTaxon(i)));
    }
    Collections.sort(tips);
    List<Node> unusedTips = new ArrayList<Node>(tips);
    double time = 0;
    List<Node> nodes = new ArrayList<Node>();
    // Add any tips with zero sampling time.
    List<Node> nodesToAdd = new ArrayList<Node>();
    for (Node tip : unusedTips) {
        if (tip.getTime() == 0.0) {
            nodesToAdd.add(tip);
            System.out.println("Time: " + time + " adding " + tip.getTaxon());
        }
    }
    nodes.addAll(nodesToAdd);
    unusedTips.removeAll(nodesToAdd);
    do {
        boolean tipsAdded;
        do {
            tipsAdded = false;
            final int lineageCount = nodes.size();
            // get time to next event...
            // create unit uniform random variate
            final double U = MathUtils.nextDouble();
            final double interval = -Math.log(U) / (lineageCount * recombinationRate);
            final double nextTime = time + interval;
            // Add any tips for which we have reached their sampling time.
            nodesToAdd.clear();
            for (Node tip : unusedTips) {
                if (tip.getTime() <= nextTime) {
                    nodesToAdd.add(tip);
                    tipsAdded = true;
                    System.out.println("Time: " + tip.getTime() + " adding " + tip.getTaxon());
                    // reset the current time (the tips are sorted in time order
                    // so this will be the oldest added tip).
                    time = tip.getTime();
                }
            }
            nodes.addAll(nodesToAdd);
            unusedTips.removeAll(nodesToAdd);
            if (!tipsAdded) {
                time = nextTime;
            }
        } while (// only continue when no tips are added
        tipsAdded);
        int r = MathUtils.nextInt(nodes.size());
        Node node = nodes.get(r);
        // create two new parent nodes
        Node parent1 = new Node(node, time);
        Node parent2 = new Node(node, time);
        // select a break point in interval [1, length - 2] on
        // a zero-indexed line.
        int breakPoint = MathUtils.nextInt(length - 2) + 1;
        // setup child node with parents and break point
        node.setParent1(parent1);
        node.setParent2(parent2);
        node.setBreakPoint(breakPoint);
        System.out.println("Time: " + time + " recombining " + (node.getTaxon() != null ? node.getTaxon() : r) + " at breakpoint " + breakPoint);
        nodes.add(parent1);
        nodes.add(parent2);
        nodes.remove(node);
    } while (unusedTips.size() > 0);
    // Construct a taxon set for coalescent simulation of deep tree
    Taxa treeTaxa = new Taxa();
    int i = 0;
    Map<Node, Taxon> nodeMap = new HashMap<Node, Taxon>();
    for (Node node : nodes) {
        Taxon taxon = new Taxon("Taxon" + i);
        treeTaxa.addTaxon(taxon);
        nodeMap.put(node, taxon);
        i++;
    }
    CoalescentSimulator coalSim = new CoalescentSimulator();
    ConstantPopulation demo = new ConstantPopulation(Units.Type.YEARS);
    demo.setN0(ancestralPopulationSize);
    Tree tree = coalSim.simulateTree(treeTaxa, demo);
    System.out.println("Tree MRCA " + tree.getNodeHeight(tree.getRoot()) + time);
    SequenceSimulator seqSim = new SequenceSimulator(tree, siteModel, branchRateModel, length);
    Alignment ancestralAlignment = seqSim.simulate();
    SimpleAlignment alignment = new SimpleAlignment();
    // generated recombinant history
    for (Node tip : tips) {
        String seqString = generateRecombinant(tip, nodeMap, ancestralAlignment);
        Sequence sequence = new Sequence(tip.getTaxon(), seqString);
        //            System.out.println(">" + tip.getTaxon() + "\r" + sequence);
        alignment.addSequence(sequence);
    }
    return alignment;
}