Examples with Chromosome uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome used on opensource projects

Example 61 with Chromosome

use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome in project SeqMonk by s-andrews.

the class BisulphiteFeaturePipeline 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.getMinCount();
    int minFeatureCountValue = optionsPanel.getMinFeatureCount();
    String combineString = optionsPanel.getCombineType();
    int combineType = 0;
    if (combineString.equals("Mean")) {
        combineType = MEAN;
    } else if (combineString.equals("Median")) {
        combineType = MEDIAN;
    } else if (combineString.equals("Max")) {
        combineType = MAX;
    } else if (combineString.equals("Min")) {
        combineType = MIN;
    } else {
        throw new IllegalStateException("Unknown combine type '" + combineString + "' found in methylation pipeline options");
    }
    Chromosome[] chrs = collection().genome().getAllChromosomes();
    if (optionsPanel.getSelectedFeatureType() != "[Existing Probes]") {
        for (int c = 0; c < chrs.length; c++) {
            if (cancel) {
                progressCancelled();
                return;
            }
            progressUpdated("Making features for chr" + chrs[c].name(), c, chrs.length * 2);
            Feature[] features = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedFeatureType());
            Arrays.sort(features);
            for (int f = 0; f < features.length; f++) {
                if (cancel) {
                    progressCancelled();
                    return;
                }
                probes.add(new Probe(chrs[c], features[f].location().start(), features[f].location().end(), features[f].location().strand(), features[f].name()));
            }
        }
        Probe[] allProbes = probes.toArray(new Probe[0]);
        collection().setProbeSet(new ProbeSet("Features over " + optionsPanel.getSelectedFeatureType(), allProbes));
    }
    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);
        Probe[] thisChrProbes = collection().probeSet().getProbesForChromosome(chrs[c]);
        for (int p = 0; p < thisChrProbes.length; p++) {
            // This data structure is going to hold the counts.  The int array will hold
            // pairs of values for the different data stores (meth and unmeth) so the length
            // of this array will be data store length * 2.
            Hashtable<Integer, int[]> counts = new Hashtable<Integer, int[]>();
            for (int d = 0; d < data.length; d++) {
                if (cancel) {
                    progressCancelled();
                    return;
                }
                long[] reads = data[d].getReadsForProbe(thisChrProbes[p]);
                for (int r = 0; r < reads.length; r++) {
                    for (int pos = SequenceRead.start(reads[r]); pos <= SequenceRead.end(reads[r]); pos++) {
                        if (pos < thisChrProbes[p].start())
                            continue;
                        if (pos > thisChrProbes[p].end())
                            break;
                        if (!counts.containsKey(pos)) {
                            counts.put(pos, new int[data.length * 2]);
                        }
                        if (SequenceRead.strand(reads[r]) == Location.FORWARD) {
                            counts.get(pos)[2 * d]++;
                        } else if (SequenceRead.strand(reads[r]) == Location.REVERSE) {
                            counts.get(pos)[(2 * d) + 1]++;
                        }
                    }
                }
            }
            if (optionsPanel.applyMinCountAcrossAllStores()) {
                Enumeration<Integer> en = counts.keys();
                OUTER: while (en.hasMoreElements()) {
                    Integer i = en.nextElement();
                    int[] values = counts.get(i);
                    for (int d = 0; d < data.length; d++) {
                        if (values[2 * d] + values[(2 * d) + 1] < minCount) {
                            counts.remove(i);
                            continue OUTER;
                        }
                    }
                }
            }
            for (int d = 0; d < data.length; d++) {
                Enumeration<Integer> en = counts.keys();
                Vector<Float> percentages = new Vector<Float>();
                while (en.hasMoreElements()) {
                    Integer pos = en.nextElement();
                    int total = counts.get(pos)[2 * d] + counts.get(pos)[(2 * d) + 1];
                    if (total < minCount)
                        continue;
                    float percent = 100f * (counts.get(pos)[2 * d] / (float) total);
                    percentages.add(percent);
                }
                // Now combine the methylation percentages we have in the way they
                // wanted.
                float[] validMeasures = new float[percentages.size()];
                float finalMeasure = Float.NaN;
                for (int i = 0; i < validMeasures.length; i++) {
                    validMeasures[i] = percentages.elementAt(i).floatValue();
                }
                if (validMeasures.length < minFeatureCountValue) {
                    finalMeasure = Float.NaN;
                } else if (validMeasures.length == 0) {
                    finalMeasure = Float.NaN;
                } else if (validMeasures.length == 1) {
                    finalMeasure = validMeasures[0];
                } else {
                    Arrays.sort(validMeasures);
                    switch(combineType) {
                        case MEAN:
                            finalMeasure = SimpleStats.mean(validMeasures);
                            break;
                        case MEDIAN:
                            finalMeasure = SimpleStats.median(validMeasures);
                            break;
                        case MIN:
                            finalMeasure = validMeasures[0];
                            break;
                        case MAX:
                            finalMeasure = validMeasures[validMeasures.length - 1];
                            break;
                    }
                }
                data[d].setValueForProbe(thisChrProbes[p], finalMeasure);
            }
        }
    }
    if (optionsPanel.getSelectedFeatureType() != "[Existing Probes]") {
        collection().probeSet().setDescription("Probes over " + optionsPanel.getSelectedFeatureType() + " features");
    }
    StringBuffer quantitationDescription = new StringBuffer();
    quantitationDescription.append("Methylation feature pipeline quantitation using ");
    quantitationDescription.append(optionsPanel.getSelectedFeatureType());
    quantitationDescription.append(" features with min count ");
    if (optionsPanel.applyMinCountAcrossAllStores()) {
        quantitationDescription.append(" in all stores ");
    }
    quantitationDescription.append("of ");
    quantitationDescription.append(optionsPanel.getMinCount());
    quantitationDescription.append(" per position, and at least ");
    quantitationDescription.append(optionsPanel.getMinFeatureCount());
    quantitationDescription.append(" observations per feature, then combining using the ");
    quantitationDescription.append(optionsPanel.getCombineType());
    collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
    quantitatonComplete();
}

