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

Example 21 with Alignment

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

the class BeagleSeqSimTest method simulateCodon.

// END: simulateAminoAcid
static void simulateCodon() {
    try {
        boolean calculateLikelihood = true;
        System.out.println("Test case 6: simulate codons");
        MathUtils.setSeed(666);
        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 site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(Utils.UNIFORM_CODON_FREQUENCIES);
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create substitution model
        Parameter alpha = new Parameter.Default(1, 10);
        Parameter beta = new Parameter.Default(1, 5);
        //			Parameter kappa = new Parameter.Default(1, 1);
        MG94CodonModel mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(mg94);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        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(simulateInPar, false);
        System.out.println(alignment.toString());
        if (calculateLikelihood) {
            // NewBeagleSequenceLikelihood nbtl = new
            // NewBeagleSequenceLikelihood(alignment, treeModel,
            // substitutionModel, (SiteModel) siteRateModel,
            // branchRateModel, null, false,
            // PartialsRescalingScheme.DEFAULT);
            ConvertAlignment convert = new ConvertAlignment(Nucleotides.INSTANCE, GeneticCode.UNIVERSAL, alignment);
            BeagleTreeLikelihood nbtl = new //
            BeagleTreeLikelihood(//
            convert, //
            treeModel, //
            substitutionModel, //
            siteRateModel, //
            branchRateModel, //
            null, //
            false, PartialsRescalingScheme.DEFAULT, true);
            System.out.println("likelihood = " + nbtl.getLogLikelihood());
        }
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch
}

Example 22 with Alignment

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

the class AlignmentScore method main.

public static void main(String[] args) throws java.io.IOException, Importer.ImportException {
    NexusImporter importer = new NexusImporter(new FileReader(args[0]));
    Alignment alignment = importer.importAlignment();
    ExtractPairs pairs = new ExtractPairs(alignment);
    Parameter muParam = new Parameter.Default(1.0);
    Parameter kappaParam = new Parameter.Default(1.0);
    kappaParam.addBounds(new Parameter.DefaultBounds(100.0, 0.0, 1));
    muParam.addBounds(new Parameter.DefaultBounds(1.0, 1.0, 1));
    Parameter freqParam = new Parameter.Default(alignment.getStateFrequencies());
    FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqParam);
    SubstitutionModel substModel = new HKY(kappaParam, freqModel);
    SiteModel siteModel = new GammaSiteModel(substModel, muParam, null, 1, null);
    ScoreMatrix scoreMatrix = new ScoreMatrix(siteModel, 0.1);
    double threshold = 0.1;
    List<PairDistance> pairDistances = new ArrayList<PairDistance>();
    Set<Integer> sequencesUsed = new HashSet<Integer>();
    List<Integer> allGaps = new ArrayList<Integer>();
    for (int i = 0; i < alignment.getSequenceCount(); i++) {
        for (int j = i + 1; j < alignment.getSequenceCount(); j++) {
            Alignment pairAlignment = pairs.getPairAlignment(i, j);
            if (pairAlignment != null) {
                SitePatterns patterns = new SitePatterns(pairAlignment);
                double distance = getGeneticDistance(scoreMatrix, patterns);
                if (distance < threshold) {
                    List gaps = new ArrayList();
                    GapUtils.getGapSizes(pairAlignment, gaps);
                    pairDistances.add(new PairDistance(i, j, distance, gaps, pairAlignment.getSiteCount()));
                    System.out.print(".");
                } else {
                    System.out.print("*");
                }
            } else {
                System.out.print("x");
            }
        }
        System.out.println();
    }
    Collections.sort(pairDistances);
    int totalLength = 0;
    for (PairDistance pairDistance : pairDistances) {
        Integer x = pairDistance.x;
        Integer y = pairDistance.y;
        if (!sequencesUsed.contains(x) && !sequencesUsed.contains(y)) {
            allGaps.addAll(pairDistance.gaps);
            sequencesUsed.add(x);
            sequencesUsed.add(y);
            System.out.println("Added pair (" + x + "," + y + ") d=" + pairDistance.distance + " L=" + pairDistance.alignmentLength);
            totalLength += pairDistance.alignmentLength;
        }
    }
    printFrequencyTable(allGaps);
    System.out.println("total length=" + totalLength);
}

Example 23 with Alignment

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

the class BeautiTesterConfig method importFromFile.

