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

Example 11 with ProbeSet

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

the class WigglePipeline method startPipeline.

protected void startPipeline() {
    // We first need to generate probes over all of the features listed in
    // the feature types.  The probes should cover the whole area of the
    // feature regardless of where it splices.
    boolean logTransform = optionsPanel.logTransform();
    boolean correctTotal = optionsPanel.correctTotal();
    int probeSize = optionsPanel.probeSize();
    int stepSize = optionsPanel.stepSize();
    super.progressUpdated("Making probes", 0, 1);
    Probe[] probes = null;
    String region = optionsPanel.getRegion();
    try {
        if (region.equals("Whole Genome")) {
            probes = makeGenomeProbes(probeSize, stepSize);
        } else if (region.equals("Current Chromosome")) {
            probes = makeChromosomeProbes(DisplayPreferences.getInstance().getCurrentChromosome(), probeSize, stepSize);
        } else if (region.equals("Currently Visible Region")) {
            probes = makeVisibleProbes(DisplayPreferences.getInstance().getCurrentChromosome(), SequenceRead.start(DisplayPreferences.getInstance().getCurrentLocation()), SequenceRead.end(DisplayPreferences.getInstance().getCurrentLocation()), probeSize, stepSize);
        }
    } catch (SeqMonkException sme) {
        progressExceptionReceived(sme);
        return;
    }
    collection().setProbeSet(new ProbeSet("Wiggle probes", probes));
    // method from inside the pipeline rather than doing this again ourselves.
    if (cancel) {
        progressCancelled();
        return;
    }
    bpq = new BasePairQuantitation(application);
    bpq.addProgressListener(this);
    bpq.quantitate(collection(), data, QuantitationStrandType.getTypeOptions()[0], correctTotal, false, false, false, logTransform, false);
}

Example 12 with ProbeSet

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

the class ProbeListProbeGenerator method run.

/* (non-Javadoc)
	 * @see java.lang.Runnable#run()
	 */
public void run() {
    Vector<Probe> newProbes = new Vector<Probe>();
    boolean reverse = reverseDirectionBox.isSelected();
    String description = collection.probeSet().description() + " then took subset of probes in " + selectedList.name();
    if (reverse) {
        description = description + " and reversed all probe directions";
    }
    Probe[] probes = selectedList.getAllProbes();
    for (int p = 0; p < probes.length; p++) {
        // See if we need to quit
        if (cancel) {
            generationCancelled();
            return;
        }
        if (p % 10000 == 0) {
            // Time for an update
            updateGenerationProgress("Processed " + p + " probes", p, probes.length);
        }
        if (reverse) {
            // We reverse the existing strand.  We don't do anything to
            // probes with unknown strand.
            int strand = probes[p].strand();
            if (strand == Location.FORWARD) {
                strand = Location.REVERSE;
            } else if (strand == Location.REVERSE) {
                strand = Location.FORWARD;
            }
            newProbes.add(new Probe(probes[p].chromosome(), probes[p].start(), probes[p].end(), strand, probes[p].name()));
        } else {
            newProbes.add(new Probe(probes[p].chromosome(), probes[p].packedPosition(), probes[p].name()));
        }
    }
    Probe[] finalList = newProbes.toArray(new Probe[0]);
    ProbeSet finalSet = new ProbeSet(description, finalList);
    generationComplete(finalSet);
}

Example 13 with ProbeSet

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

the class ReadPositionProbeGenerator method run.

/* (non-Javadoc)
	 * @see java.lang.Runnable#run()
	 */
