Examples with Pileup com.github.lindenb.jvarkit.samtools.util.Pileup used on opensource projects

Example 1 with Pileup

use of com.github.lindenb.jvarkit.samtools.util.Pileup in project jvarkit by lindenb.

the class CoverageServer method printRaster.

/**
 * print BAM for small interval, displaying reads
 */
private void printRaster(final BamInput bam, final SimpleInterval midRegion, final SimpleInterval region, final HttpServletRequest request, final HttpServletResponse response) throws IOException, ServletException {
    final IntToDoubleFunction position2pixel = X -> ((X - region.getStart()) / (double) region.getLengthOnReference()) * (double) image_width;
    final SamReaderFactory srf = SamReaderFactory.make().validationStringency(ValidationStringency.LENIENT).referenceSequence(this.faidxRef);
    final Pileup<SAMRecord> pileup = new Pileup<>((L, R) -> position2pixel.applyAsDouble(L.getUnclippedEnd() + 1) + 1 < position2pixel.applyAsDouble(R.getUnclippedStart()));
    try (SamReader sr = srf.open(bam.bamPath)) {
        try (CloseableIterator<SAMRecord> iter = sr.query(region.getContig(), // extend to get clipR
        Math.max(0, region.getStart() - this.small_region_size), region.getEnd() + this.small_region_size, false)) {
            while (iter.hasNext()) {
                final SAMRecord rec = iter.next();
                if (!acceptRead(rec))
                    continue;
                if (rec.getUnclippedEnd() < region.getStart())
                    continue;
                if (rec.getUnclippedStart() > region.getEnd())
                    continue;
                final Cigar cigar = rec.getCigar();
                if (cigar == null || cigar.isEmpty())
                    continue;
                pileup.add(rec);
            }
        }
    // end iterator
    }
    // end samreader
    ReferenceSequence refInInterval = null;
    try (ReferenceSequenceFile refseq = ReferenceSequenceFileFactory.getReferenceSequenceFile(this.faidxRef)) {
        final SAMSequenceRecord ssr = this.dictionary.getSequence(region.getContig());
        if (region.getStart() <= ssr.getSequenceLength()) {
            refInInterval = refseq.getSubsequenceAt(region.getContig(), region.getStart(), Math.min(region.getEnd(), ssr.getSequenceLength()));
        }
    }
    final BufferedImage img = new BufferedImage(image_width, image_height, BufferedImage.TYPE_INT_RGB);
    final Graphics2D g = img.createGraphics();
    g.setRenderingHint(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);
    g.setColor(new Color(240, 240, 240));
    g.fillRect(0, 0, image_width + 1, image_height + 1);
    final Stroke oldStroke = g.getStroke();
    g.setStroke(new BasicStroke(0.5f));
    g.setColor(Color.WHITE);
    final double mid_start = position2pixel.applyAsDouble(midRegion.getStart());
    final double mid_end = position2pixel.applyAsDouble(midRegion.getEnd() + 1);
    g.fill(new Rectangle2D.Double(mid_start, 0, (mid_end - mid_start), this.image_height));
    final int[] int_coverage = new int[region.getLengthOnReference()];
    final int margin_top = 12;
    final double featureHeight = Math.min(20, (this.image_height - margin_top) / Math.max(1.0, (double) pileup.getRowCount()));
    double y = image_height - featureHeight;
    for (final List<SAMRecord> row : pileup) {
        final double h2 = Math.min(featureHeight * 0.9, featureHeight - 2);
        for (final SAMRecord rec : row) {
            final Cigar cigar = rec.getCigar();
            if (cigar == null || cigar.isEmpty())
                continue;
            /* fill coverage array  */
            int ref1 = rec.getAlignmentStart();
            for (CigarElement ce : cigar.getCigarElements()) {
                if (ref1 > region.getEnd())
                    break;
                final CigarOperator op = ce.getOperator();
                if (op.consumesReferenceBases()) {
                    if (op.consumesReadBases()) {
                        for (int x = 0; x < ce.getLength(); ++x) {
                            int pos = ref1 + x;
                            if (pos < region.getStart())
                                continue;
                            if (pos > region.getEnd())
                                break;
                            int_coverage[pos - region.getStart()]++;
                        }
                    }
                    ref1 += ce.getLength();
                }
            }
            /* draw rec itself */
            final double midy = y + h2 / 2.0;
            g.setColor(Color.DARK_GRAY);
            g.draw(new Line2D.Double(position2pixel.applyAsDouble(rec.getUnclippedStart()), midy, position2pixel.applyAsDouble(rec.getUnclippedEnd()), midy));
            ref1 = rec.getUnclippedStart();
            final List<Double> insertions = new ArrayList<>();
            for (final CigarElement ce : cigar.getCigarElements()) {
                if (ref1 > region.getEnd())
                    break;
                final CigarOperator op = ce.getOperator();
                Shape shape = null;
                Color fill = null;
                Color stroke = Color.DARK_GRAY;
                switch(op) {
                    case P:
                        break;
                    // through
                    case M:
                    // through
                    case X:
                    // through
                    case EQ:
                    // through
                    case S:
                    case H:
                        final double x1 = position2pixel.applyAsDouble(ref1);
                        shape = new Rectangle2D.Double(x1, y, position2pixel.applyAsDouble(ref1 + ce.getLength()) - x1, h2);
                        ref1 += ce.getLength();
                        switch(op) {
                            case H:
                            case S:
                                fill = Color.YELLOW;
                                break;
                            case X:
                                fill = Color.RED;
                                break;
                            case EQ:
                            case M:
                                fill = Color.LIGHT_GRAY;
                                break;
                            default:
                                break;
                        }
                        break;
                    // through
                    case N:
                    case D:
                        shape = null;
                        fill = null;
                        stroke = null;
                        ref1 += ce.getLength();
                        break;
                    case I:
                        shape = null;
                        fill = null;
                        stroke = null;
                        insertions.add(position2pixel.applyAsDouble(ref1));
                        break;
                    default:
                        throw new IllegalStateException("" + op);
                }
                if (ref1 < region.getStart())
                    continue;
                if (shape != null) {
                    if (fill != null) {
                        g.setColor(fill);
                        g.fill(shape);
                    }
                    if (stroke != null && h2 > 4) {
                        g.setColor(stroke);
                        g.draw(shape);
                    }
                }
            }
            /* draw mismatched bases */
            if (refInInterval != null && rec.getReadBases() != null && rec.getReadBases() != SAMRecord.NULL_SEQUENCE) {
                final byte[] bases = rec.getReadBases();
                final IntFunction<Character> baseRead = IDX -> IDX < 0 || IDX >= bases.length || bases == SAMRecord.NULL_SEQUENCE ? 'N' : (char) Character.toUpperCase(bases[IDX]);
                int read0 = 0;
                ref1 = rec.getAlignmentStart();
                for (CigarElement ce : cigar.getCigarElements()) {
                    if (ref1 > region.getEnd())
                        break;
                    final CigarOperator op = ce.getOperator();
                    switch(op) {
                        case P:
                            break;
                        case H:
                            break;
                        case D:
                        case N:
                            ref1 += ce.getLength();
                            break;
                        case S:
                        case I:
                            read0 += ce.getLength();
                            break;
                        case EQ:
                        case M:
                        case X:
                            {
                                for (int j = 0; j < ce.getLength(); j++) {
                                    if (ref1 + j < region.getStart())
                                        continue;
                                    if (ref1 + j >= region.getStart() + refInInterval.length())
                                        break;
                                    final int ref_base_idx = ref1 - region.getStart() + j;
                                    char ctgBase = (char) (ref_base_idx < 0 || ref_base_idx >= refInInterval.length() ? 'N' : Character.toUpperCase(refInInterval.getBases()[ref_base_idx]));
                                    if (ctgBase == 'N')
                                        continue;
                                    char readBase = baseRead.apply(read0 + j);
                                    if (readBase == 'N')
                                        continue;
                                    if (readBase == ctgBase)
                                        continue;
                                    g.setColor(Color.ORANGE);
                                    final double x1 = position2pixel.applyAsDouble(ref1 + j);
                                    final double x2 = position2pixel.applyAsDouble(ref1 + j + 1);
                                    g.fill(new Rectangle2D.Double(x1, y, x2 - x1, h2));
                                }
                                read0 += ce.getLength();
                                ref1 += ce.getLength();
                                break;
                            }
                        default:
                            break;
                    }
                }
            }
            for (double px : insertions) {
                g.setColor(Color.RED);
                g.draw(new Line2D.Double(px, y - 0.5, px, y + h2 + 0.5));
            }
        }
        y -= featureHeight;
    }
    double max_cov = IntStream.of(int_coverage).max().orElse(0);
    g.setColor(Color.DARK_GRAY);
    g.drawString("Sample:" + bam.sample + " max-cov:" + (int) max_cov + " " + region.toNiceString(), 10, 10);
    /* plot coverage */
    final GeneralPath gp = new GeneralPath();
    for (int i = 0; max_cov > 0 && i < int_coverage.length; ++i) {
        final double x1 = position2pixel.applyAsDouble(region.getStart() + i);
        final double x2 = position2pixel.applyAsDouble(region.getStart() + i + 1);
        final double y1 = image_height - (int_coverage[i] / max_cov) * (image_height - margin_top);
        if (i == 0)
            gp.moveTo(x1, y1);
        else
            gp.lineTo(x1, y1);
        gp.lineTo(x2, y1);
    }
    g.setStroke(new BasicStroke(0.5f, BasicStroke.CAP_BUTT, BasicStroke.JOIN_BEVEL, 0f, new float[] { 1f, 2f, 1f }, 0f));
    g.setColor(Color.BLUE);
    g.draw(gp);
    g.setStroke(oldStroke);
    writeGenes(g, region);
    writeKnownCnv(g, region);
    g.setColor(Color.PINK);
    g.draw(new Line2D.Double(mid_start, 0, mid_start, image_height));
    g.draw(new Line2D.Double(mid_end, 0, mid_end, image_height));
    writeImage(img, bam, region, response);
}

