Search in sources :

Example 56 with SeqMonkException

use of uk.ac.babraham.SeqMonk.SeqMonkException in project SeqMonk by s-andrews.

the class SeqMonkParser method parseSamples.

/**
 * Parses the list of samples.
 *
 * @param sections The tab split initial samples line
 * @throws SeqMonkException
 * @throws IOException Signals that an I/O exception has occurred.
 */
private void parseSamples(String[] sections) throws SeqMonkException, IOException {
    if (sections.length != 2) {
        throw new SeqMonkException("Samples line didn't contain 2 sections");
    }
    if (!sections[0].equals("Samples")) {
        throw new SeqMonkException("Couldn't find expected samples line");
    }
    int n = Integer.parseInt(sections[1]);
    // We need to keep the Data Sets around to add data to later.
    dataSets = new DataSet[n];
    for (int i = 0; i < n; i++) {
        sections = br.readLine().split("\\t");
        // anything else is assumed to be a normal dataset.
        if (sections.length == 1) {
            dataSets[i] = new DataSet(sections[0], "Not known", DataSet.DUPLICATES_REMOVE_NO);
        } else if (sections.length == 2) {
            dataSets[i] = new DataSet(sections[0], sections[1], DataSet.DUPLICATES_REMOVE_NO);
        } else if (sections.length == 3) {
            if (sections[2].equals("HiC")) {
                dataSets[i] = new PairedDataSet(sections[0], sections[1], DataSet.DUPLICATES_REMOVE_NO, 0, false);
            } else {
                dataSets[i] = new DataSet(sections[0], sections[1], DataSet.DUPLICATES_REMOVE_NO);
            }
        }
        application.dataCollection().addDataSet(dataSets[i]);
    }
    // Immediately after the list of samples comes the lists of reads
    String line;
    // Iterate through the number of samples
    for (int i = 0; i < n; i++) {
        Enumeration<ProgressListener> en = listeners.elements();
        while (en.hasMoreElements()) {
            en.nextElement().progressUpdated("Reading data for " + dataSets[i].name(), i * 10, n * 10);
        }
        // The first line is
        line = br.readLine();
        sections = line.split("\t");
        if (sections.length != 2) {
            throw new SeqMonkException("Read line " + i + " didn't contain 2 sections");
        }
        int readCount = Integer.parseInt(sections[0]);
        // In versions prior to 7 we encoded everything on every line separately
        if (thisDataVersion < 7) {
            for (int r = 0; r < readCount; r++) {
                if ((r % (1 + (readCount / 10))) == 0) {
                    Enumeration<ProgressListener> en2 = listeners.elements();
                    while (en2.hasMoreElements()) {
                        en2.nextElement().progressUpdated("Reading data for " + dataSets[i].name(), i * 10 + (r / (1 + (readCount / 10))), n * 10);
                    }
                }
                line = br.readLine();
                if (line == null) {
                    throw new SeqMonkException("Ran out of data whilst parsing reads for sample " + i);
                }
                sections = line.split("\t");
                Chromosome c;
                try {
                    c = application.dataCollection().genome().getChromosome(sections[0]).chromosome();
                } catch (Exception sme) {
                    Enumeration<ProgressListener> e = listeners.elements();
                    while (e.hasMoreElements()) {
                        e.nextElement().progressWarningReceived(new SeqMonkException("Read from sample " + i + " could not be mapped to chromosome '" + sections[0] + "'"));
                    }
                    continue;
                }
                int start = Integer.parseInt(sections[1]);
                int end = Integer.parseInt(sections[2]);
                int strand = Integer.parseInt(sections[3]);
                try {
                    dataSets[i].addData(c, SequenceRead.packPosition(start, end, strand), true);
                } catch (SeqMonkException ex) {
                    Enumeration<ProgressListener> e = listeners.elements();
                    while (e.hasMoreElements()) {
                        e.nextElement().progressWarningReceived(ex);
                        continue;
                    }
                }
            }
        } else if (dataSets[i] instanceof PairedDataSet) {
            // Paired Data sets have a different format, with a packed position for
            // the first read, then a chromosome name and packed position for the second
            // read.
            // From v13 onwards we started entering HiC data pairs in both directions so
            // the data would be able to be read back in without sorting it.
            boolean hicDataIsPresorted = thisDataVersion >= 13;
            // We use this to update the progress bar every time we move to a new chromosome
            int seenChromosomeCount = 0;
            while (true) {
                // The first line should be the chromosome and a number of reads
                line = br.readLine();
                if (line == null) {
                    throw new SeqMonkException("Ran out of data whilst parsing reads for sample " + i);
                }
                // A blank line indicates the end of the sample
                if (line.length() == 0)
                    break;
                sections = line.split("\\t");
                Chromosome c = application.dataCollection().genome().getChromosome(sections[0]).chromosome();
                ++seenChromosomeCount;
                progressUpdated("Reading data for " + dataSets[i].name(), i * application.dataCollection().genome().getChromosomeCount() + seenChromosomeCount, n * application.dataCollection().genome().getChromosomeCount());
                int chrReadCount = Integer.parseInt(sections[1]);
                for (int r = 0; r < chrReadCount; r++) {
                    line = br.readLine();
                    if (line == null) {
                        throw new SeqMonkException("Ran out of data whilst parsing reads for sample " + i);
                    }
                    /*
						 * We used to have a split("\t") here, but this turned out to be the bottleneck
						 * which hugely restricted the speed at which the data could be read.  Switching
						 * for a set of index calls and substring makes this *way* faster.
						 */
                    int tabPosition = line.indexOf('\t');
                    long packedPosition = quickParseLong(line, 0);
                    int tabPosition2 = line.indexOf('\t', tabPosition + 1);
                    Chromosome c2 = application.dataCollection().genome().getChromosome(line.substring(tabPosition, tabPosition2)).chromosome();
                    long packedPosition2 = quickParseLong(line, tabPosition2 + 1);
                    try {
                        dataSets[i].addData(c, packedPosition, hicDataIsPresorted);
                        dataSets[i].addData(c2, packedPosition2, hicDataIsPresorted);
                    } catch (SeqMonkException ex) {
                        Enumeration<ProgressListener> e = listeners.elements();
                        while (e.hasMoreElements()) {
                            e.nextElement().progressWarningReceived(ex);
                            continue;
                        }
                    }
                }
            }
        } else {
            // In versions after 7 we split the section up into chromosomes
            // so we don't put the chromosome on every line, and we put out
            // the packed double value rather than the individual start, end
            // and strand
            // As of version 12 we collapse repeated reads into one line with
            // a count after it, so we need to check for this.
            // We keep count of reads processed to update the progress listeners
            int readsRead = 0;
            while (true) {
                // The first line should be the chromosome and a number of reads
                line = br.readLine();
                if (line == null) {
                    throw new SeqMonkException("Ran out of data whilst parsing reads for sample " + i);
                }
                // A blank line indicates the end of the sample
                if (line.length() == 0)
                    break;
                sections = line.split("\\t");
                // We don't try to capture this exception since we can't then process
                // any of the reads which follow.
                Chromosome c = application.dataCollection().genome().getChromosome(sections[0]).chromosome();
                int chrReadCount = Integer.parseInt(sections[1]);
                for (int r = 0; r < chrReadCount; r++) {
                    if ((readsRead % (1 + (readCount / 10))) == 0) {
                        Enumeration<ProgressListener> en2 = listeners.elements();
                        while (en2.hasMoreElements()) {
                            en2.nextElement().progressUpdated("Reading data for " + dataSets[i].name(), i * 10 + (readsRead / (1 + (readCount / 10))), n * 10);
                        }
                    }
                    line = br.readLine();
                    if (line == null) {
                        throw new SeqMonkException("Ran out of data whilst parsing reads for sample " + i);
                    }
                    // We use some custom parsing code to efficiently extract the packed position and
                    // count from the line.  This avoids having to use the generic parsing code or
                    // doing a split to determine the tab location.
                    long packedPosition = 0;
                    int count = 0;
                    // Check for a sign.
                    int sign = 1;
                    int len = line.length();
                    if (len == 0) {
                        System.err.println("Yuk - blank line found in sample data");
                    }
                    char ch = line.charAt(0);
                    if (ch == '-') {
                        sign = -1;
                    } else {
                        packedPosition = ch - '0';
                    }
                    // Build the numbers.
                    int j = 1;
                    int tabIndex = 0;
                    while (j < len) {
                        if (line.charAt(j) == '\t') {
                            tabIndex = j++;
                        } else if (tabIndex > 0) {
                            // char is ASCII so need to subtract char of '0' to get correct int value
                            count = (count * 10) + (line.charAt(j++) - '0');
                        } else {
                            packedPosition = (packedPosition * 10) + (line.charAt(j++) - '0');
                        }
                    }
                    packedPosition = sign * packedPosition;
                    if (count == 0) {
                        // There was no count so there's only one of these
                        count = 1;
                    } else {
                        // Update the true read count with the additional reads we get from this line
                        r += (count - 1);
                    }
                    // System.out.println("From line "+line+" got packed="+packedPosition+" count="+count);
                    dataSets[i].addData(c, packedPosition, count, true);
                    readsRead += count;
                }
            }
        }
        // We've now read all of the data for this sample so we can compact it
        Enumeration<ProgressListener> en2 = listeners.elements();
        while (en2.hasMoreElements()) {
            en2.nextElement().progressUpdated("Caching data for " + dataSets[i].name(), (i + 1) * 10, n * 10);
        }
        dataSets[i].finalise();
    }
}
Also used : Enumeration(java.util.Enumeration) ProgressListener(uk.ac.babraham.SeqMonk.DataTypes.ProgressListener) DataSet(uk.ac.babraham.SeqMonk.DataTypes.DataSet) PairedDataSet(uk.ac.babraham.SeqMonk.DataTypes.PairedDataSet) PairedDataSet(uk.ac.babraham.SeqMonk.DataTypes.PairedDataSet) Chromosome(uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome) SeqMonkException(uk.ac.babraham.SeqMonk.SeqMonkException) IOException(java.io.IOException) UnknownHostException(java.net.UnknownHostException) FileNotFoundException(java.io.FileNotFoundException) SeqMonkException(uk.ac.babraham.SeqMonk.SeqMonkException)