Example 62 with Chromosome

use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome in project SeqMonk by s-andrews.

the class GeneTrapQuantitationPipeline 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>();
    boolean logTransform = optionsPanel.logTransform();
    boolean rawCounts = optionsPanel.rawCountsBox.isSelected();
    if (rawCounts) {
        logTransform = false;
    }
    Chromosome[] chrs = collection().genome().getAllChromosomes();
    for (int c = 0; c < chrs.length; c++) {
        // System.err.println("Processing chr "+chrs[c].name());
        if (cancel) {
            progressCancelled();
            return;
        }
        progressUpdated("Making features for chr" + chrs[c].name(), c, chrs.length * 2);
        Feature[] features = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedFeatureType());
        Arrays.sort(features);
        FeatureGroup[] mergedTranscripts = mergeTranscripts(features);
        for (int f = 0; f < mergedTranscripts.length; f++) {
            if (cancel) {
                progressCancelled();
                return;
            }
            probes.add(new Probe(chrs[c], mergedTranscripts[f].start(), mergedTranscripts[f].end(), mergedTranscripts[f].strand(), mergedTranscripts[f].name));
        }
    }
    Probe[] allProbes = probes.toArray(new Probe[0]);
    collection().setProbeSet(new ProbeSet("Transcript features over " + optionsPanel.getSelectedFeatureType(), allProbes));
    // Having made probes we now need to quantitate them.  We'll fetch the
    // probes overlapping each sub-feature and then aggregate these together
    // to get the final quantitation.
    int currentIndex = 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());
        Arrays.sort(features);
        FeatureGroup[] mergedTranscripts = mergeTranscripts(features);
        for (int f = 0; f < mergedTranscripts.length; f++) {
            Location[] subLocations = mergedTranscripts[f].getSubLocations();
            for (int d = 0; d < data.length; d++) {
                if (cancel) {
                    progressCancelled();
                    return;
                }
                int totalCount = 0;
                long[] reads = data[d].getReadsForProbe(new Probe(chrs[c], mergedTranscripts[f].start(), mergedTranscripts[f].end()));
                READ: for (int r = 0; r < reads.length; r++) {
                    if (SequenceRead.strand(reads[r]) == mergedTranscripts[f].strand()) {
                        // TODO: Should we check if we're within the CDS?
                        ++totalCount;
                        continue;
                    }
                    // We'll still count this if it overlaps with an exon
                    for (int s = 0; s < subLocations.length; s++) {
                        if (SequenceRead.overlaps(reads[r], subLocations[s].packedPosition())) {
                            ++totalCount;
                            continue READ;
                        }
                    }
                }
                float value = totalCount;
                // If we're log transforming then we need to set zero values to 0.9
                if (logTransform && value == 0) {
                    value = 0.9f;
                }
                // We also correct by the total read count
                if (!rawCounts) {
                    value /= (data[d].getTotalReadCount() / 1000000f);
                }
                // Finally we do the log transform if we've been asked to
                if (logTransform) {
                    value = (float) Math.log(value) / log2;
                }
                data[d].setValueForProbe(allProbes[currentIndex], value);
            }
            currentIndex++;
        }
    }
    StringBuffer quantitationDescription = new StringBuffer();
    quantitationDescription.append("Gene trap quantitation on transcripts of type ");
    quantitationDescription.append(optionsPanel.getSelectedFeatureType());
    if (logTransform) {
        quantitationDescription.append(". Log transformed");
    }
    collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
    quantitatonComplete();
}

