Examples with HiCInteractionStrengthCalculator uk.ac.babraham.SeqMonk.DataTypes.Interaction.HiCInteractionStrengthCalculator used on opensource projects

Example 1 with HiCInteractionStrengthCalculator

use of uk.ac.babraham.SeqMonk.DataTypes.Interaction.HiCInteractionStrengthCalculator in project SeqMonk by s-andrews.

the class FourCEnrichmentQuantitation method run.

/* (non-Javadoc)
	 * @see java.lang.Runnable#run()
	 */
public void run() {
    // Start off by finding the right HiC data for each 4C dataset
    HiCDataStore[] parentHiC = new HiCDataStore[data.length];
    Probe[] parentProbes = new Probe[data.length];
    ProbeList[] significantLists = new ProbeList[data.length];
    for (int d = 0; d < data.length; d++) {
        String filename = data[d].fileName();
        filename = filename.replaceAll("HiC other end of ", "");
        filename = filename.replaceAll(" for region.*", "");
        System.err.println("Looking for HiC match to " + filename);
        for (int h = 0; h < hiCData.length; h++) {
            if (((DataStore) hiCData[h]).name().equals(filename)) {
                parentHiC[d] = hiCData[h];
            }
        }
        if (parentHiC[d] == null) {
            progressWarningReceived(new SeqMonkException("Failed to find HiC dataset '" + filename + "' for 4C dataset " + data[d].name()));
            continue;
        }
        significantLists[d] = new ProbeList(application.dataCollection().probeSet(), "4C p<0.05 " + data[d].name(), "4C pipeline significance < 0.05", "P-value");
        // Also make up a probe to represent the region from which
        // the dataset was made
        filename = data[d].fileName();
        filename = filename.replaceAll("^.*for region ", "");
        String[] locationSections = filename.split("[-:]");
        if (locationSections.length != 3) {
            progressWarningReceived(new SeqMonkException("Couldn't extract location from " + filename));
            continue;
        }
        try {
            parentProbes[d] = new Probe(application.dataCollection().genome().getChromosome(locationSections[0]).chromosome(), Integer.parseInt(locationSections[1]), Integer.parseInt(locationSections[2]));
        } catch (Exception e) {
            progressExceptionReceived(e);
            return;
        }
    }
    // Make strength calculators for each HiC set
    HiCInteractionStrengthCalculator[] strengthCalculators = new HiCInteractionStrengthCalculator[parentHiC.length];
    for (int i = 0; i < parentHiC.length; i++) {
        strengthCalculators[i] = new HiCInteractionStrengthCalculator(parentHiC[i], true);
    }
    // Get the cis/trans counts for the parent probes
    int[] parentProbeCisCounts = new int[parentHiC.length];
    int[] parentProbeTransCounts = new int[parentHiC.length];
    for (int p = 0; p < parentHiC.length; p++) {
        if (parentHiC[p] == null)
            continue;
        HiCHitCollection hits = parentHiC[p].getHiCReadsForProbe(parentProbes[p]);
        String[] chrNames = hits.getChromosomeNamesWithHits();
        for (int c = 0; c < chrNames.length; c++) {
            if (chrNames[c].equals(hits.getSourceChromosomeName())) {
                parentProbeCisCounts[p] = hits.getSourcePositionsForChromosome(chrNames[c]).length;
            } else {
                parentProbeTransCounts[p] += hits.getSourcePositionsForChromosome(chrNames[c]).length;
            }
        }
    }
    Probe[] probes = application.dataCollection().probeSet().getAllProbes();
    for (int p = 0; p < probes.length; p++) {
        // See if we need to quit
        if (cancel) {
            progressCancelled();
            return;
        }
        progressUpdated(p, probes.length);
        for (int d = 0; d < data.length; d++) {
            if (parentHiC[d] == null)
                continue;
            int thisProbeCisCount = 0;
            int thisProbeTransCount = 0;
            HiCHitCollection hiCHits = parentHiC[d].getHiCReadsForProbe(probes[p]);
            String[] chrNames = hiCHits.getChromosomeNamesWithHits();
            for (int c = 0; c < chrNames.length; c++) {
                if (chrNames[c].equals(hiCHits.getSourceChromosomeName())) {
                    thisProbeCisCount = hiCHits.getSourcePositionsForChromosome(chrNames[c]).length;
                } else {
                    thisProbeTransCount += hiCHits.getSourcePositionsForChromosome(chrNames[c]).length;
                }
            }
            strengthCalculators[d].calculateInteraction(data[d].getReadsForProbe(probes[p]).length, thisProbeCisCount, thisProbeTransCount, parentProbeCisCounts[d], parentProbeTransCounts[d], probes[p], parentProbes[d]);
            float obsExp = (float) strengthCalculators[d].obsExp();
            data[d].setValueForProbe(probes[p], obsExp);
            float pValue = (float) strengthCalculators[d].rawPValue() * probes.length;
            if (pValue < 0.05) {
                significantLists[d].addProbe(probes[p], pValue);
            }
        }
    }
    quantitatonComplete();
}