Examples with ProbeList uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList used on opensource projects

Example 16 with ProbeList

use of uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList in project SeqMonk by s-andrews.

the class StrandBiasPlotPanel method getFilteredProbes.

/**
 * Gets the filtered probes.
 *
 * @param probeset the probeset
 * @return the filtered probes
 */
public ProbeList getFilteredProbes(ProbeSet probeset) {
    double minDiff = Math.min(diffStart, diffEnd);
    double maxDiff = Math.max(diffStart, diffEnd);
    ProbeList list = new ProbeList(probeList, "Strand bias between " + df.format(minDiff) + " and " + df.format(maxDiff), "Strand bias difference in " + store.name() + " was between " + df.format(minDiff) + " and " + df.format(maxDiff), "Strand Bias diff");
    if (madeSelection) {
        for (int p = 0; p < allValues.length; p++) {
            if (allValues[p].proportion < minDiff)
                continue;
            if (allValues[p].proportion > maxDiff)
                continue;
            list.addProbe(allValues[p].probe, allValues[p].proportion);
        }
    }
    return list;
}

Example 17 with ProbeList

use of uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList in project SeqMonk by s-andrews.

the class IntensityDifferenceFilter method generateProbeList.

/* (non-Javadoc)
	 * @see uk.ac.babraham.SeqMonk.Filters.ProbeFilter#generateProbeList()
	 */