public void run() {
    Chromosome[] chromosomes = collection.genome().getAllChromosomes();
    if (limitWithinRegion) {
        chromosomes = new Chromosome[] { DisplayPreferences.getInstance().getCurrentChromosome() };
    }
    Vector<Probe> newProbes = new Vector<Probe>();
    for (int c = 0; c < chromosomes.length; c++) {
        // Time for an update
        updateGenerationProgress("Processed " + c + " chromosomes", c, chromosomes.length);
        // We'll merge together the reads for all of the selected DataStores and
        // compute a single set of probes which covers all of them.
        Probe[] regions = new Probe[0];
        if (limitWithinRegion) {
            if (limitRegionBox.getSelectedItem().toString().equals("Currently Visible Region")) {
                if (chromosomes[c] == DisplayPreferences.getInstance().getCurrentChromosome()) {
                    regions = new Probe[] { new Probe(chromosomes[c], DisplayPreferences.getInstance().getCurrentLocation()) };
                }
            } else if (limitRegionBox.getSelectedItem().toString().equals("Active Probe List")) {
                regions = collection.probeSet().getActiveList().getProbesForChromosome(chromosomes[c]);
            } else {
                throw new IllegalStateException("Don't know how to filter by " + limitRegionBox.getSelectedItem().toString());
            }
        } else {
            regions = new Probe[] { new Probe(chromosomes[c], 0, chromosomes[c].length()) };
        }
        for (int p = 0; p < regions.length; p++) {
            long[][] v = new long[selectedStores.length][];
            for (int s = 0; s < selectedStores.length; s++) {
                v[s] = selectedStores[s].getReadsForProbe(regions[p]);
            }
            long[] rawReads = getUsableRedundantReads(LongSetSorter.sortLongSets(v));
            v = null;
            int currentCount = 1;
            int currentStart = 0;
            int currentEnd = 0;
            int currentPositionCount = 0;
            for (int r = 1; r < rawReads.length; r++) {
                if (SequenceRead.start(rawReads[r]) == SequenceRead.start(rawReads[r - 1]) && SequenceRead.end(rawReads[r]) == SequenceRead.end(rawReads[r - 1]) && (ignoreStrand || SequenceRead.strand(rawReads[r]) == SequenceRead.strand(rawReads[r - 1]))) {
                    // It's the same
                    ++currentCount;
                } else {
                    // Check if we need to make a new probe
                    if (currentCount >= minCount) {
                        // Add this probe to the current set
                        if (currentPositionCount == 0) {
                            // Start a new position
                            currentStart = SequenceRead.start(rawReads[r - 1]);
                            currentEnd = SequenceRead.end(rawReads[r - 1]);
                        } else {
                            if (SequenceRead.end(rawReads[r - 1]) > currentEnd) {
                                currentEnd = SequenceRead.end(rawReads[r - 1]);
                            }
                        }
                        currentPositionCount++;
                        if (currentPositionCount == readsPerWindow) {
                            int strand = Probe.UNKNOWN;
                            if (!ignoreStrand) {
                                strand = SequenceRead.strand(rawReads[r - 1]);
                            }
                            newProbes.add(new Probe(chromosomes[c], currentStart, currentEnd, strand));
                            currentPositionCount = 0;
                        }
                    }
                    currentCount = 1;
                }
            }
            // See if we need to add the last read
            if (currentCount >= minCount && rawReads.length >= 1) {
                // Add this probe to the current set
                if (currentPositionCount == 0) {
                    // Start a new position
                    currentStart = SequenceRead.start(rawReads[rawReads.length - 1]);
                    currentEnd = SequenceRead.end(rawReads[rawReads.length - 1]);
                } else {
                    if (SequenceRead.end(rawReads[rawReads.length - 1]) > currentEnd) {
                        currentEnd = SequenceRead.end(rawReads[rawReads.length - 1]);
                    }
                }
                currentPositionCount++;
                // Make a probe with whatever we have left
                int strand = Probe.UNKNOWN;
                if (!ignoreStrand) {
                    strand = SequenceRead.strand(rawReads[rawReads.length - 1]);
                }
                newProbes.add(new Probe(chromosomes[c], currentStart, currentEnd, strand));
            }
        }
    }
    Probe[] finalList = newProbes.toArray(new Probe[0]);
    newProbes.clear();
    ProbeSet finalSet = new ProbeSet(getDescription(), finalList);
    generationComplete(finalSet);
}

Example 14 with ProbeSet

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

the class CodonBiasPipeline method startPipeline.

