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

Example 51 with CloseableIterator

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

the class VcfPeekVcf method doVcfToVcf.

/**
 * public for knime
 */
@Override
public int doVcfToVcf(final String inputName, final VcfIterator vcfIn, final VariantContextWriter out) {
    try {
        final Set<String> unmatchedcontigs = new HashSet<>();
        final VCFHeader h = vcfIn.getHeader();
        final VCFHeader h2 = new VCFHeader(h);
        super.addMetaData(h2);
        final Map<String, VCFInfoHeaderLine> databaseTags = new HashMap<String, VCFInfoHeaderLine>();
        final VCFHeader databaseHeader = this.indexedVcfFileReader.getFileHeader();
        final ContigNameConverter nameConverter = (h.getSequenceDictionary() != null && !h.getSequenceDictionary().isEmpty() && databaseHeader.getSequenceDictionary() != null && !databaseHeader.getSequenceDictionary().isEmpty() ? ContigNameConverter.fromDictionaries(h.getSequenceDictionary(), databaseHeader.getSequenceDictionary()) : ContigNameConverter.getIdentity()).setOnNotFound(this.onContigNotFound);
        ;
        for (final String key : this.peek_info_tags) {
            VCFInfoHeaderLine hinfo = databaseHeader.getInfoHeaderLine(key);
            if (hinfo == null) {
                final String msg = "INFO name=" + key + " missing in " + this.resourceVcfFile;
                if (this.missingIdIsError) {
                    LOG.warn(msg);
                    continue;
                } else {
                    LOG.error(msg);
                    return -1;
                }
            }
            switch(hinfo.getCountType()) {
                case G:
                    throw new JvarkitException.UserError("Cannot handle VCFHeaderLineCount.G for " + hinfo.getID());
                default:
                    databaseTags.put(hinfo.getID(), hinfo);
                    break;
            }
            hinfo = VCFUtils.renameVCFInfoHeaderLine(hinfo, this.peekTagPrefix + key);
            if (h2.getInfoHeaderLine(hinfo.getID()) != null) {
                throw new JvarkitException.UserError("key " + this.peekTagPrefix + key + " already defined in VCF header");
            }
            h2.addMetaDataLine(hinfo);
            ;
        }
        out.writeHeader(h2);
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(h).logger(LOG);
        while (vcfIn.hasNext()) {
            final VariantContext ctx = progress.watch(vcfIn.next());
            final String outContig = nameConverter.apply(ctx.getContig());
            if (outContig == null) {
                unmatchedcontigs.add(ctx.getContig());
                continue;
            }
            final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
            CloseableIterator<VariantContext> iter = this.indexedVcfFileReader.query(outContig, Math.max(0, ctx.getStart() - 1), (ctx.getEnd() + 1));
            while (iter.hasNext()) {
                final VariantContext ctx2 = iter.next();
                if (!outContig.equals(ctx2.getContig()))
                    continue;
                if (ctx.getStart() != ctx2.getStart())
                    continue;
                if (!ctx.getReference().equals(ctx2.getReference()))
                    continue;
                final boolean okAllele;
                switch(altAlleleMatcher) {
                    case all:
                        {
                            okAllele = ctx.getAlternateAlleles().stream().filter(A -> ctx2.hasAlternateAllele(A)).count() == ctx.getAlternateAlleles().size();
                            break;
                        }
                    case at_least_one:
                        {
                            okAllele = ctx.getAlternateAlleles().stream().filter(A -> ctx2.hasAlternateAllele(A)).findAny().isPresent();
                            break;
                        }
                    case none:
                        okAllele = true;
                        break;
                    default:
                        throw new IllegalStateException(altAlleleMatcher.name());
                }
                if (!okAllele)
                    continue;
                if (this.peekId && ctx2.hasID()) {
                    vcb.id(ctx2.getID());
                }
                boolean somethingWasChanged = false;
                for (final String key : databaseTags.keySet()) {
                    if (!ctx2.hasAttribute(key))
                        continue;
                    final VCFInfoHeaderLine dbHeader = databaseTags.get(key);
                    switch(dbHeader.getCountType()) {
                        case A:
                            {
                                final List<Object> newatt = new ArrayList<>();
                                final List<Object> ctx2att = ctx2.getAttributeAsList(key);
                                for (int i = 0; i < ctx.getAlternateAlleles().size(); ++i) {
                                    final Allele ctxalt = ctx.getAlternateAllele(i);
                                    int index2 = ctx2.getAlternateAlleles().indexOf(ctxalt);
                                    if (index2 == -1 || index2 >= ctx2att.size()) {
                                        newatt.add(null);
                                    } else {
                                        newatt.add(ctx2att.get(index2));
                                    }
                                }
                                if (newatt.stream().filter(Obj -> !(Obj == null || VCFConstants.EMPTY_INFO_FIELD.equals(Obj))).count() > 0) {
                                    vcb.attribute(this.peekTagPrefix + key, newatt);
                                    somethingWasChanged = true;
                                }
                                break;
                            }
                        case R:
                            {
                                final List<Object> newatt = new ArrayList<>();
                                final List<Object> ctx2att = ctx2.getAttributeAsList(key);
                                for (int i = 0; i < ctx.getAlleles().size(); ++i) {
                                    final Allele ctxalt = ctx.getAlleles().get(i);
                                    int index2 = ctx2.getAlleleIndex(ctxalt);
                                    if (index2 == -1 || index2 >= ctx2att.size()) {
                                        newatt.add(null);
                                    } else {
                                        newatt.add(ctx2att.get(index2));
                                    }
                                }
                                if (newatt.stream().filter(Obj -> !(Obj == null || VCFConstants.EMPTY_INFO_FIELD.equals(Obj))).count() > 0) {
                                    vcb.attribute(this.peekTagPrefix + key, newatt);
                                    somethingWasChanged = true;
                                }
                                break;
                            }
                        default:
                            {
                                final Object o = ctx2.getAttribute(key);
                                vcb.attribute(this.peekTagPrefix + key, o);
                                somethingWasChanged = true;
                                break;
                            }
                    }
                }
                if (somethingWasChanged)
                    break;
            }
            iter.close();
            iter = null;
            out.add(vcb.make());
            if (out.checkError())
                break;
        }
        progress.finish();
        if (!unmatchedcontigs.isEmpty()) {
            LOG.debug("Unmatched contigs: " + unmatchedcontigs.stream().collect(Collectors.joining("; ")));
        }
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    }
}

