Examples with CloseableIterator htsjdk.samtools.util.CloseableIterator used on opensource projects

Example 26 with CloseableIterator

use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.

the class VcfFilterByLiftOver method doVcfToVcf.

@Override
protected int doVcfToVcf(final String inputName, final VCFIterator in, final VariantContextWriter out) {
    final ContigNameConverter anotherVcfCtgConverter;
    if (this.anotherVcfReader != null) {
        final SAMSequenceDictionary dict = this.anotherVcfReader.getHeader().getSequenceDictionary();
        if (dict != null) {
            anotherVcfCtgConverter = ContigNameConverter.fromOneDictionary(dict);
        } else {
            anotherVcfCtgConverter = ContigNameConverter.getIdentity();
        }
    } else {
        anotherVcfCtgConverter = ContigNameConverter.getIdentity();
    }
    final LiftOver liftOver = new LiftOver(this.liftOverFile);
    liftOver.setLiftOverMinMatch(this.userMinMatch);
    final VCFHeader header = in.getHeader();
    final SAMSequenceDictionary dict = header.getSequenceDictionary();
    if (!this.disableValidation && dict != null && !dict.isEmpty()) {
        liftOver.validateToSequences(dict);
    }
    final VCFHeader header2 = new VCFHeader(header);
    final VCFFilterHeaderLine filterLiftOverFailed = new VCFFilterHeaderLine("LIFTOVER_FAILED", "liftover failed " + this.liftOverFile);
    header2.addMetaDataLine(filterLiftOverFailed);
    final VCFFilterHeaderLine filterNoSameContig = new VCFFilterHeaderLine("LIFTOVER_OTHER_CTG", "Variant is mapped to another contig after liftover with " + this.liftOverFile);
    header2.addMetaDataLine(filterNoSameContig);
    final VCFInfoHeaderLine infoLiftOverPos = new VCFInfoHeaderLine("LIFTOVER_LOC", 1, VCFHeaderLineType.String, "Position of the variant liftover-ed with " + this.liftOverFile);
    header2.addMetaDataLine(infoLiftOverPos);
    final VCFFilterHeaderLine filterDistantFromPrev = new VCFFilterHeaderLine("LIFTOVER_DISTANT", "After liftover the distance (< " + min_distance + ") with the previous variant is unusual( > " + max_distance + ") after liftover with " + this.liftOverFile);
    header2.addMetaDataLine(filterDistantFromPrev);
    /* transfert filters to new header */
    if (this.anotherVcfReader != null) {
        this.anotherVcfReader.getHeader().getFilterLines().forEach(F -> header2.addMetaDataLine(F));
    }
    JVarkitVersion.getInstance().addMetaData(this, header2);
    out.writeHeader(header2);
    Locatable prevCtx = null;
    Locatable prevLifted = null;
    final Function<Locatable, String> interval2str = R -> R.getContig() + "|" + R.getStart() + "|" + R.getEnd();
    while (in.hasNext()) {
        final VariantContext ctx = in.next();
        if (prevCtx != null && !prevCtx.getContig().equals(ctx.getContig())) {
            prevCtx = null;
            prevLifted = null;
        }
        final Interval lifted = liftOver.liftOver(new Interval(ctx.getContig(), ctx.getStart(), ctx.getEnd(), // negative strand
        false, interval2str.apply(ctx)));
        // lifover failed
        if (lifted == null) {
            final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
            vcb.filter(filterLiftOverFailed.getID());
            out.add(vcb.make());
        } else // another contig
        if (lifted != null && !sameContig(lifted, ctx)) {
            final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
            vcb.filter(filterNoSameContig.getID());
            vcb.attribute(infoLiftOverPos.getID(), interval2str.apply(lifted));
            out.add(vcb.make());
        } else // strange distance
        if (prevCtx != null && lifted != null && prevLifted != null && prevCtx.getContig().equals(ctx.getContig()) && sameContig(prevLifted, lifted) && distance(prevCtx, ctx) < this.min_distance && distance(prevLifted, lifted) > this.max_distance) {
            final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
            vcb.filter(filterDistantFromPrev.getID());
            vcb.attribute(infoLiftOverPos.getID(), interval2str.apply(lifted));
            out.add(vcb.make());
        } else // filtered in anotherVcf
        if (this.anotherVcfReader != null) {
            final Set<String> found_filtered = new HashSet<>();
            final String ctg2 = anotherVcfCtgConverter.apply(lifted.getContig());
            if (!StringUtils.isBlank(ctg2)) {
                try (CloseableIterator<VariantContext> iter2 = this.anotherVcfReader.query(ctg2, lifted.getStart(), lifted.getEnd())) {
                    while (iter2.hasNext()) {
                        final VariantContext ctx2 = iter2.next();
                        if (!ctx2.isFiltered() || ctx2.getLengthOnReference() != ctx.getLengthOnReference())
                            continue;
                        found_filtered.addAll(ctx2.getFilters());
                        break;
                    }
                }
            }
            if (found_filtered.isEmpty()) {
                out.add(ctx);
            } else {
                // add previous
                found_filtered.addAll(ctx.getFilters());
                final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
                vcb.filters(found_filtered);
                out.add(vcb.make());
            }
        } else {
            out.add(ctx);
        }
        prevCtx = ctx;
        prevLifted = lifted;
    }
    return 0;
}

