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

Example 11 with HiCDataStore

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

the class HiCPrevNextQuantitation method run.

/* (non-Javadoc)
	 * @see java.lang.Runnable#run()
	 */
public void run() {
    Probe[] probes = application.dataCollection().probeSet().getAllProbes();
    Arrays.sort(probes);
    for (int p = 0; p < probes.length; p++) {
        // See if we need to quit
        if (cancel) {
            progressCancelled();
            return;
        }
        if (probes[p].start() - distance < 0 || probes[p].end() + distance > probes[p].chromosome().length()) {
            for (int d = 0; d < data.length; d++) {
                ((DataStore) data[d]).setValueForProbe(probes[p], 0);
            }
            continue;
        }
        // Make up a pseudo probe for the previous and next regions
        Probe previousProbe = new Probe(probes[p].chromosome(), probes[p].start() - distance, probes[p].start());
        Probe nextProbe = new Probe(probes[p].chromosome(), probes[p].end(), probes[p].end() + distance);
        progressUpdated(p, probes.length);
        for (int d = 0; d < data.length; d++) {
            // Get the counts for the previous and next probes
            int previousTotalCount = data[d].getHiCReadCountForProbe(previousProbe);
            HiCHitCollection hiCHits = data[d].getHiCReadsForProbe(probes[p]);
            int nextTotalCount = data[d].getHiCReadCountForProbe(nextProbe);
            // any reads then give up on this one.
            if (previousTotalCount == 0 || nextTotalCount == 0) {
                ((DataStore) data[d]).setValueForProbe(probes[p], 0);
                continue;
            }
            int previousCount = 0;
            int nextCount = 0;
            long[] hitPositions = hiCHits.getHitPositionsForChromosome(hiCHits.getSourceChromosomeName());
            SequenceRead.sort(hitPositions);
            for (int r = 0; r < hitPositions.length; r++) {
                // Check if we can ignore this one
                if (removeDuplicates) {
                    if (r > 0 && hitPositions[r] == hitPositions[r - 1])
                        continue;
                }
                // Check if the other end maps to either the prev or next probe
                if (SequenceRead.overlaps(hitPositions[r], previousProbe.packedPosition())) {
                    ++previousCount;
                }
                if (SequenceRead.overlaps(hitPositions[r], nextProbe.packedPosition())) {
                    ++nextCount;
                }
            }
            // If both of the actual counts are zero then don't try to calculate a real value
            if (previousCount == 0 && nextCount == 0) {
                ((DataStore) data[d]).setValueForProbe(probes[p], 0);
                continue;
            }
            // Basically what happens here is that we calculate a normal X^2 value and then
            // turn it into a negative depending on the direction of the deviation from the
            // expected proportion.
            // We need to calculate expected values for the two directions based on the proportion
            // of reads falling into those two regions in total.
            // Correct the counts based on the proportions of the two totals
            double expectedProportion = (previousTotalCount / (double) (nextTotalCount + previousTotalCount));
            int totalObservedCount = previousCount + nextCount;
            double expectedPreviousCount = totalObservedCount * expectedProportion;
            double expectedNextCount = totalObservedCount - expectedPreviousCount;
            // Now we can calculate the X^2 values to compare the expected and observed values
            double chisquare = (Math.pow(previousCount - expectedPreviousCount, 2) / expectedPreviousCount) + (Math.pow(nextCount - expectedNextCount, 2) / expectedNextCount);
            // We now negate this count if the proportions bias towards the next count
            double observedProportion = (previousCount / (double) totalObservedCount);
            if (observedProportion > expectedProportion) {
                chisquare = 0 - chisquare;
            }
            // System.err.println("Raw counts "+previousCount+","+nextCount+" Totals "+previousTotalCount+","+nextTotalCount+" Expected "+expectedPreviousCount+","+expectedNextCount+" Chi "+chisquare);
            // TODO: This is icky since the inheritance between HiCDataStore and DataStore
            // isn't properly sorted out.
            ((DataStore) data[d]).setValueForProbe(probes[p], (float) chisquare);
        }
    }
    quantitatonComplete();
}

Example 12 with HiCDataStore

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

the class ChromosomeDataTrack method paint.

/* (non-Javadoc)
	 * @see javax.swing.JComponent#paint(java.awt.Graphics)
	 */
