Search in sources :

Example 1 with BufferedList

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 {
    }
}
Also used : Color(java.awt.Color) Cigar(htsjdk.samtools.Cigar) Program(com.github.lindenb.jvarkit.util.jcommander.Program) Parameter(com.beust.jcommander.Parameter) Rectangle2D(java.awt.geom.Rectangle2D) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) CigarElement(htsjdk.samtools.CigarElement) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) CigarOperator(htsjdk.samtools.CigarOperator) HashMap(java.util.HashMap) SAMRecordPartition(com.github.lindenb.jvarkit.util.samtools.SAMRecordPartition) Interval(htsjdk.samtools.util.Interval) Graphics2D(java.awt.Graphics2D) Map(java.util.Map) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) ImageIO(javax.imageio.ImageIO) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) BufferedList(com.github.lindenb.jvarkit.util.BufferedList) Hershey(com.github.lindenb.jvarkit.util.Hershey) CloserUtil(htsjdk.samtools.util.CloserUtil) Line2D(java.awt.geom.Line2D) BufferedImage(java.awt.image.BufferedImage) IntervalTreeMap(htsjdk.samtools.util.IntervalTreeMap) Logger(com.github.lindenb.jvarkit.util.log.Logger) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Set(java.util.Set) IOException(java.io.IOException) AffineTransform(java.awt.geom.AffineTransform) Collectors(java.util.stream.Collectors) File(java.io.File) SAMRecord(htsjdk.samtools.SAMRecord) SamRecordFilter(htsjdk.samtools.filter.SamRecordFilter) List(java.util.List) SamRecordJEXLFilter(com.github.lindenb.jvarkit.util.samtools.SamRecordJEXLFilter) BufferedReader(java.io.BufferedReader) BitSet(java.util.BitSet) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) SamReaderFactory(htsjdk.samtools.SamReaderFactory) BufferedList(com.github.lindenb.jvarkit.util.BufferedList) SamReaderFactory(htsjdk.samtools.SamReaderFactory) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) Line2D(java.awt.geom.Line2D) BufferedImage(java.awt.image.BufferedImage) SamReaderFactory(htsjdk.samtools.SamReaderFactory) SAMSequenceDictionaryProgress(com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress) Hershey(com.github.lindenb.jvarkit.util.Hershey) Rectangle2D(java.awt.geom.Rectangle2D) CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement) IOException(java.io.IOException) Graphics2D(java.awt.Graphics2D) Cigar(htsjdk.samtools.Cigar) SAMRecord(htsjdk.samtools.SAMRecord) AffineTransform(java.awt.geom.AffineTransform) Interval(htsjdk.samtools.util.Interval)

Aggregations

Parameter (com.beust.jcommander.Parameter)1 IOUtils (com.github.lindenb.jvarkit.io.IOUtils)1 BufferedList (com.github.lindenb.jvarkit.util.BufferedList)1 Hershey (com.github.lindenb.jvarkit.util.Hershey)1 BedLine (com.github.lindenb.jvarkit.util.bio.bed.BedLine)1 BedLineCodec (com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec)1 Launcher (com.github.lindenb.jvarkit.util.jcommander.Launcher)1 Program (com.github.lindenb.jvarkit.util.jcommander.Program)1 Logger (com.github.lindenb.jvarkit.util.log.Logger)1 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)1 SAMRecordPartition (com.github.lindenb.jvarkit.util.samtools.SAMRecordPartition)1 SamRecordJEXLFilter (com.github.lindenb.jvarkit.util.samtools.SamRecordJEXLFilter)1 Cigar (htsjdk.samtools.Cigar)1 CigarElement (htsjdk.samtools.CigarElement)1 CigarOperator (htsjdk.samtools.CigarOperator)1 SAMRecord (htsjdk.samtools.SAMRecord)1 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)1 SamReaderFactory (htsjdk.samtools.SamReaderFactory)1 SamRecordFilter (htsjdk.samtools.filter.SamRecordFilter)1 CloserUtil (htsjdk.samtools.util.CloserUtil)1