Example 27 with CloseableIterator

use of htsjdk.samtools.util.CloseableIterator 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 28 with CloseableIterator

use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.

the class DepthAnomaly method doWork.

@Override
public int doWork(final List<String> args) {
    PrintWriter pw = null;
    try {
        final IndexCovUtils indexCovUtils = new IndexCovUtils(this.treshold);
        final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(this.refPath);
        final SamReaderFactory samReaderFactory = SamReaderFactory.makeDefault().referenceSequence(DepthAnomaly.this.refPath).validationStringency(ValidationStringency.LENIENT);
        final List<Path> inputBams = IOUtils.unrollPaths(this.bamsPath);
        if (inputBams.isEmpty()) {
            LOG.error("input bam file missing.");
            return -1;
        }
        Iterator<? extends Locatable> iter;
        final String input = oneFileOrNull(args);
        if (input == null) {
            final BedLineCodec codec = new BedLineCodec();
            final LineIterator liter = new LineIterator(stdin());
            iter = new AbstractIterator<Locatable>() {

                @Override
                protected Locatable advance() {
                    while (liter.hasNext()) {
                        final String line = liter.next();
                        final BedLine bed = codec.decode(line);
                        if (bed == null) {
                            continue;
                        }
                        return bed;
                    }
                    liter.close();
                    return null;
                }
            };
        } else {
            iter = IntervalListProvider.from(input).dictionary(dict).stream().iterator();
        }
        pw = super.openPathOrStdoutAsPrintWriter(this.outputFile);
        while (iter.hasNext()) {
            final SimpleInterval locatable = new SimpleInterval(iter.next());
            boolean found_anomaly_here = false;
            if (this.min_anomaly_length * 3 >= locatable.getLengthOnReference()) {
                LOG.warning("interval " + locatable.toNiceString() + " is too short. skipping.");
                continue;
            }
            final int[] depth = new int[locatable.getLengthOnReference()];
            final int[] copy = new int[depth.length];
            for (final Path path : inputBams) {
                try (SamReader sr = samReaderFactory.open(path)) {
                    final SAMFileHeader header = sr.getFileHeader();
                    final String sample = header.getReadGroups().stream().map(RG -> RG.getSample()).filter(S -> !StringUtil.isBlank(S)).findFirst().orElse(IOUtils.getFilenameWithoutCommonSuffixes(path));
                    SequenceUtil.assertSequenceDictionariesEqual(dict, header.getSequenceDictionary());
                    Arrays.fill(depth, 0);
                    try (CloseableIterator<SAMRecord> siter = sr.queryOverlapping(locatable.getContig(), locatable.getStart(), locatable.getEnd())) {
                        while (siter.hasNext()) {
                            final SAMRecord rec = siter.next();
                            if (rec.getReadUnmappedFlag())
                                continue;
                            if (!SAMRecordDefaultFilter.accept(rec, this.min_mapq))
                                continue;
                            int ref = rec.getStart();
                            final Cigar cigar = rec.getCigar();
                            if (cigar == null)
                                continue;
                            for (CigarElement ce : cigar) {
                                final CigarOperator op = ce.getOperator();
                                final int len = ce.getLength();
                                if (op.consumesReferenceBases()) {
                                    if (op.consumesReadBases()) {
                                        for (int i = 0; i < len; i++) {
                                            final int pos = ref + i;
                                            if (pos < locatable.getStart())
                                                continue;
                                            if (pos > locatable.getEnd())
                                                break;
                                            depth[pos - locatable.getStart()]++;
                                        }
                                    }
                                    ref += len;
                                }
                            }
                        }
                    // loop cigar
                    }
                    // end samItere
                    System.arraycopy(depth, 0, copy, 0, depth.length);
                    final double median = median(copy);
                    final List<CovInterval> anomalies = new ArrayList<>();
                    // int minDepth = Arrays.stream(depth).min().orElse(0);
                    int x0 = 0;
                    while (x0 < depth.length && median > 0.0) {
                        final int xi = x0;
                        double total = 0;
                        double count = 0;
                        IndexCovUtils.SvType prevType = null;
                        while (x0 < depth.length) {
                            final IndexCovUtils.SvType type;
                            final int depthHere = depth[x0];
                            final double normDepth = depthHere / (median == 0 ? 1.0 : median);
                            if (depthHere > this.max_depth) {
                                type = null;
                            } else {
                                type = indexCovUtils.getType(normDepth);
                            }
                            x0++;
                            if (type == null || !type.isVariant())
                                break;
                            if (prevType != null && !type.equals(prevType))
                                break;
                            if (prevType == null)
                                prevType = type;
                            total += depthHere;
                            count++;
                        }
                        if (prevType != null && count >= this.min_anomaly_length) {
                            anomalies.add(new CovInterval(locatable.getContig(), locatable.getStart() + xi, locatable.getStart() + x0 - 1, prevType, Collections.singletonList(total / count)));
                        }
                    }
                    if (!anomalies.isEmpty() || force_screen) {
                        int i = 0;
                        while (i + 1 < anomalies.size() && this.merge_intervals) {
                            final CovInterval loc1 = anomalies.get(i);
                            final CovInterval loc2 = anomalies.get(i + 1);
                            if (loc1.svtype.equals(loc2.svtype) && loc1.withinDistanceOf(loc2, this.min_anomaly_length)) {
                                final List<Double> newdepths = new ArrayList<>(loc1.depths);
                                newdepths.addAll(loc2.depths);
                                anomalies.set(i, new CovInterval(loc1.getContig(), loc1.getStart(), loc2.getEnd(), loc1.svtype, newdepths));
                                anomalies.remove(i + 1);
                            } else {
                                i++;
                            }
                        }
                        if (!found_anomaly_here) {
                            pw.println(">>> " + locatable.toNiceString() + " length:" + StringUtils.niceInt(locatable.getLengthOnReference()));
                            found_anomaly_here = true;
                        }
                        if (screen_width > 0) {
                            pw.print("#");
                            pw.print(String.format("%-15s", sample));
                            pw.print("[");
                            for (i = 0; i < screen_width; i++) {
                                double t = 0;
                                double n = 0;
                                final int x1 = (int) (((i + 0) / (double) screen_width) * depth.length);
                                final int x2 = (int) (((i + 1) / (double) screen_width) * depth.length);
                                for (int x3 = x1; x3 <= x2 && x3 < depth.length; ++x3) {
                                    t += depth[x1];
                                    n++;
                                }
                                double normDepth = t /= n;
                                if (median > 0)
                                    normDepth /= median;
                                // centered
                                normDepth /= 2.0;
                                final boolean is_anomaly = anomalies.stream().anyMatch(R -> CoordMath.overlaps(R.getStart(), R.getEnd(), locatable.getStart() + x1, locatable.getStart() + x2));
                                final AnsiUtils.AnsiColor color = is_anomaly ? AnsiColor.RED : null;
                                if (color != null)
                                    pw.print(color.begin());
                                pw.print(AnsiUtils.getHistogram(normDepth));
                                if (color != null)
                                    pw.print(color.end());
                            }
                            pw.print("]");
                            pw.println();
                        }
                        for (i = 0; i < anomalies.size(); i++) {
                            final CovInterval anomalie = anomalies.get(i);
                            pw.print(anomalie.getContig());
                            pw.print("\t");
                            pw.print(anomalie.getStart() - 1);
                            pw.print("\t");
                            pw.print(anomalie.getEnd());
                            pw.print("\t");
                            pw.print(anomalie.getLengthOnReference());
                            pw.print("\t");
                            pw.print(anomalie.svtype.name());
                            pw.print("\t");
                            pw.print(sample);
                            pw.print("\t");
                            pw.print(path);
                            pw.print("\t");
                            pw.print(i + 1);
                            pw.print("\t");
                            pw.print(anomalies.size());
                            pw.print("\t");
                            pw.print(locatable.toNiceString());
                            pw.print("\t");
                            pw.print((int) median);
                            pw.print("\t");
                            pw.print((int) anomalie.depths.stream().mapToDouble(X -> X.doubleValue()).average().orElse(0));
                            pw.print("\t");
                            pw.println();
                        }
                    }
                }
            }
            if (found_anomaly_here) {
                pw.println("<<< " + locatable.toNiceString() + " length:" + StringUtils.niceInt(locatable.getLengthOnReference()));
                pw.println();
            }
        }
        // end while iter
        pw.flush();
        pw.close();
        pw = null;
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(pw);
    }
}

