Search in sources :

Example 1 with Sample

use of com.github.lindenb.jvarkit.pedigree.Sample in project jvarkit by lindenb.

the class SwingVCFGenotypesTableModel method setRows.

@Override
public synchronized void setRows(final List<Genotype> genotypes) {
    if (genotypes == null) {
        setRows(Collections.emptyList());
    }
    this.columns.clear();
    final Vector<ColumnInfo> columns = new Vector<>();
    ColumnInfo ci = new ColumnInfo();
    ci.clazz = String.class;
    ci.name = "Sample";
    ci.extractor = GT -> GT.getSampleName();
    columns.add(ci);
    if (this.pedigree != null) {
        ci = new ColumnInfo();
        ci.clazz = String.class;
        ci.name = "Status";
        ci.extractor = GT -> {
            final Sample sn = this.pedigree.getSampleById(GT.getSampleName());
            if (sn == null || !sn.isStatusSet())
                return null;
            if (sn.getStatus().equals(Status.affected))
                return "[*]";
            if (sn.getStatus().equals(Status.unaffected))
                return "[ ]";
            return null;
        };
        columns.add(ci);
    }
    if (genotypes.stream().anyMatch(G -> G.isAvailable())) {
        ci = new ColumnInfo();
        ci.clazz = String.class;
        ci.name = VCFConstants.GENOTYPE_KEY;
        ci.extractor = GT -> GT.getAlleles().stream().map(A -> A.getDisplayString()).collect(Collectors.joining(GT.isPhased() ? Genotype.PHASED_ALLELE_SEPARATOR : Genotype.UNPHASED_ALLELE_SEPARATOR));
        columns.add(ci);
    }
    ci = new ColumnInfo();
    ci.clazz = String.class;
    ci.name = "Type";
    ci.extractor = GT -> GT.getType().name();
    columns.add(ci);
    if (genotypes.stream().anyMatch(G -> G.hasDP())) {
        ci = new ColumnInfo();
        ci.clazz = String.class;
        ci.name = VCFConstants.DEPTH_KEY;
        ci.extractor = GT -> GT.hasDP() ? GT.getDP() : null;
        columns.add(ci);
    }
    if (genotypes.stream().anyMatch(G -> G.hasGQ())) {
        ci = new ColumnInfo();
        ci.clazz = String.class;
        ci.name = VCFConstants.GENOTYPE_QUALITY_KEY;
        ci.extractor = GT -> GT.hasGQ() ? GT.getGQ() : null;
        columns.add(ci);
    }
    if (genotypes.stream().anyMatch(G -> G.hasAD())) {
        ci = new ColumnInfo();
        ci.clazz = String.class;
        ci.name = VCFConstants.GENOTYPE_ALLELE_DEPTHS;
        ci.extractor = GT -> GT.hasAD() ? Arrays.stream(GT.getAD()).mapToObj(P -> String.valueOf(P)).collect(Collectors.joining(",")) : null;
        columns.add(ci);
    }
    if (genotypes.stream().anyMatch(G -> G.hasPL())) {
        ci = new ColumnInfo();
        ci.clazz = String.class;
        ci.name = VCFConstants.GENOTYPE_PL_KEY;
        ci.extractor = GT -> GT.hasPL() ? Arrays.stream(GT.getPL()).mapToObj(P -> String.valueOf(P)).collect(Collectors.joining(",")) : null;
        columns.add(ci);
    }
    for (final String att : genotypes.stream().flatMap(G -> G.getExtendedAttributes().keySet().stream()).collect(Collectors.toSet())) {
        ci = new ColumnInfo();
        ci.clazz = Object.class;
        ci.name = att;
        ci.extractor = GT -> GT.hasExtendedAttribute(att) ? GT.getExtendedAttribute(att) : null;
        columns.add(ci);
    }
    this.columns.addAll(columns);
    super.rows.clear();
    super.rows.addAll(genotypes);
    fireTableStructureChanged();
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) AbstractGenericTable(com.github.lindenb.jvarkit.util.swing.AbstractGenericTable) Arrays(java.util.Arrays) Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree) Status(com.github.lindenb.jvarkit.pedigree.Status) Function(java.util.function.Function) Collectors(java.util.stream.Collectors) List(java.util.List) Vector(java.util.Vector) VariantContext(htsjdk.variant.variantcontext.VariantContext) Collections(java.util.Collections) Sample(com.github.lindenb.jvarkit.pedigree.Sample) VCFConstants(htsjdk.variant.vcf.VCFConstants) Sample(com.github.lindenb.jvarkit.pedigree.Sample) Vector(java.util.Vector)

Example 2 with Sample

use of com.github.lindenb.jvarkit.pedigree.Sample in project jvarkit by lindenb.

the class CoverageServer method printRaster.

/**
 * print BAM for small interval, displaying reads
 */
