Examples with Hershey com.github.lindenb.jvarkit.util.hershey.Hershey used on opensource projects

Example 1 with Hershey

use of com.github.lindenb.jvarkit.util.hershey.Hershey in project jvarkit by lindenb.

the class BamMatrix method doWork.

@Override
public int doWork(final List<String> args) {
    if (pixel_size < 1) {
        LOG.error("pixel size is too small (" + this.pixel_size + ")");
        return -1;
    }
    if (StringUtils.isBlank(region2Str)) {
        this.region2Str = region1Str;
    }
    try {
        final SamReaderFactory srf = SamReaderFactory.makeDefault().validationStringency(ValidationStringency.LENIENT);
        if (this.faidx != null)
            srf.referenceSequence(this.faidx);
        final String inputX;
        final String inputY;
        if (args.size() == 1) {
            inputX = args.get(0);
            inputY = null;
        } else if (args.size() == 2) {
            inputX = args.get(0);
            inputY = args.get(1);
        } else {
            LOG.error("illegal number of arguments.");
            return -1;
        }
        this.samReaderX = srf.open(SamInputResource.of(inputX));
        if (!this.samReaderX.hasIndex()) {
            LOG.error("Input " + inputX + " is not indexed");
            return -1;
        }
        this.dict = SequenceDictionaryUtils.extractRequired(this.samReaderX.getFileHeader());
        if (inputY == null) {
            this.samReaderY = srf.open(SamInputResource.of(inputY));
            if (!this.samReaderY.hasIndex()) {
                LOG.error("Input " + inputY + " is not indexed");
                return -1;
            }
            SequenceUtil.assertSequenceDictionariesEqual(SequenceDictionaryUtils.extractRequired(this.samReaderY.getFileHeader()), this.dict);
        } else {
            this.samReaderY = this.samReaderX;
        }
        final ContigNameConverter converter = ContigNameConverter.fromOneDictionary(this.dict);
        final Function<String, Optional<SimpleInterval>> intervalParser = IntervalParserFactory.newInstance().dictionary(dict).enableWholeContig().make();
        this.userIntervalX = intervalParser.apply(this.region1Str).orElseThrow(IntervalParserFactory.exception(this.region1Str));
        this.userIntervalY = intervalParser.apply(this.region2Str).orElseThrow(IntervalParserFactory.exception(this.region2Str));
        // adjust intervals so they have the same length
        if (this.userIntervalX.getLengthOnReference() > this.userIntervalY.getLengthOnReference()) {
            final int mid = this.userIntervalY.getStart() + this.userIntervalY.getLengthOnReference() / 2;
            final int start = Math.max(1, mid - this.userIntervalX.getLengthOnReference() / 2);
            this.userIntervalY = new SimpleInterval(this.userIntervalY.getContig(), start, start + this.userIntervalX.getLengthOnReference());
            LOG.warn("Adjusting interval Y to " + this.userIntervalY + " so both intervals have the same length");
        } else if (this.userIntervalY.getLengthOnReference() > this.userIntervalX.getLengthOnReference()) {
            final int mid = this.userIntervalX.getStart() + this.userIntervalX.getLengthOnReference() / 2;
            final int start = Math.max(1, mid - this.userIntervalY.getLengthOnReference() / 2);
            this.userIntervalX = new SimpleInterval(this.userIntervalX.getContig(), start, start + this.userIntervalY.getLengthOnReference());
            LOG.warn("Adjusting interval X to " + this.userIntervalX + " so both intervals have the same length");
        }
        LOG.info("One pixel is " + (this.userIntervalX.getLengthOnReference() / (double) matrix_size) + " bases");
        final int distance = Math.max(this.userIntervalX.getLengthOnReference(), this.userIntervalY.getLengthOnReference());
        final double pixel2base = distance / (double) matrix_size;
        short max_count = 1;
        final short[] counts = new short[this.matrix_size * this.matrix_size];
        final ReadCounter counter = new MemoryReadCounter();
        /* loop over each pixel 1st axis */
        for (int pixY = 0; pixY < this.matrix_size; pixY++) {
            final int start1 = (int) (this.userIntervalY.getStart() + pixY * pixel2base);
            final int end1 = start1 + (int) pixel2base;
            final Interval qy = new Interval(this.userIntervalY.getContig(), start1, end1);
            if (!qy.overlaps(this.userIntervalY))
                continue;
            final Set<String> set1 = counter.getNamesMatching(1, qy);
            if (set1.isEmpty())
                continue;
            /* loop over each pixel 2nd axis */
            for (int pixX = 0; pixX < this.matrix_size; pixX++) {
                final int start2 = (int) (this.userIntervalX.getStart() + pixX * pixel2base);
                final int end2 = start2 + (int) pixel2base;
                final Interval qx = new Interval(this.userIntervalX.getContig(), start2, end2);
                if (!qx.overlaps(this.userIntervalX))
                    continue;
                if (!validateDisance(qy, qx))
                    continue;
                final int count_common;
                if (qx.compareTo(qy) == 0) {
                    count_common = set1.size();
                } else {
                    final HashSet<String> common = new HashSet<>(set1);
                    common.retainAll(counter.getNamesMatching(0, qx));
                    count_common = common.size();
                }
                final short count = count_common > Short.MAX_VALUE ? Short.MAX_VALUE : (short) count_common;
                max_count = (short) Math.max(count, max_count);
                counts[pixY * this.matrix_size + pixX] = count;
            }
        }
        counter.dispose();
        final int font_size = 10;
        final int cov_height = (this.hide_coverage ? 0 : 50);
        final int gene_height = 25;
        final int margin = font_size + cov_height + (this.gtfPath == null ? 0 : gene_height);
        final Insets margins = new Insets(margin, margin, 10, 10);
        final Dimension drawingAreaDim = new Dimension(this.matrix_size + margins.left + margins.right, this.matrix_size + margins.top + margins.bottom);
        final BufferedImage img = new BufferedImage(drawingAreaDim.width, drawingAreaDim.height, BufferedImage.TYPE_INT_RGB);
        final Graphics2D g = img.createGraphics();
        g.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, drawingAreaDim.width, drawingAreaDim.height);
        // draw sample
        final Hershey herschey = new Hershey();
        final String sampleX = samReaderX.getFileHeader().getReadGroups().stream().map(R -> R.getSample()).filter(S -> !StringUtils.isBlank(S)).findFirst().orElse(inputX);
        final String sampleY = (samReaderX == samReaderY ? sampleX : samReaderX.getFileHeader().getReadGroups().stream().map(R -> R.getSample()).filter(S -> !StringUtils.isBlank(S)).findFirst().orElse(inputY));
        final String sample = (sampleX.equals(sampleY) ? sampleX : String.join(" ", sampleX, sampleY));
        g.setColor(Color.DARK_GRAY);
        herschey.paint(g, sample, new Rectangle2D.Double(0, 1, margins.left - 1, font_size));
        for (int side = 0; side < 2 && !StringUtils.isBlank(this.highlightPath); ++side) {
            final int curr_side = side;
            final SimpleInterval r = (side == 0 ? this.userIntervalX : this.userIntervalY);
            final BedLineCodec bedCodec = new BedLineCodec();
            final Composite oldComposite = g.getComposite();
            g.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.3f));
            try (BufferedReader br = IOUtils.openURIForBufferedReading(this.highlightPath)) {
                br.lines().filter(L -> !(StringUtils.isBlank(L) || L.startsWith("#"))).map(L -> bedCodec.decode(L)).filter(B -> B != null).filter(K -> converter.apply(K.getContig()) != null && r.getContig().equals(converter.apply(K.getContig()))).filter(K -> CoordMath.overlaps(K.getStart(), K.getEnd(), r.getStart(), r.getEnd())).map(E -> new Interval(converter.apply(E.getContig()), E.getStart() + 1, E.getEnd())).filter(E -> CoordMath.overlaps(E.getStart(), E.getEnd(), r.getStart(), r.getEnd())).map(E -> new Interval(E.getContig(), Math.max(r.getStart(), E.getStart()), Math.min(r.getEnd(), E.getEnd()))).forEach(E -> {
                    double d = ((E.getStart() - r.getStart()) / (double) r.getLengthOnReference()) * matrix_size;
                    double dL = ((E.getLengthOnReference()) / (double) r.getLengthOnReference()) * matrix_size;
                    g.setColor(Color.YELLOW);
                    if (curr_side == 0) {
                        g.fill(new Rectangle2D.Double(d, 0, dL, margins.left));
                    } else {
                        g.fill(new Rectangle2D.Double(0, d, margins.top, dL));
                    }
                });
            }
            g.setComposite(oldComposite);
        }
        g.translate(margins.left, margins.top);
        final double logMaxV = Math.log(max_count);
        for (int pix1 = 0; pix1 < this.matrix_size; pix1++) {
            for (int pix2 = 0; pix2 < this.matrix_size; pix2++) {
                final short count = counts[pix1 * this.matrix_size + pix2];
                if (count == 0 || count < this.min_common_names)
                    continue;
                final int gray;
                switch(color_scale) {
                    case LINEAR:
                        gray = 255 - (int) (255 * (count / (double) max_count));
                        break;
                    case LOG:
                        gray = 255 - (int) (255 * ((Math.log(count)) / logMaxV));
                        break;
                    default:
                        throw new IllegalStateException(color_scale.name());
                }
                g.setColor(new Color(gray, 0, 0));
                g.fill(new Rectangle2D.Double(pix1 - pixel_size / 2.0, pix2 - pixel_size / 2.0, pixel_size, pixel_size));
            }
        }
        // draw frame
        g.setColor(Color.GRAY);
        g.drawRect(0, 0, this.matrix_size, this.matrix_size);
        g.translate(-margins.left, -margins.top);
        // used to plot depth
        final double[] coverage = new double[matrix_size];
        final List<SimpleInterval> exonsList;
        if (this.gtfPath == null) {
            exonsList = Collections.emptyList();
        } else {
            try (GtfReader gtfReader = new GtfReader(this.gtfPath)) {
                gtfReader.setContigNameConverter(converter);
                exonsList = gtfReader.getAllGenes().stream().filter(K -> K.overlaps(this.userIntervalX) || K.overlaps(this.userIntervalY)).flatMap(G -> G.getTranscripts().stream()).filter(T -> T.hasExon()).flatMap(K -> K.getExons().stream()).filter(E -> E.overlaps(this.userIntervalX) || E.overlaps(this.userIntervalY)).map(E -> new SimpleInterval(E)).collect(Collectors.toSet()).stream().collect(Collectors.toList());
            }
        }
        for (int side = 0; side < 2; ++side) {
            final SimpleInterval r = (side == 0 ? this.userIntervalX : this.userIntervalY);
            final AffineTransform oldtr = g.getTransform();
            AffineTransform tr;
            if (side == 0) {
                // horizonal axis
                tr = AffineTransform.getTranslateInstance(margins.left, 1);
            } else {
                // vertical
                tr = AffineTransform.getTranslateInstance(margins.left, margins.top);
                tr.concatenate(AffineTransform.getRotateInstance(Math.PI / 2.0));
            }
            g.setTransform(tr);
            // calculate coverage , do this only once if regionX==regionY
            if (!hide_coverage && !(side == 1 && this.userIntervalX.equals(this.userIntervalY))) {
                Arrays.fill(coverage, 0);
                final int[] count = new int[this.matrix_size];
                final IntervalList intervalList = new IntervalList(this.dict);
                intervalList.add(new Interval(r));
                try (final SamLocusIterator sli = new SamLocusIterator(this.samReaderX, intervalList, true)) {
                    while (sli.hasNext()) {
                        final LocusInfo locusInfo = sli.next();
                        final int pos = locusInfo.getPosition();
                        if (pos < r.getStart() || pos > r.getEnd())
                            continue;
                        final int depth = locusInfo.getRecordAndOffsets().size();
                        final int array_index = (int) (((pos - r.getStart()) / (double) r.getLengthOnReference()) * matrix_size);
                        coverage[array_index] += depth;
                        count[array_index]++;
                    }
                }
                for (int i = 0; i < coverage.length; ++i) {
                    if (count[i] == 0)
                        continue;
                    coverage[i] /= count[i];
                }
            }
            // draw ruler
            int y = 0;
            if (!this.hide_coverage) {
                final double max_cov = Arrays.stream(coverage).max().orElse(1);
                final GeneralPath gp = new GeneralPath();
                gp.moveTo(0, cov_height);
                for (int x = 0; x < coverage.length; ++x) {
                    gp.lineTo(x, y + cov_height - (coverage[x] / max_cov) * cov_height);
                }
                gp.lineTo(coverage.length, cov_height);
                gp.closePath();
                g.setColor(Color.GRAY);
                g.fill(gp);
                // string for max cov
                String label = StringUtils.niceInt((int) Arrays.stream(coverage).max().orElse(9));
                g.setColor(Color.DARK_GRAY);
                herschey.paint(g, label, new Rectangle2D.Double(matrix_size - label.length() * font_size, y, label.length() * font_size, font_size));
                y += cov_height;
            }
            // draw label
            g.setColor(Color.DARK_GRAY);
            // label is 'start position'
            String label = StringUtils.niceInt(r.getStart());
            herschey.paint(g, label, new Rectangle2D.Double(0, y, label.length() * font_size, font_size));
            // label is 'end position'
            label = StringUtils.niceInt(r.getEnd());
            herschey.paint(g, label, new Rectangle2D.Double(matrix_size - (label.length() * font_size), y, label.length() * font_size, font_size));
            // label is 'chromosome and length'
            label = r.getContig() + " ( " + StringUtils.niceInt(r.getLengthOnReference()) + " bp )";
            herschey.paint(g, label, new Rectangle2D.Double(matrix_size / 2.0 - (label.length() * font_size) / 2.0, y, label.length() * font_size, font_size));
            y += font_size;
            // draw genes
            if (this.gtfPath != null) {
                final double curr_y = y;
                double midy = y + gene_height / 2.0;
                g.setColor(Color.CYAN);
                g.draw(new Line2D.Double(0, midy, matrix_size, midy));
                exonsList.stream().filter(E -> E.overlaps(r)).map(E -> new SimpleInterval(E.getContig(), Math.max(r.getStart(), E.getStart()), Math.min(r.getEnd(), E.getEnd()))).forEach(E -> {
                    final double x = ((E.getStart() - r.getStart()) / (double) r.getLengthOnReference()) * matrix_size;
                    final double width = ((E.getLengthOnReference()) / (double) r.getLengthOnReference()) * matrix_size;
                    g.setColor(Color.BLUE);
                    g.fill(new Rectangle2D.Double(x, curr_y, width, gene_height));
                });
            }
            g.setTransform(oldtr);
        }
        g.dispose();
        try {
            if (this.outputFile == null) {
                ImageIO.write(img, "PNG", stdout());
            } else {
                ImageIO.write(img, this.outputFile.getName().endsWith(".png") ? "PNG" : "JPG", this.outputFile);
            }
        } catch (final IOException err) {
            throw new RuntimeIOException(err);
        }
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(this.samReaderX);
        CloserUtil.close(this.samReaderY);
    }
}