protected void startPipeline() {
    // We first need to generate probes over all of the features listed in
    // the feature types.  The probes should cover the whole area of the
    // feature regardless of where it splices.
    Vector<Probe> probes = new Vector<Probe>();
    double pValue = optionsPanel.pValue();
    String libraryType = optionsPanel.libraryType();
    Chromosome[] chrs = collection().genome().getAllChromosomes();
    for (int c = 0; c < chrs.length; c++) {
        if (cancel) {
            progressCancelled();
            return;
        }
        progressUpdated("Making probes for chr" + chrs[c].name(), c, chrs.length * 2);
        Feature[] features = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedFeatureType());
        for (int f = 0; f < features.length; f++) {
            if (cancel) {
                progressCancelled();
                return;
            }
            Probe p = new Probe(chrs[c], features[f].location().start(), features[f].location().end(), features[f].location().strand(), features[f].name());
            probes.add(p);
        }
    }
    allProbes = probes.toArray(new Probe[0]);
    collection().setProbeSet(new ProbeSet("Features over " + optionsPanel.getSelectedFeatureType(), allProbes));
    // Now we can quantitate each individual feature and test for whether it is significantly
    // showing codon bias
    ArrayList<Vector<ProbeTTestValue>> significantProbes = new ArrayList<Vector<ProbeTTestValue>>();
    // data contains the data stores that this pipeline is going to use. We need to test each data store.
    for (int d = 0; d < data.length; d++) {
        significantProbes.add(new Vector<ProbeTTestValue>());
    }
    int probeCounter = 0;
    for (int c = 0; c < chrs.length; c++) {
        if (cancel) {
            progressCancelled();
            return;
        }
        progressUpdated("Quantitating features on chr" + chrs[c].name(), chrs.length + c, chrs.length * 2);
        Feature[] features = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedFeatureType());
        for (int p = 0; p < features.length; p++) {
            // Get the corresponding feature and work out the mapping between genomic position and codon sub position.
            int[] mappingArray = createGenomeMappingArray(features[p]);
            DATASTORE_LOOP: for (int d = 0; d < data.length; d++) {
                if (cancel) {
                    progressCancelled();
                    return;
                }
                long[] reads = data[d].getReadsForProbe(allProbes[probeCounter]);
                // TODO: make this configurable
                if (reads.length < 5) {
                    data[d].setValueForProbe(allProbes[probeCounter], Float.NaN);
                    continue DATASTORE_LOOP;
                }
                int pos1Count = 0;
                int pos2Count = 0;
                int pos3Count = 0;
                READ_LOOP: for (int r = 0; r < reads.length; r++) {
                    int genomicReadStart = SequenceRead.start(reads[r]);
                    int genomicReadEnd = SequenceRead.end(reads[r]);
                    int readStrand = SequenceRead.strand(reads[r]);
                    int relativeReadStart = -1;
                    // forward reads
                    if (readStrand == 1) {
                        if (libraryType == "Same strand specific") {
                            if (features[p].location().strand() == Location.FORWARD) {
                                // The start of the read needs to be within the feature
                                if (genomicReadStart - features[p].location().start() < 0) {
                                    continue READ_LOOP;
                                } else {
                                    // look up the read start pos in the mapping array
                                    relativeReadStart = mappingArray[genomicReadStart - features[p].location().start()];
                                }
                            }
                        } else if (libraryType == "Opposing strand specific") {
                            if (features[p].location().strand() == Location.REVERSE) {
                                // The "start" of a reverse read/probe is actually the end
                                if (features[p].location().end() - genomicReadEnd < 0) {
                                    continue READ_LOOP;
                                } else {
                                    relativeReadStart = mappingArray[features[p].location().end() - genomicReadEnd];
                                }
                            }
                        }
                    }
                    // reverse reads
                    if (readStrand == -1) {
                        if (libraryType == "Same strand specific") {
                            if (features[p].location().strand() == Location.REVERSE) {
                                if (features[p].location().end() - genomicReadEnd < 0) {
                                    continue READ_LOOP;
                                } else {
                                    // System.out.println("features[p].location().end() is " + features[p].location().end() + ", genomicReadEnd is " + genomicReadEnd);
                                    // System.out.println("mapping array[0] is " + mappingArray[0]);
                                    relativeReadStart = mappingArray[features[p].location().end() - genomicReadEnd];
                                }
                            }
                        } else if (libraryType == "Opposing strand specific") {
                            if (features[p].location().strand() == Location.FORWARD) {
                                // The start of the read needs to be within the feature
                                if (genomicReadStart - features[p].location().start() < 0) {
                                    continue READ_LOOP;
                                } else {
                                    // look up the read start position in the mapping array
                                    relativeReadStart = mappingArray[genomicReadStart - features[p].location().start()];
                                }
                            }
                        }
                    }
                    // find out which position the read is in
                    if (relativeReadStart == -1) {
                        continue READ_LOOP;
                    } else if (relativeReadStart % 3 == 0) {
                        pos3Count++;
                        continue READ_LOOP;
                    } else if ((relativeReadStart + 1) % 3 == 0) {
                        pos2Count++;
                        continue READ_LOOP;
                    } else if ((relativeReadStart + 2) % 3 == 0) {
                        pos1Count++;
                    }
                }
                // closing bracket for read loop
                // System.out.println("pos1Count for "+ features[p].name() + " is " + pos1Count);
                // System.out.println("pos2Count for "+ features[p].name() + " is " + pos2Count);
                // System.out.println("pos3Count for "+ features[p].name() + " is " + pos3Count);
                int interestingCodonCount = 0;
                int otherCodonCount = 0;
                if (optionsPanel.codonSubPosition() == 1) {
                    interestingCodonCount = pos1Count;
                    otherCodonCount = pos2Count + pos3Count;
                } else if (optionsPanel.codonSubPosition() == 2) {
                    interestingCodonCount = pos2Count;
                    otherCodonCount = pos1Count + pos3Count;
                } else if (optionsPanel.codonSubPosition() == 3) {
                    interestingCodonCount = pos3Count;
                    otherCodonCount = pos1Count + pos2Count;
                }
                int totalCount = interestingCodonCount + otherCodonCount;
                // BinomialDistribution bd = new BinomialDistribution(interestingCodonCount+otherCodonCount, 1/3d);
                BinomialDistribution bd = new BinomialDistribution(totalCount, 1 / 3d);
                // Since the binomial distribution gives the probability of getting a value higher than
                // this we need to subtract one so we get the probability of this or higher.
                double thisPValue = 1 - bd.cumulativeProbability(interestingCodonCount - 1);
                if (interestingCodonCount == 0)
                    thisPValue = 1;
                // We have to add all results at this stage so we don't mess up the multiple
                // testing correction later on.
                significantProbes.get(d).add(new ProbeTTestValue(allProbes[probeCounter], thisPValue));
                float percentageCount;
                if (totalCount == 0) {
                    percentageCount = 0;
                } else {
                    percentageCount = ((float) interestingCodonCount / (float) totalCount) * 100;
                }
                data[d].setValueForProbe(allProbes[probeCounter], percentageCount);
            // System.out.println("totalCount = " + totalCount);
            // System.out.println("interestingCodonCount " + interestingCodonCount);
            // System.out.println("pValue = " + thisPValue);
            // System.out.println("percentageCount = " + percentageCount);
            // System.out.println("");
            }
            probeCounter++;
        }
    }
    // filtering those which pass our p-value cutoff
    for (int d = 0; d < data.length; d++) {
        ProbeTTestValue[] ttestResults = significantProbes.get(d).toArray(new ProbeTTestValue[0]);
        BenjHochFDR.calculateQValues(ttestResults);
        ProbeList newList = new ProbeList(collection().probeSet(), "Codon bias < " + pValue + " in " + data[d].name(), "Probes showing significant codon bias for position " + optionsPanel.codonSubPosition() + " with a cutoff of " + pValue, "FDR");
        for (int i = 0; i < ttestResults.length; i++) {
            if (ttestResults[i].q < pValue) {
                newList.addProbe(ttestResults[i].probe, (float) ttestResults[i].q);
            }
        }
    }
    StringBuffer quantitationDescription = new StringBuffer();
    quantitationDescription.append("Codon bias pipeline using codon position " + optionsPanel.codonSubPosition() + " for " + optionsPanel.libraryType() + " library.");
    collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
    quantitatonComplete();
}

