use of beast.util.ClusterTree in project bacter by tgvaughan.
the class SimulatedAlignmentTest method test.
@Test
public void test() throws Exception {
// Randomizer.setSeed(26);
Randomizer.setSeed(7);
Locus locus = new Locus("locus", 100000);
TaxonSet taxonSet = getTaxonSet(10);
ConstantPopulation popFunc = new ConstantPopulation();
popFunc.initByName("popSize", new RealParameter("1.0"));
ConversionGraph acg = new SimulatedACG();
acg.initByName("rho", 1.0 / 100000, "delta", 1000.0, "populationModel", popFunc, "locus", locus, "taxonset", taxonSet);
System.out.println(acg);
// Site model:
JukesCantor jc = new JukesCantor();
jc.initByName();
SiteModel siteModel = new SiteModel();
siteModel.initByName("mutationRate", new RealParameter("1.0"), "substModel", jc);
// Simulate alignment:
SimulatedAlignment alignment = new SimulatedAlignment();
alignment.initByName("acg", acg, "siteModel", siteModel, "outputFileName", "simulated_alignment.nexus", "useNexus", true);
for (Region region : acg.getRegions(locus)) System.out.println(new MarginalTree(acg, region));
// (Should be enough info here for precise agreement)
for (Region region : acg.getRegions(locus)) {
Alignment margAlign = createMarginalAlignment(alignment, region);
ClusterTree upgmaTree = new ClusterTree();
upgmaTree.initByName("clusterType", "upgma", "taxa", margAlign);
MarginalTree marginalTree = new MarginalTree(acg, region);
// outfMarg.println(marginalTree.getRoot() + ";");
// outfMarg.flush();
//
// outfUPGMA.println(upgmaTree.getRoot() + ";");
// outfUPGMA.flush();
assertTrue(topologiesEquivalent(marginalTree.getRoot(), upgmaTree.getRoot()));
}
// outfMarg.close();
// outfUPGMA.close();
}
use of beast.util.ClusterTree in project bacter by tgvaughan.
the class SerializationDeserializationTest method testExtendedNewick.
@Test
public void testExtendedNewick() throws Exception {
Alignment alignment = getAlignment();
alignment.setID("alignment");
Locus locus = new Locus("locus", alignment.getSiteCount());
// ConversionGraph
ConversionGraph acg = new ConversionGraph();
ClusterTree tree = new ClusterTree();
tree.initByName("clusterType", "upgma", "taxa", alignment);
acg.assignFrom(tree);
acg.initByName("locus", locus);
// Add recombination event 1
Node node1 = acg.getExternalNodes().get(0);
Node node2 = node1.getParent();
double height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
double height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
int startLocus = 100;
int endLocus = 200;
Conversion conv1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
acg.addConversion(conv1);
// Add recombination event 2
node1 = acg.getExternalNodes().get(8);
node2 = acg.getRoot();
height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
height2 = node2.getHeight() + 1.0;
startLocus = 300;
endLocus = 400;
Conversion conv2 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
acg.addConversion(conv2);
String newickString = acg.getExtendedNewick();
System.out.println(newickString);
ConversionGraph acgNew = new ConversionGraph();
acgNew.initByName("locus", locus);
acgNew.fromExtendedNewick(newickString);
// Check that new ACG matches old
Conversion newConv1 = acgNew.getConversions(locus).get(0);
assertEquals(newConv1.getNode1().getNr(), conv1.getNode1().getNr());
assertEquals(newConv1.getNode2().getNr(), conv1.getNode2().getNr());
assertEquals(newConv1.getHeight1(), conv1.getHeight1(), 1e-15);
assertEquals(newConv1.getHeight2(), conv1.getHeight2(), 1e-15);
assertEquals(newConv1.getStartSite(), conv1.getStartSite());
assertEquals(newConv1.getEndSite(), conv1.getEndSite());
Conversion newConv2 = acgNew.getConversions(locus).get(1);
assertEquals(newConv2.getNode1().getNr(), conv2.getNode1().getNr());
assertEquals(newConv2.getNode2().getNr(), conv2.getNode2().getNr());
assertEquals(newConv2.getHeight1(), conv2.getHeight1(), 1e-15);
assertEquals(newConv2.getHeight2(), conv2.getHeight2(), 1e-15);
assertEquals(newConv2.getStartSite(), conv2.getStartSite());
assertEquals(newConv2.getEndSite(), conv2.getEndSite());
}
use of beast.util.ClusterTree in project bacter by tgvaughan.
the class ACGLikelihoodApproxTest method testPairwiseDistances.
@Test
public void testPairwiseDistances() throws Exception {
List<Sequence> sequences = new ArrayList<>();
// 01234567890123456789
sequences.add(new Sequence("t1", "AAAAAAAAAAAAAAAAAAAA"));
sequences.add(new Sequence("t2", "AAAACAAAAAAGAAAAAAAA"));
sequences.add(new Sequence("t3", "CCCCCCCCCCCCCCCCCCCC"));
Alignment alignment = new Alignment(sequences, "nucleotide");
Locus locus = new Locus("locus", alignment);
ConversionGraph acg = new ConversionGraph();
ClusterTree tree = new ClusterTree();
tree.initByName("clusterType", "upgma", "taxa", locus.getAlignment());
acg.assignFrom(tree);
acg.initByName("locus", locus);
ACGLikelihoodApprox likelihoodApprox = new ACGLikelihoodApprox();
likelihoodApprox.initByName("acg", acg, "substitutionRate", "1.0", "alignment", alignment, "locus", locus);
Assert.assertEquals(0, likelihoodApprox.getPairwiseDistance(0, 1, 0, 4));
Assert.assertEquals(1, likelihoodApprox.getPairwiseDistance(0, 1, 0, 5));
Assert.assertEquals(1, likelihoodApprox.getPairwiseDistance(0, 1, 0, 11));
Assert.assertEquals(2, likelihoodApprox.getPairwiseDistance(0, 1, 0, 12));
Assert.assertEquals(2, likelihoodApprox.getPairwiseDistance(0, 1, 0, 20));
Assert.assertEquals(5, likelihoodApprox.getPairwiseDistance(1, 2, 2, 8));
Assert.assertEquals(19, likelihoodApprox.getPairwiseDistance(1, 2, 0, 20));
Assert.assertEquals(6, likelihoodApprox.getPairwiseDistance(0, 2, 2, 8));
Assert.assertEquals(20, likelihoodApprox.getPairwiseDistance(0, 2, 0, 20));
}
use of beast.util.ClusterTree in project bacter by tgvaughan.
the class ACGLikelihoodTest method testBeagleLikelihood.
@Test
public void testBeagleLikelihood() throws Exception {
Locus locus = new Locus("locus", getAlignment());
// ConversionGraph
ConversionGraph acg = new ConversionGraph();
ClusterTree tree = new ClusterTree();
tree.initByName("clusterType", "upgma", "taxa", locus.getAlignment());
acg.assignFrom(tree);
acg.initByName("locus", locus);
// Site model:
JukesCantor jc = new JukesCantor();
jc.initByName();
SiteModel siteModel = new SiteModel();
siteModel.initByName("substModel", jc);
// Likelihood
ACGLikelihoodBeagle argLikelihood = new ACGLikelihoodBeagle();
argLikelihood.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
acg.setEverythingDirty(true);
try {
double logP = argLikelihood.calculateLogP();
double logPtrue = argLikelihoodSlow.calculateLogP();
double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
assertTrue(relativeDiff < 1e-14);
// Add a single recombination event
Node node1 = acg.getExternalNodes().get(0);
Node node2 = node1.getParent();
double height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
double height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
int startLocus = 100;
int endLocus = 200;
Conversion recomb1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
acg.addConversion(recomb1);
logP = argLikelihood.calculateLogP();
logPtrue = argLikelihoodSlow.calculateLogP();
relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
assertTrue(relativeDiff < 1e-14);
// Add another recombination event
node1 = acg.getExternalNodes().get(0);
node2 = acg.getNode(20);
height1 = 0.75 * (node1.getHeight() + node1.getParent().getHeight());
height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
startLocus = 250;
endLocus = 300;
Conversion recomb2 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
acg.addConversion(recomb2);
logP = argLikelihood.calculateLogP();
logPtrue = argLikelihoodSlow.calculateLogP();
relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
assertTrue(relativeDiff < 1e-14);
} catch (RuntimeException ex) {
System.err.println("Beagle library not found: skipping beagle likelihood test.");
}
}
use of beast.util.ClusterTree in project bacter by tgvaughan.
the class ACGLikelihoodTest method testLikelihoodFixedData.
@Test
public void testLikelihoodFixedData() throws Exception {
Locus locus = new Locus("locus", getAlignment());
// ConversionGraph
ConversionGraph acg = new ConversionGraph();
ClusterTree tree = new ClusterTree();
tree.initByName("clusterType", "upgma", "taxa", locus.getAlignment());
acg.assignFrom(tree);
acg.initByName("locus", locus);
// Site model:
JukesCantor jc = new JukesCantor();
jc.initByName();
SiteModel siteModel = new SiteModel();
siteModel.initByName("substModel", jc);
// Likelihood
ACGLikelihood argLikelihood = new ACGLikelihood();
argLikelihood.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
ACGLikelihoodSlow argLikelihoodSlow = new ACGLikelihoodSlow();
argLikelihoodSlow.initByName("locus", locus, "tree", acg, "siteModel", siteModel);
acg.setEverythingDirty(true);
double logP = argLikelihood.calculateLogP();
double logPtrue = argLikelihoodSlow.calculateLogP();
double relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
assertTrue(relativeDiff < 1e-14);
// Add a single recombination event
Node node1 = acg.getExternalNodes().get(0);
Node node2 = node1.getParent();
double height1 = 0.5 * (node1.getHeight() + node1.getParent().getHeight());
double height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
int startLocus = 100;
int endLocus = 200;
Conversion recomb1 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
acg.addConversion(recomb1);
logP = argLikelihood.calculateLogP();
logPtrue = argLikelihoodSlow.calculateLogP();
relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
assertTrue(relativeDiff < 1e-14);
// Add another recombination event
node1 = acg.getExternalNodes().get(0);
node2 = acg.getNode(20);
height1 = 0.75 * (node1.getHeight() + node1.getParent().getHeight());
height2 = 0.5 * (node2.getHeight() + node2.getParent().getHeight());
startLocus = 250;
endLocus = 300;
Conversion recomb2 = new Conversion(node1, height1, node2, height2, startLocus, endLocus, acg, locus);
acg.addConversion(recomb2);
logP = argLikelihood.calculateLogP();
logPtrue = argLikelihoodSlow.calculateLogP();
relativeDiff = Math.abs(2.0 * (logPtrue - logP) / (logPtrue + logP));
assertTrue(relativeDiff < 1e-14);
}
Aggregations