private void printRaster(final BamInput bam, final SimpleInterval midRegion, final SimpleInterval region, final HttpServletRequest request, final HttpServletResponse response) throws IOException, ServletException {
    final IntToDoubleFunction position2pixel = X -> ((X - region.getStart()) / (double) region.getLengthOnReference()) * (double) image_width;
    final SamReaderFactory srf = SamReaderFactory.make().validationStringency(ValidationStringency.LENIENT).referenceSequence(this.faidxRef);
    final Pileup<SAMRecord> pileup = new Pileup<>((L, R) -> position2pixel.applyAsDouble(L.getUnclippedEnd() + 1) + 1 < position2pixel.applyAsDouble(R.getUnclippedStart()));
    try (SamReader sr = srf.open(bam.bamPath)) {
        try (CloseableIterator<SAMRecord> iter = sr.query(region.getContig(), // extend to get clipR
        Math.max(0, region.getStart() - this.small_region_size), region.getEnd() + this.small_region_size, false)) {
            while (iter.hasNext()) {
                final SAMRecord rec = iter.next();
                if (!acceptRead(rec))
                    continue;
                if (rec.getUnclippedEnd() < region.getStart())
                    continue;
                if (rec.getUnclippedStart() > region.getEnd())
                    continue;
                final Cigar cigar = rec.getCigar();
                if (cigar == null || cigar.isEmpty())
                    continue;
                pileup.add(rec);
            }
        }
    // end iterator
    }
    // end samreader
    ReferenceSequence refInInterval = null;
    try (ReferenceSequenceFile refseq = ReferenceSequenceFileFactory.getReferenceSequenceFile(this.faidxRef)) {
        final SAMSequenceRecord ssr = this.dictionary.getSequence(region.getContig());
        if (region.getStart() <= ssr.getSequenceLength()) {
            refInInterval = refseq.getSubsequenceAt(region.getContig(), region.getStart(), Math.min(region.getEnd(), ssr.getSequenceLength()));
        }
    }
    final BufferedImage img = new BufferedImage(image_width, image_height, BufferedImage.TYPE_INT_RGB);
    final Graphics2D g = img.createGraphics();
    g.setRenderingHint(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);
    g.setColor(new Color(240, 240, 240));
    g.fillRect(0, 0, image_width + 1, image_height + 1);
    final Stroke oldStroke = g.getStroke();
    g.setStroke(new BasicStroke(0.5f));
    g.setColor(Color.WHITE);
    final double mid_start = position2pixel.applyAsDouble(midRegion.getStart());
    final double mid_end = position2pixel.applyAsDouble(midRegion.getEnd() + 1);
    g.fill(new Rectangle2D.Double(mid_start, 0, (mid_end - mid_start), this.image_height));
    final int[] int_coverage = new int[region.getLengthOnReference()];
    final int margin_top = 12;
    final double featureHeight = Math.min(20, (this.image_height - margin_top) / Math.max(1.0, (double) pileup.getRowCount()));
    double y = image_height - featureHeight;
    for (final List<SAMRecord> row : pileup) {
        final double h2 = Math.min(featureHeight * 0.9, featureHeight - 2);
        for (final SAMRecord rec : row) {
            final Cigar cigar = rec.getCigar();
            if (cigar == null || cigar.isEmpty())
                continue;
            /* fill coverage array  */
            int ref1 = rec.getAlignmentStart();
            for (CigarElement ce : cigar.getCigarElements()) {
                if (ref1 > region.getEnd())
                    break;
                final CigarOperator op = ce.getOperator();
                if (op.consumesReferenceBases()) {
                    if (op.consumesReadBases()) {
                        for (int x = 0; x < ce.getLength(); ++x) {
                            int pos = ref1 + x;
                            if (pos < region.getStart())
                                continue;
                            if (pos > region.getEnd())
                                break;
                            int_coverage[pos - region.getStart()]++;
                        }
                    }
                    ref1 += ce.getLength();
                }
            }
            /* draw rec itself */
            final double midy = y + h2 / 2.0;
            g.setColor(Color.DARK_GRAY);
            g.draw(new Line2D.Double(position2pixel.applyAsDouble(rec.getUnclippedStart()), midy, position2pixel.applyAsDouble(rec.getUnclippedEnd()), midy));
            ref1 = rec.getUnclippedStart();
            final List<Double> insertions = new ArrayList<>();
            for (final CigarElement ce : cigar.getCigarElements()) {
                if (ref1 > region.getEnd())
                    break;
                final CigarOperator op = ce.getOperator();
                Shape shape = null;
                Color fill = null;
                Color stroke = Color.DARK_GRAY;
                switch(op) {
                    case P:
                        break;
                    // through
                    case M:
                    // through
                    case X:
                    // through
                    case EQ:
                    // through
                    case S:
                    case H:
                        final double x1 = position2pixel.applyAsDouble(ref1);
                        shape = new Rectangle2D.Double(x1, y, position2pixel.applyAsDouble(ref1 + ce.getLength()) - x1, h2);
                        ref1 += ce.getLength();
                        switch(op) {
                            case H:
                            case S:
                                fill = Color.YELLOW;
                                break;
                            case X:
                                fill = Color.RED;
                                break;
                            case EQ:
                            case M:
                                fill = Color.LIGHT_GRAY;
                                break;
                            default:
                                break;
                        }
                        break;
                    // through
                    case N:
                    case D:
                        shape = null;
                        fill = null;
                        stroke = null;
                        ref1 += ce.getLength();
                        break;
                    case I:
                        shape = null;
                        fill = null;
                        stroke = null;
                        insertions.add(position2pixel.applyAsDouble(ref1));
                        break;
                    default:
                        throw new IllegalStateException("" + op);
                }
                if (ref1 < region.getStart())
                    continue;
                if (shape != null) {
                    if (fill != null) {
                        g.setColor(fill);
                        g.fill(shape);
                    }
                    if (stroke != null && h2 > 4) {
                        g.setColor(stroke);
                        g.draw(shape);
                    }
                }
            }
            /* draw mismatched bases */
            if (refInInterval != null && rec.getReadBases() != null && rec.getReadBases() != SAMRecord.NULL_SEQUENCE) {
                final byte[] bases = rec.getReadBases();
                final IntFunction<Character> baseRead = IDX -> IDX < 0 || IDX >= bases.length || bases == SAMRecord.NULL_SEQUENCE ? 'N' : (char) Character.toUpperCase(bases[IDX]);
                int read0 = 0;
                ref1 = rec.getAlignmentStart();
                for (CigarElement ce : cigar.getCigarElements()) {
                    if (ref1 > region.getEnd())
                        break;
                    final CigarOperator op = ce.getOperator();
                    switch(op) {
                        case P:
                            break;
                        case H:
                            break;
                        case D:
                        case N:
                            ref1 += ce.getLength();
                            break;
                        case S:
                        case I:
                            read0 += ce.getLength();
                            break;
                        case EQ:
                        case M:
                        case X:
                            {
                                for (int j = 0; j < ce.getLength(); j++) {
                                    if (ref1 + j < region.getStart())
                                        continue;
                                    if (ref1 + j >= region.getStart() + refInInterval.length())
                                        break;
                                    final int ref_base_idx = ref1 - region.getStart() + j;
                                    char ctgBase = (char) (ref_base_idx < 0 || ref_base_idx >= refInInterval.length() ? 'N' : Character.toUpperCase(refInInterval.getBases()[ref_base_idx]));
                                    if (ctgBase == 'N')
                                        continue;
                                    char readBase = baseRead.apply(read0 + j);
                                    if (readBase == 'N')
                                        continue;
                                    if (readBase == ctgBase)
                                        continue;
                                    g.setColor(Color.ORANGE);
                                    final double x1 = position2pixel.applyAsDouble(ref1 + j);
                                    final double x2 = position2pixel.applyAsDouble(ref1 + j + 1);
                                    g.fill(new Rectangle2D.Double(x1, y, x2 - x1, h2));
                                }
                                read0 += ce.getLength();
                                ref1 += ce.getLength();
                                break;
                            }
                        default:
                            break;
                    }
                }
            }
            for (double px : insertions) {
                g.setColor(Color.RED);
                g.draw(new Line2D.Double(px, y - 0.5, px, y + h2 + 0.5));
            }
        }
        y -= featureHeight;
    }
    double max_cov = IntStream.of(int_coverage).max().orElse(0);
    g.setColor(Color.DARK_GRAY);
    g.drawString("Sample:" + bam.sample + " max-cov:" + (int) max_cov + " " + region.toNiceString(), 10, 10);
    /* plot coverage */
    final GeneralPath gp = new GeneralPath();
    for (int i = 0; max_cov > 0 && i < int_coverage.length; ++i) {
        final double x1 = position2pixel.applyAsDouble(region.getStart() + i);
        final double x2 = position2pixel.applyAsDouble(region.getStart() + i + 1);
        final double y1 = image_height - (int_coverage[i] / max_cov) * (image_height - margin_top);
        if (i == 0)
            gp.moveTo(x1, y1);
        else
            gp.lineTo(x1, y1);
        gp.lineTo(x2, y1);
    }
    g.setStroke(new BasicStroke(0.5f, BasicStroke.CAP_BUTT, BasicStroke.JOIN_BEVEL, 0f, new float[] { 1f, 2f, 1f }, 0f));
    g.setColor(Color.BLUE);
    g.draw(gp);
    g.setStroke(oldStroke);
    writeGenes(g, region);
    writeKnownCnv(g, region);
    g.setColor(Color.PINK);
    g.draw(new Line2D.Double(mid_start, 0, mid_start, image_height));
    g.draw(new Line2D.Double(mid_end, 0, mid_end, image_height));
    writeImage(img, bam, region, response);
}
Also used : Color(java.awt.Color) ServletContextHandler(org.eclipse.jetty.servlet.ServletContextHandler) Arrays(java.util.Arrays) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) Program(com.github.lindenb.jvarkit.util.jcommander.Program) ServletException(javax.servlet.ServletException) Rectangle2D(java.awt.geom.Rectangle2D) CigarElement(htsjdk.samtools.CigarElement) CigarOperator(htsjdk.samtools.CigarOperator) RenderingHints(java.awt.RenderingHints) IntervalListProvider(com.github.lindenb.jvarkit.samtools.util.IntervalListProvider) IntervalParserFactory(com.github.lindenb.jvarkit.samtools.util.IntervalParserFactory) SAMFileHeader(htsjdk.samtools.SAMFileHeader) ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile) Hyperlink(com.github.lindenb.jvarkit.net.Hyperlink) DistanceParser(com.github.lindenb.jvarkit.util.bio.DistanceParser) Vector(java.util.Vector) XMLStreamException(javax.xml.stream.XMLStreamException) ImageIO(javax.imageio.ImageIO) XMLStreamWriter(javax.xml.stream.XMLStreamWriter) HttpStatus(org.eclipse.jetty.http.HttpStatus) Path(java.nio.file.Path) Server(org.eclipse.jetty.server.Server) CloserUtil(htsjdk.samtools.util.CloserUtil) Shape(java.awt.Shape) PrintWriter(java.io.PrintWriter) SimpleInterval(com.github.lindenb.jvarkit.samtools.util.SimpleInterval) SequenceDictionaryUtils(com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils) HttpServlet(javax.servlet.http.HttpServlet) Composite(java.awt.Composite) BufferedImage(java.awt.image.BufferedImage) Logger(com.github.lindenb.jvarkit.util.log.Logger) StandardOpenOption(java.nio.file.StandardOpenOption) ServletHolder(org.eclipse.jetty.servlet.ServletHolder) Set(java.util.Set) Collectors(java.util.stream.Collectors) GTFCodec(com.github.lindenb.jvarkit.util.bio.gtf.GTFCodec) SAMRecord(htsjdk.samtools.SAMRecord) DoubleStream(java.util.stream.DoubleStream) ReferenceSequenceFileFactory(htsjdk.samtools.reference.ReferenceSequenceFileFactory) List(java.util.List) Stream(java.util.stream.Stream) StringUtils(com.github.lindenb.jvarkit.lang.StringUtils) FileExtensions(htsjdk.samtools.util.FileExtensions) ToDoubleFunction(java.util.function.ToDoubleFunction) BasicStroke(java.awt.BasicStroke) ReferenceSequence(htsjdk.samtools.reference.ReferenceSequence) GeneralPath(java.awt.geom.GeneralPath) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) SamReaderFactory(htsjdk.samtools.SamReaderFactory) IntStream(java.util.stream.IntStream) Cigar(htsjdk.samtools.Cigar) CloseableIterator(htsjdk.samtools.util.CloseableIterator) SequenceUtil(htsjdk.samtools.util.SequenceUtil) ContigNameConverter(com.github.lindenb.jvarkit.util.bio.fasta.ContigNameConverter) Parameter(com.beust.jcommander.Parameter) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) IterableAdapter(htsjdk.samtools.util.IterableAdapter) ValidationStringency(htsjdk.samtools.ValidationStringency) ArrayList(java.util.ArrayList) GTFLine(com.github.lindenb.jvarkit.util.bio.gtf.GTFLine) Interval(htsjdk.samtools.util.Interval) AlphaComposite(java.awt.AlphaComposite) HttpServletRequest(javax.servlet.http.HttpServletRequest) NoSplitter(com.github.lindenb.jvarkit.util.jcommander.NoSplitter) Graphics2D(java.awt.Graphics2D) PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) TabixReader(htsjdk.tribble.readers.TabixReader) VCFReaderFactory(com.github.lindenb.jvarkit.variant.vcf.VCFReaderFactory) StreamSupport(java.util.stream.StreamSupport) IntFunction(java.util.function.IntFunction) VCFConstants(htsjdk.variant.vcf.VCFConstants) Stroke(java.awt.Stroke) Line2D(java.awt.geom.Line2D) Counter(com.github.lindenb.jvarkit.util.Counter) AbstractIterator(htsjdk.samtools.util.AbstractIterator) Locatable(htsjdk.samtools.util.Locatable) IntToDoubleFunction(java.util.function.IntToDoubleFunction) Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree) Files(java.nio.file.Files) BufferedWriter(java.io.BufferedWriter) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) HttpServletResponse(javax.servlet.http.HttpServletResponse) VCFReader(htsjdk.variant.vcf.VCFReader) IOException(java.io.IOException) JVarkitVersion(com.github.lindenb.jvarkit.util.JVarkitVersion) SamReader(htsjdk.samtools.SamReader) XMLOutputFactory(javax.xml.stream.XMLOutputFactory) Pileup(com.github.lindenb.jvarkit.samtools.util.Pileup) BufferedReader(java.io.BufferedReader) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) Sample(com.github.lindenb.jvarkit.pedigree.Sample) BasicStroke(java.awt.BasicStroke) Shape(java.awt.Shape) GeneralPath(java.awt.geom.GeneralPath) ArrayList(java.util.ArrayList) IntToDoubleFunction(java.util.function.IntToDoubleFunction) ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) ReferenceSequence(htsjdk.samtools.reference.ReferenceSequence) Line2D(java.awt.geom.Line2D) BufferedImage(java.awt.image.BufferedImage) Pileup(com.github.lindenb.jvarkit.samtools.util.Pileup) SamReader(htsjdk.samtools.SamReader) BasicStroke(java.awt.BasicStroke) Stroke(java.awt.Stroke) SamReaderFactory(htsjdk.samtools.SamReaderFactory) Color(java.awt.Color) Rectangle2D(java.awt.geom.Rectangle2D) CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement) Graphics2D(java.awt.Graphics2D) Cigar(htsjdk.samtools.Cigar) SAMRecord(htsjdk.samtools.SAMRecord)