Example 57 with SeqMonkException

use of uk.ac.babraham.SeqMonk.SeqMonkException in project SeqMonk by s-andrews.

the class VisibleStoresParser method processNormalDataStore.

private DataSet processNormalDataStore(DataStore store) {
    int extendBy = prefs.extendReads();
    boolean reverse = prefs.reverseReads();
    boolean removeStrand = prefs.removeStrandInfo();
    DataSet newData = new DataSet(store.name() + "_reimport", "Reimported from " + store.name(), prefs.removeDuplicates());
    // Now process the data
    Chromosome[] chrs = dataCollection().genome().getAllChromosomes();
    for (int c = 0; c < chrs.length; c++) {
        progressUpdated("Processing " + store.name() + " chr " + chrs[c].name(), c, chrs.length);
        ReadsWithCounts reads = store.getReadsForChromosome(chrs[c]);
        Feature[] features = null;
        if (filterByFeature) {
            features = collection.genome().annotationCollection().getFeaturesForType(chrs[c], featureType);
            Arrays.sort(features);
        }
        int currentFeaturePostion = 0;
        for (int r = 0; r < reads.reads.length; r++) {
            for (int ct = 0; ct < reads.counts[r]; ct++) {
                long thisRead = reads.reads[r];
                if (cancel) {
                    progressCancelled();
                    return null;
                }
                if (downsample && downsampleProbabilty < 1) {
                    if (Math.random() > downsampleProbabilty) {
                        continue;
                    }
                }
                long read;
                int start = SequenceRead.start(thisRead);
                int end = SequenceRead.end(thisRead);
                int strand = SequenceRead.strand(thisRead);
                if (filterByStrand) {
                    if (strand == Location.FORWARD && !keepForward)
                        continue;
                    if (strand == Location.REVERSE && !keepReverse)
                        continue;
                    if (strand == Location.UNKNOWN && !keepUnknown)
                        continue;
                }
                if (filterByLength) {
                    int length = SequenceRead.length(thisRead);
                    if (minLength != null && length < minLength)
                        continue;
                    if (maxLength != null && length > maxLength)
                        continue;
                }
                if (strand == Location.FORWARD) {
                    start += forwardOffset;
                    end += forwardOffset;
                }
                if (strand == Location.REVERSE) {
                    start -= reverseOffset;
                    end -= reverseOffset;
                }
                if (filterByFeature && features.length == 0 && !excludeFeature)
                    continue;
                if (filterByFeature && features.length > 0) {
                    // See if we're comparing against the right feature
                    while (SequenceRead.start(thisRead) > features[currentFeaturePostion].location().end() && currentFeaturePostion < (features.length - 1)) {
                        currentFeaturePostion++;
                    }
                    // Test to see if we overlap
                    if (SequenceRead.overlaps(thisRead, features[currentFeaturePostion].location().packedPosition())) {
                        if (excludeFeature)
                            continue;
                    } else {
                        if (!excludeFeature)
                            continue;
                    }
                }
                if (reverse) {
                    if (strand == Location.FORWARD) {
                        strand = Location.REVERSE;
                    } else if (strand == Location.REVERSE) {
                        strand = Location.FORWARD;
                    }
                }
                if (removeStrand) {
                    strand = Location.UNKNOWN;
                }
                if (extractCentres) {
                    int centre = start + ((end - start) / 2);
                    start = centre - centreExtractContext;
                    end = centre + centreExtractContext;
                }
                if (extendBy != 0) {
                    // We now allow negative extensions to shorten reads
                    if (strand == Location.FORWARD || strand == Location.UNKNOWN) {
                        end += extendBy;
                        if (end < start)
                            end = start;
                    } else if (strand == Location.REVERSE) {
                        start -= extendBy;
                        if (start > end)
                            start = end;
                    }
                }
                // We don't allow reads before the start of the chromosome
                if (start < 1) {
                    int overrun = (0 - start) + 1;
                    progressWarningReceived(new SeqMonkException("Reading position " + start + " was " + overrun + "bp before the start of chr" + chrs[c].name() + " (" + chrs[c].length() + ")"));
                    continue;
                }
                // We also don't allow readings which are beyond the end of the chromosome
                if (end > chrs[c].length()) {
                    int overrun = end - chrs[c].length();
                    progressWarningReceived(new SeqMonkException("Reading position " + end + " was " + overrun + "bp beyond the end of chr" + chrs[c].name() + " (" + chrs[c].length() + ")"));
                    continue;
                }
                // We can now make the new reading
                try {
                    read = SequenceRead.packPosition(start, end, strand);
                    if (!prefs.isHiC()) {
                        // HiC additions are deferred until we know the other end is OK too.
                        newData.addData(chrs[c], read);
                    }
                } catch (SeqMonkException e) {
                    progressWarningReceived(e);
                    continue;
                }
            }
        }
    }
    return newData;
}
Also used : DataSet(uk.ac.babraham.SeqMonk.DataTypes.DataSet) PairedDataSet(uk.ac.babraham.SeqMonk.DataTypes.PairedDataSet) Chromosome(uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome) ReadsWithCounts(uk.ac.babraham.SeqMonk.DataTypes.Sequence.ReadsWithCounts) SeqMonkException(uk.ac.babraham.SeqMonk.SeqMonkException) Feature(uk.ac.babraham.SeqMonk.DataTypes.Genome.Feature)

