Examples with Locus bacter.Locus used on opensource projects

Example 26 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ACGAnnotator method summarizeConversions.

/**
 * Add summarized conversions to given ACG.
 *
 * @param cladeSystem information summarizing ACG posterior
 * @param acg conversion graph
 * @param threshold significance threshold
 * @param summaryStrategy strategy used when summarizing event ages/heights
 */
protected void summarizeConversions(ACGCladeSystem cladeSystem, ConversionGraph acg, int nACGs, double threshold, SummaryStrategy summaryStrategy) {
    BitSet[] bitSets = cladeSystem.getBitSets(acg);
    for (int fromNr = 0; fromNr < acg.getNodeCount(); fromNr++) {
        BitSet from = bitSets[fromNr];
        for (int toNr = 0; toNr < acg.getNodeCount(); toNr++) {
            BitSet to = bitSets[toNr];
            for (Locus locus : acg.getLoci()) {
                List<ACGCladeSystem.ConversionSummary> conversionSummaries = cladeSystem.getConversionSummaries(from, to, locus, nACGs, threshold);
                for (ACGCladeSystem.ConversionSummary conversionSummary : conversionSummaries) {
                    Conversion conv = new Conversion();
                    conv.setLocus(locus);
                    conv.setNode1(acg.getNode(fromNr));
                    conv.setNode2(acg.getNode(toNr));
                    double posteriorSupport = conversionSummary.nIncludedACGs / (double) nACGs;
                    double[] height1s = new double[conversionSummary.summarizedConvCount()];
                    double[] height2s = new double[conversionSummary.summarizedConvCount()];
                    double[] startSites = new double[conversionSummary.summarizedConvCount()];
                    double[] endSites = new double[conversionSummary.summarizedConvCount()];
                    for (int i = 0; i < conversionSummary.summarizedConvCount(); i++) {
                        height1s[i] = conversionSummary.height1s.get(i);
                        height2s[i] = conversionSummary.height2s.get(i);
                        startSites[i] = conversionSummary.startSites.get(i);
                        endSites[i] = conversionSummary.ends.get(i);
                    }
                    if (summaryStrategy == SummaryStrategy.MEAN) {
                        conv.setHeight1(DiscreteStatistics.mean(height1s));
                        conv.setHeight2(DiscreteStatistics.mean(height2s));
                        conv.setStartSite((int) Math.round(DiscreteStatistics.mean(startSites)));
                        conv.setEndSite((int) Math.round(DiscreteStatistics.mean(endSites)));
                    } else {
                        conv.setHeight1(DiscreteStatistics.median(height1s));
                        conv.setHeight2(DiscreteStatistics.median(height2s));
                        conv.setStartSite((int) Math.round(DiscreteStatistics.median(startSites)));
                        conv.setEndSite((int) Math.round(DiscreteStatistics.median(endSites)));
                    }
                    Arrays.sort(height1s);
                    double minHeight1HPD = height1s[(int) (0.025 * height1s.length)];
                    double maxHeight1HPD = height1s[(int) (0.975 * height1s.length)];
                    Arrays.sort(height2s);
                    double minHeight2HPD = height2s[(int) (0.025 * height2s.length)];
                    double maxHeight2HPD = height2s[(int) (0.975 * height2s.length)];
                    Arrays.sort(startSites);
                    int minStartHPD = (int) startSites[(int) (0.025 * startSites.length)];
                    int maxStartHPD = (int) startSites[(int) (0.975 * startSites.length)];
                    Arrays.sort(endSites);
                    int minEndHPD = (int) endSites[(int) (0.025 * endSites.length)];
                    int maxEndHPD = (int) endSites[(int) (0.975 * endSites.length)];
                    conv.newickMetaDataBottom = "height_95%_HPD={" + minHeight1HPD + "," + maxHeight1HPD + "}";
                    conv.newickMetaDataMiddle = "posterior=" + String.format("%2.2f", posteriorSupport) + ", startSite_95%_HPD={" + minStartHPD + "," + maxStartHPD + "}" + ", endSite_95%_HPD={" + minEndHPD + "," + maxEndHPD + "}";
                    conv.newickMetaDataTop = "height_95%_HPD={" + minHeight2HPD + "," + maxHeight2HPD + "}";
                    acg.addConversion(conv);
                }
            }
        }
    }
}

Example 27 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class ACGCladeSystem method collectConversions.

/**
 * Add conversions described on provided acg to the internal list
 * for later summary.
 *
 * @param acg conversion graph from which to extract conversions
 */