Example 2 with Hershey

use of com.github.lindenb.jvarkit.util.hershey.Hershey in project jvarkit by lindenb.

the class BamCmpCoverage method doWork.

@Override
public int doWork(final List<String> args) {
    if (outputFile == null) {
        LOG.error("output image file not defined");
        return -1;
    }
    if (this.imgageSize < 1) {
        LOG.error("Bad image size:" + this.imgageSize);
        return -1;
    }
    if (this.minDepth < 0) {
        LOG.error("Bad min depth : " + this.minDepth);
        return -1;
    }
    if (this.minDepth >= this.maxDepth) {
        LOG.error("Bad min<max depth : " + this.minDepth + "<" + this.maxDepth);
        return 1;
    }
    if (this.bedFile != null) {
        readBedFile(this.bedFile);
    }
    if (regionStr != null && this.intervals != null) {
        LOG.error("bed and interval both defined.");
        return -1;
    }
    try {
        final ConcatSam.Factory concatSamFactory = new ConcatSam.Factory();
        final SamReaderFactory srf = concatSamFactory.getSamReaderFactory();
        srf.disable(SamReaderFactory.Option.EAGERLY_DECODE);
        srf.disable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS);
        srf.disable(SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS);
        if (this.regionStr != null) {
            concatSamFactory.addInterval(this.regionStr);
        }
        ConcatSam.ConcatSamIterator concatIter = concatSamFactory.open(args);
        final SAMSequenceDictionary dict = concatIter.getFileHeader().getSequenceDictionary();
        final Set<String> samples = concatIter.getFileHeader().getReadGroups().stream().map(RG -> this.samRecordPartition.apply(RG, "N/A")).collect(Collectors.toSet());
        LOG.info("Samples:" + samples.size());
        for (String sample : samples) {
            this.sample2column.put(sample, this.sample2column.size());
        }
        // create image
        LOG.info("Creating image " + this.imgageSize + "x" + this.imgageSize);
        this.image = new BufferedImage(this.imgageSize, this.imgageSize, BufferedImage.TYPE_INT_RGB);
        Graphics2D g = this.image.createGraphics();
        this.marginWidth = this.imgageSize * 0.05;
        double drawingWidth = (this.imgageSize - 1) - marginWidth;
        this.sampleWidth = drawingWidth / samples.size();
        // g.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, this.imgageSize, this.imgageSize);
        g.setColor(Color.BLACK);
        Hershey hershey = new Hershey();
        for (final String sample_x : samples) {
            double labelHeight = marginWidth;
            if (labelHeight > 50)
                labelHeight = 50;
            g.setColor(Color.BLACK);
            hershey.paint(g, sample_x, marginWidth + sample2column.get(sample_x) * sampleWidth, marginWidth - labelHeight, sampleWidth * 0.9, labelHeight * 0.9);
            AffineTransform old = g.getTransform();
            AffineTransform tr = AffineTransform.getTranslateInstance(marginWidth, marginWidth + sample2column.get(sample_x) * sampleWidth);
            tr.rotate(Math.PI / 2);
            g.setTransform(tr);
            hershey.paint(g, sample_x, 0.0, 0.0, sampleWidth * 0.9, labelHeight * 0.9);
            // g.drawString(this.tabixFile.getFile().getName(),0,0);
            g.setTransform(old);
            for (String sample_y : samples) {
                Rectangle2D rect = new Rectangle2D.Double(marginWidth + sample2column.get(sample_x) * sampleWidth, marginWidth + sample2column.get(sample_y) * sampleWidth, sampleWidth, sampleWidth);
                g.setColor(Color.BLUE);
                g.draw(new Line2D.Double(rect.getMinX(), rect.getMinY(), rect.getMaxX(), rect.getMaxY()));
                g.setColor(Color.BLACK);
                g.draw(rect);
            }
        }
        // ceate bit-array
        BitSampleMatrix bitMatrix = new BitSampleMatrix(samples.size());
        // preivous chrom
        // int prev_tid=-1;
        BufferedList<Depth> depthList = new BufferedList<Depth>();
        g.setColor(Color.BLACK);
        SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(dict).logger(LOG);
        LOG.info("Scanning bams...");
        while (concatIter.hasNext()) {
            final SAMRecord rec = progress.watch(concatIter.next());
            if (this.samRecordFilter.filterOut(rec))
                continue;
            final String sample = this.samRecordPartition.getPartion(rec, "N/A");
            final int sample_id = this.sample2column.get(sample);
            final Cigar cigar = rec.getCigar();
            if (cigar == null || cigar.isEmpty())
                continue;
            int refPos = rec.getAlignmentStart();
            /* cleanup front pos */
            while (!depthList.isEmpty()) {
                final Depth front = depthList.getFirst();
                if (front.tid != rec.getReferenceIndex().intValue() || front.pos < refPos) {
                    paint(bitMatrix, front);
                    depthList.removeFirst();
                    continue;
                } else {
                    break;
                }
            }
            for (final CigarElement ce : cigar.getCigarElements()) {
                final CigarOperator op = ce.getOperator();
                if (!op.consumesReferenceBases())
                    continue;
                if (op.consumesReadBases()) {
                    for (int i = 0; i < ce.getLength(); ++i) {
                        Depth depth = null;
                        int pos = refPos + i;
                        // ignore non-overlapping BED
                        if (this.intervals != null && !this.intervals.containsOverlapping(new Interval(rec.getReferenceName(), pos, pos))) {
                            continue;
                        } else if (depthList.isEmpty()) {
                            depth = new Depth();
                            depth.pos = pos;
                            depth.tid = rec.getReferenceIndex();
                            depthList.add(depth);
                        } else if (depthList.getLast().pos < pos) {
                            Depth prev = depthList.getLast();
                            while (prev.pos < pos) {
                                depth = new Depth();
                                depth.pos = prev.pos + 1;
                                depth.tid = rec.getReferenceIndex();
                                depthList.add(depth);
                                prev = depth;
                            }
                            depth = prev;
                        } else {
                            int lastPos = depthList.get(depthList.size() - 1).pos;
                            int distance = lastPos - pos;
                            int indexInList = (depthList.size() - 1) - (distance);
                            if (indexInList < 0) {
                                // can appen when BED declared and partially overlap the read
                                continue;
                            }
                            depth = depthList.get((depthList.size() - 1) - (distance));
                            if (depth.pos != pos) {
                                LOG.error(" " + pos + " vs " + depth.pos + " " + lastPos);
                                return -1;
                            }
                        }
                        depth.depths[sample_id]++;
                    }
                }
                refPos += ce.getLength();
            }
        }
        while (!depthList.isEmpty()) {
            // paint(g,depthList.remove(0));
            paint(bitMatrix, depthList.remove(0));
        }
        progress.finish();
        concatIter.close();
        concatIter = null;
        for (int x = 0; x < bitMatrix.n_samples; ++x) {
            for (int y = 0; y < bitMatrix.n_samples; ++y) {
                LOG.info("Painting...(" + x + "/" + y + ")");
                paint(g, bitMatrix.get(x, y));
            }
        }
        g.dispose();
        // save file
        LOG.info("saving " + this.outputFile);
        if (this.outputFile.getName().toLowerCase().endsWith(".png")) {
            ImageIO.write(this.image, "PNG", this.outputFile);
        } else {
            ImageIO.write(this.image, "JPG", this.outputFile);
        }
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
    }
}