Example 58 with SeqMonkException

use of uk.ac.babraham.SeqMonk.SeqMonkException in project SeqMonk by s-andrews.

the class SimilarProbeListsDialog method mouseClicked.

public void mouseClicked(MouseEvent me) {
    // We're only interested in double clicks
    if (me.getClickCount() != 2)
        return;
    // If the table hasn't been drawn yet then bail out
    if (table == null)
        return;
    // Get the selected row
    int selectedRow = table.getSelectedRow();
    if (selectedRow < 0)
        return;
    ProbeList newActiveList = (ProbeList) table.getModel().getValueAt(selectedRow, 0);
    try {
        collection.probeSet().setActiveList(newActiveList);
    } catch (SeqMonkException e) {
        throw new IllegalStateException(e);
    }
}
Also used : ProbeList(uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList) SeqMonkException(uk.ac.babraham.SeqMonk.SeqMonkException)

Example 59 with SeqMonkException

use of uk.ac.babraham.SeqMonk.SeqMonkException in project SeqMonk by s-andrews.

the class HeatmapGenomePanel method paint.

public void paint(Graphics g) {
    super.paint(g);
    if (drawnPixels.length != getWidth() || drawnPixels[0].length != getHeight()) {
        drawnPixels = new boolean[getWidth()][getHeight()];
    } else {
        for (int i = 0; i < getWidth(); i++) {
            for (int j = 0; j < getHeight(); j++) {
                drawnPixels[i][j] = false;
            }
        }
    }
    g.setColor(Color.WHITE);
    g.fillRect(0, 0, getWidth(), getHeight());
    g.setColor(Color.BLACK);
    // Add a label at the top to signify current filters
    StringBuffer topLabel = new StringBuffer();
    // First include the position
    Chromosome xStartChr = getChrForPosition(currentXStartBp);
    if (xStartChr != null) {
        topLabel.append("X=Chr");
        topLabel.append(xStartChr.name());
        topLabel.append(":");
        topLabel.append(PositionFormat.formatLength(currentXStartBp - chromosomeBaseOffsets.get(xStartChr)));
        topLabel.append("-");
        Chromosome xEndChr = getChrForPosition(currentXEndBp);
        if (xStartChr != xEndChr) {
            topLabel.append("Chr");
            topLabel.append(xEndChr.name());
            topLabel.append(":");
        }
        topLabel.append(PositionFormat.formatLength(currentXEndBp - chromosomeBaseOffsets.get(xEndChr)));
        if (probeSortingOrder == null) {
            topLabel.append(" Y=Chr");
            Chromosome yStartChr = getChrForPosition(currentYStartBp);
            topLabel.append(yStartChr.name());
            topLabel.append(":");
            topLabel.append(PositionFormat.formatLength(currentYStartBp - chromosomeBaseOffsets.get(yStartChr)));
            topLabel.append("-");
            Chromosome yEndChr = getChrForPosition(currentYEndBp);
            if (yStartChr != yEndChr) {
                topLabel.append("Chr");
                topLabel.append(yEndChr.name());
                topLabel.append(":");
            }
            topLabel.append(PositionFormat.formatLength(currentYEndBp - chromosomeBaseOffsets.get(yEndChr)));
        }
        topLabel.append(" ");
    }
    // Now append any current limits
    if (matrix.currentMinStrength() > 0) {
        topLabel.append("Strength > ");
        topLabel.append(df.format(matrix.currentMinStrength()));
        topLabel.append(" ");
    }
    if (matrix.minDifference() > 0) {
        topLabel.append("Difference > ");
        topLabel.append(df.format(matrix.minDifference()));
        topLabel.append(" ");
    }
    if (matrix.currentMaxSignficance() < 1) {
        topLabel.append("P-value < ");
        topLabel.append(df.format(matrix.currentMaxSignficance()));
        topLabel.append(" ");
    }
    if (topLabel.length() == 0) {
        topLabel.append("No filters");
    }
    g.drawString(topLabel.toString(), getWidth() / 2 - (g.getFontMetrics().stringWidth(topLabel.toString()) / 2), 15 + (g.getFontMetrics().getAscent() / 2));
    Chromosome[] chrs = genome.getAllChromosomes();
    Arrays.sort(chrs);
    // Find the max height and width of the chromosome names in this genome
    if (maxNameWidth == 0) {
        nameHeight = g.getFontMetrics().getHeight();
        maxNameWidth = 0;
        for (int c = 0; c < chrs.length; c++) {
            int thisWidth = g.getFontMetrics().stringWidth(chrs[c].name());
            if (thisWidth > maxNameWidth)
                maxNameWidth = thisWidth;
        }
        // Give both the width and height a bit of breathing space
        nameHeight += 6;
        maxNameWidth += 6;
    }
    // Make the background of the plot black
    g.setColor(Color.WHITE);
    g.fillRect(maxNameWidth, 30, getWidth() - (maxNameWidth + 10), getHeight() - (nameHeight + 30));
    // Draw the actual data
    InteractionProbePair[] interactions = matrix.filteredInteractions();
    // Cache some values for use with the quantitation colouring
    double minQuantitatedValue;
    double maxQuantitatedValue;
    if (DisplayPreferences.getInstance().getScaleType() == DisplayPreferences.SCALE_TYPE_POSITIVE) {
        minQuantitatedValue = 0;
        maxQuantitatedValue = DisplayPreferences.getInstance().getMaxDataValue();
    } else {
        maxQuantitatedValue = DisplayPreferences.getInstance().getMaxDataValue();
        minQuantitatedValue = 0 - maxQuantitatedValue;
    }
    for (int i = 0; i < interactions.length; i++) {
        // the plot symmetrical
        for (int forRev = 0; forRev <= 1; forRev++) {
            int yIndex;
            Probe probe1;
            Probe probe2;
            if (forRev == 0) {
                yIndex = interactions[i].probe2Index();
                probe1 = interactions[i].probe1();
                probe2 = interactions[i].probe2();
            } else {
                yIndex = interactions[i].probe1Index();
                probe2 = interactions[i].probe1();
                probe1 = interactions[i].probe2();
            }
            if (probeSortingOrder != null) {
                yIndex = probeSortingOrder[yIndex];
                if (yIndex < currentYStartIndex || yIndex > currentYEndIndex) {
                    // System.err.println("Skipping for y zoom");
                    continue;
                }
            }
            int xStart = getXForPosition(chromosomeBaseOffsets.get(probe1.chromosome()) + probe1.start());
            if (xStart < maxNameWidth)
                continue;
            if (chromosomeBaseOffsets.get(probe1.chromosome()) + probe1.start() > currentXEndBp) {
                // System.err.println("Skipping for x end");
                continue;
            }
            int xEnd = getXForPosition(chromosomeBaseOffsets.get(probe1.chromosome()) + probe1.end());
            if (xEnd > getWidth() - 10)
                continue;
            if (chromosomeBaseOffsets.get(probe1.chromosome()) + probe1.end() < currentXStartBp) {
                // System.err.println("Skipping for x start");
                continue;
            }
            int yStart;
            int yEnd;
            if (probeSortingOrder == null) {
                if (chromosomeBaseOffsets.get(probe2.chromosome()) + probe2.start() > currentYEndBp)
                    continue;
                if (chromosomeBaseOffsets.get(probe2.chromosome()) + probe2.end() < currentYStartBp)
                    continue;
                yStart = getYForPosition(chromosomeBaseOffsets.get(probe2.chromosome()) + probe2.start());
                yEnd = getYForPosition(chromosomeBaseOffsets.get(probe2.chromosome()) + probe2.end());
            } else {
                yStart = getYForIndex(probeSortingOrder[yIndex]);
                yEnd = getYForIndex(probeSortingOrder[yIndex] + 1);
            }
            if (yStart > getHeight() - nameHeight) {
                continue;
            }
            if (yEnd < 30)
                continue;
            if (xEnd - xStart < 3) {
                xEnd += 1;
                xStart -= 1;
            }
            if (yStart - yEnd < 3) {
                // System.err.println("Expanding y selection");
                yEnd -= 1;
                yStart += 1;
            }
            // To be skipped there has to be colour at two corners and the middle.
            if (drawnPixels[xStart][yEnd] && drawnPixels[xEnd][yStart] && drawnPixels[xStart + ((xEnd - xStart) / 2)][yEnd + ((yStart - yEnd) / 2)]) {
                // System.err.println("Skipping for overlap with existing");
                continue;
            }
            switch(matrix.currentColourSetting()) {
                case HeatmapMatrix.COLOUR_BY_OBS_EXP:
                    if (matrix.initialMinStrength() < 1) {
                        // They're interested in depletion as well as enrichment.
                        // Make a symmetrical gradient around 0 and the max strength
                        g.setColor(matrix.colourGradient().getColor(Math.log10(interactions[i].strength()), Math.log10(1 / matrix.maxValue()), Math.log10(matrix.maxValue())));
                    } else {
                        g.setColor(matrix.colourGradient().getColor(Math.log10(interactions[i].strength() - matrix.initialMinStrength()), Math.log10(matrix.initialMinStrength()), Math.log10(matrix.maxValue() - matrix.initialMinStrength())));
                    }
                    break;
                case HeatmapMatrix.COLOUR_BY_INTERACTIONS:
                    g.setColor(matrix.colourGradient().getColor(interactions[i].absolute(), matrix.initialMinAbsolute(), 50));
                    break;
                case HeatmapMatrix.COLOUR_BY_P_VALUE:
                    g.setColor(matrix.colourGradient().getColor(Math.log10(interactions[i].signficance()) * -10, Math.log10(matrix.initialMaxSignificance()) * -10, 50));
                    break;
                case HeatmapMatrix.COLOUR_BY_QUANTITATION:
                    Probe probeForQuantitation;
                    if (forRev == 0) {
                        probeForQuantitation = interactions[i].lowestProbe();
                    } else {
                        probeForQuantitation = interactions[i].highestProbe();
                    }
                    try {
                        g.setColor(matrix.colourGradient().getColor(((DataStore) dataSet).getValueForProbe(probeForQuantitation), minQuantitatedValue, maxQuantitatedValue));
                    } catch (SeqMonkException e) {
                    }
                    break;
            }
            g.fillRect(xStart, yEnd, xEnd - xStart, yStart - yEnd);
            // If we're looking for selected probes check this now.
            if (displayXProbe || displayYProbe) {
                if (xStart <= displayX && xEnd >= displayX && yEnd <= displayY && yStart >= displayY) {
                    if (displayXProbe) {
                        DisplayPreferences.getInstance().setLocation(interactions[i].probe1().chromosome(), interactions[i].probe1().packedPosition());
                        displayXProbe = false;
                    }
                    if (displayYProbe) {
                        DisplayPreferences.getInstance().setLocation(interactions[i].probe2().chromosome(), interactions[i].probe2().packedPosition());
                        displayYProbe = false;
                    }
                }
            }
            // them again
            for (int x = Math.min(xStart, xEnd); x <= Math.min(xStart, xEnd) + Math.abs(xStart - xEnd); x++) {
                for (int y = Math.min(yStart, yEnd); y <= Math.min(yStart, yEnd) + Math.abs(yStart - yEnd); y++) {
                    drawnPixels[x][y] = true;
                }
            }
        }
    }
    // System.err.println("Skipped "+skipped+" and drew "+drawn+" out of "+(interactions.length*2)+" interactions");
    // Draw the chromosome lines
    g.setColor(Color.GRAY);
    // Draw Chr Lines on X axis
    long runningGenomeLength = 0;
    for (int c = 0; c < chrs.length; c++) {
        int startPos = getXForPosition(runningGenomeLength);
        int endPos = getXForPosition(runningGenomeLength + chrs[c].length());
        if (c > 0) {
            if (startPos >= maxNameWidth && startPos <= getWidth() - 10) {
                g.drawLine(startPos, 30, startPos, getHeight() - nameHeight);
            }
        }
        if (c + 1 == chrs.length) {
            if (endPos >= maxNameWidth && endPos <= getWidth() - 10) {
                g.drawLine(endPos, 30, endPos, getHeight() - nameHeight);
            }
        }
        int nameWidth = g.getFontMetrics().stringWidth(chrs[c].name());
        g.drawString(chrs[c].name(), (startPos + ((endPos - startPos) / 2)) - (nameWidth / 2), getHeight() - 3);
        runningGenomeLength += chrs[c].length();
    }
    // Draw Chr Lines on Y axis
    if (probeSortingOrder == null) {
        runningGenomeLength = 0;
        for (int c = 0; c < chrs.length; c++) {
            int startPos = getYForPosition(runningGenomeLength);
            int endPos = getYForPosition(runningGenomeLength + chrs[c].length());
            if (c > 0) {
                if (startPos <= getHeight() - nameHeight && startPos >= 30) {
                    g.drawLine(maxNameWidth, startPos, getWidth() - 10, startPos);
                }
            }
            if (c + 1 == chrs.length) {
                if (endPos <= getHeight() - nameHeight && endPos >= 30) {
                    g.drawLine(maxNameWidth, endPos, getWidth() - 10, endPos);
                }
            }
            int nameWidth = g.getFontMetrics().stringWidth(chrs[c].name());
            g.drawString(chrs[c].name(), (maxNameWidth / 2) - (nameWidth / 2), (endPos + ((startPos - endPos) / 2)) + (g.getFontMetrics().getAscent() / 2));
            runningGenomeLength += chrs[c].length();
        }
    } else {
        int runningListPosition = 0;
        for (int l = 0; l < clusterIntervals.length; l++) {
            runningListPosition += clusterIntervals[l];
            if (runningListPosition < currentYStartIndex)
                continue;
            if (runningListPosition > currentYEndIndex)
                break;
            int pos = getYForIndex(runningListPosition);
            g.drawLine(maxNameWidth, pos, getWidth() - 10, pos);
        }
    }
    // Draw the axes
    g.drawLine(maxNameWidth, getHeight() - nameHeight, getWidth() - 10, getHeight() - nameHeight);
    g.drawLine(maxNameWidth, getHeight() - nameHeight, maxNameWidth, 30);
    // Draw a selection if we're making one
    if (makingSelection) {
        g.setColor(ColourScheme.DRAGGED_SELECTION);
        g.drawRect(Math.min(selectionEndX, selectionStartX), Math.min(selectionEndY, selectionStartY), Math.abs(selectionEndX - selectionStartX), Math.abs(selectionEndY - selectionStartY));
    }
    displayXProbe = false;
    displayYProbe = false;
}
Also used : InteractionProbePair(uk.ac.babraham.SeqMonk.DataTypes.Interaction.InteractionProbePair) DataStore(uk.ac.babraham.SeqMonk.DataTypes.DataStore) HiCDataStore(uk.ac.babraham.SeqMonk.DataTypes.HiCDataStore) Chromosome(uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome) SeqMonkException(uk.ac.babraham.SeqMonk.SeqMonkException) Probe(uk.ac.babraham.SeqMonk.DataTypes.Probes.Probe)