protected void generateProbeList() {
    applyMultipleTestingCorrection = optionsPanel.multipleTestingBox.isSelected();
    Probe[] probes = startingList.getAllProbes();
    // We'll pull the number of probes to sample from the preferences if they've changed it
    Integer updatedProbesPerSet = optionsPanel.probesPerSet();
    if (updatedProbesPerSet != null)
        probesPerSet = updatedProbesPerSet;
    ProbeList newList = new ProbeList(startingList, "Filtered Probes", "", "Diff p-value");
    // We'll build up a set of p-values as we go along
    float[] lowestPValues = new float[probes.length];
    for (int p = 0; p < lowestPValues.length; p++) {
        lowestPValues[p] = 1;
    }
    // This is going to be the temporary array we populate with the set of
    // differences we are going to analyse.
    double[] currentDiffSet = new double[probesPerSet];
    // First work out the set of comparisons we're going to make
    Vector<SingleComparison> comparisons = new Vector<IntensityDifferenceFilter.SingleComparison>();
    for (int fromIndex = 0; fromIndex < fromStores.length; fromIndex++) {
        for (int toIndex = 0; toIndex < toStores.length; toIndex++) {
            if (fromStores[fromIndex] == toStores[toIndex])
                continue;
            // If we can find the fromStore in the toStores we've already done and the
            // toStore anywhere in the fromStores then we can skip this.
            boolean canSkip = false;
            for (int i = 0; i < fromIndex; i++) {
                if (fromStores[i] == toStores[toIndex]) {
                    for (int j = 0; j < toStores.length; j++) {
                        if (toStores[j] == fromStores[fromIndex]) {
                            canSkip = true;
                            break;
                        }
                    }
                    break;
                }
            }
            if (canSkip)
                continue;
            comparisons.add(new SingleComparison(fromIndex, toIndex));
        }
    }
    // Put something in the progress whilst we're ordering the probe values to make
    // the comparison.
    progressUpdated("Generating background model", 0, 1);
    for (int comparisonIndex = 0; comparisonIndex < comparisons.size(); comparisonIndex++) {
        int fromIndex = comparisons.elementAt(comparisonIndex).fromIndex;
        int toIndex = comparisons.elementAt(comparisonIndex).toIndex;
        // We need to generate a set of probe indices ordered by their average intensity
        Integer[] indices = new Integer[probes.length];
        for (int i = 0; i < probes.length; i++) {
            indices[i] = i;
        }
        Comparator<Integer> comp = new AverageIntensityComparator(fromStores[fromIndex], toStores[toIndex], probes);
        Arrays.sort(indices, comp);
        progressUpdated("Made " + comparisonIndex + " out of " + comparisons.size() + " comparisons", comparisonIndex, comparisons.size());
        IndexTTestValue[] currentPValues = new IndexTTestValue[indices.length];
        for (int i = 0; i < indices.length; i++) {
            if (cancel) {
                cancel = false;
                progressCancelled();
                return;
            }
            if (i % 1000 == 0) {
                int progress = (i * 100) / indices.length;
                progress += 100 * comparisonIndex;
                progressUpdated("Made " + comparisonIndex + " out of " + comparisons.size() + " comparisons", progress, comparisons.size() * 100);
            }
            // We need to make up the set of differences to represent this probe
            int startingIndex = i - (probesPerSet / 2);
            if (startingIndex < 0)
                startingIndex = 0;
            if (startingIndex + (probesPerSet + 1) >= probes.length)
                startingIndex = probes.length - (probesPerSet + 1);
            try {
                for (int j = startingIndex; j < startingIndex + (probesPerSet + 1); j++) {
                    if (// Don't include the point being tested in the background model
                    j == startingIndex)
                        // Don't include the point being tested in the background model
                        continue;
                    else if (j < startingIndex) {
                        currentDiffSet[j - startingIndex] = fromStores[fromIndex].getValueForProbe(probes[indices[j]]) - toStores[toIndex].getValueForProbe(probes[indices[j]]);
                    } else {
                        currentDiffSet[(j - startingIndex) - 1] = fromStores[fromIndex].getValueForProbe(probes[indices[j]]) - toStores[toIndex].getValueForProbe(probes[indices[j]]);
                    }
                }
                // Should we fix the mean at 0?
                double mean = 0;
                // double mean = SimpleStats.mean(currentDiffSet);
                double stdev = SimpleStats.stdev(currentDiffSet, mean);
                if (stdev == 0) {
                    currentPValues[indices[i]] = new IndexTTestValue(indices[i], 1);
                    continue;
                }
                // Get the difference for this point
                double diff = fromStores[fromIndex].getValueForProbe(probes[indices[i]]) - toStores[toIndex].getValueForProbe(probes[indices[i]]);
                NormalDistribution nd = new NormalDistribution(mean, stdev);
                double significance;
                if (diff < mean) {
                    significance = nd.cumulativeProbability(diff);
                } else {
                    significance = 1 - nd.cumulativeProbability(diff);
                }
                currentPValues[indices[i]] = new IndexTTestValue(indices[i], significance);
            } catch (SeqMonkException sme) {
                progressExceptionReceived(sme);
                return;
            }
        }
        // We now need to correct the set of pValues
        if (applyMultipleTestingCorrection) {
            BenjHochFDR.calculateQValues(currentPValues);
        }
        // the combined set
        if (applyMultipleTestingCorrection) {
            for (int i = 0; i < currentPValues.length; i++) {
                if (currentPValues[i].q < lowestPValues[currentPValues[i].index]) {
                    lowestPValues[currentPValues[i].index] = (float) currentPValues[i].q;
                }
            }
        } else {
            for (int i = 0; i < currentPValues.length; i++) {
                if (currentPValues[i].p < lowestPValues[currentPValues[i].index]) {
                    lowestPValues[currentPValues[i].index] = (float) currentPValues[i].p;
                }
            }
        }
    }
    // pass the filter.
    for (int i = 0; i < lowestPValues.length; i++) {
        if (lowestPValues[i] < pValueLimit) {
            newList.addProbe(probes[i], lowestPValues[i]);
        }
    }
    filterFinished(newList);
}

Example 18 with ProbeList

use of uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList in project SeqMonk by s-andrews.

the class LimmaFilter method generateProbeList.

/* (non-Javadoc)
	 * @see uk.ac.babraham.SeqMonk.Filters.ProbeFilter#generateProbeList()
	 */