Example 63 with Chromosome

use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome in project SeqMonk by s-andrews.

the class SplicingEfficiencyPipeline 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>();
    boolean correctReadLength = optionsPanel.correctForReadLength();
    boolean logTransform = optionsPanel.logTransform();
    boolean applyTranscriptLengthCorrection = optionsPanel.applyTranscriptLengthCorrection();
    // We need to make a lookup for chromosomes from their names for later on.
    Chromosome[] chrs = collection().genome().getAllChromosomes();
    Hashtable<String, Chromosome> chrLookup = new Hashtable<String, Chromosome>();
    for (int c = 0; c < chrs.length; c++) {
        chrLookup.put(chrs[c].name(), chrs[c]);
    }
    // We're going to cache the gene to transcript mappings so we don't have to re-do them
    // later on.
    Hashtable<String, Vector<Feature>> transcriptLookup = new Hashtable<String, Vector<Feature>>();
    Vector<Feature> usedGeneFeatures = new Vector<Feature>();
    // We'll record the number of failures so we can report it
    int noTranscriptCount = 0;
    int noIntronCount = 0;
    for (int c = 0; c < chrs.length; c++) {
        if (cancel) {
            progressCancelled();
            return;
        }
        progressUpdated("Making features for chr" + chrs[c].name(), c, chrs.length * 2);
        Feature[] geneFeatures = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedGeneFeatureType());
        Arrays.sort(geneFeatures);
        Feature[] transcriptFeatures = collection().genome().annotationCollection().getFeaturesForType(chrs[c], optionsPanel.getSelectedTranscriptFeatureType());
        Arrays.sort(transcriptFeatures);
        // We need to figure out whether we can match based on name, or whether we need to do it based on location.
        boolean matchOnNames = true;
        for (int t = 0; t < transcriptFeatures.length; t++) {
            if (!transcriptFeatures[t].name().matches("-\\d\\d\\d$")) {
                matchOnNames = false;
                break;
            }
        }
        if (!matchOnNames) {
            progressWarningReceived(new SeqMonkException("Non-Ensembl transcript names found - matching based on position"));
        }
        if (matchOnNames) {
            for (int i = 0; i < transcriptFeatures.length; i++) {
                String name = transcriptFeatures[i].name();
                name = name.replaceAll("-\\d\\d\\d$", "");
                if (!transcriptLookup.containsKey(name)) {
                    transcriptLookup.put(name, new Vector<Feature>());
                }
                transcriptLookup.get(name).add(transcriptFeatures[i]);
            }
        } else {
            // We need to go through the genes and transcripts in parallel and match them up based on containment
            // and direction.
            int lastGeneIndex = 0;
            for (int t = 0; t < transcriptFeatures.length; t++) {
                for (int g = lastGeneIndex; g < geneFeatures.length; g++) {
                    if (transcriptFeatures[t].location().start() > geneFeatures[g].location().end()) {
                        lastGeneIndex = g;
                        continue;
                    }
                    // If the gene is beyond the end of the transcript we can stop looking
                    if (geneFeatures[g].location().start() > transcriptFeatures[t].location().end()) {
                        break;
                    }
                    // If we're on the same strand and contained within the gene then we have a match
                    if (geneFeatures[g].location().strand() == transcriptFeatures[t].location().strand() && transcriptFeatures[t].location().start() >= geneFeatures[g].location().start() && transcriptFeatures[t].location().end() <= geneFeatures[g].location().end()) {
                        String name = geneFeatures[g].name();
                        if (!transcriptLookup.containsKey(name)) {
                            transcriptLookup.put(name, new Vector<Feature>());
                        }
                        transcriptLookup.get(name).add(transcriptFeatures[t]);
                    }
                }
            }
        }
        for (int f = 0; f < geneFeatures.length; f++) {
            if (cancel) {
                progressCancelled();
                return;
            }
            if (!transcriptLookup.containsKey(geneFeatures[f].name())) {
                ++noTranscriptCount;
                // No good making features for genes with no transcript.
                continue;
            }
            Feature[] localTranscripts = transcriptLookup.get(geneFeatures[f].name()).toArray(new Feature[0]);
            boolean[] exonPositions = new boolean[geneFeatures[f].location().length()];
            for (int t = 0; t < localTranscripts.length; t++) {
                if (!SequenceRead.overlaps(geneFeatures[f].location().packedPosition(), localTranscripts[t].location().packedPosition())) {
                    // It's not a match for this gene
                    continue;
                }
                if (!(localTranscripts[t].location() instanceof SplitLocation)) {
                    for (int p = localTranscripts[t].location().start() - geneFeatures[f].location().start(); p <= localTranscripts[t].location().end() - geneFeatures[f].location().start(); p++) {
                        if (p < 0)
                            continue;
                        if (p >= exonPositions.length)
                            continue;
                        exonPositions[p] = true;
                    }
                    continue;
                }
                Location[] exons = ((SplitLocation) localTranscripts[t].location()).subLocations();
                for (int e = 0; e < exons.length; e++) {
                    for (int p = exons[e].start() - geneFeatures[f].location().start(); p <= exons[e].end() - geneFeatures[f].location().start(); p++) {
                        if (p < 0)
                            continue;
                        if (p >= exonPositions.length)
                            continue;
                        exonPositions[p] = true;
                    }
                }
            }
            int exonPositionCount = 0;
            int intronPositionCount = 0;
            for (int p = 0; p < exonPositions.length; p++) {
                if (exonPositions[p])
                    exonPositionCount++;
                else
                    intronPositionCount++;
            }
            if (exonPositionCount == 0) {
                ++noTranscriptCount;
                continue;
            }
            if (intronPositionCount == 0) {
                ++noIntronCount;
                continue;
            }
            probes.add(new Probe(chrs[c], geneFeatures[f].location().start(), geneFeatures[f].location().end(), geneFeatures[f].location().strand(), geneFeatures[f].name()));
            usedGeneFeatures.add(geneFeatures[f]);
        }
    }
    if (noTranscriptCount > 0) {
        progressWarningReceived(new SeqMonkException("" + noTranscriptCount + " genes had no associated transcripts"));
    }
    if (noIntronCount > 0) {
        progressWarningReceived(new SeqMonkException("" + noIntronCount + " genes had no intron only areas"));
    }
    Probe[] allProbes = probes.toArray(new Probe[0]);
    collection().setProbeSet(new ProbeSet("Gene features over " + optionsPanel.getSelectedGeneFeatureType(), allProbes));
    // Having made probes we now need to quantitate them.  We'll get the full set of reads for the full gene
    // and then work out a combined set of exons.  We can then quantitate the bases which fall into introns
    // and exons separately and then work out a ratio from there.
    QuantitationStrandType readFilter = optionsPanel.readFilter();
    Feature[] geneFeatures = usedGeneFeatures.toArray(new Feature[0]);
    for (int g = 0; g < geneFeatures.length; g++) {
        if (cancel) {
            progressCancelled();
            return;
        }
        if (g % 100 == 0) {
            progressUpdated("Quantitating features", g, geneFeatures.length);
        }
        int[] readLengths = new int[data.length];
        if (correctReadLength) {
            for (int d = 0; d < data.length; d++) {
                readLengths[d] = data[d].getMaxReadLength();
            }
        }
        // Find the set of transcripts which relate to this gene.
        Feature[] transcripts = transcriptLookup.get(geneFeatures[g].name()).toArray(new Feature[0]);
        Vector<Feature> validTranscripts = new Vector<Feature>();
        for (int t = 0; t < transcripts.length; t++) {
            if (SequenceRead.overlaps(geneFeatures[g].location().packedPosition(), transcripts[t].location().packedPosition())) {
                validTranscripts.add(transcripts[t]);
            }
        }
        transcripts = validTranscripts.toArray(new Feature[0]);
        // before so if we do then something has gone wrong.
        if (transcripts.length == 0) {
            throw new IllegalStateException("No transcripts for gene " + geneFeatures[g] + " on chr " + geneFeatures[g].chromosomeName());
        }
        Vector<Location> allExonsVector = new Vector<Location>();
        for (int t = 0; t < transcripts.length; t++) {
            if (transcripts[t].location() instanceof SplitLocation) {
                Location[] sublocs = ((SplitLocation) transcripts[t].location()).subLocations();
                for (int i = 0; i < sublocs.length; i++) allExonsVector.add(sublocs[i]);
            } else {
                allExonsVector.add(transcripts[t].location());
            }
        }
        // if (geneFeatures[f].name().equals("Cpa6")) {
        // System.err.println("Cpa6 had "+allExonsVector.size()+" total exons");
        // }
        Collections.sort(allExonsVector);
        // Now go through and make a merged set of exons to remove any redundancy
        Vector<Location> mergedExonsVector = new Vector<Location>();
        Location lastLocation = null;
        Enumeration<Location> en = allExonsVector.elements();
        while (en.hasMoreElements()) {
            Location l = en.nextElement();
            if (lastLocation == null) {
                // if (geneFeatures[f].name().equals("Cpa6")) {
                // System.err.println("Setting as first location");
                // }
                lastLocation = l;
                continue;
            }
            // Check if it's the same as the last, which is likely
            if (l.start() == lastLocation.start() && l.end() == lastLocation.end()) {
                continue;
            }
            // Check if they overlap and can be merged
            if (l.start() <= lastLocation.end() && l.end() >= lastLocation.start()) {
                // if (geneFeatures[f].name().equals("Cpa6")) {
                // System.err.println("Overlaps with last location");
                // }
                // It overlaps with the last location so merge them
                lastLocation = new Location(Math.min(l.start(), lastLocation.start()), Math.max(l.end(), lastLocation.end()), geneFeatures[g].location().strand());
                // }
                continue;
            }
            // Start a new location
            // if (geneFeatures[f].name().equals("Cpa6")) {
            // System.err.println("Doesn't overlap - adding last location and creating new one");
            // }
            mergedExonsVector.add(lastLocation);
            lastLocation = l;
        }
        if (lastLocation != null) {
            mergedExonsVector.add(lastLocation);
        }
        // if (geneFeatures[f].name().equals("Cpa6")) {
        // System.err.println("Cpa6 had "+mergedExonsVector.size()+" merged exons");
        // }
        // Now we can start the quantitation.
        int intronLength = geneFeatures[g].location().length();
        int exonLength = 0;
        Location[] subLocs = mergedExonsVector.toArray(new Location[0]);
        for (int l = 0; l < subLocs.length; l++) {
            exonLength += subLocs[l].length();
        }
        // if (geneFeatures[f].name().equals("Cpa6")) {
        // System.err.println("Cpa6 total intron length="+intronLength+" exon length="+exonLength);
        // }
        intronLength -= exonLength;
        if (intronLength <= 0) {
            progressWarningReceived(new IllegalStateException("Intron length of " + intronLength + " for gene " + geneFeatures[g]));
            for (int d = 0; d < data.length; d++) {
                data[d].setValueForProbe(allProbes[g], Float.NaN);
            }
            continue;
        }
        if (exonLength <= 0) {
            progressWarningReceived(new IllegalStateException("Exon length of " + exonLength + " for gene " + geneFeatures[g]));
            for (int d = 0; d < data.length; d++) {
                data[d].setValueForProbe(allProbes[g], Float.NaN);
            }
            continue;
        }
        for (int d = 0; d < data.length; d++) {
            if (cancel) {
                progressCancelled();
                return;
            }
            int totalIntronCount = 0;
            int totalExonCount = 0;
            long[] reads = data[d].getReadsForProbe(new Probe(chrLookup.get(geneFeatures[g].chromosomeName()), geneFeatures[g].location().start(), geneFeatures[g].location().end()));
            for (int r = 0; r < reads.length; r++) {
                if (!readFilter.useRead(geneFeatures[g].location(), reads[r])) {
                    continue;
                }
                int overlap = (Math.min(geneFeatures[g].location().end(), SequenceRead.end(reads[r])) - Math.max(geneFeatures[g].location().start(), SequenceRead.start(reads[r]))) + 1;
                totalIntronCount += overlap;
                // Now we see if we overlap with any of the exons
                for (int s = 0; s < subLocs.length; s++) {
                    if (subLocs[s].start() <= SequenceRead.end(reads[r]) && subLocs[s].end() >= SequenceRead.start(reads[r])) {
                        int exonOverlap = (Math.min(subLocs[s].end(), SequenceRead.end(reads[r])) - Math.max(subLocs[s].start(), SequenceRead.start(reads[r]))) + 1;
                        if (exonOverlap > 0) {
                            totalExonCount += exonOverlap;
                            totalIntronCount -= exonOverlap;
                        }
                    }
                }
            }
            if (totalIntronCount < 0) {
                progressWarningReceived(new SeqMonkException("Intron count of " + totalIntronCount + " for " + geneFeatures[g].name()));
                continue;
            }
            // reads which this count could comprise, rounding down to a whole number.
            if (correctReadLength) {
                totalIntronCount /= readLengths[d];
                totalExonCount /= readLengths[d];
            }
            float intronValue = totalIntronCount;
            float exonValue = totalExonCount;
            // If we're log transforming then we need to set zero values to 0.9
            if (logTransform && intronValue == 0) {
                intronValue = 0.9f;
            }
            if (logTransform && exonValue == 0) {
                exonValue = 0.9f;
            }
            // been asked to.
            if (applyTranscriptLengthCorrection) {
                intronValue /= (intronLength / 1000f);
                exonValue /= (exonLength / 1000f);
            }
            // We also correct by the total read count, or length
            if (correctReadLength) {
                intronValue /= (data[d].getTotalReadCount() / 1000000f);
                exonValue /= (data[d].getTotalReadCount() / 1000000f);
            } else {
                intronValue /= (data[d].getTotalReadLength() / 1000000f);
                exonValue /= (data[d].getTotalReadLength() / 1000000f);
            }
            // Finally we do the log transform if we've been asked to
            if (logTransform) {
                if (intronValue == 0) {
                    intronValue = 0.0001f;
                }
                if (exonValue == 0) {
                    exonValue = 0.0001f;
                }
                intronValue = (float) Math.log(intronValue) / log2;
                exonValue = (float) Math.log(exonValue) / log2;
            }
            // Now we check what value they actually wanted to record
            switch(optionsPanel.quantitationType()) {
                case (EXONS):
                    {
                        data[d].setValueForProbe(allProbes[g], exonValue);
                        break;
                    }
                case (INTRONS):
                    {
                        data[d].setValueForProbe(allProbes[g], intronValue);
                        break;
                    }
                case (RATIO):
                    {
                        if (logTransform) {
                            data[d].setValueForProbe(allProbes[g], exonValue - intronValue);
                        } else {
                            if (intronValue == 0) {
                                intronValue = 0.0001f;
                            }
                            data[d].setValueForProbe(allProbes[g], exonValue / intronValue);
                        }
                        break;
                    }
                default:
                    throw new IllegalStateException("Unknonwn quantitation type " + optionsPanel.quantitationType());
            }
        // if (geneFeatures[f].name().equals("Cpa6")) {
        // try {
        // System.err.println("Stored value was "+data[d].getValueForProbe(allProbes[currentIndex]));
        // } catch (SeqMonkException e) {
        // e.printStackTrace();
        // }
        // }
        }
    }
    StringBuffer quantitationDescription = new StringBuffer();
    quantitationDescription.append("Splicing efficiency quantitation");
    if (correctReadLength) {
        quantitationDescription.append(" counting reads");
    } else {
        quantitationDescription.append(" counting bases");
    }
    if (optionsPanel.logTransform()) {
        quantitationDescription.append(" log transformed");
    }
    if (applyTranscriptLengthCorrection) {
        quantitationDescription.append(" correcting for feature length");
    }
    // TODO: Add more description
    collection().probeSet().setCurrentQuantitation(quantitationDescription.toString());
    quantitatonComplete();
}

