Examples with Sequence dr.evolution.sequence.Sequence used on opensource projects

Example 16 with Sequence

use of dr.evolution.sequence.Sequence in project beast-mcmc by beast-dev.

the class Coevolve method makePairAlignment.

/**
     * @param alignment
     * @param dist the distance |i-j| between the site pairs (j = i + dist).
     * @param start the first i value considered
     * @param step the distance between successive i values (i += step)
     * @return an alignment of sites, where each site corresponds to a pair (i, j) in the old alignment.
     */
private Alignment makePairAlignment(Alignment alignment, int dist, int start, int step) {
    DataType dataType = alignment.getDataType();
    int stateCount = dataType.getStateCount();
    DataType newDataType = makePairDataType(alignment);
    SimpleAlignment pairAlignment = new SimpleAlignment();
    pairAlignment.setDataType(newDataType);
    for (int k = 0; k < alignment.getSequenceCount(); k++) {
        StringBuffer sequence = new StringBuffer();
        for (int i = start; (i + dist) < alignment.getSiteCount(); i += step) {
            int j = i + dist;
            int state = alignment.getState(k, i) * stateCount + alignment.getState(k, j);
            if (state < newDataType.getStateCount()) {
                sequence.append(newDataType.getChar(state));
            } else {
                sequence.append('?');
            }
        }
        pairAlignment.addSequence(new Sequence(alignment.getTaxon(k), sequence.toString()));
    }
    return pairAlignment;
}

Example 17 with Sequence

use of dr.evolution.sequence.Sequence in project beast-mcmc by beast-dev.

the class NexusExporter method exportAlignment.

public String exportAlignment(Alignment alignment) throws IOException, IllegalArgumentException {
    StringBuffer buffer = new StringBuffer();
    DataType dataType = null;
    int seqLength = 0;
    for (int i = 0; i < alignment.getSequenceCount(); i++) {
        Sequence sequence = alignment.getSequence(i);
        if (sequence.getLength() > seqLength) {
            seqLength = sequence.getLength();
        }
        if (dataType == null) {
            dataType = sequence.getDataType();
        } else if (dataType != sequence.getDataType()) {
            throw new RuntimeException("Sequences must have the same data type.");
        }
    // END: dataType check
    }
    // END: sequences loop
    buffer.append("#NEXUS\n");
    buffer.append("begin data;\n");
    buffer.append("\tdimensions" + " " + "ntax=" + alignment.getTaxonCount() + " " + "nchar=" + seqLength + ";\n");
    buffer.append("\tformat datatype=" + dataType.getDescription() + " missing=" + DataType.UNKNOWN_CHARACTER + " gap=" + DataType.GAP_CHARACTER + ";\n");
    buffer.append("\tmatrix\n");
    int maxRowLength = seqLength;
    for (int n = 0; n < Math.ceil((double) seqLength / maxRowLength); n++) {
        for (int i = 0; i < alignment.getSequenceCount(); i++) {
            Sequence sequence = alignment.getSequence(i);
            StringBuilder builder = new StringBuilder("\t");
            appendTaxonName(sequence.getTaxon(), builder);
            String sequenceString = sequence.getSequenceString();
            builder.append("\t").append(sequenceString.subSequence(n * maxRowLength, Math.min((n + 1) * maxRowLength, sequenceString.length())));
            int shortBy = Math.min(Math.min(n * maxRowLength, seqLength) - sequence.getLength(), maxRowLength);
            if (shortBy > 0) {
                for (int j = 0; j < shortBy; j++) {
                    builder.append(DataType.GAP_CHARACTER);
                }
            }
            buffer.append(builder + "\n");
        }
    // END: sequences loop
    }
    buffer.append(";\nend;");
    return buffer.toString();
}

Example 18 with Sequence

use of dr.evolution.sequence.Sequence 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
}

Example 19 with Sequence

use of dr.evolution.sequence.Sequence 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;
}

Example 20 with Sequence

use of dr.evolution.sequence.Sequence in project beast-mcmc by beast-dev.

the class SequenceSimulator method main.

// getDefaultSiteModel
public static void main(String[] args) {
    try {
        int nReplications = 10;
        // create tree
        NewickImporter importer = new NewickImporter("((A:1.0,B:1.0)AB:1.0,(C:1.0,D:1.0)CD:1.0)ABCD;");
        Tree tree = importer.importTree(null);
        // create site model
        SiteModel siteModel = getDefaultSiteModel();
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // feed to sequence simulator and generate leaves
        SequenceSimulator treeSimulator = new SequenceSimulator(tree, siteModel, branchRateModel, nReplications);
        Sequence ancestralSequence = new Sequence();
        ancestralSequence.appendSequenceString("TCAGGTCAAG");
        treeSimulator.setAncestralSequence(ancestralSequence);
        System.out.println(treeSimulator.simulate().toString());
    } catch (Exception e) {
        e.printStackTrace();
    }
//END: try-catch block
}