public void importFromFile(String fileName, BeautiOptions options, boolean translate) {
    TaxonList taxa = null;
    Alignment alignment = null;
    Tree tree = null;
    PartitionSubstitutionModel model = null;
    java.util.List<NexusApplicationImporter.CharSet> charSets = new ArrayList<NexusApplicationImporter.CharSet>();
    try {
        FileReader reader = new FileReader(fileName);
        NexusApplicationImporter importer = new NexusApplicationImporter(reader);
        boolean done = false;
        while (!done) {
            try {
                NexusImporter.NexusBlock block = importer.findNextBlock();
                if (block == NexusImporter.TAXA_BLOCK) {
                    if (taxa != null) {
                        throw new NexusImporter.MissingBlockException("TAXA block already defined");
                    }
                    taxa = importer.parseTaxaBlock();
                } else if (block == NexusImporter.CALIBRATION_BLOCK) {
                    if (taxa == null) {
                        throw new NexusImporter.MissingBlockException("TAXA or DATA block must be defined before a CALIBRATION block");
                    }
                    importer.parseCalibrationBlock(taxa);
                } else if (block == NexusImporter.CHARACTERS_BLOCK) {
                    if (taxa == null) {
                        throw new NexusImporter.MissingBlockException("TAXA block must be defined before a CHARACTERS block");
                    }
                    if (alignment != null) {
                        throw new NexusImporter.MissingBlockException("CHARACTERS or DATA block already defined");
                    }
                    alignment = importer.parseCharactersBlock(options.taxonList);
                } else if (block == NexusImporter.DATA_BLOCK) {
                    if (alignment != null) {
                        throw new NexusImporter.MissingBlockException("CHARACTERS or DATA block already defined");
                    }
                    // A data block doesn't need a taxon block before it
                    // but if one exists then it will use it.
                    alignment = importer.parseDataBlock(options.taxonList);
                    if (taxa == null) {
                        taxa = alignment;
                    }
                } else if (block == NexusImporter.TREES_BLOCK) {
                    if (taxa == null) {
                        throw new NexusImporter.MissingBlockException("TAXA or DATA block must be defined before a TREES block");
                    }
                    if (tree != null) {
                        throw new NexusImporter.MissingBlockException("TREES block already defined");
                    }
                    Tree[] trees = importer.parseTreesBlock(taxa);
                    if (trees.length > 0) {
                        tree = trees[0];
                    }
                } else if (block == NexusApplicationImporter.PAUP_BLOCK) {
                    model = importer.parsePAUPBlock(options, charSets);
                } else if (block == NexusApplicationImporter.MRBAYES_BLOCK) {
                    model = importer.parseMrBayesBlock(options, charSets);
                } else if (block == NexusApplicationImporter.ASSUMPTIONS_BLOCK) {
                    importer.parseAssumptionsBlock(charSets);
                } else {
                // Ignore the block..
                }
            } catch (EOFException ex) {
                done = true;
            }
        }
        // Allow the user to load taxa only (perhaps from a tree file) so that they can sample from a prior...
        if (alignment == null && taxa == null) {
            throw new NexusImporter.MissingBlockException("TAXON, DATA or CHARACTERS block is missing");
        }
    } catch (FileNotFoundException fnfe) {
        System.err.println("File not found: " + fnfe);
        System.exit(1);
    } catch (Importer.ImportException ime) {
        System.err.println("Error parsing imported file: " + ime);
        System.exit(1);
    } catch (IOException ioex) {
        System.err.println("File I/O Error: " + ioex);
        System.exit(1);
    } catch (Exception ex) {
        System.err.println("Fatal exception: " + ex);
        System.exit(1);
    }
    if (options.taxonList == null) {
        // This is the first partition to be loaded...
        options.taxonList = new Taxa(taxa);
        // check the taxon names for invalid characters
        boolean foundAmp = false;
        for (int i = 0; i < taxa.getTaxonCount(); i++) {
            String name = taxa.getTaxon(i).getId();
            if (name.indexOf('&') >= 0) {
                foundAmp = true;
            }
        }
        if (foundAmp) {
            System.err.println("One or more taxon names include an illegal character ('&').");
            return;
        }
        // make sure they all have dates...
        for (int i = 0; i < taxa.getTaxonCount(); i++) {
            if (taxa.getTaxonAttribute(i, "date") == null) {
                java.util.Date origin = new java.util.Date(0);
                dr.evolution.util.Date date = dr.evolution.util.Date.createTimeSinceOrigin(0.0, Units.Type.YEARS, origin);
                taxa.getTaxon(i).setAttribute("date", date);
            }
        }
    } else {
        // This is an additional partition so check it uses the same taxa
        java.util.List<String> oldTaxa = new ArrayList<String>();
        for (int i = 0; i < options.taxonList.getTaxonCount(); i++) {
            oldTaxa.add(options.taxonList.getTaxon(i).getId());
        }
        java.util.List<String> newTaxa = new ArrayList<String>();
        for (int i = 0; i < taxa.getTaxonCount(); i++) {
            newTaxa.add(taxa.getTaxon(i).getId());
        }
        if (!(oldTaxa.containsAll(newTaxa) && oldTaxa.size() == newTaxa.size())) {
            System.err.println("This file contains different taxa from the previously loaded");
            return;
        }
    }
    String fileNameStem = dr.app.util.Utils.trimExtensions(fileName, new String[] { "NEX", "NEXUS", "TRE", "TREE" });
    if (alignment != null) {
        if (translate) {
            alignment = new ConvertAlignment(AminoAcids.INSTANCE, GeneticCode.UNIVERSAL, alignment);
        }
        java.util.List<PartitionData> partitions = new ArrayList<PartitionData>();
        if (charSets != null && charSets.size() > 0) {
            for (NexusApplicationImporter.CharSet charSet : charSets) {
                partitions.add(new PartitionData(options, charSet.getName(), fileName, new CharSetAlignment(charSet, alignment)));
            }
        } else {
            partitions.add(new PartitionData(options, fileNameStem, fileName, alignment));
        }
        for (PartitionData partition : partitions) {
            if (model != null) {
                partition.setPartitionSubstitutionModel(model);
            //                    options.addPartitionSubstitutionModel(model);
            } else {
                for (PartitionSubstitutionModel pm : options.getPartitionSubstitutionModels()) {
                    if (pm.getDataType() == alignment.getDataType()) {
                        partition.setPartitionSubstitutionModel(pm);
                    }
                }
                if (partition.getPartitionSubstitutionModel() == null) {
                    PartitionSubstitutionModel pm = new PartitionSubstitutionModel(options, partition);
                    partition.setPartitionSubstitutionModel(pm);
                //                        options.addPartitionSubstitutionModel(pm);
                }
            }
            options.dataPartitions.add(partition);
        }
    }
    calculateHeights(options);
}