Example 64 with Chromosome

use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome in project SeqMonk by s-andrews.

the class WindowedReplicateStatsFilter method generateProbeList.

/* (non-Javadoc)
	 * @see uk.ac.babraham.SeqMonk.Filters.ProbeFilter#generateProbeList()
	 */
@Override
protected void generateProbeList() {
    Chromosome[] chromosomes = collection.genome().getAllChromosomes();
    Vector<ProbeGroupTTestValue> newListProbesVector = new Vector<ProbeGroupTTestValue>();
    for (int c = 0; c < chromosomes.length; c++) {
        progressUpdated("Processing windows on Chr" + chromosomes[c].name(), c, chromosomes.length);
        Probe[] probes = startingList.getProbesForChromosome(chromosomes[c]);
        ProbeGroupGenerator gen = null;
        if (windowType == DISTANCE_WINDOW) {
            gen = new ProbeWindowGenerator(probes, windowSize);
        } else if (windowType == CONSECUTIVE_WINDOW) {
            gen = new ConsecutiveProbeGenerator(probes, windowSize);
        } else if (windowType == FEATURE_WINDOW) {
            gen = new FeatureProbeGroupGenerator(probes, collection.genome().annotationCollection().getFeaturesForType(optionsPanel.featureTypeBox.getSelectedItem().toString()));
        }
        while (true) {
            if (cancel) {
                cancel = false;
                progressCancelled();
                return;
            }
            Probe[] theseProbes = gen.nextSet();
            if (theseProbes == null) {
                // System.err.println("List of probes was null");
                break;
            }
            // We need at least 3 probes in a set to calculate a p-value
            if (theseProbes.length < 3) {
                // System.err.println("Only "+theseProbes.length+" probes in the set");
                continue;
            }
            double[][] values = new double[stores.length][];
            for (int j = 0; j < stores.length; j++) {
                if (splitReplicateSets & stores[j] instanceof ReplicateSet) {
                    values[j] = new double[theseProbes.length * ((ReplicateSet) stores[j]).dataStores().length];
                } else {
                    values[j] = new double[theseProbes.length];
                }
            }
            for (int j = 0; j < stores.length; j++) {
                int index = 0;
                for (int i = 0; i < theseProbes.length; i++) {
                    try {
                        if (splitReplicateSets & stores[j] instanceof ReplicateSet) {
                            DataStore[] localStores = ((ReplicateSet) stores[j]).dataStores();
                            for (int l = 0; l < localStores.length; l++) {
                                values[j][index] = localStores[l].getValueForProbe(theseProbes[i]);
                                index++;
                            }
                        } else {
                            values[j][index] = stores[j].getValueForProbe(theseProbes[i]);
                            index++;
                        }
                    } catch (SeqMonkException e) {
                    }
                }
                if (index != values[j].length) {
                    throw new IllegalStateException("Didn't fill all values total=" + values[j].length + " index=" + index);
                }
            }
            double pValue = 0;
            try {
                if (stores.length == 1) {
                    pValue = TTest.calculatePValue(values[0], 0);
                } else if (stores.length == 2) {
                    pValue = TTest.calculatePValue(values[0], values[1]);
                } else {
                    pValue = AnovaTest.calculatePValue(values);
                }
            } catch (SeqMonkException e) {
                throw new IllegalStateException(e);
            }
            newListProbesVector.add(new ProbeGroupTTestValue(theseProbes, pValue));
        }
    }
    ProbeGroupTTestValue[] newListProbes = newListProbesVector.toArray(new ProbeGroupTTestValue[0]);
    // Do the multi-testing correction if necessary
    if (multiTest) {
        BenjHochFDR.calculateQValues(newListProbes);
    }
    ProbeList newList;
    // We need to handle duplicate hits internally since probe lists can't do
    // this themselves any more.
    Hashtable<Probe, Float> newListTemp = new Hashtable<Probe, Float>();
    if (multiTest) {
        newList = new ProbeList(startingList, "", "", "Q-value");
        for (int i = 0; i < newListProbes.length; i++) {
            if (newListProbes[i].q <= cutoff) {
                Probe[] passedProbes = newListProbes[i].probes;
                for (int p = 0; p < passedProbes.length; p++) {
                    if (newListTemp.containsKey(passedProbes[p])) {
                        // We always give a probe the lowest possible q-value
                        if (newListTemp.get(passedProbes[p]) <= newListProbes[i].q) {
                            continue;
                        }
                    }
                    newListTemp.put(passedProbes[p], (float) newListProbes[i].q);
                }
            }
        }
    } else {
        newList = new ProbeList(startingList, "", "", "P-value");
        for (int i = 0; i < newListProbes.length; i++) {
            if (newListProbes[i].p <= cutoff) {
                Probe[] passedProbes = newListProbes[i].probes;
                for (int p = 0; p < passedProbes.length; p++) {
                    if (newListTemp.containsKey(passedProbes[p])) {
                        // We always give a probe the lowest possible p-value
                        if (newListTemp.get(passedProbes[p]) <= newListProbes[i].p) {
                            continue;
                        }
                    }
                    newListTemp.put(passedProbes[p], (float) newListProbes[i].p);
                }
            }
        }
    }
    // Add the cached hits to the new list
    Enumeration<Probe> en = newListTemp.keys();
    while (en.hasMoreElements()) {
        Probe p = en.nextElement();
        newList.addProbe(p, newListTemp.get(p));
    }
    filterFinished(newList);
}