Example 29 with CloseableIterator

use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.

the class SamScanSplitReads method processInput.

@Override
protected int processInput(final SAMFileHeader samHeader, final CloseableIterator<SAMRecord> iter) {
    final IntervalTreeMap<List<Arc>> database = new IntervalTreeMap<>();
    final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(samHeader);
    final Set<String> samples = samHeader.getReadGroups().stream().map(RG -> this.partition.apply(RG)).filter(S -> !StringUtil.isBlank(S)).collect(Collectors.toCollection(TreeSet::new));
    if (samples.isEmpty()) {
        iter.close();
        LOG.error("No samples defined");
        return -1;
    }
    final Set<VCFHeaderLine> meta = new HashSet<>();
    VCFStandardHeaderLines.addStandardFormatLines(meta, false, VCFConstants.GENOTYPE_KEY, VCFConstants.DEPTH_KEY);
    VCFStandardHeaderLines.addStandardInfoLines(meta, false, VCFConstants.DEPTH_KEY, VCFConstants.END_KEY);
    meta.add(new VCFFormatHeaderLine("N5", 1, VCFHeaderLineType.Integer, "Number of clipped reads in 5'"));
    meta.add(new VCFFormatHeaderLine("N3", 1, VCFHeaderLineType.Integer, "Number of clipped reads in 3'"));
    meta.add(new VCFFormatHeaderLine("M5", 1, VCFHeaderLineType.Float, "Median size of the clip in 5'"));
    meta.add(new VCFFormatHeaderLine("M3", 1, VCFHeaderLineType.Float, "Median size of the clip in 3'"));
    meta.add(new VCFInfoHeaderLine("SVLEN", 1, VCFHeaderLineType.Integer, "SV length"));
    meta.add(new VCFInfoHeaderLine("NSAMPLES", 1, VCFHeaderLineType.Integer, "Number of sample having some split reads"));
    final VCFHeader header = new VCFHeader(meta, samples);
    JVarkitVersion.getInstance().addMetaData(this, header);
    header.setSequenceDictionary(dict);
    try {
        try (VariantContextWriter vcw = this.writingVcfConfig.dictionary(dict).open(outputFile)) {
            vcw.writeHeader(header);
            final ProgressFactory.Watcher<SAMRecord> progress = ProgressFactory.newInstance().dictionary(dict).logger(LOG).build();
            String prevContig = null;
            while (iter.hasNext()) {
                final SAMRecord rec = progress.apply(iter.next());
                if (!SAMRecordDefaultFilter.accept(rec))
                    continue;
                final Cigar cigar = rec.getCigar();
                if (cigar == null || cigar.isEmpty() || !cigar.isClipped())
                    continue;
                final String sample = this.partition.getPartion(rec, null);
                if (StringUtil.isBlank(sample))
                    continue;
                if (!rec.getContig().equals(prevContig)) {
                    dump(dict, database, vcw, samples, null);
                    database.clear();
                    prevContig = rec.getContig();
                } else {
                    final int before = rec.getUnclippedStart() - this.buffer_dump_size;
                    dump(dict, database, vcw, samples, before);
                    database.entrySet().removeIf(entries -> entries.getKey().getEnd() < before);
                }
                for (int side = 0; side < 2; ++side) {
                    final int chromStart;
                    final int chromEnd;
                    final byte type;
                    if (side == 0) {
                        if (!cigar.isLeftClipped()) {
                            continue;
                        }
                        chromStart = rec.getUnclippedStart();
                        chromEnd = rec.getStart() - 1;
                        type = VOID_TO_LEFT;
                    } else {
                        if (!cigar.isRightClipped()) {
                            continue;
                        }
                        chromStart = rec.getStart() + 1;
                        chromEnd = rec.getUnclippedEnd();
                        type = RIGHT_TO_VOID;
                    }
                    // final int length = chromEnd  - chromStart;
                    // NON peux augmenter la puissance if(length < 2) continue;
                    final Arc arc = new Arc();
                    arc.sample = sample;
                    arc.tid = rec.getReferenceIndex();
                    arc.chromStart = chromStart;
                    arc.chromEnd = chromEnd;
                    arc.type = type;
                    final Interval rgn = new Interval(rec.getReferenceName(), chromStart, chromEnd);
                    List<Arc> list = database.get(rgn);
                    if (list == null) {
                        list = new ArrayList<>();
                        database.put(rgn, list);
                    }
                    list.add(arc);
                }
            }
            // end iter
            dump(dict, database, vcw, samples, null);
        }
        // end writer
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    }
}