public void paint(Graphics g) {
    super.paint(g);
    if (useAntiAliasing && g instanceof Graphics2D) {
        ((Graphics2D) g).setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
    }
    int displayMode = DisplayPreferences.getInstance().getDisplayMode();
    int scaleMode = DisplayPreferences.getInstance().getScaleType();
    dataColours = DisplayPreferences.getInstance().getColourType();
    thisReadDensity = DisplayPreferences.getInstance().getReadDensity();
    drawProbes = true;
    drawReads = true;
    showNegative = true;
    if (displayMode == DisplayPreferences.DISPLAY_MODE_QUANTITATION_ONLY) {
        drawReads = false;
    }
    if (displayMode == DisplayPreferences.DISPLAY_MODE_READS_ONLY) {
        drawProbes = false;
    }
    if (scaleMode == DisplayPreferences.SCALE_TYPE_POSITIVE) {
        showNegative = false;
    }
    if (thisReadDensity == DisplayPreferences.READ_DENSITY_LOW) {
        readHeight = 3;
        readSpace = 2;
    } else if (thisReadDensity == DisplayPreferences.READ_DENSITY_MEDIUM) {
        readHeight = 2;
        readSpace = 1;
    } else if (thisReadDensity == DisplayPreferences.READ_DENSITY_HIGH) {
        readHeight = 1;
        readSpace = 1;
    }
    if (drawReads) {
        assignSlots();
    }
    // Cache a copy of the current reads array so we don't get it changed from under us
    long[] reads = this.reads;
    drawnReads.removeAllElements();
    drawnProbes.removeAllElements();
    height = getHeight();
    maxValue = (float) DisplayPreferences.getInstance().getMaxDataValue();
    if (showNegative) {
        minValue = 0 - maxValue;
    } else {
        minValue = 0;
    }
    width = getWidth();
    // between tracks.
    if (viewer.getIndex(this) % 2 == 0) {
        g.setColor(ColourScheme.DATA_BACKGROUND_EVEN);
    } else {
        g.setColor(ColourScheme.DATA_BACKGROUND_ODD);
    }
    g.fillRect(x, 0, width, height);
    if (viewer.makingSelection()) {
        int selStart = viewer.selectionStart();
        int selEnd = viewer.selectionEnd();
        int useStart = (selEnd > selStart) ? selStart : selEnd;
        int selWidth = selEnd - selStart;
        if (selWidth < 0)
            selWidth = 0 - selWidth;
        g.setColor(ColourScheme.DRAGGED_SELECTION);
        g.fillRect(useStart, 0, selWidth, height);
    }
    int startBp = viewer.currentStart();
    int endBp = viewer.currentEnd();
    if (drawProbes) {
        if (data.isQuantitated()) {
            // the origin
            if (showNegative) {
                g.setColor(Color.LIGHT_GRAY);
                if (drawReads) {
                    g.drawLine(x, (height * 3) / 4, x + width, (height * 3) / 4);
                } else {
                    g.drawLine(x, height / 2, x + width, height / 2);
                }
            }
            // Now go through all the probes figuring out whether they
            // need to be displayed
            // Reset the values used to optimise drawing
            lastProbeXEnd = 0;
            lastProbeXMid = 0;
            lastProbeY = 0;
            lastProbeValue = 0;
            for (int i = 0; i < probes.length; i++) {
                // System.out.println("Looking at read "+i+" Start: "+reads[i].start()+" End: "+reads[i].end()+" Global Start:"+startBp+" End:"+endBp);
                if (probes[i].end() > startBp && probes[i].start() < endBp) {
                    drawProbe(probes[i], g);
                }
            }
            // Always draw the active probe last so it stays on top, unless we're doing a line graph
            if (activeProbe != null && DisplayPreferences.getInstance().getGraphType() != DisplayPreferences.GRAPH_TYPE_LINE) {
                drawProbe(activeProbe, g);
            }
        } else {
            // This DataStore isn't quantitated
            g.setColor(Color.GRAY);
            if (displayMode == DisplayPreferences.DISPLAY_MODE_QUANTITATION_ONLY) {
                g.drawString("DataStore Not Quantitated", x + ((width / 2) - 20), (getHeight() / 2) - 2);
            } else {
                g.drawString("DataStore Not Quantitated", x + ((width / 2) - 20), ((getHeight() * 3) / 4) - 2);
            }
        }
    }
    if (drawReads && reads.length > 0) {
        // Wipe the cached end values
        for (int i = 0; i < lastReadXEnds.length; i++) {
            lastReadXEnds[i] = 0;
        }
        // Now go through all the reads figuring out whether they
        // need to be displayed
        // Find where we need to start
        // System.err.println("Starting looking at "+lastReadIndexStart+" at pos "+SequenceRead.start(reads[lastReadIndexStart])+" startBp is "+startBp);
        // Go back until we're sure we're not going to see any more reads
        boolean hitLimit = false;
        if (lastReadIndexStart > reads.length - 1) {
            lastReadIndexStart = reads.length - 1;
        }
        for (int i = lastReadIndexStart; i >= 0; i--) {
            if (SequenceRead.start(reads[i]) < startBp - data.getMaxReadLength()) {
                lastReadIndexStart = i;
                hitLimit = true;
                break;
            }
        }
        if (!hitLimit) {
            // System.err.println("Didn't hit the limit so bailing out");
            lastReadIndexStart = 0;
        }
        // Go forward until we hit our first read
        for (int i = lastReadIndexStart; i < reads.length; i++) {
            if (SequenceRead.start(reads[i]) >= startBp - data.getMaxReadLength()) {
                lastReadIndexStart = i;
                break;
            }
        }
        for (int i = lastReadIndexStart; i < reads.length; i++) {
            // System.out.println("Looking at read "+i+" Start: "+reads[i].start()+" End: "+reads[i].end()+" Global Start:"+startBp+" End:"+endBp);
            if (SequenceRead.end(reads[i]) > startBp && SequenceRead.start(reads[i]) < endBp) {
                drawRead(reads[i], i, g);
            }
            if (SequenceRead.start(reads[i]) > endBp)
                break;
        }
        // Always draw the active read last
        if (activeRead != 0) {
            drawRead(activeRead, activeReadIndex, g);
        }
    }
    if (drawReads && reads.length > 0 && !drawProbes && isHiC) {
        int[][] hicPixelCounts = getHiCPixelCounts();
        boolean foundSomething = false;
        boolean drawnSomthing = false;
        int drawnCount = 0;
        INTERACTIONS: for (int length = width; length >= 1; length--) {
            for (int x1 = 0; x1 < width - length; x1++) {
                int x2 = x1 + length;
                if (hicPixelCounts[x1][x2] != 0)
                    foundSomething = true;
                if (hicPixelCounts[x1][x2] > maxValue) {
                    ++drawnCount;
                    // anything else on screen and everything gets super slow.
                    if (drawnCount >= 10000)
                        break INTERACTIONS;
                    // Draw this interaction.
                    drawnSomthing = true;
                    // Work out the height of the crossover.  This is proportional to the distance between the pairs
                    // and the total width of the window.
                    double proportion = (x2 - x1) / (double) width;
                    int interactionHeight = (int) ((height / 2) * proportion);
                    // Scale the colours by distance so you can tell where things go
                    g.setColor(DisplayPreferences.getInstance().getGradient().getColor((height / 2) - interactionHeight, 0, height / 2));
                    // Draw the interaction
                    int midPoint = x1 + ((x2 - x1) / 2);
                    g.drawLine(x1, height / 2, midPoint, (height / 2) + interactionHeight);
                    g.drawLine(midPoint, (height / 2) + interactionHeight, x2, height / 2);
                }
            }
        }
        if (foundSomething && !drawnSomthing) {
            g.setColor(Color.GRAY);
            g.drawString("Lower the data zoom level to see HiC interactions", x + ((width / 2) - 100), ((getHeight() * 3) / 4) - 2);
        }
        if (drawnCount >= 10000 && maxValue < 2) {
            g.setColor(Color.DARK_GRAY);
            g.drawString("Increase the data zoom level to filter HiC interactions", x + ((width / 2) - 100), ((getHeight() * 3) / 4) - 2);
        }
    }
    // System.out.println("Drew "+drawnReads.size()+" reads");
    // Draw a line across the bottom of the display
    g.setColor(Color.LIGHT_GRAY);
    g.drawLine(x, height - 1, x + width, height - 1);
    // so catch this
    try {
        if (viewer.application().dataCollection().getActiveDataStore() == data || (enclosingSet != null && viewer.application().dataCollection().getActiveDataStore() == enclosingSet)) {
            g.setColor(Color.RED);
            g.drawLine(x, height - 2, x + width, height - 2);
            g.drawLine(x, height - 1, x + width, height - 1);
            g.drawLine(x, 0, x + width, 0);
            g.drawLine(x, 1, x + width, 1);
        }
    } catch (NullPointerException npe) {
    }
    String name = data.name();
    if (data instanceof HiCDataStore && ((HiCDataStore) data).isValidHiC()) {
        name = "[HiC] " + name;
    }
    if (enclosingSet != null) {
        name = "[" + enclosingSet.name() + "] " + name;
    }
    int fontSize = 12;
    int nameHeight = height;
    while (fontSize > 1) {
        g.setFont(new Font("Default", Font.PLAIN, fontSize));
        nameHeight = g.getFontMetrics().getAscent() + g.getFontMetrics().getDescent();
        if (nameHeight < getHeight())
            break;
        fontSize--;
    }
    // Draw a box into which we'll put the track name so it's not obscured
    // by the data
    int nameWidth = g.getFontMetrics().stringWidth(name);
    if (viewer.getIndex(this) % 2 == 0) {
        g.setColor(ColourScheme.DATA_BACKGROUND_EVEN);
    } else {
        g.setColor(ColourScheme.DATA_BACKGROUND_ODD);
    }
    g.fillRect(x, (height / 2) - (nameHeight / 2), nameWidth + 5, nameHeight);
    // Finally draw the name of the data track
    g.setColor(Color.GRAY);
    g.drawString(name, x + 2, (height / 2) + ((nameHeight / 2) - g.getFontMetrics().getDescent()));
}