Example 15 with ProbeSet

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

the class TranscriptionTerminationPipeline method startPipeline.

protected void startPipeline() {
    // We first need to generate probes over all of the features listed in
    // the feature types.  The probes should cover the whole area of the
    // feature regardless of where it splices.
    Vector<Probe> probes = new Vector<Probe>();
    int minCount = optionsPanel.minCount();
    int measurementLength = optionsPanel.measurementLength();
    QuantitationStrandType readFilter = optionsPanel.readFilter();
    Chromosome[] chrs = collection().genome().getAllChromosomes();
    for (int c = 0; c < chrs.length; c++) {
        if (cancel) {
            progressCancelled();
            return;
        }
        progressUpdated("Creating Probes" + chrs[c].name(), c, chrs.length * 2);
        Feature[] features = getValidFeatures(chrs[c], measurementLength);
        for (int f = 0; f < features.length; f++) {
            if (cancel) {
                progressCancelled();
                return;
            }
            if (features[f].location().strand() == Location.REVERSE) {
                Probe p = new Probe(chrs[c], features[f].location().start() - measurementLength, features[f].location().start() + measurementLength, features[f].location().strand(), features[f].name());
                probes.add(p);
            } else {
                Probe p = new Probe(chrs[c], features[f].location().end() - measurementLength, features[f].location().end() + measurementLength, features[f].location().strand(), features[f].name());
                probes.add(p);
            }
        }
    }
    Probe[] allProbes = probes.toArray(new Probe[0]);
    collection().setProbeSet(new ProbeSet("Features " + measurementLength + "bp around the end of " + optionsPanel.getSelectedFeatureType(), allProbes));
    int probeIndex = 0;
    for (int c = 0; c < chrs.length; c++) {
        if (cancel) {
            progressCancelled();
            return;
        }
        progressUpdated("Quantitating features on chr" + chrs[c].name(), chrs.length + c, chrs.length * 2);
        Feature[] features = getValidFeatures(chrs[c], measurementLength);
        for (int f = 0; f < features.length; f++) {
            if (cancel) {
                progressCancelled();
                return;
            }
            for (int d = 0; d < data.length; d++) {
                if (allProbes[probeIndex].strand() == Location.REVERSE) {
                    Probe downstreamProbe = new Probe(chrs[c], features[f].location().start() - measurementLength, features[f].location().start(), features[f].location().strand(), features[f].name());
                    Probe upstreamProbe = new Probe(chrs[c], features[f].location().start(), features[f].location().start() + measurementLength, features[f].location().strand(), features[f].name());
                    long[] upstreamReads = data[d].getReadsForProbe(upstreamProbe);
                    long[] downstreamReads = data[d].getReadsForProbe(downstreamProbe);
                    int upstreamCount = 0;
                    for (int i = 0; i < upstreamReads.length; i++) {
                        if (readFilter.useRead(allProbes[probeIndex], upstreamReads[i]))
                            ++upstreamCount;
                    }
                    int downstreamCount = 0;
                    for (int i = 0; i < downstreamReads.length; i++) {
                        if (readFilter.useRead(allProbes[probeIndex], downstreamReads[i]))
                            ++downstreamCount;
                    }
                    float percentage = ((upstreamCount - downstreamCount) / (float) upstreamCount) * 100f;
                    if (upstreamCount >= minCount) {
                        data[d].setValueForProbe(allProbes[probeIndex], percentage);
                    } else {
                        data[d].setValueForProbe(allProbes[probeIndex], Float.NaN);
                    }
                } else {
                    Probe upstreamProbe = new Probe(chrs[c], features[f].location().end() - measurementLength, features[f].location().end(), features[f].location().strand(), features[f].name());
                    Probe downstreamProbe = new Probe(chrs[c], features[f].location().end(), features[f].location().end() + measurementLength, features[f].location().strand(), features[f].name());
                    long[] upstreamReads = data[d].getReadsForProbe(upstreamProbe);
                    long[] downstreamReads = data[d].getReadsForProbe(downstreamProbe);
                    int upstreamCount = 0;
                    for (int i = 0; i < upstreamReads.length; i++) {
                        if (readFilter.useRead(allProbes[probeIndex], upstreamReads[i]))
                            ++upstreamCount;
                    }
                    int downstreamCount = 0;
                    for (int i = 0; i < downstreamReads.length; i++) {
                        if (readFilter.useRead(allProbes[probeIndex], downstreamReads[i]))
                            ++downstreamCount;
                    }
                    float percentage = ((upstreamCount - downstreamCount) / (float) upstreamCount) * 100f;
                    if (upstreamCount >= minCount) {
                        data[d].setValueForProbe(allProbes[probeIndex], percentage);
                    } else {
                        data[d].setValueForProbe(allProbes[probeIndex], Float.NaN);
                    }
                }
            }
            ++probeIndex;
        }
    }
    StringBuffer quantitationDescription = new StringBuffer();
    quantitationDescription.append("Transcription termination pipeline quantitation ");
    quantitationDescription.append(". Directionality was ");
    quantitationDescription.append(optionsPanel.libraryTypeBox.getSelectedItem());
    quantitationDescription.append(". Measurement length was ");
    quantitationDescription.append(optionsPanel.measurementLength());
    quantitationDescription.append(". Min count was ");
    quantitationDescription.append(optionsPanel.minCount());
    collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
    quantitatonComplete();
}