use of dr.evolution.alignment.ConvertAlignment in project beast-mcmc by beast-dev.
the class AncestralSequenceAnnotator method processTree.
private Tree processTree(Tree tree) {
// Remake tree to fix node ordering - Marc
GammaSiteRateModel siteModel = loadSiteModel(tree);
SimpleAlignment alignment = new SimpleAlignment();
alignment.setDataType(siteModel.getSubstitutionModel().getDataType());
if (siteModel.getSubstitutionModel().getDataType().getClass().equals(Codons.class)) {
//System.out.println("trololo");
alignment.setDataType(Nucleotides.INSTANCE);
}
//System.out.println("BOO BOO " + siteModel.getSubstitutionModel().getDataType().getClass().getName()+"\t" + Codons.UNIVERSAL.getClass().getName() + "\t" + alignment.getDataType().getClass().getName());
// Get sequences
String[] sequence = new String[tree.getNodeCount()];
for (int i = 0; i < tree.getNodeCount(); i++) {
NodeRef node = tree.getNode(i);
sequence[i] = (String) tree.getNodeAttribute(node, SEQ_STRING);
if (tree.isExternal(node)) {
Taxon taxon = tree.getNodeTaxon(node);
alignment.addSequence(new Sequence(taxon, sequence[i]));
//System.out.println("seq " + sequence[i]);
}
}
// Make evolutionary model
BranchRateModel rateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
FlexibleTree flexTree;
if (siteModel.getSubstitutionModel().getDataType().getClass().equals(Codons.class)) {
ConvertAlignment convertAlignment = new ConvertAlignment(siteModel.getSubstitutionModel().getDataType(), ((Codons) siteModel.getSubstitutionModel().getDataType()).getGeneticCode(), alignment);
flexTree = sampleTree(tree, convertAlignment, siteModel, rateModel);
//flexTree = sampleTree(tree, alignment, siteModel, rateModel);
} else {
flexTree = sampleTree(tree, alignment, siteModel, rateModel);
}
introduceGaps(flexTree, tree);
return flexTree;
}
use of dr.evolution.alignment.ConvertAlignment 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.alignment.ConvertAlignment 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);
}
use of dr.evolution.alignment.ConvertAlignment 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();
}
}
use of dr.evolution.alignment.ConvertAlignment in project beast-mcmc by beast-dev.
the class ConvertAlignmentParser method parseXMLObject.
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
Alignment alignment = (Alignment) xo.getChild(Alignment.class);
// Old parser always returned UNIVERSAL type for codon conversion
DataType dataType = DataTypeUtils.getDataType(xo);
GeneticCode geneticCode = GeneticCode.UNIVERSAL;
if (dataType instanceof Codons) {
geneticCode = ((Codons) dataType).getGeneticCode();
}
ConvertAlignment convert = new ConvertAlignment(dataType, geneticCode, alignment);
Logger.getLogger("dr.evoxml").info("Converted alignment, '" + xo.getId() + "', from " + alignment.getDataType().getDescription() + " to " + dataType.getDescription());
return convert;
}
Aggregations