Example 24 with Alignment

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

the class BeagleBranchLikelihood method main.

// END: LikelihoodColumn class
// ////////////
// ---TEST---//
// ////////////
public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("((SimSeq1:22.0,SimSeq2:22.0):12.0,(SimSeq3:23.1,SimSeq4:23.1):10.899999999999999);");
        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);
        HKY hky1 = new HKY(kappa1, freqModel);
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 1.000);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogeneousBranchModel, //
        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);
        System.out.println(alignment);
        BeagleTreeLikelihood btl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("BTL(homogeneous) = " + btl.getLogLikelihood());
        BeagleBranchLikelihood bbl = new BeagleBranchLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, freqModel, branchRateModel);
        int branchIndex = 4;
        System.out.println(bbl.getBranchLogLikelihood(branchIndex));
        bbl.finalizeBeagle();
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    } catch (Throwable e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 25 with Alignment

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

the class LineageSpecificBranchModel method main.

// END: acceptState
public static void main(String[] args) {
    try {
        // the seed of the BEAST
        MathUtils.setSeed(666);
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        TreeModel tree = new TreeModel(importer.importTree(null));
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { //
        0.0163936, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344 });
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create substitution model
        Parameter alpha = new Parameter.Default(1, 10);
        Parameter beta = new Parameter.Default(1, 5);
        MG94CodonModel mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(mg94);
        // create partition
        Partition partition1 = new //
        Partition(//
        tree, //
        substitutionModel, //
        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);
        ConvertAlignment convert = new ConvertAlignment(Nucleotides.INSTANCE, GeneticCode.UNIVERSAL, alignment);
        List<SubstitutionModel> substModels = new ArrayList<SubstitutionModel>();
        for (int i = 0; i < 2; i++) {
            //				alpha = new Parameter.Default(1, 10 );
            //				beta = new Parameter.Default(1, 5 );
            //				mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta,
            //						freqModel);
            substModels.add(mg94);
        }
        Parameter uCategories = new Parameter.Default(2, 0);
        //            CountableBranchCategoryProvider provider = new CountableBranchCategoryProvider.IndependentBranchCategoryModel(tree, uCategories);
        LineageSpecificBranchModel branchSpecific = new //provider, 
        LineageSpecificBranchModel(//provider, 
        tree, //provider, 
        freqModel, //provider, 
        substModels, uCategories);
        BeagleTreeLikelihood like = new //
        BeagleTreeLikelihood(//
        convert, //
        tree, //
        branchSpecific, //
        siteRateModel, //
        branchRateModel, //
        null, //
        false, PartialsRescalingScheme.DEFAULT, true);
        BeagleTreeLikelihood gold = new //
        BeagleTreeLikelihood(//
        convert, //
        tree, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        null, //
        false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("likelihood (gold) = " + gold.getLogLikelihood());
        System.out.println("likelihood = " + like.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
    }
}