Examples with LineIterator com.github.lindenb.jvarkit.util.iterator.LineIterator used on opensource projects

Example 1 with LineIterator

use of com.github.lindenb.jvarkit.util.iterator.LineIterator 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 2 with LineIterator

use of com.github.lindenb.jvarkit.util.iterator.LineIterator in project jvarkit by lindenb.

the class VcfUcscGdb method readRemoteResources.

private List<RemoteBigFile> readRemoteResources(final Path path) throws IOException {
    final List<RemoteBigFile> remoteBigFiles = new ArrayList<>();
    IOUtil.assertFileIsReadable(path);
    try (BufferedReader br = IOUtil.openFileForBufferedReading(path)) {
        final HashMap<String, String> hash = new HashMap<>();
        final Function<String, String> required = (K) -> {
            if (!hash.containsKey(K))
                throw new RuntimeIOException("Key \"" + K + "\" missing. Found: " + hash.keySet());
            final String v = hash.get(K).trim();
            if (StringUtils.isBlank(v))
                throw new RuntimeIOException("Key \"" + K + "\" is empty");
            return v;
        };
        try (LineIterator iter = new LineIterator(br)) {
            for (; ; ) {
                final String line = (iter.hasNext() ? iter.next() : null);
                if (StringUtils.isBlank(line)) {
                    if (hash.getOrDefault("enabled", "true").equals("false")) {
                        hash.clear();
                    }
                    if (!hash.isEmpty()) {
                        final RemoteBigFile bf = new RemoteBigFile();
                        bf.url = required.apply("url");
                        if (hash.containsKey("name")) {
                            bf.name = hash.get("name");
                        } else {
                            bf.name = bf.url;
                            int slah = bf.name.lastIndexOf('/');
                            bf.name = bf.name.substring(slah + 1);
                            int dot = bf.name.lastIndexOf('.');
                            bf.name = bf.name.substring(0, dot).replace('.', '_').replace('-', '_').replace(',', '_');
                        }
                        if (remoteBigFiles.stream().anyMatch(R -> R.name.equals(bf.name))) {
                            bf.close();
                            throw new RuntimeIOException("Duplicate remote resource: " + hash);
                        }
                        if (hash.containsKey("accept")) {
                            bf.accept = new JexlPredicate(hash.get("accept"));
                        }
                        if (hash.containsKey("tostring")) {
                            bf.converter = new JexlToString(hash.get("tostring"));
                        }
                        if (hash.containsKey("desc")) {
                            bf.description = hash.get("desc");
                        } else if (hash.containsKey("description")) {
                            bf.description = hash.get("description");
                        } else {
                            bf.description = "Data from " + bf.url;
                        }
                        if (hash.containsKey("limit")) {
                            bf.limit = Integer.parseInt(hash.get("limit"));
                        }
                        if (hash.containsKey("fractV")) {
                            bf.fractionOfVariant = Double.parseDouble(hash.get("fractV"));
                        }
                        if (hash.containsKey("fractF")) {
                            bf.fractionOfVariant = Double.parseDouble(hash.get("fractF"));
                        }
                        if (hash.containsKey("aggregate")) {
                            bf.wigAggregate = hash.get("aggregate");
                            if (!(bf.wigAggregate.equals("min") || bf.wigAggregate.equals("max"))) {
                                bf.close();
                                throw new RuntimeIOException("Bad value for aggregate accepted:(min/max))");
                            }
                        }
                        remoteBigFiles.add(bf);
                    }
                    if (line == null)
                        break;
                    hash.clear();
                    continue;
                }
                if (line.startsWith("#"))
                    continue;
                int sep = line.indexOf(':');
                if (sep == -1)
                    sep = line.indexOf('=');
                if (sep == -1)
                    throw new RuntimeIOException("Cannot find ':' or '=' in  " + line);
                final String key = line.substring(0, sep).toLowerCase().trim();
                if (hash.containsKey(key))
                    throw new RuntimeIOException("Duplicate key " + key + " in resource: " + hash);
                final String value = line.substring(sep + 1).trim();
                hash.put(key, value);
            }
        }
    }
    return remoteBigFiles;
}