Example 60 with SeqMonkException

use of uk.ac.babraham.SeqMonk.SeqMonkException in project SeqMonk by s-andrews.

the class EdgeRFilter 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/EdgeRFilter/edger_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", "FDR");
        } else {
            template.setValue("CORRECTED", "PValue");
        }
        StringBuffer sb = new StringBuffer();
        for (int i = 0; i < fromStores.length; i++) {
            if (i > 0)
                sb.append(",");
            sb.append("1");
        }
        for (int i = 0; i < toStores.length; i++) {
            sb.append(",");
            sb.append("2");
        }
        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]);
                if (value != (int) value) {
                    progressExceptionReceived(new IllegalArgumentException("Inputs to the EdgeR filter MUST be raw, incorrected counts, not things like " + value));
                    pr.close();
                    return;
                }
                sb.append(value);
            }
            for (int i = 0; i < toStores.length; i++) {
                sb.append("\t");
                value = toStores[i].getValueForProbe(probes[p]);
                if (value != (int) value) {
                    progressExceptionReceived(new IllegalArgumentException("Inputs to the EdgeR filter MUST be raw, incorrected counts, not things like " + value));
                    pr.close();
                    return;
                }
                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(new SeqMonkException("R Script failed"));
                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 = Float.parseFloat(sections[sections.length - 1]);
            newList.addProbe(probes[probeIndex], pValue);
        }
        br.close();
        runner.cleanUp();
        filterFinished(newList);
    } catch (Exception ioe) {
        progressExceptionReceived(ioe);
        return;
    }
