use of com.github.lindenb.jvarkit.util.BufferedList in project jvarkit by lindenb.
the class BamCmpCoverage method doWork.
@Override
public int doWork(final List<String> args) {
if (outputFile == null) {
LOG.error("output image file not defined");
return -1;
}
if (this.imgageSize < 1) {
LOG.error("Bad image size:" + this.imgageSize);
return -1;
}
if (this.minDepth < 0) {
LOG.error("Bad min depth : " + this.minDepth);
return -1;
}
if (this.minDepth >= this.maxDepth) {
LOG.error("Bad min<max depth : " + this.minDepth + "<" + this.maxDepth);
return 1;
}
if (this.bedFile != null) {
readBedFile(this.bedFile);
}
if (regionStr != null && this.intervals != null) {
LOG.error("bed and interval both defined.");
return -1;
}
try {
final ConcatSam.Factory concatSamFactory = new ConcatSam.Factory();
final SamReaderFactory srf = concatSamFactory.getSamReaderFactory();
srf.disable(SamReaderFactory.Option.EAGERLY_DECODE);
srf.disable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS);
srf.disable(SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS);
if (this.regionStr != null) {
concatSamFactory.addInterval(this.regionStr);
}
ConcatSam.ConcatSamIterator concatIter = concatSamFactory.open(args);
final SAMSequenceDictionary dict = concatIter.getFileHeader().getSequenceDictionary();
final Set<String> samples = concatIter.getFileHeader().getReadGroups().stream().map(RG -> this.samRecordPartition.apply(RG, "N/A")).collect(Collectors.toSet());
LOG.info("Samples:" + samples.size());
for (String sample : samples) {
this.sample2column.put(sample, this.sample2column.size());
}
// create image
LOG.info("Creating image " + this.imgageSize + "x" + this.imgageSize);
this.image = new BufferedImage(this.imgageSize, this.imgageSize, BufferedImage.TYPE_INT_RGB);
Graphics2D g = this.image.createGraphics();
this.marginWidth = this.imgageSize * 0.05;
double drawingWidth = (this.imgageSize - 1) - marginWidth;
this.sampleWidth = drawingWidth / samples.size();
// g.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
g.setColor(Color.WHITE);
g.fillRect(0, 0, this.imgageSize, this.imgageSize);
g.setColor(Color.BLACK);
Hershey hershey = new Hershey();
for (final String sample_x : samples) {
double labelHeight = marginWidth;
if (labelHeight > 50)
labelHeight = 50;
g.setColor(Color.BLACK);
hershey.paint(g, sample_x, marginWidth + sample2column.get(sample_x) * sampleWidth, marginWidth - labelHeight, sampleWidth * 0.9, labelHeight * 0.9);
AffineTransform old = g.getTransform();
AffineTransform tr = AffineTransform.getTranslateInstance(marginWidth, marginWidth + sample2column.get(sample_x) * sampleWidth);
tr.rotate(Math.PI / 2);
g.setTransform(tr);
hershey.paint(g, sample_x, 0.0, 0.0, sampleWidth * 0.9, labelHeight * 0.9);
// g.drawString(this.tabixFile.getFile().getName(),0,0);
g.setTransform(old);
for (String sample_y : samples) {
Rectangle2D rect = new Rectangle2D.Double(marginWidth + sample2column.get(sample_x) * sampleWidth, marginWidth + sample2column.get(sample_y) * sampleWidth, sampleWidth, sampleWidth);
g.setColor(Color.BLUE);
g.draw(new Line2D.Double(rect.getMinX(), rect.getMinY(), rect.getMaxX(), rect.getMaxY()));
g.setColor(Color.BLACK);
g.draw(rect);
}
}
// ceate bit-array
BitSampleMatrix bitMatrix = new BitSampleMatrix(samples.size());
// preivous chrom
// int prev_tid=-1;
BufferedList<Depth> depthList = new BufferedList<Depth>();
g.setColor(Color.BLACK);
SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(dict).logger(LOG);
LOG.info("Scanning bams...");
while (concatIter.hasNext()) {
final SAMRecord rec = progress.watch(concatIter.next());
if (this.samRecordFilter.filterOut(rec))
continue;
final String sample = this.samRecordPartition.getPartion(rec, "N/A");
final int sample_id = this.sample2column.get(sample);
final Cigar cigar = rec.getCigar();
if (cigar == null || cigar.isEmpty())
continue;
int refPos = rec.getAlignmentStart();
/* cleanup front pos */
while (!depthList.isEmpty()) {
final Depth front = depthList.getFirst();
if (front.tid != rec.getReferenceIndex().intValue() || front.pos < refPos) {
paint(bitMatrix, front);
depthList.removeFirst();
continue;
} else {
break;
}
}
for (final CigarElement ce : cigar.getCigarElements()) {
final CigarOperator op = ce.getOperator();
if (!op.consumesReferenceBases())
continue;
if (op.consumesReadBases()) {
for (int i = 0; i < ce.getLength(); ++i) {
Depth depth = null;
int pos = refPos + i;
// ignore non-overlapping BED
if (this.intervals != null && !this.intervals.containsOverlapping(new Interval(rec.getReferenceName(), pos, pos))) {
continue;
} else if (depthList.isEmpty()) {
depth = new Depth();
depth.pos = pos;
depth.tid = rec.getReferenceIndex();
depthList.add(depth);
} else if (depthList.getLast().pos < pos) {
Depth prev = depthList.getLast();
while (prev.pos < pos) {
depth = new Depth();
depth.pos = prev.pos + 1;
depth.tid = rec.getReferenceIndex();
depthList.add(depth);
prev = depth;
}
depth = prev;
} else {
int lastPos = depthList.get(depthList.size() - 1).pos;
int distance = lastPos - pos;
int indexInList = (depthList.size() - 1) - (distance);
if (indexInList < 0) {
// can appen when BED declared and partially overlap the read
continue;
}
depth = depthList.get((depthList.size() - 1) - (distance));
if (depth.pos != pos) {
LOG.error(" " + pos + " vs " + depth.pos + " " + lastPos);
return -1;
}
}
depth.depths[sample_id]++;
}
}
refPos += ce.getLength();
}
}
while (!depthList.isEmpty()) {
// paint(g,depthList.remove(0));
paint(bitMatrix, depthList.remove(0));
}
progress.finish();
concatIter.close();
concatIter = null;
for (int x = 0; x < bitMatrix.n_samples; ++x) {
for (int y = 0; y < bitMatrix.n_samples; ++y) {
LOG.info("Painting...(" + x + "/" + y + ")");
paint(g, bitMatrix.get(x, y));
}
}
g.dispose();
// save file
LOG.info("saving " + this.outputFile);
if (this.outputFile.getName().toLowerCase().endsWith(".png")) {
ImageIO.write(this.image, "PNG", this.outputFile);
} else {
ImageIO.write(this.image, "JPG", this.outputFile);
}
return 0;
} catch (final Exception err) {
LOG.error(err);
return -1;
} finally {
}
}
Aggregations