Example 3 with LineIterator

use of com.github.lindenb.jvarkit.util.iterator.LineIterator in project jvarkit by lindenb.

the class BedMergeCnv method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.fraction <= 0) {
        LOG.error("Bad fraction : " + this.fraction);
        return -1;
    }
    PrintWriter w = null;
    try {
        if (args.isEmpty()) {
            parseBed(new LineIterator(new InputStreamReader(stdin())));
        } else {
            for (final String fname : args) {
                final IntervalListProvider ivp = IntervalListProvider.from(fname);
                final Iterator<? extends Locatable> iter = ivp.stream().iterator();
                parseLocatables(iter);
            }
        }
        w = super.openPathOrStdoutAsPrintWriter(this.outputFile);
        final List<Cluster> ordered = new ArrayList<>(this.clusters.values());
        Collections.sort(ordered, (A, B) -> {
            int i = A.getContig().compareTo(B.getContig());
            if (i != 0)
                return i;
            i = Integer.compare(A.getStart(), B.getStart());
            if (i != 0)
                return i;
            return Integer.compare(A.getEnd(), B.getEnd());
        });
        Cluster prev = null;
        for (final Cluster bc : ordered) {
            w.print(bc.getContig());
            w.print("\t");
            w.print(bc.getStart() - 1);
            w.print("\t");
            w.print(bc.getEnd());
            w.print("\t");
            w.print(bc.intervals.size());
            w.print("\t");
            w.print(prev != null && prev.contigsMatch(bc) ? String.format("%.2f", this.fractionOverlap(prev, bc)) : "0.0");
            w.print("\t");
            w.print(bc.getLabel());
            w.println();
            prev = bc;
        }
        w.flush();
        w.close();
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(w);
    }
}

Example 4 with LineIterator

use of com.github.lindenb.jvarkit.util.iterator.LineIterator in project jvarkit by lindenb.

the class GffLiftOver method doWork.

@Override
public int doWork(final List<String> args) {
    PrintWriter fail = null;
    try {
        fail = this.failedFile == null ? new PrintWriter(new NullOuputStream()) : super.openPathOrStdoutAsPrintWriter(this.failedFile);
        this.liftOver = new LiftOver(this.liftOverFile);
        try (BufferedReader in = openBufferedReader(oneFileOrNull(args))) {
            final Gff3Codec codec = new Gff3Codec();
            final LineIterator liter = new LineIterator(in);
            final FeatureCodecHeader header = codec.readHeader(liter);
            try (Gff3Writer p = new Gff3Writer(super.openPathOrStdoutAsStream(this.outputFile))) {
                while (!codec.isDone(liter)) {
                    final Gff3Feature feature = codec.decode(liter);
                    feature.getChildren().stream().filter(F -> F.getType().equals("gene")).forEach(GENE -> {
                        if (!canLiftOver(GENE)) {
                            continue;
                        }
                        p.addFeature(GENE);
                        GENE.getChildren().stream().forEach(TRANSCRIPT -> {
                            final List<Gff3Feature> features = new ArrayList<>();
                            features.add(TRANSCRIPT);
                            TRANSCRIPT.getChildren().forEach(FEAT -> features.add(FEAT));
                            if (!features.stream().allMatch(F2 -> canLiftOver(F2))) {
                                continue;
                            }
                            for (Gff3Feature feat : features) {
                                final Interval lifted = applyLiftOver(feat);
                                final Gff3FeatureImpl feat2 = new Gff3FeatureImpl(lifted.getContig(), feat.getSource(), feat.getType(), lifted.getStart(), lifted.getEnd(), feat.getScore(), lifted.getStrand(), feat.getPhase(), feat.getAttributes());
                                p.addFeature(feat2);
                            }
                        });
                    });
                }
                codec.close(liter);
                p.flush();
            }
        }
        fail.flush();
        fail.close();
        fail = null;
        return RETURN_OK;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(fail);
        CloserUtil.close(p);
    }
}

Example 5 with LineIterator