// filterFinished(newList);
}
Also used : ProbeList(uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList) Probe(uk.ac.babraham.SeqMonk.DataTypes.Probes.Probe) RProgressListener(uk.ac.babraham.SeqMonk.R.RProgressListener) SeqMonkException(uk.ac.babraham.SeqMonk.SeqMonkException) Template(uk.ac.babraham.SeqMonk.Utilities.Templates.Template) ProgressRecordDialog(uk.ac.babraham.SeqMonk.Dialogs.ProgressRecordDialog) DataStore(uk.ac.babraham.SeqMonk.DataTypes.DataStore) BufferedReader(java.io.BufferedReader) FileReader(java.io.FileReader) SeqMonkException(uk.ac.babraham.SeqMonk.SeqMonkException) File(java.io.File) RScriptRunner(uk.ac.babraham.SeqMonk.R.RScriptRunner) PrintWriter(java.io.PrintWriter)

Aggregations

SeqMonkException (uk.ac.babraham.SeqMonk.SeqMonkException)91 Probe (uk.ac.babraham.SeqMonk.DataTypes.Probes.Probe)49 ProbeList (uk.ac.babraham.SeqMonk.DataTypes.Probes.ProbeList)30 Chromosome (uk.ac.babraham.SeqMonk.DataTypes.Genome.Chromosome)22 Vector (java.util.Vector)21 DataSet (uk.ac.babraham.SeqMonk.DataTypes.DataSet)20 File (java.io.File)19 DataStore (uk.ac.babraham.SeqMonk.DataTypes.DataStore)17 BufferedReader (java.io.BufferedReader)16 FileReader (java.io.FileReader)16 ChromosomeWithOffset (uk.ac.babraham.SeqMonk.Utilities.ChromosomeWithOffset)14 PairedDataSet (uk.ac.babraham.SeqMonk.DataTypes.PairedDataSet)13 FileInputStream (java.io.FileInputStream)11 IOException (java.io.IOException)11 InputStreamReader (java.io.InputStreamReader)11 GZIPInputStream (java.util.zip.GZIPInputStream)11 HiCDataStore (uk.ac.babraham.SeqMonk.DataTypes.HiCDataStore)8 ProgressListener (uk.ac.babraham.SeqMonk.DataTypes.ProgressListener)8 FileNotFoundException (java.io.FileNotFoundException)7 SequenceReadWithChromosome (uk.ac.babraham.SeqMonk.DataTypes.Sequence.SequenceReadWithChromosome)7