use of dr.evolution.io.NewickImporter 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
}
use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.
the class BeagleSeqSimTest method simulateThreePartitions.
// END: simulateTwoPartitions
static void simulateThreePartitions(int i, int N) {
try {
MathUtils.setSeed(666);
System.out.println("Test case 3: simulateThreePartitions");
int sequenceLength = 100000;
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
3);
// create partition
Partition Partition = new //
Partition(//
treeModel, //
substitutionModel, //
siteRateModel, //
branchRateModel, //
freqModel, // from
1, // to
sequenceLength - 1, // every
3);
// create partition
Partition partition3 = new //
Partition(//
treeModel, //
substitutionModel, //
siteRateModel, //
branchRateModel, //
freqModel, // from
2, // to
sequenceLength - 1, // every
3);
partitionsList.add(partition1);
partitionsList.add(Partition);
partitionsList.add(partition3);
// feed to sequence simulator and generate data
BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
if (i == (N - 1)) {
System.out.println(simulator.simulate(simulateInPar, false).toString());
} else {
simulator.simulate(simulateInPar, false);
}
} catch (Exception e) {
e.printStackTrace();
System.exit(-1);
}
// END: try-catch block
}
use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.
the class BeagleSeqSimTest method ancestralSequenceTree.
//END: simulateRandomBranchAssignment
static void ancestralSequenceTree() {
try {
LinkedHashMap<NodeRef, int[]> sequenceMap = new LinkedHashMap<NodeRef, int[]>();
DataType dataType = Nucleotides.INSTANCE;
// 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);
for (NodeRef node : treeModel.getNodes()) {
if (treeModel.isExternal(node)) {
int[] seq = new int[] { 1, 1, 1 };
sequenceMap.put(node, seq);
} else {
int[] seq = new int[] { 2, 2, 2 };
sequenceMap.put(node, seq);
}
}
// END: nodes loop
AncestralSequenceTrait ancestralSequence = new AncestralSequenceTrait(sequenceMap, dataType);
TreeTraitProvider[] treeTraitProviders = new TreeTraitProvider[] { ancestralSequence };
StringBuffer buffer = new StringBuffer();
NumberFormat format = NumberFormat.getNumberInstance(Locale.ENGLISH);
boolean useTipLabels = true;
//
TreeUtils.newick(//
treeModel, //
treeModel.getRoot(), //
useTipLabels, //
TreeUtils.BranchLengthType.LENGTHS_AS_TIME, //
format, //
null, //
treeTraitProviders, null, buffer);
System.out.println(buffer);
} catch (IOException e) {
e.printStackTrace();
} catch (ImportException e) {
e.printStackTrace();
}
// END: try-catch
}
use of dr.evolution.io.NewickImporter in project beast-mcmc by beast-dev.
the class BeagleSeqSimTest method simulateTopology.
// END: annotateTree
static void simulateTopology() {
try {
System.out.println("Test case 1: simulateTopology");
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);
// set demographic function
ExponentialGrowth exponentialGrowth = new ExponentialGrowth(Units.Type.YEARS);
exponentialGrowth.setN0(10);
exponentialGrowth.setGrowthRate(0.5);
Taxa taxa = new Taxa();
for (Taxon taxon : tree.asList()) {
double absoluteHeight = Utils.getAbsoluteTaxonHeight(taxon, tree);
taxon.setAttribute(Utils.ABSOLUTE_HEIGHT, absoluteHeight);
// taxon.setAttribute("date", new Date(absoluteHeight,
// Units.Type.YEARS, true));
taxa.addTaxon(taxon);
}
// END: taxon loop
CoalescentSimulator topologySimulator = new CoalescentSimulator();
TreeModel treeModel = new TreeModel(topologySimulator.simulateTree(taxa, exponentialGrowth));
System.out.println(treeModel.toString());
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);
partitionsList.add(partition1);
// feed to sequence simulator and generate data
BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
System.out.println(simulator.simulate(simulateInPar, false).toString());
} catch (Exception e) {
e.printStackTrace();
System.exit(-1);
}
// END: try-catch block
}
use of dr.evolution.io.NewickImporter 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
}
Aggregations