@Override
protected void generateProbeList() {
    // We need to make a temporary directory, save the data into it, write out the R script
    // and then run it an collect the list of results, then clean up.
    // Make up the list of DataStores in each replicate set
    DataStore[] fromStores = replicateSets[0].dataStores();
    DataStore[] toStores = replicateSets[1].dataStores();
    File tempDir;
    try {
        progressUpdated("Creating temp directory", 0, 1);
        tempDir = TempDirectory.createTempDirectory();
        // System.err.println("Temp dir is "+tempDir.getAbsolutePath());
        progressUpdated("Writing R script", 0, 1);
        // Get the template script
        Template template = new Template(ClassLoader.getSystemResource("uk/ac/babraham/SeqMonk/Filters/LimmaFilter/limma_template.r"));
        // Substitute in the variables we need to change
        template.setValue("WORKING", tempDir.getAbsolutePath().replace("\\", "/"));
        // Say which p value column we're filtering on
        if (multiTest) {
            template.setValue("CORRECTED", "adj.P.Val");
        } else {
            template.setValue("CORRECTED", "P.Value");
        }
        StringBuffer sb = new StringBuffer();
        for (int i = 0; i < fromStores.length; i++) {
            if (i > 0)
                sb.append(",");
            sb.append("\"from\"");
        }
        for (int i = 0; i < toStores.length; i++) {
            sb.append(",");
            sb.append("\"to\"");
        }
        template.setValue("CONDITIONS", sb.toString());
        template.setValue("PVALUE", "" + cutoff);
        // Write the script file
        File scriptFile = new File(tempDir.getAbsoluteFile() + "/script.r");
        PrintWriter pr = new PrintWriter(scriptFile);
        pr.print(template.toString());
        pr.close();
        // Write the count data
        File countFile = new File(tempDir.getAbsoluteFile() + "/counts.txt");
        pr = new PrintWriter(countFile);
        sb = new StringBuffer();
        sb.append("probe");
        for (int i = 0; i < fromStores.length; i++) {
            sb.append("\t");
            sb.append("from");
            sb.append(i);
        }
        for (int i = 0; i < toStores.length; i++) {
            sb.append("\t");
            sb.append("to");
            sb.append(i);
        }
        pr.println(sb.toString());
        progressUpdated("Writing count data", 0, 1);
        Probe[] probes = startingList.getAllProbes();
        float value;
        for (int p = 0; p < probes.length; p++) {
            if (p % 1000 == 0) {
                progressUpdated("Writing count data", p, probes.length);
            }
            sb = new StringBuffer();
            sb.append(p);
            for (int i = 0; i < fromStores.length; i++) {
                sb.append("\t");
                value = fromStores[i].getValueForProbe(probes[p]);
                sb.append(value);
            }
            for (int i = 0; i < toStores.length; i++) {
                sb.append("\t");
                value = toStores[i].getValueForProbe(probes[p]);
                sb.append(value);
            }
            pr.println(sb.toString());
        }
        pr.close();
        progressUpdated("Running R Script", 0, 1);
        RScriptRunner runner = new RScriptRunner(tempDir);
        RProgressListener listener = new RProgressListener(runner);
        runner.addProgressListener(new ProgressRecordDialog("R Session", runner));
        runner.runScript();
        while (true) {
            if (listener.cancelled()) {
                progressCancelled();
                return;
            }
            if (listener.exceptionReceived()) {
                progressExceptionReceived(listener.exception());
                return;
            }
            if (listener.complete())
                break;
            Thread.sleep(500);
        }
        // We can now parse the results and put the hits into a new probe list
        ProbeList newList;
        newList = new ProbeList(startingList, "", "", "FDR");
        File hitsFile = new File(tempDir.getAbsolutePath() + "/hits.txt");
        BufferedReader br = new BufferedReader(new FileReader(hitsFile));
        String line = br.readLine();
        while ((line = br.readLine()) != null) {
            String[] sections = line.split("\t");
            int probeIndex = Integer.parseInt(sections[0]);
            float pValue;
            if (multiTest) {
                pValue = Float.parseFloat(sections[sections.length - 2]);
            } else {
                pValue = Float.parseFloat(sections[sections.length - 3]);
            }
            newList.addProbe(probes[probeIndex], pValue);
        }
        br.close();
        runner.cleanUp();
        // TODO: Show log of R session?
        filterFinished(newList);
    } catch (Exception ioe) {
        progressExceptionReceived(ioe);
        return;
    }
// filterFinished(newList);
}

Example 19 with ProbeList

use of uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList in project SeqMonk by s-andrews.

the class LogisticRegressionFilter method generateProbeList.

/* (non-Javadoc)
	 * @see uk.ac.babraham.SeqMonk.Filters.ProbeFilter#generateProbeList()
	 */
