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

Example 1 with Alignment

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

the class BeagleTreeLikelihood method main.

public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 1: simulateOnePartition");
        int sequenceLength = 1000;
        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 branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        Parameter kappa2 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        HKY hky2 = new HKY(kappa2, freqModel);
        HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        substitutionModels.add(hky2);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        Parameter epochTimes = new Parameter.Default(1, 20);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 0.001);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogenousBranchSubstitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        Alignment alignment = simulator.simulate(false, false);
        BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
        nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 2 with Alignment

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

the class DataModelImporter method importFastaFile.

// FASTA
private void importFastaFile(File file, DateGuesser guesser, Map dataModel) throws IOException, ImportException {
    try {
        FileReader reader = new FileReader(file);
        FastaImporter importer = new FastaImporter(reader, Nucleotides.INSTANCE);
        Alignment alignment = importer.importAlignment();
        reader.close();
        setData(dataModel, guesser, file.getName(), alignment, null, alignment, null, null, null);
    } catch (ImportException e) {
        throw new ImportException(e.getMessage());
    } catch (IOException e) {
        throw new IOException(e.getMessage());
    }
}

Example 3 with Alignment

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

the class MicrosatelliteSimulator method simulateMsatPattern.

// intArray2Sequence
/**
     * Convert an alignment to a pattern
     */
public Patterns simulateMsatPattern() {
    Alignment align = simulate();
    int[] pattern = new int[align.getTaxonCount()];
    for (int i = 0; i < pattern.length; i++) {
        String taxonName = align.getSequence(i).getTaxon().getId();
        int index = taxa.getTaxonIndex(taxonName);
        pattern[index] = Integer.parseInt(align.getSequence(i).getSequenceString());
    }
    Patterns patterns = new Patterns(dataType, taxa);
    patterns.addPattern(pattern);
    for (int i = 0; i < pattern.length; i++) {
        System.out.print(pattern[i] + ",");
    }
    System.out.println();
    return patterns;
}

Example 4 with Alignment

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

the class HiddenLinkageModel method getPartialsForGroupSizeOne.

private void getPartialsForGroupSizeOne(Taxon tax, double[] tipPartials) {
    Alignment aln = data.getAlignment();
    int sc = aln.getStateCount();
    int index = aln.getTaxonIndex(tax);
    int j = 0;
    for (int i = 0; i < aln.getSiteCount(); i++) {
        int s = aln.getState(index, i);
        if (s >= sc) {
            for (int k = 0; k < sc; k++) tipPartials[j + k] = 1.0;
        } else
            System.arraycopy(internalMatrix, s * sc, tipPartials, j, sc);
        j += sc;
    }
}

Example 5 with Alignment

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

the class LinkageGroupSwap method doOperation.

@Override
public double doOperation() {
    if (MathUtils.nextBoolean()) {
        // swap all taxa in one group to a new group
        int A = MathUtils.nextInt(groupCount);
        int B = MathUtils.nextInt(groupCount);
        HashSet<Taxon> aTaxa = new HashSet<Taxon>(hlm.getGroup(A));
        HashSet<Taxon> bTaxa = new HashSet<Taxon>(hlm.getGroup(B));
        for (Taxon taxon : aTaxa) {
            hlm.moveReadGroup(taxon, A, B);
        }
        for (Taxon taxon : bTaxa) {
            hlm.moveReadGroup(taxon, B, A);
        }
    } else {
        // pick two linkage groups uniformly A, B
        // pick an alignment column uniformly X
        // move all reads in group A that are defined in column X to group B, and
        // vice versa from B to A.
        // pick a read uniformly at random, add it to a linkage group uniformly at random
        ArrayList<Taxon> aTaxa = new ArrayList<Taxon>();
        ArrayList<Taxon> bTaxa = new ArrayList<Taxon>();
        int A = 0, B = 0;
        // choice of X and Y to groups that have something in col X.
        while (aTaxa.size() == 0 && bTaxa.size() == 0) {
            int X = MathUtils.nextInt(columnCount);
            A = MathUtils.nextInt(groupCount);
            B = MathUtils.nextInt(groupCount);
            if (A == B)
                // nothing to do.
                continue;
            // find all reads intersecting column X
            Alignment aln = hlm.getData().getAlignment();
            for (int i = 0; i < aln.getTaxonCount(); i++) {
                int state = aln.getPatternState(i, X);
                if (state == hlm.getDataType().getGapState() || state == hlm.getDataType().getUnknownState())
                    // seq undefined in this column
                    continue;
                if (hlm.getGroup(A).contains(aln.getTaxon(i))) {
                    aTaxa.add(aln.getTaxon(i));
                }
                if (hlm.getGroup(B).contains(aln.getTaxon(i))) {
                    bTaxa.add(aln.getTaxon(i));
                }
            }
        }
        // move taxa from A to B and from B to A
        for (Taxon taxon : aTaxa) {
            hlm.moveReadGroup(taxon, A, B);
        }
        for (Taxon taxon : bTaxa) {
            hlm.moveReadGroup(taxon, B, A);
        }
    }
    return 0;
}