use of dr.evolution.alignment.SimpleAlignment 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
}
use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.
the class IstvanOperator method doOperation.
public double doOperation() {
Tree tree = likelihood.getTreeModel();
alignment = likelihood.getAlignment();
//System.out.println("Incoming alignment");
//System.out.println(alignment);
//System.out.println();
SubstitutionModel substModel = likelihood.getSiteModel().getSubstitutionModel();
// initialize the iParent and iTau arrays based on the given tree.
initTree(tree, likelihood.getSiteModel().getMu());
int[] treeIndex = new int[tree.getTaxonCount()];
for (int i = 0; i < treeIndex.length; i++) {
treeIndex[i] = tree.getTaxonIndex(alignment.getTaxonId(i));
}
// initialize the iAlignment array from the given alignment.
initAlignment(alignment, treeIndex);
// initialize the iProbs array from the substitution model -- must be called after populating tree!
initSubstitutionModel(substModel);
DataType dataType = substModel.getDataType();
proposal.setGapSymbol(dataType.getGapState());
int[][] returnedAlignment = new int[iAlignment.length][];
//System.out.println("Initialization done, starting proposal proper...");
double logq = proposal.propose(iAlignment, iProbs, iBaseFreqs, iParent, iTau, returnedAlignment, tuning, exponent, gapPenalty);
//System.out.println("Proposal finished, logq=" + logq);
//create new alignment object
SimpleAlignment newAlignment = new SimpleAlignment();
for (int i = 0; i < alignment.getTaxonCount(); i++) {
StringBuffer seqBuffer = new StringBuffer();
for (int j = 0; j < returnedAlignment[i].length; j++) {
seqBuffer.append(dataType.getChar(returnedAlignment[treeIndex[i]][j]));
}
// add sequences in order of tree
String seqString = seqBuffer.toString();
Sequence sequence = new Sequence(alignment.getTaxon(i), seqString);
newAlignment.addSequence(sequence);
String oldunaligned = alignment.getUnalignedSequenceString(i);
String unaligned = newAlignment.getUnalignedSequenceString(i);
if (!unaligned.equals(oldunaligned)) {
System.err.println("Sequence changed from:");
System.err.println("old:'" + oldunaligned + "'");
System.err.println("new:'" + unaligned + "'");
throw new RuntimeException();
}
}
//System.out.println("Outgoing alignment");
//System.out.println(newAlignment);
//System.out.println();
likelihood.setAlignment(newAlignment);
return logq;
}
use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.
the class AncestralStateBeagleTreeLikelihoodTest method testJointLikelihood.
public void testJointLikelihood() {
TreeModel treeModel = new TreeModel("treeModel", tree);
Sequence[] sequence = new Sequence[3];
sequence[0] = new Sequence(new Taxon("0"), "A");
sequence[1] = new Sequence(new Taxon("1"), "C");
sequence[2] = new Sequence(new Taxon("2"), "C");
Taxa taxa = new Taxa();
for (Sequence s : sequence) {
taxa.addTaxon(s.getTaxon());
}
SimpleAlignment alignment = new SimpleAlignment();
for (Sequence s : sequence) {
alignment.addSequence(s);
}
Parameter mu = new Parameter.Default(1, 1.0);
Parameter kappa = new Parameter.Default(1, 1.0);
double[] pi = { 0.25, 0.25, 0.25, 0.25 };
Parameter freqs = new Parameter.Default(pi);
FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
HKY hky = new HKY(kappa, f);
GammaSiteRateModel siteRateModel = new GammaSiteRateModel("gammaModel", mu, null, -1, null);
siteRateModel.setSubstitutionModel(hky);
BranchModel branchModel = new HomogeneousBranchModel(siteRateModel.getSubstitutionModel());
BranchRateModel branchRateModel = null;
AncestralStateBeagleTreeLikelihood treeLikelihood = new AncestralStateBeagleTreeLikelihood(alignment, treeModel, branchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true, null, hky.getDataType(), "stateTag", // useMap = true
true, false);
double logLike = treeLikelihood.getLogLikelihood();
StringBuffer buffer = new StringBuffer();
// Tree.Utils.newick(treeModel, treeModel.getRoot(), false, Tree.BranchLengthType.LENGTHS_AS_TIME,
// null, null, new NodeAttributeProvider[]{treeLikelihood}, null, null, buffer);
TreeUtils.newick(treeModel, treeModel.getRoot(), false, TreeUtils.BranchLengthType.LENGTHS_AS_TIME, null, null, new TreeTraitProvider[] { treeLikelihood }, null, buffer);
System.out.println(buffer);
System.out.println("t_CA(2) = " + t(false, 2.0));
System.out.println("t_CC(1) = " + t(true, 1.0));
double trueValue = 0.25 * t(false, 2.0) * Math.pow(t(true, 1.0), 3.0);
assertEquals(logLike, Math.log(trueValue), 1e-6);
}
use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.
the class NexusImporter method readCharactersBlock.
/**
* Reads a 'CHARACTERS' block.
*/
private Alignment readCharactersBlock(TaxonList taxonList) throws ImportException, IOException {
siteCount = 0;
dataType = null;
readDataBlockHeader("MATRIX", CHARACTERS_BLOCK);
SimpleAlignment alignment = new SimpleAlignment();
readSequenceData(alignment, taxonList);
alignment.updateSiteCount();
findEndBlock();
return alignment;
}
use of dr.evolution.alignment.SimpleAlignment in project beast-mcmc by beast-dev.
the class PhylipSequentialImporter method importAlignment.
/**
* importAlignment.
*/
public Alignment importAlignment() throws IOException, Importer.ImportException {
SimpleAlignment alignment = null;
try {
int taxonCount = readInteger();
int siteCount = readInteger();
String firstSeq = null;
for (int i = 0; i < taxonCount; i++) {
StringBuffer name = new StringBuffer();
char ch = read();
int n = 0;
while (!Character.isWhitespace(ch) && (maxNameLength < 1 || n < maxNameLength)) {
name.append(ch);
ch = read();
n++;
}
StringBuffer seq = new StringBuffer(siteCount);
readSequence(seq, dataType, "", siteCount, "-", "?", ".", firstSeq);
if (firstSeq == null) {
firstSeq = seq.toString();
}
if (alignment == null) {
alignment = new SimpleAlignment();
}
alignment.addSequence(new Sequence(new Taxon(name.toString()), seq.toString()));
}
} catch (EOFException e) {
}
return alignment;
}
Aggregations