use of com.github.lindenb.jvarkit.util.iterator.LineIterator in project jvarkit by lindenb.

the class CoveragePlotter method doWork.

@Override
public int doWork(final List<String> args) {
    ArchiveFactory archive = null;
    PrintWriter manifest = null;
    try {
        if (extend < 1.0) {
            LOG.error("extend is lower than 1 :" + this.extend);
            return -1;
        }
        final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(this.refPath);
        final SamReaderFactory samReaderFactory = SamReaderFactory.makeDefault().referenceSequence(CoveragePlotter.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();
        }
        final BufferedImage image = new BufferedImage(this.dimension.width, this.dimension.height, BufferedImage.TYPE_INT_ARGB);
        final BufferedImage offscreen = new BufferedImage(this.dimension.width, this.dimension.height, BufferedImage.TYPE_INT_ARGB);
        final double y_mid = this.dimension.getHeight() / 2.0;
        final ToDoubleFunction<Double> normToPixelY = NORM -> this.dimension.getHeight() - NORM * y_mid;
        final DiscreteMedian<Integer> discreteMedian = new DiscreteMedian<>();
        manifest = (this.manifestPath == null ? new PrintWriter(new NullOuputStream()) : IOUtils.openPathForPrintWriter(this.manifestPath));
        archive = ArchiveFactory.open(this.outputFile);
        while (iter.hasNext()) {
            final Locatable the_locatable = iter.next();
            String label = "";
            if (the_locatable instanceof BedLine) {
                final BedLine bedline = BedLine.class.cast(the_locatable);
                label = bedline.getOrDefault(3, label);
            }
            final SimpleInterval rawRegion = new SimpleInterval(the_locatable);
            final SimpleInterval extendedRegion;
            /* extend interval */
            if (this.extend > 1) {
                final int L1 = rawRegion.getLengthOnReference();
                final int L2 = (int) Math.ceil(L1 * this.extend);
                final int mid = rawRegion.getCenter().getPosition();
                int x0 = mid - L2 / 2;
                if (x0 < 0)
                    x0 = 1;
                int x1 = mid + L2 / 2;
                final SAMSequenceRecord ssr = dict.getSequence(rawRegion.getContig());
                if (ssr != null)
                    x1 = Math.min(ssr.getSequenceLength(), x1);
                if (x0 > x1)
                    continue;
                extendedRegion = new SimpleInterval(rawRegion.getContig(), x0, x1);
            } else {
                extendedRegion = rawRegion;
            }
            final ToDoubleFunction<Integer> pos2pixel = POS -> (POS - extendedRegion.getStart()) / (double) extendedRegion.getLengthOnReference() * this.dimension.getWidth();
            final String theGeneName = this.getGenes(extendedRegion).filter(G -> G.getType().equals("gene")).map(G -> G.getAttribute("gene_name")).filter(S -> !StringUtils.isBlank(S)).findFirst().orElse(null);
            final String outputFilename = this.prefix + extendedRegion.getContig() + "_" + extendedRegion.getStart() + "_" + extendedRegion.getEnd() + (StringUtils.isBlank(label) ? "" : "." + label.replaceAll("[^A-Za-z\\-\\.0-9]+", "_")) + (StringUtils.isBlank(theGeneName) || !StringUtils.isBlank(label) ? "" : "." + theGeneName.replaceAll("[^A-Za-z\\-\\.0-9]+", "_")) + ".png";
            final Graphics2D g2 = offscreen.createGraphics();
            g2.setColor(Color.BLACK);
            g2.setComposite(AlphaComposite.getInstance(AlphaComposite.CLEAR));
            g2.fillRect(0, 0, this.dimension.width, this.dimension.height);
            g2.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER));
            g2.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
            g2.setRenderingHint(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);
            final Graphics2D g = image.createGraphics();
            g.setColor(Color.WHITE);
            g.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
            g.setRenderingHint(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);
            g.fillRect(0, 0, this.dimension.width, this.dimension.height);
            // draw exons
            getGenes(extendedRegion).filter(G -> G.getType().equals("exon")).forEach(EX -> {
                g.setColor(new Color(240, 240, 240));
                double x1 = pos2pixel.applyAsDouble(EX.getStart());
                double x2 = pos2pixel.applyAsDouble(EX.getEnd());
                g.fill(new Rectangle2D.Double(x1, 0, (x2 - x1), this.dimension.height));
            });
            int y = (int) (this.dimension.height / 2.0);
            g.setColor(Color.BLUE);
            g.drawLine(0, y, image.getWidth(), y);
            y = (int) (this.dimension.height / 4.0);
            g.setColor(Color.CYAN);
            g.drawLine(0, y, image.getWidth(), y);
            y = (int) (3.0 * this.dimension.height / 4.0);
            g.drawLine(0, y, image.getWidth(), y);
            g.setColor(Color.DARK_GRAY);
            g.drawRect(0, 0, this.dimension.width - 1, this.dimension.height - 1);
            drawGenes(g, new Rectangle(0, 0, image.getWidth(), image.getHeight()), extendedRegion);
            drawKnownCnv(g, new Rectangle(0, 0, image.getWidth(), image.getHeight()), extendedRegion);
            if (this.extend > 1) {
                g.setColor(Color.GREEN);
                int x = (int) pos2pixel.applyAsDouble(rawRegion.getStart());
                g.drawLine(x, 0, x, image.getHeight());
                x = (int) pos2pixel.applyAsDouble(rawRegion.getEnd());
                g.drawLine(x, 0, x, image.getHeight());
            }
            final int[] depth = new int[extendedRegion.getLengthOnReference()];
            final int[] copy = new int[depth.length];
            final BitSet blackListedPositions = new BitSet(depth.length);
            final Map<String, Point2D> sample2maxPoint = new HashMap<>(inputBams.size());
            boolean drawAbove = false;
            // fill black listed regions
            if (this.blackListedPath != null) {
                try (TabixReader tbr = new TabixReader(this.blackListedPath.toString())) {
                    final ContigNameConverter cvt = ContigNameConverter.fromContigSet(tbr.getChromosomes());
                    final String ctg = cvt.apply(extendedRegion.getContig());
                    if (!StringUtils.isBlank(ctg)) {
                        final BedLineCodec codec = new BedLineCodec();
                        final TabixReader.Iterator tbxr = tbr.query(ctg, extendedRegion.getStart(), extendedRegion.getEnd());
                        for (; ; ) {
                            final String line = tbxr.next();
                            if (line == null)
                                break;
                            final BedLine bed = codec.decode(line);
                            if (bed == null)
                                continue;
                            int p1 = Math.max(bed.getStart(), extendedRegion.getStart());
                            while (p1 <= extendedRegion.getEnd() && p1 <= bed.getEnd()) {
                                blackListedPositions.set(p1 - extendedRegion.getStart());
                                ++p1;
                            }
                        }
                    }
                } catch (Throwable err) {
                    LOG.warn(err);
                }
            }
            if (this.skip_original_interval_for_median && this.extend > 1) {
                for (int i = the_locatable.getStart(); i <= the_locatable.getEnd(); i++) {
                    blackListedPositions.set(i - extendedRegion.getStart());
                }
            }
            final Set<String> sampleNames = new TreeSet<>();
            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));
                    sampleNames.add(sample);
                    /* sample has deletion */
                    int count_has_dp_le_0_5 = 0;
                    /* sample has dup */
                    int count_has_dp_ge_1_5 = 0;
                    /* longest arc */
                    int longest_arc = 0;
                    SequenceUtil.assertSequenceDictionariesEqual(dict, header.getSequenceDictionary());
                    Arrays.fill(depth, 0);
                    try (CloseableIterator<SAMRecord> siter = sr.queryOverlapping(extendedRegion.getContig(), extendedRegion.getStart(), extendedRegion.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;
                            if (this.alpha_arc > 0.0 && rec.getReadPairedFlag() && !rec.getMateUnmappedFlag() && !rec.getProperPairFlag() && rec.getReferenceIndex().equals(rec.getMateReferenceIndex()) && (!this.arc_only_rr_ff || (rec.getReadNegativeStrandFlag() == rec.getMateNegativeStrandFlag()))) {
                                longest_arc = Math.max(longest_arc, Math.abs(rec.getMateAlignmentStart() - rec.getAlignmentStart()));
                                final double xstart = pos2pixel.applyAsDouble(rec.getAlignmentStart());
                                final double xend = pos2pixel.applyAsDouble(rec.getMateAlignmentStart());
                                final double len = (xend - xstart);
                                if (Math.abs(len) > this.min_invert && Math.abs(len) < this.max_invert) {
                                    final double y2 = y_mid + (drawAbove ? -1 : 1) * Math.min(y_mid, Math.abs(len / 2.0));
                                    final Composite oldComposite = g2.getComposite();
                                    g2.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, (float) this.alpha_arc));
                                    g2.setColor(getColor("arc", Color.ORANGE));
                                    final GeneralPath curve = new GeneralPath();
                                    curve.moveTo(xstart, y_mid);
                                    curve.quadTo(xstart + len / 2.0, y2, xend, y_mid);
                                    g2.draw(curve);
                                    g2.setComposite(oldComposite);
                                    drawAbove = !drawAbove;
                                }
                            }
                            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 < extendedRegion.getStart())
                                                continue;
                                            if (pos > extendedRegion.getEnd())
                                                break;
                                            depth[pos - extendedRegion.getStart()]++;
                                        }
                                    }
                                    ref += len;
                                }
                            }
                        }
                    // loop cigar
                    }
                    if (extendedRegion.getLengthOnReference() > image.getWidth() && extendedRegion.getLengthOnReference() > this.window_smooth_size) {
                        // smooth
                        final int bases_per_pixel = this.window_smooth_size;
                        System.arraycopy(depth, 0, copy, 0, depth.length);
                        for (int i = 0; i < depth.length && bases_per_pixel > 1; i++) {
                            double t = 0;
                            int count = 0;
                            for (int j = Math.max(0, i - bases_per_pixel); j <= i + bases_per_pixel && j < depth.length; j++) {
                                t += copy[j];
                                count++;
                            }
                            if (count == 0)
                                continue;
                            depth[i] = (int) (t / count);
                        }
                    }
                    discreteMedian.clear();
                    for (int i = 0; i < depth.length; i++) {
                        if (depth[i] > this.max_depth)
                            continue;
                        if (blackListedPositions.get(i))
                            continue;
                        discreteMedian.add(depth[i]);
                    }
                    final double median = discreteMedian.getMedian().orElse(0);
                    if (median <= 0) {
                        final String msg = "Skipping " + sample + " " + extendedRegion + " because median is 0";
                        LOG.warning(msg);
                        manifest.println("#" + msg);
                        continue;
                    }
                    Point2D max_position = null;
                    double max_distance_to_1 = 0.0;
                    final GeneralPath line = new GeneralPath();
                    for (int x = 0; x < image.getWidth(); x++) {
                        discreteMedian.clear();
                        int pos1 = (int) Math.floor(((x + 0) / (double) image.getWidth()) * depth.length);
                        final int pos2 = (int) Math.ceil(((x + 0) / (double) image.getWidth()) * depth.length);
                        while (pos1 <= pos2 && pos1 < depth.length) {
                            discreteMedian.add(depth[pos1]);
                            pos1++;
                        }
                        final double average = discreteMedian.getMedian().orElse(0);
                        final double normDepth = (average / median);
                        if (normDepth <= 0.5)
                            count_has_dp_le_0_5++;
                        if (normDepth >= 1.5)
                            count_has_dp_ge_1_5++;
                        final double y2 = normToPixelY.applyAsDouble(normDepth);
                        double distance_to_1 = Math.abs(normDepth - 1.0);
                        if (distance_to_1 > 0.3 && (max_position == null || distance_to_1 > max_distance_to_1)) {
                            max_distance_to_1 = distance_to_1;
                            max_position = new Point2D.Double(x, y2);
                        }
                        if (x == 0) {
                            line.moveTo(x, y2);
                        } else {
                            line.lineTo(x, y2);
                        }
                    }
                    g.setColor(max_distance_to_1 < 0.1 ? Color.lightGray : Color.GRAY);
                    final Stroke oldStroke = g.getStroke();
                    final Composite oldComposite = g.getComposite();
                    g.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, (float) this.alpha));
                    g.setStroke(new BasicStroke(0.5f, BasicStroke.CAP_BUTT, BasicStroke.JOIN_ROUND));
                    g.draw(line);
                    g.setStroke(oldStroke);
                    g.setComposite(oldComposite);
                    if (max_position != null)
                        sample2maxPoint.put(sample, max_position);
                    // fill manifest
                    manifest.print(rawRegion.getContig());
                    manifest.print("\t");
                    manifest.print(rawRegion.getStart() - 1);
                    manifest.print("\t");
                    manifest.print(rawRegion.getEnd());
                    manifest.print("\t");
                    manifest.print(extendedRegion.getStart() - 1);
                    manifest.print("\t");
                    manifest.print(extendedRegion.getEnd());
                    manifest.print("\t");
                    manifest.print(outputFilename);
                    manifest.print("\t");
                    manifest.print(StringUtils.isBlank(theGeneName) ? "." : theGeneName);
                    manifest.print("\t");
                    manifest.print(sample);
                    manifest.print("\t");
                    manifest.print(count_has_dp_le_0_5);
                    manifest.print("\t");
                    manifest.print((int) ((count_has_dp_le_0_5 / (double) extendedRegion.getLengthOnReference()) * 100.0));
                    manifest.print("\t");
                    manifest.print(count_has_dp_ge_1_5);
                    manifest.print("\t");
                    manifest.print((int) ((count_has_dp_ge_1_5 / (double) extendedRegion.getLengthOnReference()) * 100.0));
                    manifest.print("\t");
                    manifest.print(median);
                    manifest.print("\t");
                    manifest.print(discreteMedian.getStandardDeviation().orElse(-99999));
                    manifest.print("\t");
                    manifest.print(longest_arc);
                    manifest.println();
                }
            // end loop over samples
            }
            g2.dispose();
            g.drawImage(offscreen, 0, 0, null);
            g.setColor(Color.BLACK);
            final String label_samples = sampleNames.size() > 10 ? "N=" + StringUtils.niceInt(sampleNames.size()) : String.join(";", sampleNames);
            final Set<String> all_genes = getGenes(extendedRegion).filter(G -> G.getType().equals("gene")).map(G -> G.getAttribute("gene_name")).filter(F -> !StringUtils.isBlank(F)).collect(Collectors.toCollection(TreeSet::new));
            g.drawString(extendedRegion.toNiceString() + " Length:" + StringUtils.niceInt(extendedRegion.getLengthOnReference()) + " Sample(s):" + label_samples + " " + label + " Gene(s):" + (all_genes.size() > 20 ? "N=" + StringUtils.niceInt(all_genes.size()) : String.join(";", all_genes)), 10, 10);
            if (!sample2maxPoint.isEmpty()) {
                /**
                 * draw sample names
                 */
                g.setColor(getColor("sample", Color.BLUE));
                final int sampleFontSize = 7;
                final Font oldFont = g.getFont();
                g.setFont(new Font(oldFont.getName(), Font.PLAIN, sampleFontSize));
                for (final String sample : sample2maxPoint.keySet()) {
                    final Point2D pt = sample2maxPoint.get(sample);
                    double sny = pt.getY();
                    if (sny > y_mid)
                        sny += sampleFontSize;
                    g.drawString(sample, (int) pt.getX(), (int) Math.min(this.dimension.height - sampleFontSize, Math.max(sampleFontSize, sny)));
                }
                g.setFont(oldFont);
            }
            g.dispose();
            try (OutputStream out = archive.openOuputStream(outputFilename)) {
                ImageIO.write(image, "PNG", out);
                out.flush();
            }
        }
        // end while iter
        archive.close();
        archive = null;
        manifest.flush();
        manifest.close();
        manifest = null;
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(manifest);
        CloserUtil.close(archive);
    }
}