Example 3 with Sample

use of com.github.lindenb.jvarkit.pedigree.Sample in project jvarkit by lindenb.

the class CoverageServer method printPage.

/**
 * print HTML page
 */
private void printPage(final HttpServletRequest request, final HttpServletResponse response) throws IOException, ServletException {
    String message = "";
    final String charset = StringUtils.ifBlank(request.getCharacterEncoding(), "UTF-8");
    response.setContentType("text/html; charset=" + charset.toLowerCase());
    response.setCharacterEncoding(charset);
    PrintWriter pw = response.getWriter();
    SimpleInterval interval = null;
    int columns_count = this.images_per_row;
    if (!StringUtils.isBlank(request.getParameter("columns"))) {
        columns_count = Math.max(1, StringUtils.parseInt(request.getParameter("columns")).orElse(this.images_per_row));
    }
    final boolean normalize = request.getParameter("normalize") != null;
    if (!StringUtils.isBlank(request.getParameter("custominterval"))) {
        interval = parseInterval(request.getParameter("custominterval"));
        if (interval == null) {
            message += " Cannot parse user interval '" + request.getParameter("custominterval") + "'.";
        }
    }
    if (interval == null && !StringUtils.isBlank(request.getParameter("interval"))) {
        interval = parseInterval(request.getParameter("interval"));
        if (interval == null) {
            message += " Cannot parse interval. using default: " + interval + ".";
        }
    }
    if (interval == null && !this.named_intervals.isEmpty()) {
        /* first non reviewed */
        interval = this.named_intervals.stream().filter(R -> !R.reviewed).findFirst().map(R -> new SimpleInterval(R)).orElse(null);
        /* all reviewed ? */
        if (interval == null)
            interval = new SimpleInterval(this.named_intervals.get(0));
    }
    /* still no interval ?*/
    if (interval == null) {
        final SAMSequenceRecord ssr = this.dictionary.getSequence(0);
        interval = new SimpleInterval(ssr.getSequenceName(), 1, Math.min(ssr.getSequenceLength(), 100));
    }
    if (!StringUtils.isBlank(request.getParameter("move"))) {
        final String value = request.getParameter("move");
        double factor;
        switch(value.length()) {
            case 0:
                factor = 0;
                break;
            case 1:
                factor = 0.1;
                break;
            case 2:
                factor = 0.475;
                break;
            default:
                factor = 0.95;
                break;
        }
        int length0 = interval.getLengthOnReference();
        int length1 = (int) (length0 * factor);
        int shift = (value.startsWith(">") ? 1 : -1) * length1;
        int start = interval.getStart() + shift;
        int end = interval.getEnd() + shift;
        if (start < 1)
            start = 1;
        final SAMSequenceRecord ssr = this.dictionary.getSequence(interval.getContig());
        if (ssr != null && end >= ssr.getSequenceLength()) {
            end = ssr.getSequenceLength();
        }
        if (ssr != null && start >= ssr.getSequenceLength()) {
            start = ssr.getSequenceLength();
        }
        interval = new SimpleInterval(interval.getContig(), start, end);
    }
    /* ZOOM buttons */
    for (int side = 0; side < 2; ++side) {
        final String param = "zoom" + (side == 0 ? "in" : "out");
        String value = request.getParameter(param);
        if (StringUtils.isBlank(value))
            continue;
        double factor = Double.parseDouble(value);
        if (side == 0)
            factor = 1.0 / factor;
        int length0 = interval.getLengthOnReference();
        int length1 = (int) (length0 * factor);
        if (length1 < 1)
            length1 = 1;
        if (length1 > this.max_window_size)
            length1 = this.max_window_size;
        int mid = interval.getStart() + length0 / 2;
        int start = mid - length1 / 2;
        if (start < 1)
            start = 1;
        int end = mid + length1 / 2;
        final SAMSequenceRecord ssr = this.dictionary.getSequence(interval.getContig());
        if (ssr != null && end >= ssr.getSequenceLength()) {
            end = ssr.getSequenceLength();
        }
        interval = new SimpleInterval(interval.getContig(), start, end);
        break;
    }
    final String title = interval.toNiceString() + " (" + StringUtils.niceInt(interval.getLengthOnReference()) + " bp.)";
    try {
        final XMLStreamWriter w = XMLOutputFactory.newFactory().createXMLStreamWriter(pw);
        w.writeStartElement("html");
        w.writeStartElement("head");
        w.writeEmptyElement("meta");
        w.writeAttribute("charset", charset);
        w.writeStartElement("title");
        w.writeCharacters(title);
        // title
        w.writeEndElement();
        w.writeStartElement("style");
        w.writeCharacters("body {background-color:#f0f0f0;color:#070707;font-size:18px;}" + "h1 {text-align:center;color:#070707;}" + ".grid-container {display: grid; grid-template-columns: " + IntStream.range(0, Math.max(1, columns_count)).mapToObj(X -> "auto").collect(Collectors.joining(" ")) + ";grid-gap: 10px;  padding: 10px;}" + ".span1 {border:1px dotted blue;}" + ".lbl {font-weight: bold;}" + ".bampath {font-family:monospace;font-size:12px; font-style: normal;color:gray;}" + ".headerform {background-color:lavender;text-align:center;font-size:14px;}" + ".comment {background-color:khaki;text-align:center;font-size:14px;}" + ".highlight {background-color:#DB7093;}" + ".parents {font-size:75%;color:gray;}" + ".children {font-size:75%;color:gray;}" + ".allsamples {font-size:125%;}" + ".message {color:red;}" + ".affected {background-color:#e6cccc;}" + ".gtf {background-color:moccasin;text-align:center;}" + ".known {background-color:wheat;text-align:center;}");
        // title
        w.writeEndElement();
        w.writeStartElement("script");
        {
            w.writeCharacters("");
            w.flush();
            /* messy with < and > characters + xmlstream */
            pw.write("function highlight(name) {" + "var e = document.getElementById(name+\"_div\"); if(e==null) return;" + "e.classList.add(\"highlight\");" + "setTimeout(function () {e.classList.remove(\"highlight\");},2000);" + "}" + "function sendComment() {" + "var commentstr=document.getElementById(\"comment\").value;" + "console.log(\"comment is\"+commentstr);" + "if(commentstr.trim().length==0) return;" + "var xhttp = new XMLHttpRequest();" + "xhttp.onreadystatechange = function() {if (this.readyState == 4) if(this.status == 200) { alert(this.responseText);} else alert(\"Cannot send Comment.\")};" + "xhttp.open(\"GET\",\"/comment?interval=" + StringUtils.escapeHttp(interval.toString()) + "&comment=\"+encodeURI(commentstr), true);" + "xhttp.send();" + "}" + "function loadImage(idx) {" + "if(idx>=" + this.bamInput.size() + ") return;" + "var img = document.getElementById(\"bamid\"+idx);" + "img.addEventListener('load',(event) => {img.width=" + image_width + ";img.height=" + image_height + ";loadImage(idx+1);});" + "img.setAttribute(\"src\",\"/getimage?id=\"+idx+\"&interval=" + StringUtils.escapeHttp(interval.toString()) + (normalize ? "&normalize=1" : "") + "\");" + "img.setAttribute(\"alt\",\"bam idx\"+idx);" + "}" + "function init() {" + "var span=document.getElementById(\"spaninterval\");" + "if(span!=null) span.addEventListener('click',(evt)=>{document.getElementById(\"custominterval\").value = span.textContent; });" + "var sel=document.getElementById(\"selinterval\");" + "if(sel!=null) sel.addEventListener('change',(evt)=>{document.getElementById(\"custominterval\").value = evt.target.value; });" + "var cbox=document.getElementById(\"norminput\");" + "if(cbox!=null) cbox.addEventListener('change',(evt)=>{cbox.form.submit(); });" + "var comment=document.getElementById(\"commentbtn\");" + "if(comment!=null) comment.addEventListener('click',(evt)=>{ console.log(\"send comment\");sendComment(); });" + "var shortcuts=document.getElementById(\"shortcuts\");" + "if(shortcuts!=null) shortcuts.addEventListener('change',(evt)=>{document.getElementById(\"comment\").value += evt.target.value; });" + "loadImage(0);" + "}" + "window.addEventListener('load', (event) => {init();});");
            pw.flush();
        }
        // script
        w.writeEndElement();
        // head
        w.writeEndElement();
        w.writeStartElement("body");
        w.writeStartElement("h1");
        w.writeCharacters(SequenceDictionaryUtils.getBuildName(this.dictionary).orElse("") + " ");
        final String outlinkurl = hyperlink.apply(interval);
        if (StringUtils.isBlank(outlinkurl)) {
            w.writeCharacters(title);
        } else {
            w.writeStartElement("a");
            w.writeAttribute("href", outlinkurl);
            w.writeAttribute("target", "_blank");
            w.writeAttribute("title", title);
            w.writeCharacters(title);
            w.writeEndElement();
        }
        // h1
        w.writeEndElement();
        if (!StringUtils.isBlank(message)) {
            w.writeStartElement("h2");
            w.writeAttribute("class", "message");
            w.writeCharacters(message);
            // h2
            w.writeEndElement();
        }
        w.writeEmptyElement("a");
        w.writeAttribute("name", "top");
        w.writeComment("BEGIN FORM");
        w.writeStartElement("div");
        w.writeAttribute("class", "headerform");
        w.writeStartElement("form");
        w.writeAttribute("method", "GET");
        w.writeAttribute("action", "/page");
        w.writeEmptyElement("input");
        w.writeAttribute("name", "interval");
        w.writeAttribute("value", interval.toString());
        w.writeAttribute("type", "hidden");
        /* number of images per row */
        if (!StringUtils.isBlank(request.getParameter("columns"))) {
            w.writeEmptyElement("input");
            w.writeAttribute("name", "columns");
            w.writeAttribute("value", String.valueOf(columns_count));
            w.writeAttribute("type", "hidden");
        }
        /* write select box with predefined interval */
        if (!this.named_intervals.isEmpty()) {
            w.writeStartElement("select");
            w.writeAttribute("class", "btn");
            w.writeAttribute("id", "selinterval");
            w.writeEmptyElement("option");
            for (int side = 0; side < 2; ++side) {
                for (final ReviewedInterval r : this.named_intervals) {
                    if (side == 0 && r.reviewed)
                        continue;
                    if (side == 1 && !r.reviewed)
                        continue;
                    final SimpleInterval simple = new SimpleInterval(r);
                    w.writeStartElement("option");
                    w.writeAttribute("value", simple.toString());
                    if (simple.equals(interval)) {
                        r.reviewed = true;
                        w.writeAttribute("selected", "true");
                    }
                    if (r.reviewed) {
                        w.writeEntityRef("#x2713");
                        w.writeCharacters(" ");
                    }
                    w.writeCharacters(simple.toNiceString() + (StringUtils.isBlank(r.getName()) ? "" : " [" + r.getName() + "]"));
                    w.writeEndElement();
                    w.writeCharacters("\n");
                }
            }
            // select
            w.writeEndElement();
        }
        w.writeEmptyElement("br");
        // 
        /* move */
        w.writeStartElement("label");
        w.writeAttribute("class", "lbl");
        w.writeCharacters("move");
        w.writeEndElement();
        for (final String mv : new String[] { "<<<", "<<", "<", ">", ">>", ">>>" }) {
            w.writeEmptyElement("input");
            w.writeAttribute("class", "btn");
            w.writeAttribute("type", "submit");
            w.writeAttribute("name", "move");
            w.writeAttribute("value", mv);
        }
        /* zoom in */
        w.writeStartElement("label");
        w.writeAttribute("class", "lbl");
        w.writeCharacters("zoom in");
        w.writeEndElement();
        for (final String zoom : new String[] { "1.5", "3", "10" }) {
            w.writeEmptyElement("input");
            w.writeAttribute("class", "btn");
            w.writeAttribute("type", "submit");
            w.writeAttribute("name", "zoomin");
            w.writeAttribute("value", zoom);
        }
        /* zoom in */
        w.writeStartElement("label");
        w.writeAttribute("class", "lbl");
        w.writeCharacters("zoom out");
        w.writeEndElement();
        for (final String zoom : new String[] { "1.5", "3", "10", "100" }) {
            w.writeEmptyElement("input");
            w.writeAttribute("type", "submit");
            w.writeAttribute("class", "btn");
            w.writeAttribute("name", "zoomout");
            w.writeAttribute("value", zoom);
        }
        /* checkbox normalize */
        w.writeEmptyElement("input");
        w.writeAttribute("id", "norminput");
        w.writeAttribute("name", "normalize");
        w.writeAttribute("type", "checkbox");
        if (normalize)
            w.writeAttribute("checked", "true");
        w.writeStartElement("label");
        w.writeAttribute("class", "lbl");
        w.writeAttribute("for", "norminput");
        w.writeCharacters("normalize");
        w.writeEndElement();
        w.writeEmptyElement("br");
        w.writeStartElement("span");
        w.writeAttribute("class", "span1");
        w.writeAttribute("id", "spaninterval");
        w.writeCharacters(interval.toNiceString());
        // span
        w.writeEndElement();
        w.writeEmptyElement("input");
        w.writeAttribute("name", "custominterval");
        w.writeAttribute("id", "custominterval");
        w.writeAttribute("type", "text");
        w.writeAttribute("value", "");
        w.writeStartElement("button");
        w.writeAttribute("class", "btn");
        w.writeAttribute("name", "parseinterval");
        w.writeCharacters("GO");
        // button
        w.writeEndElement();
        // form
        w.writeEndElement();
        // div
        w.writeEndElement();
        if (this.commentPath != null) {
            w.writeStartElement("div");
            w.writeAttribute("class", "comment");
            w.writeStartElement("label");
            w.writeAttribute("for", "comment");
            w.writeCharacters("Comment:");
            w.writeEndElement();
            w.writeEmptyElement("input");
            w.writeAttribute("id", "comment");
            w.writeAttribute("type", "text");
            w.writeAttribute("placeholder", "Comment about " + interval.toNiceString());
            w.writeStartElement("button");
            w.writeAttribute("class", "btn");
            w.writeAttribute("id", "commentbtn");
            w.writeCharacters("Send comment");
            // button
            w.writeEndElement();
            w.writeCharacters(" ");
            w.writeStartElement("select");
            w.writeAttribute("class", "btn");
            w.writeAttribute("id", "shortcuts");
            w.writeEmptyElement("option");
            for (final String opt : new String[] { "OK", "BAD", "Noisy", "False positive", "Deletion", "Later" }) {
                w.writeStartElement("option");
                w.writeAttribute("value", " " + opt + ".");
                w.writeCharacters(opt);
                w.writeEndElement();
            }
            // select
            w.writeEndElement();
            // div
            w.writeEndElement();
        }
        w.writeComment("END FORM");
        if (this.gtfFile != null) {
            w.writeStartElement("div");
            w.writeAttribute("class", "gtf");
            getGenes(interval).filter(G -> G.getType().equals("gene")).forEach(gtfLine -> {
                try {
                    String key = null;
                    if (gtfLine.getAttributes().containsKey("gene_id")) {
                        key = gtfLine.getAttribute("gene_id");
                    }
                    if (gtfLine.getAttributes().containsKey("gene_name")) {
                        key = gtfLine.getAttribute("gene_name");
                    }
                    if (StringUtils.isBlank(key))
                        return;
                    final String url = "https://www.ncbi.nlm.nih.gov/gene/?term=" + StringUtils.escapeHttp(key);
                    w.writeStartElement("a");
                    w.writeAttribute("href", url);
                    w.writeAttribute("target", "_blank");
                    w.writeAttribute("title", key);
                    w.writeCharacters(key);
                    w.writeEndElement();
                    w.writeCharacters(" ");
                } catch (final XMLStreamException err) {
                }
            });
            // div
            w.writeEndElement();
        }
        if (this.knownCnvFile != null) {
            final Interval final_interval = new Interval(interval);
            final ToDoubleFunction<Interval> getOvelap = R -> {
                if (!R.overlaps(final_interval))
                    return 0;
                final double L1 = R.getIntersectionLength(final_interval);
                final double r1 = L1 / R.getLengthOnReference();
                final double r2 = L1 / final_interval.getLengthOnReference();
                return Double.min(r1, r2);
            };
            w.writeStartElement("div");
            w.writeAttribute("class", "known");
            w.writeStartElement("label");
            w.writeCharacters("Known: ");
            w.writeEndElement();
            getKnownCnv(interval).sorted((A, B) -> Double.compare(getOvelap.applyAsDouble(B), (getOvelap.applyAsDouble(A)))).limit(100).forEach(R -> {
                try {
                    w.writeStartElement("span");
                    w.writeCharacters(new SimpleInterval(R).toNiceString());
                    if (!StringUtils.isBlank(R.getName())) {
                        w.writeCharacters(" ");
                        w.writeCharacters(R.getName());
                    }
                    w.writeEndElement();
                    w.writeCharacters("; ");
                    // span
                    w.writeEndElement();
                } catch (final Throwable err) {
                    LOG.warn(err);
                }
            });
            // div
            w.writeEndElement();
        }
        /* write anchors to samples */
        w.writeStartElement("div");
        w.writeAttribute("class", "allsamples");
        w.writeCharacters("Samples: ");
        for (final BamInput bi : this.bamInput.stream().sorted((A, B) -> A.sample.compareTo(B.sample)).collect(Collectors.toList())) {
            w.writeStartElement("a");
            w.writeAttribute("title", bi.sample);
            w.writeAttribute("href", "#" + bi.sample);
            w.writeAttribute("onclick", "highlight('" + bi.sample + "');");
            w.writeCharacters("[" + bi.sample);
            if (pedigree != null) {
                final Sample sn = this.pedigree.getSampleById(bi.sample);
                if (sn != null && sn.isAffected()) {
                    w.writeCharacters(" ");
                    w.writeEntityRef("#128577");
                }
            }
            w.writeCharacters("] ");
            w.writeEndElement();
        }
        // div
        w.writeEndElement();
        w.writeCharacters("\n");
        w.writeStartElement("div");
        w.writeAttribute("class", "grid-container");
        for (int i = 0; i < this.bamInput.size(); i++) {
            final BamInput bamInput = this.bamInput.get(i);
            final Path bam = bamInput.bamPath;
            final Sample sample = this.pedigree == null ? null : this.pedigree.getSampleById(bamInput.sample);
            w.writeStartElement("div");
            w.writeAttribute("id", bamInput.sample + "_div");
            w.writeAttribute("class", "sample" + (sample != null && sample.isAffected() ? " affected" : ""));
            w.writeEmptyElement("a");
            w.writeAttribute("name", bamInput.sample);
            w.writeStartElement("h3");
            w.writeCharacters(bamInput.sample);
            w.writeCharacters(" ");
            w.writeStartElement("a");
            w.writeAttribute("href", "#top");
            w.writeAttribute("title", "top");
            w.writeCharacters("[^]");
            w.writeEndElement();
            {
                if (sample != null) {
                    writeSample(w, sample);
                    w.writeStartElement("span");
                    w.writeAttribute("class", "parents");
                    for (int p = 0; p < 2; ++p) {
                        if (p == 0 && !sample.hasFather())
                            continue;
                        if (p == 1 && !sample.hasMother())
                            continue;
                        final Sample parent = (p == 0 ? sample.getFather() : sample.getMother());
                        boolean has_bam = this.bamInput.stream().anyMatch(B -> B.sample.equals(parent.getId()));
                        w.writeCharacters(" ");
                        w.writeCharacters(p == 0 ? "Father " : "Mother ");
                        if (has_bam) {
                            w.writeStartElement("a");
                            w.writeAttribute("href", "#" + parent.getId());
                            w.writeAttribute("title", parent.getId());
                            w.writeAttribute("onclick", "highlight('" + parent.getId() + "');");
                            w.writeCharacters("[" + parent.getId() + "].");
                            w.writeEndElement();
                        } else {
                            w.writeCharacters(parent.getId());
                        }
                        writeSample(w, parent);
                    }
                    w.writeEndElement();
                    final Set<Sample> children = sample.getChildren();
                    if (!children.isEmpty()) {
                        w.writeStartElement("span");
                        w.writeAttribute("class", "children");
                        w.writeCharacters(" Has Children :");
                        for (final Sample child : children) {
                            boolean has_bam = this.bamInput.stream().anyMatch(B -> B.sample.equals(child.getId()));
                            w.writeCharacters(" ");
                            if (has_bam) {
                                w.writeStartElement("a");
                                w.writeAttribute("href", "#" + child.getId());
                                w.writeAttribute("title", child.getId());
                                w.writeAttribute("onclick", "highlight('" + child.getId() + "');");
                                w.writeCharacters("[" + child.getId() + "].");
                                w.writeEndElement();
                            } else {
                                w.writeCharacters(child.getId());
                            }
                            writeSample(w, child);
                        }
                        w.writeEndElement();
                    }
                // end has children
                }
            }
            w.writeEndElement();
            w.writeStartElement("div");
            w.writeAttribute("class", "bampath");
            w.writeCharacters(bam.toString());
            w.writeEndElement();
            w.writeEmptyElement("img");
            w.writeAttribute("id", "bamid" + i);
            w.writeAttribute("src", "https://upload.wikimedia.org/wikipedia/commons/d/de/Ajax-loader.gif");
            w.writeAttribute("width", "32");
            w.writeAttribute("height", "32");
            // div
            w.writeEndElement();
            w.writeCharacters("\n");
        }
        w.writeEndElement();
        w.writeEmptyElement("hr");
        w.writeStartElement("div");
        w.writeCharacters("Author: Pierre Lindenbaum. ");
        w.writeCharacters(JVarkitVersion.getInstance().getLabel());
        w.writeEndElement();
        // body
        w.writeEndElement();
        // html
        w.writeEndElement();
        w.flush();
        w.close();
    } catch (XMLStreamException err) {
        LOG.warn(err);
        throw new IOException(err);
    } finally {
        pw.close();
    }
}
Also used : Color(java.awt.Color) ServletContextHandler(org.eclipse.jetty.servlet.ServletContextHandler) Arrays(java.util.Arrays) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) Program(com.github.lindenb.jvarkit.util.jcommander.Program) ServletException(javax.servlet.ServletException) Rectangle2D(java.awt.geom.Rectangle2D) CigarElement(htsjdk.samtools.CigarElement) CigarOperator(htsjdk.samtools.CigarOperator) RenderingHints(java.awt.RenderingHints) IntervalListProvider(com.github.lindenb.jvarkit.samtools.util.IntervalListProvider) IntervalParserFactory(com.github.lindenb.jvarkit.samtools.util.IntervalParserFactory) SAMFileHeader(htsjdk.samtools.SAMFileHeader) ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile) Hyperlink(com.github.lindenb.jvarkit.net.Hyperlink) DistanceParser(com.github.lindenb.jvarkit.util.bio.DistanceParser) Vector(java.util.Vector) XMLStreamException(javax.xml.stream.XMLStreamException) ImageIO(javax.imageio.ImageIO) XMLStreamWriter(javax.xml.stream.XMLStreamWriter) HttpStatus(org.eclipse.jetty.http.HttpStatus) Path(java.nio.file.Path) Server(org.eclipse.jetty.server.Server) CloserUtil(htsjdk.samtools.util.CloserUtil) Shape(java.awt.Shape) PrintWriter(java.io.PrintWriter) SimpleInterval(com.github.lindenb.jvarkit.samtools.util.SimpleInterval) SequenceDictionaryUtils(com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils) HttpServlet(javax.servlet.http.HttpServlet) Composite(java.awt.Composite) BufferedImage(java.awt.image.BufferedImage) Logger(com.github.lindenb.jvarkit.util.log.Logger) StandardOpenOption(java.nio.file.StandardOpenOption) ServletHolder(org.eclipse.jetty.servlet.ServletHolder) Set(java.util.Set) Collectors(java.util.stream.Collectors) GTFCodec(com.github.lindenb.jvarkit.util.bio.gtf.GTFCodec) SAMRecord(htsjdk.samtools.SAMRecord) DoubleStream(java.util.stream.DoubleStream) ReferenceSequenceFileFactory(htsjdk.samtools.reference.ReferenceSequenceFileFactory) List(java.util.List) Stream(java.util.stream.Stream) StringUtils(com.github.lindenb.jvarkit.lang.StringUtils) FileExtensions(htsjdk.samtools.util.FileExtensions) ToDoubleFunction(java.util.function.ToDoubleFunction) BasicStroke(java.awt.BasicStroke) ReferenceSequence(htsjdk.samtools.reference.ReferenceSequence) GeneralPath(java.awt.geom.GeneralPath) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) SamReaderFactory(htsjdk.samtools.SamReaderFactory) IntStream(java.util.stream.IntStream) Cigar(htsjdk.samtools.Cigar) CloseableIterator(htsjdk.samtools.util.CloseableIterator) SequenceUtil(htsjdk.samtools.util.SequenceUtil) ContigNameConverter(com.github.lindenb.jvarkit.util.bio.fasta.ContigNameConverter) Parameter(com.beust.jcommander.Parameter) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) IterableAdapter(htsjdk.samtools.util.IterableAdapter) ValidationStringency(htsjdk.samtools.ValidationStringency) ArrayList(java.util.ArrayList) GTFLine(com.github.lindenb.jvarkit.util.bio.gtf.GTFLine) Interval(htsjdk.samtools.util.Interval) AlphaComposite(java.awt.AlphaComposite) HttpServletRequest(javax.servlet.http.HttpServletRequest) NoSplitter(com.github.lindenb.jvarkit.util.jcommander.NoSplitter) Graphics2D(java.awt.Graphics2D) PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) TabixReader(htsjdk.tribble.readers.TabixReader) VCFReaderFactory(com.github.lindenb.jvarkit.variant.vcf.VCFReaderFactory) StreamSupport(java.util.stream.StreamSupport) IntFunction(java.util.function.IntFunction) VCFConstants(htsjdk.variant.vcf.VCFConstants) Stroke(java.awt.Stroke) Line2D(java.awt.geom.Line2D) Counter(com.github.lindenb.jvarkit.util.Counter) AbstractIterator(htsjdk.samtools.util.AbstractIterator) Locatable(htsjdk.samtools.util.Locatable) IntToDoubleFunction(java.util.function.IntToDoubleFunction) Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree) Files(java.nio.file.Files) BufferedWriter(java.io.BufferedWriter) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) HttpServletResponse(javax.servlet.http.HttpServletResponse) VCFReader(htsjdk.variant.vcf.VCFReader) IOException(java.io.IOException) JVarkitVersion(com.github.lindenb.jvarkit.util.JVarkitVersion) SamReader(htsjdk.samtools.SamReader) XMLOutputFactory(javax.xml.stream.XMLOutputFactory) Pileup(com.github.lindenb.jvarkit.samtools.util.Pileup) BufferedReader(java.io.BufferedReader) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) Sample(com.github.lindenb.jvarkit.pedigree.Sample) Path(java.nio.file.Path) GeneralPath(java.awt.geom.GeneralPath) Set(java.util.Set) Sample(com.github.lindenb.jvarkit.pedigree.Sample) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) IOException(java.io.IOException) XMLStreamException(javax.xml.stream.XMLStreamException) XMLStreamWriter(javax.xml.stream.XMLStreamWriter) SimpleInterval(com.github.lindenb.jvarkit.samtools.util.SimpleInterval) PrintWriter(java.io.PrintWriter) SimpleInterval(com.github.lindenb.jvarkit.samtools.util.SimpleInterval) Interval(htsjdk.samtools.util.Interval)

