use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.
the class BeagleSequenceSimulatorParser method parseXMLObject.
// END: getSyntaxRules
@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
String msg = "";
boolean parallel = false;
boolean outputAncestralSequences = false;
if (xo.hasAttribute(PARALLEL)) {
parallel = xo.getBooleanAttribute(PARALLEL);
}
if (xo.hasAttribute(OUTPUT_ANCESTRAL_SEQUENCES)) {
outputAncestralSequences = xo.getBooleanAttribute(OUTPUT_ANCESTRAL_SEQUENCES);
}
SimpleAlignment.OutputType output = SimpleAlignment.OutputType.FASTA;
if (xo.hasAttribute(OUTPUT)) {
output = SimpleAlignment.OutputType.parseFromString(xo.getStringAttribute(OUTPUT));
}
int siteCount = 0;
int to = 0;
for (int i = 0; i < xo.getChildCount(); i++) {
Partition partition = (Partition) xo.getChild(i);
to = partition.to + 1;
if (to > siteCount) {
siteCount = to;
}
}
// END: partitions loop
ArrayList<Partition> partitionsList = new ArrayList<Partition>();
for (int i = 0; i < xo.getChildCount(); i++) {
Partition partition = (Partition) xo.getChild(i);
if (partition.from > siteCount) {
throw new XMLParseException("Illegal 'from' attribute in " + PartitionParser.PARTITION + " element");
}
if (partition.to > siteCount) {
throw new XMLParseException("Illegal 'to' attribute in " + PartitionParser.PARTITION + " element");
}
if (partition.to == -1) {
partition.to = siteCount - 1;
}
if (partition.getRootSequence() != null) {
// TODO: what about 'every'?
int partitionSiteCount = (partition.to - partition.from) + 1;
if (partition.getRootSequence().getLength() != 3 * partitionSiteCount && partition.getFreqModel().getDataType() instanceof Codons) {
throw new RuntimeException("Root codon sequence " + "for partition " + (i + 1) + " has " + partition.getRootSequence().getLength() + " characters " + "expecting " + 3 * partitionSiteCount + " characters");
} else if (partition.getRootSequence().getLength() != partitionSiteCount && partition.getFreqModel().getDataType() instanceof Nucleotides) {
throw new RuntimeException("Root nuleotide sequence " + "for partition " + (i + 1) + " has " + partition.getRootSequence().getLength() + " characters " + "expecting " + partitionSiteCount + " characters");
}
// END: dataType check
// System.exit(-1);
}
// END: ancestralSequence check
partitionsList.add(partition);
}
// END: partitions loop
msg += "\n\t" + siteCount + ((siteCount > 1) ? " replications " : " replication");
if (msg.length() > 0) {
Logger.getLogger("dr.app.beagle.tools").info("Using Beagle Sequence Simulator: " + msg);
}
BeagleSequenceSimulator s = new BeagleSequenceSimulator(partitionsList);
SimpleAlignment alignment = s.simulate(parallel, outputAncestralSequences);
alignment.setOutputType(output);
return alignment;
}
use of dr.app.beagle.tools.Partition in project beast-mcmc by beast-dev.
the class PartitionParser method parseXMLObject.
@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
int from = 0;
int to = -1;
int every = xo.getAttribute(EVERY, 1);
if (xo.hasAttribute(FROM)) {
from = xo.getIntegerAttribute(FROM) - 1;
if (from < 0) {
throw new XMLParseException("Illegal 'from' attribute in patterns element");
}
}
if (xo.hasAttribute(TO)) {
to = xo.getIntegerAttribute(TO) - 1;
if (to < 0 || to < from) {
throw new XMLParseException("Illegal 'to' attribute in patterns element");
}
}
if (every <= 0) {
throw new XMLParseException("Illegal 'every' attribute in patterns element");
}
// END: every check
TreeModel tree = (TreeModel) xo.getChild(TreeModel.class);
GammaSiteRateModel siteModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
FrequencyModel freqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
Sequence rootSequence = (Sequence) xo.getChild(Sequence.class);
BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
if (rateModel == null) {
rateModel = new DefaultBranchRateModel();
}
BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
if (branchModel == null) {
SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
branchModel = new HomogeneousBranchModel(substitutionModel);
}
Partition partition = new Partition(tree, branchModel, siteModel, rateModel, freqModel, from, to, every);
if (rootSequence != null) {
partition.setRootSequence(rootSequence);
}
return partition;
}
use of dr.app.beagle.tools.Partition 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.app.beagle.tools.Partition 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.app.beagle.tools.Partition 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
}
Aggregations