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();
}
}
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;
}
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);
}
}
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;
}
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);
}
Aggregations