Example 4 with Sample

use of com.github.lindenb.jvarkit.pedigree.Sample in project jvarkit by lindenb.

the class CoverageServer method printImage.

private void printImage(final HttpServletRequest request, final HttpServletResponse response) throws IOException, ServletException {
    int bam_id;
    try {
        bam_id = Integer.parseInt(StringUtils.ifBlank(request.getParameter("id"), "-1"));
    } catch (Exception err) {
        bam_id = -1;
    }
    final SimpleInterval midRegion = parseInterval(request.getParameter("interval"));
    if (midRegion == null || bam_id < 0 || bam_id >= this.bamInput.size()) {
        response.reset();
        response.sendError(HttpStatus.BAD_REQUEST_400, "id:" + bam_id);
        response.flushBuffer();
        return;
    }
    final int extend = (int) (midRegion.getLengthOnReference() * this.extend_factor);
    int xstart = Math.max(midRegion.getStart() - extend, 0);
    int xend = midRegion.getEnd() + extend;
    final SAMSequenceRecord ssr = this.dictionary.getSequence(midRegion.getContig());
    if (ssr != null) {
        xend = Math.min(xend, ssr.getSequenceLength());
    }
    final SimpleInterval region = new SimpleInterval(midRegion.getContig(), xstart, xend);
    if (region.getLengthOnReference() > this.max_window_size) {
        response.reset();
        response.sendError(HttpStatus.BAD_REQUEST_400, "contig:" + midRegion);
        response.flushBuffer();
        return;
    }
    final Counter<Arc> sashimiArcs = new Counter<>();
    final BamInput bam = this.bamInput.get(bam_id);
    if (region.length() <= this.small_region_size) {
        printRaster(bam, midRegion, region, request, response);
        return;
    }
    final boolean normalize = request.getParameter("normalize") != null;
    final SamReaderFactory srf = SamReaderFactory.make().validationStringency(ValidationStringency.LENIENT).referenceSequence(this.faidxRef);
    try (SamReader sr = srf.open(bam.bamPath)) {
        final int[] int_coverage = new int[region.getLengthOnReference()];
        Arrays.fill(int_coverage, 0);
        try (CloseableIterator<SAMRecord> iter = sr.query(region.getContig(), region.getStart(), region.getEnd(), false)) {
            while (iter.hasNext()) {
                final SAMRecord rec = iter.next();
                if (!acceptRead(rec))
                    continue;
                final Cigar cigar = rec.getCigar();
                if (cigar == null || cigar.isEmpty())
                    continue;
                int ref = rec.getAlignmentStart();
                for (final CigarElement ce : cigar) {
                    final CigarOperator op = ce.getOperator();
                    if (op.consumesReferenceBases()) {
                        if (this.enable_sashimi && op.equals(CigarOperator.N)) {
                            sashimiArcs.incr(new Arc(ref, ref + ce.getLength()));
                        }
                        if (op.consumesReadBases()) {
                            for (int x = 0; x < ce.getLength(); ++x) {
                                int pos = ref + x;
                                if (pos < region.getStart())
                                    continue;
                                if (pos > region.getEnd())
                                    break;
                                int_coverage[pos - region.getStart()]++;
                            }
                        }
                        ref += ce.getLength();
                    }
                }
            }
        }
        /* smooth coverage */
        if (int_coverage.length > image_width) {
            final int[] copy = Arrays.copyOf(int_coverage, int_coverage.length);
            final int len = Math.max(1, int_coverage.length / 100);
            for (int i = 0; i < int_coverage.length; i++) {
                int j = Math.max(0, i - len);
                double sum = 0;
                int count = 0;
                while (j < i + len && j < copy.length) {
                    sum += copy[j];
                    j++;
                    count++;
                }
                int_coverage[i] = (int) (sum / count);
            }
        }
        final double[] norm_coverage = new double[int_coverage.length];
        final double median;
        /* normalize on median */
        if (normalize) {
            final Coverage leftrightcov = new Coverage(extend * 2);
            for (int x = region.getStart(); x < midRegion.getStart(); x++) {
                final int idx = x - region.getStart();
                leftrightcov.add(int_coverage[idx]);
            }
            for (int x = midRegion.getEnd() + 1; x <= region.getEnd(); x++) {
                final int idx = x - region.getStart();
                leftrightcov.add(int_coverage[idx]);
            }
            median = Math.max(1.0, leftrightcov.median());
            // LOG.info("median is "+median+" "+leftrightcov.median());
            for (int x = 0; x < int_coverage.length; ++x) {
                norm_coverage[x] = int_coverage[x] / median;
            }
        } else /* no normalisation */
        {
            /* won't be used */
            median = Double.NaN;
            for (int x = 0; x < int_coverage.length; ++x) {
                norm_coverage[x] = int_coverage[x];
            }
        }
        final double real_max_cov = DoubleStream.of(norm_coverage).max().orElse(1.0);
        final double max_cov = Math.max((normalize ? 2 : 10), real_max_cov);
        final double pixelperbase = image_width / (double) norm_coverage.length;
        final IntFunction<Double> pos2pixel = POS -> ((POS - region.getStart()) / (double) region.getLengthOnReference()) * image_width;
        final BufferedImage img = new BufferedImage(image_width, image_height, BufferedImage.TYPE_INT_RGB);
        final Graphics2D g = img.createGraphics();
        g.setRenderingHint(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, image_width + 1, image_height + 1);
        for (int x = 0; x < norm_coverage.length; ++x) {
            final double height = image_height * (norm_coverage[x] / max_cov);
            if (normalize)
                g.setColor(Color.DARK_GRAY);
            else if (max_cov < 10)
                g.setColor(Color.RED);
            else if (max_cov < 20)
                g.setColor(Color.BLUE);
            else
                g.setColor(Color.DARK_GRAY);
            g.fill(new Rectangle2D.Double(x * pixelperbase, image_height - height, pixelperbase, height));
        }
        g.setColor(Color.DARK_GRAY);
        g.drawString("max-cov:" + IntStream.of(int_coverage).max().orElse(0) + (normalize ? " normalized on median (" + median + ")" : "") + " sample:" + bam.sample + " " + region.toNiceString(), 10, 10);
        /* ticks for vertical axis */
        g.setColor(Color.MAGENTA);
        for (int i = 1; i < 10; i++) {
            double cov = max_cov / 10.0 * i;
            if (!normalize)
                cov = Math.ceil(cov);
            final double y = image_height - image_height / 10.0 * i;
            if (!normalize && i > 0 && (int) cov == Math.ceil(max_cov / 10.0 * (i - 1)))
                continue;
            g.drawLine(0, (int) y, 5, (int) y);
            g.drawString(normalize ? String.format("%.2f", cov) : String.valueOf((int) cov), 7, (int) y);
        }
        /* vertical line for original view */
        g.setColor(Color.PINK);
        double vertical = ((midRegion.getStart() - region.getStart()) / (double) region.getLengthOnReference()) * image_width;
        g.draw(new Line2D.Double(vertical, 0, vertical, image_height));
        vertical = ((midRegion.getEnd() - region.getStart()) / (double) region.getLengthOnReference()) * image_width;
        g.draw(new Line2D.Double(vertical, 0, vertical, image_height));
        if (normalize) {
            /* horizontal line for median 0.5 / 1 / 1.5 */
            for (int t = 1; t < 4; ++t) {
                g.setColor(t == 2 ? Color.ORANGE : Color.PINK);
                final double mediany = image_height - ((0.5 * t) / max_cov) * image_height;
                g.draw(new Line2D.Double(0, mediany, image_width, mediany));
            }
        }
        if (this.enable_sashimi && !sashimiArcs.isEmpty()) {
            final double max_count = sashimiArcs.getMaxCount().orElse(1L);
            g.setColor(Color.GREEN);
            for (final Arc arc : sashimiArcs.keySet()) {
                final double x1 = pos2pixel.apply(arc.start);
                final double x2 = pos2pixel.apply(arc.end);
                final double distance = x2 - x1;
                final GeneralPath curve = new GeneralPath();
                curve.moveTo(x1, image_height);
                curve.curveTo(x1, image_height, x1 + distance / 2.0, image_height - Math.min(distance, image_height * 0.75), x2, image_height);
                final double weight = (sashimiArcs.count(arc) / max_count) * 5;
                final Stroke oldStroke = g.getStroke();
                final Composite oldComposite = g.getComposite();
                g.setComposite(AlphaComposite.getInstance(AlphaComposite.SRC_OVER, 0.5f));
                g.setStroke(new BasicStroke((float) weight, BasicStroke.CAP_ROUND, BasicStroke.JOIN_ROUND));
                g.draw(curve);
                g.setStroke(oldStroke);
                g.setComposite(oldComposite);
            }
        }
        writeGenes(g, region);
        writeKnownCnv(g, region);
        g.setColor(Color.GRAY);
        g.drawRect(0, 0, img.getWidth(), img.getHeight());
        writeImage(img, bam, region, response);
    }
}
Also used : Color(java.awt.Color) ServletContextHandler(org.eclipse.jetty.servlet.ServletContextHandler) Arrays(java.util.Arrays) CharSplitter(com.github.lindenb.jvarkit.lang.CharSplitter) Program(com.github.lindenb.jvarkit.util.jcommander.Program) ServletException(javax.servlet.ServletException) Rectangle2D(java.awt.geom.Rectangle2D) CigarElement(htsjdk.samtools.CigarElement) CigarOperator(htsjdk.samtools.CigarOperator) RenderingHints(java.awt.RenderingHints) IntervalListProvider(com.github.lindenb.jvarkit.samtools.util.IntervalListProvider) IntervalParserFactory(com.github.lindenb.jvarkit.samtools.util.IntervalParserFactory) SAMFileHeader(htsjdk.samtools.SAMFileHeader) ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile) Hyperlink(com.github.lindenb.jvarkit.net.Hyperlink) DistanceParser(com.github.lindenb.jvarkit.util.bio.DistanceParser) Vector(java.util.Vector) XMLStreamException(javax.xml.stream.XMLStreamException) ImageIO(javax.imageio.ImageIO) XMLStreamWriter(javax.xml.stream.XMLStreamWriter) HttpStatus(org.eclipse.jetty.http.HttpStatus) Path(java.nio.file.Path) Server(org.eclipse.jetty.server.Server) CloserUtil(htsjdk.samtools.util.CloserUtil) Shape(java.awt.Shape) PrintWriter(java.io.PrintWriter) SimpleInterval(com.github.lindenb.jvarkit.samtools.util.SimpleInterval) SequenceDictionaryUtils(com.github.lindenb.jvarkit.util.bio.SequenceDictionaryUtils) HttpServlet(javax.servlet.http.HttpServlet) Composite(java.awt.Composite) BufferedImage(java.awt.image.BufferedImage) Logger(com.github.lindenb.jvarkit.util.log.Logger) StandardOpenOption(java.nio.file.StandardOpenOption) ServletHolder(org.eclipse.jetty.servlet.ServletHolder) Set(java.util.Set) Collectors(java.util.stream.Collectors) GTFCodec(com.github.lindenb.jvarkit.util.bio.gtf.GTFCodec) SAMRecord(htsjdk.samtools.SAMRecord) DoubleStream(java.util.stream.DoubleStream) ReferenceSequenceFileFactory(htsjdk.samtools.reference.ReferenceSequenceFileFactory) List(java.util.List) Stream(java.util.stream.Stream) StringUtils(com.github.lindenb.jvarkit.lang.StringUtils) FileExtensions(htsjdk.samtools.util.FileExtensions) ToDoubleFunction(java.util.function.ToDoubleFunction) BasicStroke(java.awt.BasicStroke) ReferenceSequence(htsjdk.samtools.reference.ReferenceSequence) GeneralPath(java.awt.geom.GeneralPath) BedLine(com.github.lindenb.jvarkit.util.bio.bed.BedLine) SamReaderFactory(htsjdk.samtools.SamReaderFactory) IntStream(java.util.stream.IntStream) Cigar(htsjdk.samtools.Cigar) CloseableIterator(htsjdk.samtools.util.CloseableIterator) SequenceUtil(htsjdk.samtools.util.SequenceUtil) ContigNameConverter(com.github.lindenb.jvarkit.util.bio.fasta.ContigNameConverter) Parameter(com.beust.jcommander.Parameter) BedLineCodec(com.github.lindenb.jvarkit.util.bio.bed.BedLineCodec) IterableAdapter(htsjdk.samtools.util.IterableAdapter) ValidationStringency(htsjdk.samtools.ValidationStringency) ArrayList(java.util.ArrayList) GTFLine(com.github.lindenb.jvarkit.util.bio.gtf.GTFLine) Interval(htsjdk.samtools.util.Interval) AlphaComposite(java.awt.AlphaComposite) HttpServletRequest(javax.servlet.http.HttpServletRequest) NoSplitter(com.github.lindenb.jvarkit.util.jcommander.NoSplitter) Graphics2D(java.awt.Graphics2D) PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser) IOUtils(com.github.lindenb.jvarkit.io.IOUtils) Launcher(com.github.lindenb.jvarkit.util.jcommander.Launcher) TabixReader(htsjdk.tribble.readers.TabixReader) VCFReaderFactory(com.github.lindenb.jvarkit.variant.vcf.VCFReaderFactory) StreamSupport(java.util.stream.StreamSupport) IntFunction(java.util.function.IntFunction) VCFConstants(htsjdk.variant.vcf.VCFConstants) Stroke(java.awt.Stroke) Line2D(java.awt.geom.Line2D) Counter(com.github.lindenb.jvarkit.util.Counter) AbstractIterator(htsjdk.samtools.util.AbstractIterator) Locatable(htsjdk.samtools.util.Locatable) IntToDoubleFunction(java.util.function.IntToDoubleFunction) Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree) Files(java.nio.file.Files) BufferedWriter(java.io.BufferedWriter) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) HttpServletResponse(javax.servlet.http.HttpServletResponse) VCFReader(htsjdk.variant.vcf.VCFReader) IOException(java.io.IOException) JVarkitVersion(com.github.lindenb.jvarkit.util.JVarkitVersion) SamReader(htsjdk.samtools.SamReader) XMLOutputFactory(javax.xml.stream.XMLOutputFactory) Pileup(com.github.lindenb.jvarkit.samtools.util.Pileup) BufferedReader(java.io.BufferedReader) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) Sample(com.github.lindenb.jvarkit.pedigree.Sample) BasicStroke(java.awt.BasicStroke) GeneralPath(java.awt.geom.GeneralPath) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) Line2D(java.awt.geom.Line2D) BufferedImage(java.awt.image.BufferedImage) SamReader(htsjdk.samtools.SamReader) Counter(com.github.lindenb.jvarkit.util.Counter) SimpleInterval(com.github.lindenb.jvarkit.samtools.util.SimpleInterval) BasicStroke(java.awt.BasicStroke) Stroke(java.awt.Stroke) SamReaderFactory(htsjdk.samtools.SamReaderFactory) Composite(java.awt.Composite) AlphaComposite(java.awt.AlphaComposite) Rectangle2D(java.awt.geom.Rectangle2D) CigarOperator(htsjdk.samtools.CigarOperator) CigarElement(htsjdk.samtools.CigarElement) ServletException(javax.servlet.ServletException) XMLStreamException(javax.xml.stream.XMLStreamException) IOException(java.io.IOException) Graphics2D(java.awt.Graphics2D) Cigar(htsjdk.samtools.Cigar) SAMRecord(htsjdk.samtools.SAMRecord)