@Override
protected void generateProbeList() {
    // We need to make a temporary directory, save the data into it, write out the R script
    // and then run it an collect the list of results, then clean up.
    // Make up the list of DataStores in each replicate set
    DataStore[] fromStores = replicateSets[0].dataStores();
    DataStore[] toStores = replicateSets[1].dataStores();
    if (minObservations < 3)
        minObservations = 3;
    File tempDir;
    try {
        Probe[] probes = startingList.getAllProbes();
        if (resample) {
            // We need to check that the data stores are quantitated
            for (int i = 0; i < fromStores.length; i++) {
                if (!fromStores[i].isQuantitated()) {
                    progressExceptionReceived(new SeqMonkException("Data Store " + fromStores[i].name() + " wasn't quantitated"));
                    return;
                }
                for (int p = 0; p < probes.length; p++) {
                    float value = fromStores[i].getValueForProbe(probes[p]);
                    if ((!Float.isNaN(value)) && (value < 0 || value > 100)) {
                        progressExceptionReceived(new SeqMonkException("Data Store " + fromStores[i].name() + " had a value outside the range 0-100 (" + value + ")"));
                        return;
                    }
                }
            }
            for (int i = 0; i < toStores.length; i++) {
                if (!toStores[i].isQuantitated()) {
                    progressExceptionReceived(new SeqMonkException("Data Store " + toStores[i].name() + " wasn't quantitated"));
                    return;
                }
                for (int p = 0; p < probes.length; p++) {
                    float value = toStores[i].getValueForProbe(probes[p]);
                    if ((!Float.isNaN(value)) && (value < 0 || value > 100)) {
                        progressExceptionReceived(new SeqMonkException("Data Store " + toStores[i].name() + " had a value outside the range 0-100 (" + value + ")"));
                        return;
                    }
                }
            }
        }
        progressUpdated("Creating temp directory", 0, 1);
        tempDir = TempDirectory.createTempDirectory();
        System.err.println("Temp dir is " + tempDir.getAbsolutePath());
        progressUpdated("Writing R script", 0, 1);
        // Get the template script
        Template template = new Template(ClassLoader.getSystemResource("uk/ac/babraham/SeqMonk/Filters/LogisticRegressionFilter/logistic_regression_template.r"));
        // Write the count data
        // Sort these so we can get probes from the same chromosome together
        Arrays.sort(probes);
        PrintWriter pr = null;
        String lastChr = "";
        // Rather than not testing probes based on their absolute difference
        // we should just post-filter them.  The easiest way to do this will
        // be to not test (as we do now) but explicity pass in the number of
        // tests we should have performed to the multiple testing correction.
        int numberOfTestsToCorrectBy = 0;
        PROBE: for (int p = 0; p < probes.length; p++) {
            if (!probes[p].chromosome().name().equals(lastChr)) {
                if (pr != null)
                    pr.close();
                File outFile = new File(tempDir.getAbsoluteFile() + "/data_chr" + probes[p].chromosome().name() + ".txt");
                pr = new PrintWriter(outFile);
                lastChr = probes[p].chromosome().name();
                pr.println("id\tgroup\treplicate\tstate\tcount");
            }
            if (p % 1000 == 0) {
                progressUpdated("Writing data for chr" + lastChr, p, probes.length);
            }
            int[] fromMethCounts = new int[fromStores.length];
            int[] fromUnmethCounts = new int[fromStores.length];
            int[] toMethCounts = new int[toStores.length];
            int[] toUnmethCounts = new int[toStores.length];
            for (int i = 0; i < fromStores.length; i++) {
                long[] reads = fromStores[i].getReadsForProbe(probes[p]);
                int totalCount = 0;
                int methCount = 0;
                if (resample) {
                    float value = fromStores[i].getValueForProbe(probes[p]);
                    if (Float.isNaN(value)) {
                        continue PROBE;
                    }
                    totalCount = reads.length;
                    methCount = Math.round((totalCount * value) / 100f);
                } else {
                    for (int r = 0; r < reads.length; r++) {
                        totalCount++;
                        if (SequenceRead.strand(reads[r]) == Location.FORWARD) {
                            ++methCount;
                        }
                    }
                }
                fromMethCounts[i] = methCount;
                fromUnmethCounts[i] = totalCount - methCount;
            }
            for (int i = 0; i < toStores.length; i++) {
                long[] reads = toStores[i].getReadsForProbe(probes[p]);
                int totalCount = 0;
                int methCount = 0;
                if (resample) {
                    float value = toStores[i].getValueForProbe(probes[p]);
                    if (Float.isNaN(value)) {
                        continue PROBE;
                    }
                    totalCount = reads.length;
                    methCount = Math.round((totalCount * value) / 100f);
                } else {
                    for (int r = 0; r < reads.length; r++) {
                        totalCount++;
                        if (SequenceRead.strand(reads[r]) == Location.FORWARD) {
                            ++methCount;
                        }
                    }
                }
                toMethCounts[i] = methCount;
                toUnmethCounts[i] = totalCount - methCount;
            }
            // Check to see we meet the requirements for the min amount of information
            // and the min diff.
            int totalFromMeth = 0;
            int totalFrom = 0;
            int totalToMeth = 0;
            int totalTo = 0;
            int validFrom = 0;
            for (int i = 0; i < fromStores.length; i++) {
                totalFromMeth += fromMethCounts[i];
                totalFrom += fromMethCounts[i];
                totalFrom += fromUnmethCounts[i];
                if (fromMethCounts[i] + fromUnmethCounts[i] >= minObservations) {
                    ++validFrom;
                }
            }
            int validTo = 0;
            for (int i = 0; i < toStores.length; i++) {
                totalToMeth += toMethCounts[i];
                totalTo += toMethCounts[i];
                totalTo += toUnmethCounts[i];
                if (toMethCounts[i] + toUnmethCounts[i] >= minObservations) {
                    ++validTo;
                }
            }
            // have enough data in all stores to go ahead and do the test.
            if (validFrom < fromStores.length || validTo < toStores.length) {
                // We don't have enough data to measure this one
                continue;
            }
            // At this point we have to count this probe as valid for the
            // purposes of multiple testing correction
            ++numberOfTestsToCorrectBy;
            float[] fromPercentages = new float[validFrom];
            float[] toPercentages = new float[validTo];
            int lastFromIndex = 0;
            int lastToIndex = 0;
            for (int i = 0; i < fromMethCounts.length; i++) {
                if (fromMethCounts[i] + fromUnmethCounts[i] == 0)
                    continue;
                fromPercentages[lastFromIndex] = fromMethCounts[i] * 100f / (fromMethCounts[i] + fromUnmethCounts[i]);
                ++lastFromIndex;
            }
            for (int i = 0; i < toMethCounts.length; i++) {
                if (toMethCounts[i] + toUnmethCounts[i] == 0)
                    continue;
                toPercentages[lastToIndex] = toMethCounts[i] * 100f / (toMethCounts[i] + toUnmethCounts[i]);
                ++lastToIndex;
            }
            for (int i = 0; i < fromMethCounts.length; i++) {
                pr.println(p + "\tfrom\tfrom" + i + "\tmeth\t" + fromMethCounts[i]);
                pr.println(p + "\tfrom\tfrom" + i + "\tunmeth\t" + fromUnmethCounts[i]);
            }
            for (int i = 0; i < toMethCounts.length; i++) {
                pr.println(p + "\tto\tto" + i + "\tmeth\t" + toMethCounts[i]);
                pr.println(p + "\tto\tto" + i + "\tunmeth\t" + toUnmethCounts[i]);
            }
        }
        pr.close();
        // Sanity check to make sure we have something to work with.
        if (numberOfTestsToCorrectBy == 0) {
            progressExceptionReceived(new IllegalStateException("No probes had enough data to test."));
        }
        // Now we can complete the template
        // Substitute in the variables we need to change
        template.setValue("WORKING", tempDir.getAbsolutePath().replace("\\", "/"));
        template.setValue("CORRECTCOUNT", "" + numberOfTestsToCorrectBy);
        template.setValue("PVALUE", "" + pValueCutoff);
        if (multiTest) {
            template.setValue("MULTITEST", "TRUE");
        } else {
            template.setValue("MULTITEST", "FALSE");
        }
        // Write the script file
        File scriptFile = new File(tempDir.getAbsoluteFile() + "/script.r");
        pr = new PrintWriter(scriptFile);
        pr.print(template.toString());
        pr.close();
        progressUpdated("Running R Script", 0, 1);
        RScriptRunner runner = new RScriptRunner(tempDir);
        RProgressListener listener = new RProgressListener(runner);
        runner.addProgressListener(new ProgressRecordDialog("R Session", runner));
        runner.runScript();
        while (true) {
            if (listener.cancelled()) {
                progressCancelled();
                pr.close();
                return;
            }
            if (listener.exceptionReceived()) {
                progressExceptionReceived(new SeqMonkException("R Script failed"));
                pr.close();
                return;
            }
            if (listener.complete())
                break;
            Thread.sleep(500);
        }
        // We can now parse the results and put the hits into a new probe list
        ProbeList newList;
        if (multiTest) {
            newList = new ProbeList(startingList, "", "", "FDR");
        } else {
            newList = new ProbeList(startingList, "", "", "p-value");
        }
        File hitsFile = new File(tempDir.getAbsolutePath() + "/hits.txt");
        BufferedReader br = new BufferedReader(new FileReader(hitsFile));
        String line = br.readLine();
        while ((line = br.readLine()) != null) {
            String[] sections = line.split("\t");
            String[] indexSections = sections[0].split("\\.");
            int probeIndex = Integer.parseInt(indexSections[indexSections.length - 1]);
            float pValue = Float.parseFloat(sections[sections.length - 1]);
            newList.addProbe(probes[probeIndex], pValue);
        }
        br.close();
        runner.cleanUp();
        filterFinished(newList);
    } catch (Exception ioe) {
        progressExceptionReceived(ioe);
        return;
    }
}