Example 30 with CloseableIterator

use of htsjdk.samtools.util.CloseableIterator in project jvarkit by lindenb.

the class CoverageServer method printImage.

private void printImage(final HttpServletRequest request, final HttpServletResponse response) throws IOException, ServletException {
    int bam_id;
    try {
        bam_id = Integer.parseInt(StringUtils.ifBlank(request.getParameter("id"), "-1"));
    } catch (Exception err) {
        bam_id = -1;
    }
    final SimpleInterval midRegion = parseInterval(request.getParameter("interval"));
    if (midRegion == null || bam_id < 0 || bam_id >= this.bamInput.size()) {
        response.reset();
        response.sendError(HttpStatus.BAD_REQUEST_400, "id:" + bam_id);
        response.flushBuffer();
        return;
    }
    final int extend = (int) (midRegion.getLengthOnReference() * this.extend_factor);
    int xstart = Math.max(midRegion.getStart() - extend, 0);
    int xend = midRegion.getEnd() + extend;
    final SAMSequenceRecord ssr = this.dictionary.getSequence(midRegion.getContig());
    if (ssr != null) {
        xend = Math.min(xend, ssr.getSequenceLength());
    }
    final SimpleInterval region = new SimpleInterval(midRegion.getContig(), xstart, xend);
    if (region.getLengthOnReference() > this.max_window_size) {
        response.reset();
        response.sendError(HttpStatus.BAD_REQUEST_400, "contig:" + midRegion);
        response.flushBuffer();
        return;
    }
    final Counter<Arc> sashimiArcs = new Counter<>();
    final BamInput bam = this.bamInput.get(bam_id);
    if (region.length() <= this.small_region_size) {
        printRaster(bam, midRegion, region, request, response);
        return;
    }
    final boolean normalize = request.getParameter("normalize") != null;
    final SamReaderFactory srf = SamReaderFactory.make().validationStringency(ValidationStringency.LENIENT).referenceSequence(this.faidxRef);
    try (SamReader sr = srf.open(bam.bamPath)) {
        final int[] int_coverage = new int[region.getLengthOnReference()];
        Arrays.fill(int_coverage, 0);
        try (CloseableIterator<SAMRecord> iter = sr.query(region.getContig(), region.getStart(), region.getEnd(), false)) {
            while (iter.hasNext()) {
                final SAMRecord rec = iter.next();
                if (!acceptRead(rec))
                    continue;
                final Cigar cigar = rec.getCigar();
                if (cigar == null || cigar.isEmpty())
                    continue;
                int ref = rec.getAlignmentStart();
                for (final CigarElement ce : cigar) {
                    final CigarOperator op = ce.getOperator();
                    if (op.consumesReferenceBases()) {
                        if (this.enable_sashimi && op.equals(CigarOperator.N)) {
                            sashimiArcs.incr(new Arc(ref, ref + ce.getLength()));
                        }
                        if (op.consumesReadBases()) {
                            for (int x = 0; x < ce.getLength(); ++x) {
                                int pos = ref + x;
                                if (pos < region.getStart())
                                    continue;
                                if (pos > region.getEnd())
                                    break;
                                int_coverage[pos - region.getStart()]++;
                            }
                        }
                        ref += ce.getLength();
                    }
                }
            }
        }
        /* smooth coverage */
        if (int_coverage.length > image_width) {
            final int[] copy = Arrays.copyOf(int_coverage, int_coverage.length);
            final int len = Math.max(1, int_coverage.length / 100);
            for (int i = 0; i < int_coverage.length; i++) {
                int j = Math.max(0, i - len);
                double sum = 0;
                int count = 0;
                while (j < i + len && j < copy.length) {
                    sum += copy[j];
                    j++;
                    count++;
                }
                int_coverage[i] = (int) (sum / count);
            }
        }
        final double[] norm_coverage = new double[int_coverage.length];
        final double median;
        /* normalize on median */
        if (normalize) {
            final Coverage leftrightcov = new Coverage(extend * 2);
            for (int x = region.getStart(); x < midRegion.getStart(); x++) {
                final int idx = x - region.getStart();
                leftrightcov.add(int_coverage[idx]);
            }
            for (int x = midRegion.getEnd() + 1; x <= region.getEnd(); x++) {
                final int idx = x - region.getStart();
                leftrightcov.add(int_coverage[idx]);
            }
            median = Math.max(1.0, leftrightcov.median());
            // LOG.info("median is "+median+" "+leftrightcov.median());
            for (int x = 0; x < int_coverage.length; ++x) {
                norm_coverage[x] = int_coverage[x] / median;
            }
        } else /* no normalisation */
        {
            /* won't be used */
            median = Double.NaN;
            for (int x = 0; x < int_coverage.length; ++x) {
                norm_coverage[x] = int_coverage[x];
            }
        }
        final double real_max_cov = DoubleStream.of(norm_coverage).max().orElse(1.0);
        final double max_cov = Math.max((normalize ? 2 : 10), real_max_cov);
        final double pixelperbase = image_width / (double) norm_coverage.length;
        final IntFunction<Double> pos2pixel = POS -> ((POS - region.getStart()) / (double) region.getLengthOnReference()) * image_width;
        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(Color.WHITE);
        g.fillRect(0, 0, image_width + 1, image_height + 1);
        for (int x = 0; x < norm_coverage.length; ++x) {
            final double height = image_height * (norm_coverage[x] / max_cov);
            if (normalize)
                g.setColor(Color.DARK_GRAY);
            else if (max_cov < 10)
                g.setColor(Color.RED);
            else if (max_cov < 20)
                g.setColor(Color.BLUE);
            else
                g.setColor(Color.DARK_GRAY);
            g.fill(new Rectangle2D.Double(x * pixelperbase, image_height - height, pixelperbase, height));
        }
        g.setColor(Color.DARK_GRAY);
        g.drawString("max-cov:" + IntStream.of(int_coverage).max().orElse(0) + (normalize ? " normalized on median (" + median + ")" : "") + " sample:" + bam.sample + " " + region.toNiceString(), 10, 10);
        /* ticks for vertical axis */
        g.setColor(Color.MAGENTA);
        for (int i = 1; i < 10; i++) {
            double cov = max_cov / 10.0 * i;
            if (!normalize)
                cov = Math.ceil(cov);
            final double y = image_height - image_height / 10.0 * i;
            if (!normalize && i > 0 && (int) cov == Math.ceil(max_cov / 10.0 * (i - 1)))
                continue;
            g.drawLine(0, (int) y, 5, (int) y);
            g.drawString(normalize ? String.format("%.2f", cov) : String.valueOf((int) cov), 7, (int) y);
        }
        /* vertical line for original view */
        g.setColor(Color.PINK);
        double vertical = ((midRegion.getStart() - region.getStart()) / (double) region.getLengthOnReference()) * image_width;
        g.draw(new Line2D.Double(vertical, 0, vertical, image_height));
        vertical = ((midRegion.getEnd() - region.getStart()) / (double) region.getLengthOnReference()) * image_width;
        g.draw(new Line2D.Double(vertical, 0, vertical, image_height));
        if (normalize) {
            /* horizontal line for median 0.5 / 1 / 1.5 */
            for (int t = 1; t < 4; ++t) {
                g.setColor(t == 2 ? Color.ORANGE : Color.PINK);
                final double mediany = image_height - ((0.5 * t) / max_cov) * image_height;
                g.draw(new Line2D.Double(0, mediany, image_width, mediany));
            }
        }
        if (this.enable_sashimi && !sashimiArcs.isEmpty()) {
            final double max_count = sashimiArcs.getMaxCount().orElse(1L);
            g.setColor(Color.GREEN);
            for (final Arc arc : sashimiArcs.keySet()) {
                final double x1 = pos2pixel.apply(arc.start);
                final double x2 = pos2pixel.apply(arc.end);
                final double distance = x2 - x1;
                final GeneralPath curve = new GeneralPath();
                curve.moveTo(x1, image_height);
                curve.curveTo(x1, image_height, x1 + distance / 2.0, image_height - Math.min(distance, image_height * 0.75), x2, image_height);
                final double weight = (sashimiArcs.count(arc) / max_count) * 5;
                final Stroke oldStroke = g.getStroke();
                final Composite oldComposite = g.getComposite();
                g.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.5f));
                g.setStroke(new BasicStroke((float) weight, BasicStroke.CAP_ROUND, BasicStroke.JOIN_ROUND));
                g.draw(curve);
                g.setStroke(oldStroke);
                g.setComposite(oldComposite);
            }
        }
        writeGenes(g, region);
        writeKnownCnv(g, region);
        g.setColor(Color.GRAY);
        g.drawRect(0, 0, img.getWidth(), img.getHeight());
        writeImage(img, bam, region, response);
    }
}