Example 5 with Sample

use of com.github.lindenb.jvarkit.pedigree.Sample in project jvarkit by lindenb.

the class VcfFilterNotInPedigree method doVcfToVcf.

@Override
protected int doVcfToVcf(String inputName, VCFIterator iterin, VariantContextWriter out) {
    final int IGNORE_SINGLETON = -1;
    final Predicate<Genotype> acceptFilteredGenotype = G -> G != null && (this.ignoreFilteredGT == false || !G.isFiltered());
    final Pedigree pedigree;
    final VCFHeader header = iterin.getHeader();
    try {
        pedigree = new PedigreeParser().parse(this.pedigreeFile);
    } catch (final IOException err) {
        LOG.error("Cannot read pedigree in file: " + this.pedigreeFile, err);
        return -1;
    }
    if (pedigree.isEmpty()) {
        throw new JvarkitException.PedigreeError("No pedigree found in header. use VcfInjectPedigree to add it");
    }
    pedigree.checkUniqIds();
    final Set<String> samplesNames = new HashSet<>(header.getSampleNamesInOrder());
    final Set<Sample> individuals = new HashSet<>(pedigree.getSamples());
    final Iterator<Sample> iter = individuals.iterator();
    while (iter.hasNext()) {
        final Sample person = iter.next();
        if (!(samplesNames.contains(person.getId()))) {
            LOG.warn("Ignoring " + person + " because not in VCF header.");
            iter.remove();
        }
    }
    final VCFFilterHeaderLine filter = new VCFFilterHeaderLine(this.filterName, "Will be set for variant where the only genotypes non-homref are NOT in the pedigree ");
    final VCFFilterHeaderLine singletonFilter = new VCFFilterHeaderLine(this.singletonfilterName, "The ALT allele is found in less or equals than " + this.singleton + " individuals in the cases/controls");
    final VCFHeader h2 = new VCFHeader(header);
    JVarkitVersion.getInstance().addMetaData(this, h2);
    if (this.singleton != IGNORE_SINGLETON) {
        h2.addMetaDataLine(singletonFilter);
    }
    out.writeHeader(h2);
    while (iterin.hasNext()) {
        final VariantContext ctx = iterin.next();
        final boolean in_pedigree = individuals.stream().map(P -> ctx.getGenotype(P.getId())).filter(acceptFilteredGenotype).anyMatch(g -> (!(g == null || !g.isCalled() || !g.isAvailable() || g.isNoCall() || g.isHomRef())));
        if (!in_pedigree) {
            if (this.dicardVariant)
                continue;
            final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
            vcb.filter(filter.getID());
            out.add(vcb.make());
        } else {
            boolean is_singleton;
            if (this.singleton != IGNORE_SINGLETON) {
                is_singleton = true;
                for (final Allele alt : ctx.getAlternateAlleles()) {
                    if (individuals.stream().map(P -> ctx.getGenotype(P.getId())).filter(acceptFilteredGenotype).filter(g -> g != null && g.isCalled() && g.getAlleles().contains(alt)).count() > this.singleton) {
                        is_singleton = false;
                        break;
                    }
                }
            } else {
                is_singleton = false;
            }
            if (is_singleton) {
                if (this.dicardVariant)
                    continue;
                final VariantContextBuilder vcb = new VariantContextBuilder(ctx);
                vcb.filter(singletonFilter.getID());
                out.add(vcb.make());
            } else {
                out.add(ctx);
            }
        }
    }
    return 0;
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) Allele(htsjdk.variant.variantcontext.Allele) VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) Iterator(java.util.Iterator) Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree) Program(com.github.lindenb.jvarkit.util.jcommander.Program) VCFIterator(htsjdk.variant.vcf.VCFIterator) Parameter(com.beust.jcommander.Parameter) Predicate(java.util.function.Predicate) VCFHeader(htsjdk.variant.vcf.VCFHeader) Logger(com.github.lindenb.jvarkit.util.log.Logger) Set(java.util.Set) IOException(java.io.IOException) JVarkitVersion(com.github.lindenb.jvarkit.util.JVarkitVersion) JvarkitException(com.github.lindenb.jvarkit.lang.JvarkitException) HashSet(java.util.HashSet) PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser) VariantContextWriter(htsjdk.variant.variantcontext.writer.VariantContextWriter) VariantContext(htsjdk.variant.variantcontext.VariantContext) OnePassVcfLauncher(com.github.lindenb.jvarkit.jcommander.OnePassVcfLauncher) Path(java.nio.file.Path) Sample(com.github.lindenb.jvarkit.pedigree.Sample) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) Sample(com.github.lindenb.jvarkit.pedigree.Sample) Genotype(htsjdk.variant.variantcontext.Genotype) VariantContext(htsjdk.variant.variantcontext.VariantContext) IOException(java.io.IOException) PedigreeParser(com.github.lindenb.jvarkit.pedigree.PedigreeParser) Allele(htsjdk.variant.variantcontext.Allele) Pedigree(com.github.lindenb.jvarkit.pedigree.Pedigree) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VCFFilterHeaderLine(htsjdk.variant.vcf.VCFFilterHeaderLine) VCFHeader(htsjdk.variant.vcf.VCFHeader) HashSet(java.util.HashSet)

Aggregations

Sample (com.github.lindenb.jvarkit.pedigree.Sample)16 PedigreeParser (com.github.lindenb.jvarkit.pedigree.PedigreeParser)13 Parameter (com.beust.jcommander.Parameter)12 Pedigree (com.github.lindenb.jvarkit.pedigree.Pedigree)12 Program (com.github.lindenb.jvarkit.util.jcommander.Program)12 Logger (com.github.lindenb.jvarkit.util.log.Logger)12 Genotype (htsjdk.variant.variantcontext.Genotype)12 VariantContext (htsjdk.variant.variantcontext.VariantContext)12 Path (java.nio.file.Path)12 List (java.util.List)12 ArrayList (java.util.ArrayList)11 Set (java.util.Set)11 Collectors (java.util.stream.Collectors)11 IOException (java.io.IOException)10 StringUtils (com.github.lindenb.jvarkit.lang.StringUtils)9 Launcher (com.github.lindenb.jvarkit.util.jcommander.Launcher)9 JVarkitVersion (com.github.lindenb.jvarkit.util.JVarkitVersion)8 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)8 VCFHeader (htsjdk.variant.vcf.VCFHeader)8 PrintWriter (java.io.PrintWriter)8