Example 20 with ProbeList

use of uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList in project SeqMonk by s-andrews.

the class ManualCorrelationFilter method generateProbeList.

/* (non-Javadoc)
	 * @see uk.ac.babraham.SeqMonk.Filters.ProbeFilter#generateProbeList()
	 */
@Override
protected void generateProbeList() {
    // System.out.println("Data store size="+stores.length+" lower="+lowerLimit+" upper="+upperLimit+" type="+limitType+" chosen="+chosenNumber);
    profiles = new double[optionsPanel.graphs.length][];
    for (int g = 0; g < optionsPanel.graphs.length; g++) {
        profiles[g] = optionsPanel.graphs[g].profile();
    }
    Probe[] probes = startingList.getAllProbes();
    ProbeList parentList = null;
    ProbeList[] newLists = new ProbeList[profiles.length];
    if (newLists.length > 1) {
        parentList = new ProbeList(startingList, "Filtered Probes", "", null);
        for (int p = 0; p < newLists.length; p++) {
            newLists[p] = new ProbeList(parentList, listName() + " " + (p + 1), listDescription(), "R-value");
        }
    } else {
        // We won't have a parent list
        newLists[0] = new ProbeList(startingList, listName(), listDescription(), "R-value");
    }
    double[] currentProfile = new double[stores.length];
    for (int p = 0; p < probes.length; p++) {
        progressUpdated(p, probes.length);
        if (cancel) {
            cancel = false;
            progressCancelled();
            return;
        }
        try {
            for (int s = 0; s < stores.length; s++) {
                currentProfile[s] = stores[s].getValueForProbe(probes[p]);
            }
            float bestR = 0;
            int bestRIndex = -1;
            for (int pr = 0; pr < profiles.length; pr++) {
                float r = PearsonCorrelation.calculateCorrelation(profiles[pr], currentProfile);
                if (correlationCutoff > 0) {
                    if (r >= correlationCutoff && r > bestR) {
                        bestR = r;
                        bestRIndex = pr;
                    }
                } else {
                    if (r <= correlationCutoff && r < bestR) {
                        bestR = r;
                        bestRIndex = pr;
                    }
                }
            }
            if (bestRIndex >= 0) {
                newLists[bestRIndex].addProbe(probes[p], bestR);
                if (parentList != null) {
                    parentList.addProbe(probes[p], null);
                }
            }
        } catch (SeqMonkException ex) {
            continue;
        }
    }
    // If we only have one cluster then we just return that list
    if (newLists.length == 1) {
        filterFinished(newLists[0]);
    } else {
        filterFinished(parentList);
    }
}