Example 3 with Hershey

use of com.github.lindenb.jvarkit.util.hershey.Hershey in project jvarkit by lindenb.

the class VcfToHilbert method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.imgOut == null) {
        LOG.error("output image file not defined");
        return -1;
    }
    if (this.imageWidth < 1) {
        LOG.error("Bad image size:" + this.imageWidth);
        return -1;
    }
    VCFIterator iter = null;
    try {
        iter = this.openVCFIterator(oneFileOrNull(args));
        final VCFHeader header = iter.getHeader();
        this.dict = header.getSequenceDictionary();
        if (this.dict == null) {
            throw new JvarkitException.FastaDictionaryMissing("no dict in input");
        }
        final List<String> samples = header.getSampleNamesInOrder();
        if (samples.isEmpty()) {
            throw new JvarkitException.SampleMissing("no.sample.in.vcf");
        }
        LOG.info("N-Samples:" + samples.size());
        double marginWidth = (this.imageWidth - 2) * 0.05;
        this.sampleWidth = ((this.imageWidth - 2) - marginWidth) / samples.size();
        LOG.info("sample Width:" + sampleWidth);
        BufferedImage img = new BufferedImage(this.imageWidth, this.imageWidth, BufferedImage.TYPE_INT_RGB);
        this.g = (Graphics2D) img.getGraphics();
        this.g.setColor(Color.WHITE);
        this.g.fillRect(0, 0, imageWidth, imageWidth);
        g.setColor(Color.BLACK);
        final Hershey hershey = new Hershey();
        EvalCurve evalCurve = new EvalCurve();
        evalCurve.run();
        this.genomicSizePerCurveUnit = ((double) dict.getReferenceLength() / (double) (evalCurve.count));
        if (this.genomicSizePerCurveUnit < 1)
            this.genomicSizePerCurveUnit = 1;
        LOG.info("genomicSizePerCurveUnit:" + genomicSizePerCurveUnit);
        for (int x = 0; x < samples.size(); ++x) {
            String samplex = samples.get(x);
            double labelHeight = marginWidth;
            if (labelHeight > 50)
                labelHeight = 50;
            g.setColor(Color.BLACK);
            hershey.paint(g, samplex, marginWidth + x * sampleWidth, marginWidth - labelHeight, sampleWidth * 0.9, labelHeight * 0.9);
            AffineTransform old = g.getTransform();
            AffineTransform tr = AffineTransform.getTranslateInstance(marginWidth, marginWidth + x * sampleWidth);
            tr.rotate(Math.PI / 2);
            g.setTransform(tr);
            hershey.paint(g, samplex, 0.0, 0.0, sampleWidth * 0.9, labelHeight * 0.9);
            // g.drawString(this.tabixFile.getFile().getName(),0,0);
            g.setTransform(old);
            double tx = marginWidth + x * sampleWidth;
            for (int y = 0; y < samples.size(); ++y) {
                double ty = marginWidth + y * sampleWidth;
                g.translate(tx, ty);
                g.setColor(Color.BLUE);
                g.draw(new Rectangle2D.Double(0, 0, sampleWidth, sampleWidth));
                // paint each chromosome
                paintReference();
                g.translate(-tx, -ty);
            }
        }
        LOG.info("genomicSizePerCurveUnit:" + (long) genomicSizePerCurveUnit * evalCurve.count + " " + dict.getReferenceLength() + " count=" + evalCurve.count);
        LOG.info("Scanning variants");
        final SAMSequenceDictionaryProgress progress = new SAMSequenceDictionaryProgress(this.dict);
        while (iter.hasNext()) {
            final VariantContext var = progress.watch(iter.next());
            for (int x = 0; x < samples.size(); ++x) {
                final String samplex = samples.get(x);
                final Genotype gx = var.getGenotype(samplex);
                if (!gx.isCalled())
                    continue;
                double tx = marginWidth + x * sampleWidth;
                for (int y = 0; y < samples.size(); ++y) {
                    final String sampley = samples.get(y);
                    final Genotype gy = var.getGenotype(sampley);
                    if (!gy.isCalled())
                        continue;
                    if (gx.isHomRef() && gy.isHomRef())
                        continue;
                    double ty = marginWidth + y * sampleWidth;
                    g.translate(tx, ty);
                    final PaintVariant paint = new PaintVariant(var, x, y);
                    paint.run();
                    g.translate(-tx, -ty);
                }
            }
        }
        progress.finish();
        this.g.dispose();
        // save file
        LOG.info("saving " + imgOut);
        if (imgOut != null) {
            ImageIO.write(img, imgOut.getName().toLowerCase().endsWith(".png") ? "PNG" : "JPG", imgOut);
        } else {
            ImageIO.write(img, "PNG", stdout());
        }
        return 0;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(iter);
    }
}