Example 65 with Chromosome

use of uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome in project SeqMonk by s-andrews.

the class WindowedValuesFilter method generateProbeList.

/* (non-Javadoc)
	 * @see uk.ac.babraham.SeqMonk.Filters.ProbeFilter#generateProbeList()
	 */
@Override
protected void generateProbeList() {
    ProbeList newList = new ProbeList(startingList, "Filtered Probes", "", null);
    Chromosome[] chromosomes = collection.genome().getAllChromosomes();
    for (int c = 0; c < chromosomes.length; c++) {
        HashSet<Probe> passedProbes = new HashSet<Probe>();
        progressUpdated("Processing windows on Chr" + chromosomes[c].name(), c, chromosomes.length);
        Probe[] probes = startingList.getProbesForChromosome(chromosomes[c]);
        ProbeGroupGenerator gen = null;
        if (windowType == DISTANCE_WINDOW) {
            gen = new ProbeWindowGenerator(probes, windowSize);
        } else if (windowType == CONSECUTIVE_WINDOW) {
            gen = new ConsecutiveProbeGenerator(probes, windowSize);
        } else if (windowType == FEATURE_WINDOW) {
            gen = new FeatureProbeGroupGenerator(probes, collection.genome().annotationCollection().getFeaturesForType(optionsPanel.featureTypeBox.getSelectedItem().toString()));
        } else {
            progressExceptionReceived(new SeqMonkException("No window type known with number " + windowType));
            return;
        }
        while (true) {
            if (cancel) {
                cancel = false;
                progressCancelled();
                return;
            }
            Probe[] theseProbes = gen.nextSet();
            if (theseProbes == null) {
                break;
            }
            int count = 0;
            for (int s = 0; s < stores.length; s++) {
                double totalValue = 0;
                for (int i = 0; i < theseProbes.length; i++) {
                    // Get the values for the probes in this set
                    try {
                        totalValue += stores[s].getValueForProbe(theseProbes[i]);
                    } catch (SeqMonkException e) {
                    }
                }
                totalValue /= theseProbes.length;
                // Now we have the value we need to know if it passes the test
                if (upperLimit != null)
                    if (totalValue > upperLimit.doubleValue())
                        continue;
                if (lowerLimit != null)
                    if (totalValue < lowerLimit.doubleValue())
                        continue;
                // This one passes, we can add it to the count
                ++count;
            }
            // probe to the probe set.
            switch(limitType) {
                case EXACTLY:
                    if (count == storesLimit)
                        for (int i = 0; i < theseProbes.length; i++) {
                            if (passedProbes.contains(theseProbes[i]))
                                continue;
                            newList.addProbe(theseProbes[i], null);
                            passedProbes.add(theseProbes[i]);
                        }
                    break;
                case AT_LEAST:
                    if (count >= storesLimit)
                        for (int i = 0; i < theseProbes.length; i++) {
                            if (passedProbes.contains(theseProbes[i]))
                                continue;
                            newList.addProbe(theseProbes[i], null);
                            passedProbes.add(theseProbes[i]);
                        }
                    break;
                case NO_MORE_THAN:
                    if (count <= storesLimit)
                        for (int i = 0; i < theseProbes.length; i++) {
                            if (passedProbes.contains(theseProbes[i]))
                                continue;
                            newList.addProbe(theseProbes[i], null);
                            passedProbes.add(theseProbes[i]);
                        }
                    break;
            }
        }
    }
    filterFinished(newList);
}