Example 2 with Pileup

use of com.github.lindenb.jvarkit.samtools.util.Pileup in project jvarkit by lindenb.

the class CoveragePlotter method drawKnownCnv.

private void drawKnownCnv(final Graphics2D g, final Rectangle rectangle, final Locatable region) {
    if (this.knownCnvFile == null)
        return;
    final String fname = this.knownCnvFile.getFileName().toString();
    final Pileup<Interval> pileup = new Pileup<>();
    final Predicate<Interval> rejectCnv = cnv -> (this.ignore_cnv_overlapping && cnv.getStart() < region.getStart() && cnv.getEnd() > region.getEnd());
    if (fname.endsWith(".bed.gz")) {
        try (TabixReader tbr = new TabixReader(this.knownCnvFile.toString())) {
            final ContigNameConverter cvt = ContigNameConverter.fromContigSet(tbr.getChromosomes());
            final String ctg = cvt.apply(region.getContig());
            if (!StringUtils.isBlank(ctg)) {
                final BedLineCodec codec = new BedLineCodec();
                final TabixReader.Iterator iter = tbr.query(ctg, region.getStart(), region.getEnd());
                for (; ; ) {
                    final String line = iter.next();
                    if (line == null)
                        break;
                    final BedLine bed = codec.decode(line);
                    if (bed == null)
                        continue;
                    final Interval rgn = new Interval(region.getContig(), bed.getStart(), bed.getEnd(), false, bed.getOrDefault(3, ""));
                    if (rejectCnv.test(rgn))
                        return;
                    pileup.add(rgn);
                }
            }
        } catch (final Throwable err) {
            LOG.error(err);
        }
    } else if (FileExtensions.VCF_LIST.stream().anyMatch(X -> fname.endsWith(X))) {
        try (VCFReader vcfFileReader = VCFReaderFactory.makeDefault().open(this.knownCnvFile, true)) {
            final ContigNameConverter cvt = ContigNameConverter.fromOneDictionary(SequenceDictionaryUtils.extractRequired(vcfFileReader.getHeader()));
            final String ctg = cvt.apply(region.getContig());
            if (!StringUtils.isBlank(ctg)) {
                vcfFileReader.query(ctg, region.getStart(), region.getEnd()).stream().filter(VC -> !VC.isSNP()).forEach(VC -> {
                    final List<String> list = new ArrayList<>();
                    if (VC.hasID())
                        list.add(VC.getID());
                    if (VC.hasAttribute(VCFConstants.SVTYPE))
                        list.add(VC.getAttributeAsString(VCFConstants.SVTYPE, "."));
                    final Interval rgn = new Interval(region.getContig(), VC.getStart(), VC.getEnd(), false, String.join(";", list));
                    if (rejectCnv.test(rgn))
                        return;
                    pileup.add(rgn);
                });
            }
        } catch (final Throwable err) {
            LOG.error(err);
        }
    } else {
        LOG.warn("not a vcf of bed.gz file " + this.knownCnvFile);
    }
    if (!pileup.isEmpty()) {
        final Composite oldComposite = g.getComposite();
        final Stroke oldStroke = g.getStroke();
        g.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.4f));
        g.setColor(getColor("known", Color.MAGENTA));
        final IntToDoubleFunction position2pixel = X -> ((X - region.getStart()) / (double) region.getLengthOnReference()) * rectangle.getWidth();
        final double featureHeight = 4.0 / pileup.getRowCount();
        for (int row = 0; row < pileup.getRowCount(); ++row) {
            for (final Interval cnv : pileup.getRow(row)) {
                final double y = rectangle.getHeight() - 8.0 + row * featureHeight;
                final double x1 = position2pixel.applyAsDouble(cnv.getStart());
                final double x2 = position2pixel.applyAsDouble(cnv.getEnd());
                g.draw(new Rectangle2D.Double(x1, y - 1, Math.max(1.0, x2 - x1), featureHeight * 0.9));
            }
        }
        g.setComposite(oldComposite);
        g.setStroke(oldStroke);
    }
}