public void collectConversions(ConversionGraph acg) {
    getBitSets(acg);
    Map<BitSet, Map<BitSet, Long>> geneFlowTemp = new HashMap<>();
    // Assemble list of conversions for each pair of clades on each locus
    for (Locus locus : acg.getLoci()) {
        conversionListsTemp.clear();
        for (Conversion conv : acg.getConversions(locus)) {
            conv.acgIndex = acgIndex;
            BitSetPair bsPair = new BitSetPair(conv);
            if (!conversionListsTemp.containsKey(bsPair))
                conversionListsTemp.put(bsPair, new ArrayList<>());
            conversionListsTemp.get(bsPair).add(conv);
            // Record gene flow
            if (!geneFlowTemp.containsKey(bsPair.from))
                geneFlowTemp.put(bsPair.from, new HashMap<>());
            long oldFlow = 0;
            if (geneFlowTemp.get(bsPair.from).containsKey(bsPair.to))
                oldFlow = geneFlowTemp.get(bsPair.from).get(bsPair.to);
            geneFlowTemp.get(bsPair.from).put(bsPair.to, oldFlow + conv.getSiteCount());
        }
        // Merge overlapping conversions:
        for (BitSetPair bsPair : conversionListsTemp.keySet()) {
            List<Conversion> merged = mergeOverlappingConvs(conversionListsTemp.get(bsPair));
            if (!conversionLists.containsKey(bsPair))
                conversionLists.put(bsPair, new HashMap<>());
            if (!conversionLists.get(bsPair).containsKey(locus))
                conversionLists.get(bsPair).put(locus, new ArrayList<>());
            conversionLists.get(bsPair).get(locus).addAll(merged);
        }
    }
    geneFlow.add(geneFlowTemp);
    acgIndex += 1;
}

Example 28 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class CFConvSwapExperiment method main.

public static void main(String[] args) throws Exception {
    // Load model
    XMLParser parser = new XMLParser();
    MCMC mcmc = (MCMC) parser.parseFile(new File("inferencePreSimulatedData.xml"));
    State state = mcmc.startStateInput.get();
    state.setStateFileName("problem.state");
    state.restoreFromFile();
    double oldPosterior = state.robustlyCalcPosterior(mcmc.posteriorInput.get());
    ConversionGraph acg = (ConversionGraph) state.getStateNode(0);
    PrintStream ps = new PrintStream("proposal.trees");
    ps.println(acg);
    Node srcNode = acg.getNode(3);
    Node srcNodeS = getSibling(srcNode);
    Node destNode = acg.getNode(6);
    double t_srcNodeP = srcNode.getParent().getHeight();
    disconnectEdge(acg, srcNode);
    Locus locus = acg.getLoci().get(0);
    Conversion convToReplace = acg.getConversions(locus).get(27);
    acg.deleteConversion(convToReplace);
    Node srcNodeP = srcNode.getParent();
    connectEdge(acg, srcNode, destNode, convToReplace.getHeight2());
    // Move Conversions
    for (Conversion conv : acg.getConversions(locus)) {
        boolean moved = false;
        if (conv.getNode1() == srcNode && conv.getHeight1() > srcNodeP.getHeight()) {
            conv.setNode1(destNode);
            moved = true;
        }
        if (conv.getNode2() == srcNode && conv.getHeight2() > srcNodeP.getHeight()) {
            conv.setNode2(destNode);
            moved = true;
        }
        if (moved) {
            while (conv.getHeight1() > conv.getNode1().getParent().getHeight()) conv.setNode1(conv.getNode1().getParent());
            while (conv.getHeight2() > conv.getNode2().getParent().getHeight()) conv.setNode2(conv.getNode2().getParent());
        }
    }
    Conversion convNew = new Conversion();
    convNew.setLocus(locus);
    convNew.setStartSite(0);
    convNew.setEndSite(4000);
    convNew.setNode1(srcNode);
    convNew.setNode2(srcNodeS);
    convNew.setHeight1(convToReplace.getHeight1());
    convNew.setHeight2(t_srcNodeP);
    acg.addConversion(convNew);
    ps.println(acg);
    double newPosterior = state.robustlyCalcPosterior(mcmc.posteriorInput.get());
    System.out.println(newPosterior - oldPosterior);
// state.setStateFileName("proposal.state");
// state.storeToFile(1);
// Open state file
}

Example 29 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class CFConvSwapExperiment method disconnectEdge.

public static void disconnectEdge(ConversionGraph acg, Node node) {
    if (node.isRoot())
        throw new IllegalArgumentException("Programmer error: " + "root argument passed to disconnectEdge().");
    Node parent = node.getParent();
    Node sister = getSibling(node);
    if (parent.isRoot()) {
        parent.removeChild(sister);
        sister.setParent(null);
    } else {
        Node grandParent = parent.getParent();
        grandParent.removeChild(parent);
        parent.setParent(null);
        parent.removeChild(sister);
        grandParent.addChild(sister);
    }
    for (Locus locus : acg.getLoci()) {
        for (Conversion conv : acg.getConversions(locus)) {
            if (conv.getNode1() == parent)
                conv.setNode1(sister);
            if (conv.getNode2() == parent)
                conv.setNode2(sister);
        }
    }
}

Example 30 with Locus

use of bacter.Locus in project bacter by tgvaughan.

the class BacterACGLogReader method extractLoci.

/**
 * Retrieve list of loci from preamble or postamble.
 */
private void extractLoci() {
    List<String> prepost = new ArrayList<>();
    prepost.addAll(preamble);
    prepost.addAll(postamble);
    for (String line : prepost) {
        line = line.trim();
        if (line.startsWith("loci ") && line.endsWith(";")) {
            for (String locusEntry : line.substring(5, line.length() - 1).split(" ")) {
                String[] locusPair = locusEntry.split(":");
                loci.add(new Locus(locusPair[0], Integer.parseInt(locusPair[1])));
            }
        }
    }
}