Examples with Conversion bacter.Conversion used on opensource projects

Example 26 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ConversionGraphStatsLogger method getMeanStartSite.

/**
 * Obtain average position of conversion startSites.
 *
 * @param acg
 * @param locus
 * @return average site index or NaN if no conversions in ACG.
 */
public static double getMeanStartSite(ConversionGraph acg, Locus locus) {
    if (acg.getConvCount(locus) < 1)
        return Double.NaN;
    double mean = 0.0;
    for (Conversion conv : acg.getConversions(locus)) mean += conv.getStartSite();
    mean /= acg.getConvCount(locus);
    return mean;
}

Example 27 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class ConvertedRegionLogger method log.

@Override
public void log(long nSample, PrintStream out) {
    for (Locus locus : acgInput.get().getLoci()) {
        if (acgInput.get().getConvCount(locus) == 0) {
            out.print("NA\t");
            return;
        }
        for (int r = 0; r < acgInput.get().getConvCount(locus); r++) {
            if (r > 0)
                out.print(",");
            Conversion recomb = acgInput.get().getConversions(locus).get(r);
            out.print(recomb.getStartSite() + ":" + recomb.getEndSite());
        }
        out.print("\t");
    }
}

Example 28 with Conversion

use of bacter.Conversion 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.");
    }
}

Example 29 with Conversion

use of bacter.Conversion 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);
}

Example 30 with Conversion

use of bacter.Conversion in project bacter by tgvaughan.

the class AddRemoveConversionTest method testProbability.

/**
 * Tests whether probability of proposing a conversion lines up with
 * conversion probability found in ACGCoalescent.
 * @throws java.lang.Exception
 */
@Test
public void testProbability() throws Exception {
    ConstantPopulation popFunc = new ConstantPopulation();
    popFunc.initByName("popSize", new RealParameter("1.0"));
    Locus locus = new Locus("locus", 10000);
    TaxonSet taxonSet = getTaxonSet(10);
    SimulatedACG acg = new SimulatedACG();
    acg.initByName("rho", 1.0 / locus.getSiteCount(), "delta", 50.0, "locus", locus, "taxonset", taxonSet, "populationModel", popFunc);
    RealParameter rho = new RealParameter(Double.toString(1.0 / locus.getSiteCount()));
    RealParameter delta = new RealParameter("50.0");
    AddRemoveConversion operator = new AddRemoveConversion();
    operator.initByName("weight", 1.0, "acg", acg, "delta", delta, "populationModel", popFunc);
    ACGCoalescent coal = new ACGCoalescent();
    coal.initByName("tree", acg, "populationModel", popFunc, "rho", rho, "delta", delta);
    double logP1 = 0.0;
    double logP2 = 0.0;
    for (Conversion conv : acg.getConversions(locus)) {
        logP1 += operator.getConversionProb(conv);
        logP2 += coal.calculateConversionLogP(conv);
    }
    System.out.println("logP1 = " + logP1);
    System.out.println("logP2 = " + logP2);
    assertTrue(Math.abs(logP1 - logP2) < 1e-10);
}