Example 52 with CloseableIterator

use of htsjdk.samtools.util.CloseableIterator in project gridss by PapenfussLab.

the class GeneratePonBedpe method filteredMerge.

private Iterator<Pair<BreakendSummary, Integer>> filteredMerge(GenomicProcessingContext pc, List<File> file, List<Integer> ordinals) {
    List<CloseableIterator<Pair<BreakendSummary, Integer>>> fileIt = new ArrayList<>();
    for (File f : file) {
        fileIt.add(getFilteredIterator(pc, f, ordinals));
    }
    // allocate VCF files to worker threads
    final AtomicInteger counter = new AtomicInteger(0);
    final int size = WORKER_THREADS;
    final List<Iterator<Pair<BreakendSummary, Integer>>> partitioned = fileIt.stream().collect(Collectors.groupingBy(it -> counter.getAndIncrement() % size)).values().stream().map(list -> list.size() == 1 ? list.get(0) : DeterministicIterators.mergeSorted(list, ByBreakendStartEnd)).map(it -> new AsyncBufferedIterator<>(it, "AsyncVCF")).collect(Collectors.toList());
    Iterator<Pair<BreakendSummary, Integer>> mergedIt = new AsyncBufferedIterator<>(DeterministicIterators.mergeSorted(partitioned, ByBreakendStartEnd), "Merged VCF reader");
    return mergedIt;
}

Example 53 with CloseableIterator

use of htsjdk.samtools.util.CloseableIterator in project gridss by PapenfussLab.

the class AllocateEvidence method annotateAssembly.

private CloseableIterator<DirectedEvidence> annotateAssembly(CloseableIterator<DirectedEvidence> it) {
    List<Closeable> assToClose = new ArrayList<>();
    List<Iterator<SAMRecord>> rawAssemblies = new ArrayList<>();
    int windowSize = 0;
    for (AssemblyEvidenceSource aes : getAssemblySource()) {
        // need to use the raw breakend assembly file (prior to realignment) so we annotate correctly
        File assemblyFile = aes.getFile();
        if (assemblyFile == null || !assemblyFile.exists()) {
            if (getContext().getConfig().getVariantCalling().ignoreMissingAssemblyFile) {
                log.error("Missing assembly file. BAN* annotations will be incorrect.");
                return it;
            } else {
                throw new RuntimeException("Missing assembly file " + (assemblyFile == null ? "" : assemblyFile.getName()));
            }
        }
        windowSize = Math.max(windowSize, aes.getMaxAssemblyLength() + 2 * aes.getMaxConcordantFragmentSize());
        // defensive over-eager loading
        windowSize *= 2;
        SamReader reader = getContext().getSamReader(assemblyFile);
        SAMRecordIterator assit = reader.iterator();
        rawAssemblies.add(assit);
        assToClose.add(assit);
        assToClose.add(reader);
    }
    AutoClosingMergedIterator mergedAssemblies = new AutoClosingMergedIterator(rawAssemblies, new SAMRecordCoordinateOnlyComparator());
    return new AutoClosingIterator<>(new AssemblyAssociator(it, mergedAssemblies, windowSize), assToClose.toArray(new Closeable[0]));
}

Example 54 with CloseableIterator

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

the class TestSupport method wc.

public long wc(final Path f) throws IOException {
    if (f.getFileName().toString().endsWith(".bam") || f.getFileName().toString().endsWith(".sam")) {
        final SamReader sr = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.SILENT).open(f);
        final SAMRecordIterator iter = sr.iterator();
        final long n = iter.stream().count();
        iter.close();
        sr.close();
        return n;
    } else if (f.getFileName().toString().endsWith(".vcf") || f.getFileName().toString().endsWith(".vcf.gz")) {
        final VCFReader sr = VCFReaderFactory.makeDefault().open(f, false);
        final CloseableIterator<VariantContext> iter = sr.iterator();
        final long n = iter.stream().count();
        iter.close();
        sr.close();
        return n;
    }
    final BufferedReader br = IOUtils.openPathForBufferedReading(f);
    final long n = br.lines().count();
    br.close();
    return n;
}

Example 55 with CloseableIterator

use of htsjdk.samtools.util.CloseableIterator 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);
    }
}