Example 3 with Pileup

use of com.github.lindenb.jvarkit.samtools.util.Pileup in project jvarkit by lindenb.

the class BamToSVG method plotSVG.

private void plotSVG(final ArchiveFactory archive, final Path path) throws IOException, XMLStreamException {
    XMLStreamWriter w = null;
    final SamReaderFactory sfrf = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.LENIENT).referenceSequence(this.referenceFile);
    final Context context = new Context();
    try (SamReader samReader = sfrf.open(path)) {
        final SAMFileHeader header = samReader.getFileHeader();
        final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(header);
        SequenceUtil.assertSequenceDictionariesEqual(dict, this.referenceDict);
        context.sampleName = header.getReadGroups().stream().map(RG -> RG.getSample()).filter(S -> !StringUtils.isBlank(S)).findFirst().orElse(IOUtils.getFilenameWithoutCommonSuffixes(path));
        // get a chance to get clipping close to bounds
        final int extend_for_clip = 100;
        // read SamRecord
        try (CloseableIterator<SAMRecord> iter = samReader.query(this.interval.getContig(), Math.max(1, this.interval.getStart() - extend_for_clip), this.interval.getEnd() + extend_for_clip, false)) {
            final Pileup<SAMRecord> pileup = new Pileup<>((A, B) -> !CoordMath.overlaps(left(A) - 1, right(A) + 1, left(B), right(B)));
            while (iter.hasNext()) {
                final SAMRecord rec = iter.next();
                if (rec.getReadUnmappedFlag())
                    continue;
                // if( this.samRecordFilter.filterOut(rec)) continue;
                if (!SAMRecordDefaultFilter.accept(rec, this.mappingQuality))
                    continue;
                if (!rec.getReferenceName().equals(this.interval.getContig()))
                    continue;
                if (right(rec) < this.interval.getStart())
                    continue;
                if (left(rec) > this.interval.getEnd())
                    continue;
                pileup.add(rec);
            }
            context.lines = pileup.getRows();
        }
        if (this.vcf != null) {
            context.readVariantFile(vcf);
        }
        context.featureWidth = this.drawinAreaWidth / (double) ((this.interval.getEnd() - this.interval.getStart()) + 1);
        context.featureHeight = Math.min(Math.max(5.0, context.featureWidth), 30);
        context.HEIGHT_RULER = (int) (StringUtils.niceInt(this.interval.getEnd()).length() * context.featureHeight + 5);
        LOG.info("Feature height:" + context.featureHeight);
        final String buildName = SequenceDictionaryUtils.getBuildName(this.referenceDict).orElse("");
        final String filename = this.prefix + buildName + (buildName.isEmpty() ? "" : ".") + this.interval.getContig() + "_" + this.interval.getStart() + "_" + this.interval.getEnd() + "." + context.sampleName + ".svg";
        LOG.info("writing " + filename);
        final XMLOutputFactory xof = XMLOutputFactory.newFactory();
        try (OutputStream fout = archive.openOuputStream(filename)) {
            w = xof.createXMLStreamWriter(fout, "UTF-8");
            w.writeStartDocument("UTF-8", "1.0");
            printDocument(w, context);
            w.writeEndDocument();
            w.flush();
            w.close();
            fout.flush();
        }
    }
}