Examples with CigarElement htsjdk.samtools.CigarElement used on opensource projects

Example 71 with CigarElement

use of htsjdk.samtools.CigarElement in project jvarkit by lindenb.

the class ExtendReferenceWithReads method scanRegion.

/**
 *scanRegion
 * @param contig    Reference sequence of interest.
 * @param start     0-based, inclusive start of interval of interest. Zero implies start of the reference sequence.
 * @param end       0-based, exclusive end of interval of interest. Zero implies end of the reference sequence.
 */
private Map<Integer, Counter<Byte>> scanRegion(final SAMSequenceRecord contig, final int chromStart, final int chromEnd, final Rescue rescue) {
    final Map<Integer, Counter<Byte>> pos2bases = new HashMap<>(1 + chromEnd - chromStart);
    LOG.info("Scanning :" + contig.getSequenceName() + ":" + chromStart + "-" + chromEnd);
    for (int side = 0; side < 2; ++side) {
        if (// 5' or 3'
        !rescue.equals(Rescue.CENTER) && side > 0) {
            // already done
            break;
        }
        for (final SamReader samReader : samReaders) {
            final SAMSequenceDictionary dict2 = samReader.getFileHeader().getSequenceDictionary();
            final SAMSequenceRecord ssr2 = dict2.getSequence(contig.getSequenceName());
            if (ssr2 == null || ssr2.getSequenceLength() != contig.getSequenceLength()) {
                LOG.info("No contig " + contig.getSequenceName() + " with L=" + contig.getSequenceLength() + " bases in " + samReader.getResourceDescription());
                continue;
            }
            int mappedPos = -1;
            switch(rescue) {
                case LEFT:
                    mappedPos = 1;
                    break;
                case RIGHT:
                    mappedPos = contig.getSequenceLength();
                    break;
                case CENTER:
                    mappedPos = (side == 0 ? chromStart + 1 : chromEnd);
                    break;
                default:
                    throw new IllegalStateException();
            }
            final SAMRecordIterator iter = samReader.query(contig.getSequenceName(), mappedPos, mappedPos, false);
            while (iter.hasNext()) {
                final SAMRecord rec = iter.next();
                if (rec.getReadUnmappedFlag())
                    continue;
                if (this.filter.filterOut(rec))
                    continue;
                final Cigar cigar = rec.getCigar();
                if (cigar == null)
                    continue;
                final byte[] bases = rec.getReadBases();
                if (bases == null || bases.length == 0)
                    continue;
                int refPos1 = rec.getUnclippedStart();
                int readpos = 0;
                for (final CigarElement ce : cigar) {
                    final CigarOperator op = ce.getOperator();
                    for (int L = 0; L < ce.getLength(); ++L) {
                        if (op.consumesReadBases()) {
                            Counter<Byte> count = pos2bases.get(refPos1 - 1);
                            if (count == null) {
                                count = new Counter<>();
                                pos2bases.put(refPos1 - 1, count);
                            }
                            count.incr((byte) Character.toLowerCase(bases[readpos]));
                            readpos++;
                        }
                        if (op.consumesReferenceBases()) {
                            refPos1++;
                        }
                    }
                }
            }
            iter.close();
        }
    }
    return pos2bases;
}

Example 72 with CigarElement

use of htsjdk.samtools.CigarElement in project jvarkit by lindenb.

the class BamStage method reloadData.

@Override
void reloadData() {
    updateStatusBar(AlertType.NONE, "");
    final int max_items = super.maxItemsLimitSpinner.getValue();
    final List<SAMRecord> L = new ArrayList<SAMRecord>(max_items);
    final String location = this.gotoField.getText().trim();
    final CloseableIterator<SAMRecord> iter;
    final java.util.function.Predicate<SAMRecord> recFilter = makeFlagPredicate();
    Long canvasFirstRecordGenomicIndex = null;
    Long canvasLastRecordGenomicIndex = null;
    try {
        if (location.isEmpty()) {
            iter = this.getBamFile().iterator();
        } else if (location.equalsIgnoreCase("unmapped")) {
            iter = this.getBamFile().queryUnmapped();
        } else {
            final Interval interval = parseInterval(location);
            if (interval == null) {
                iter = null;
            } else {
                iter = this.getBamFile().iterator(interval.getContig(), interval.getStart(), interval.getEnd());
            }
        }
    } catch (final IOException err) {
        err.printStackTrace();
        JfxNgs.showExceptionDialog(this, err);
        return;
    }
    Optional<BamJavascripFilter> bamjsfilter = Optional.empty();
    if (this.owner.javascriptCompiler.isPresent() && !this.javascriptArea.getText().trim().isEmpty()) {
        try {
            bamjsfilter = Optional.of(new BamJavascripFilter(this.getBamFile().getHeader(), Optional.of(this.owner.javascriptCompiler.get().compile(this.javascriptArea.getText()))));
        } catch (final Exception err) {
            LOG.warning(err.getMessage());
            updateStatusBar(AlertType.ERROR, err);
            bamjsfilter = Optional.empty();
        }
    }
    final Map<ContigPos, Pileup> pos2pileup = new TreeMap<>();
    final Function<ContigPos, Pileup> getpileup = new Function<ContigPos, Pileup>() {

        @Override
        public Pileup apply(ContigPos t) {
            Pileup p = pos2pileup.get(t);
            if (p == null) {
                p = new Pileup(t.contig, t.position);
                pos2pileup.put(t, p);
            }
            return p;
        }
    };
    int count_items = 0;
    while (iter != null && iter.hasNext() && count_items < max_items) {
        final SAMRecord rec = iter.next();
        ++count_items;
        if (bamjsfilter.isPresent()) {
            if (bamjsfilter.get().eval(rec) == null)
                continue;
        }
        if (!recFilter.test(rec))
            continue;
        L.add(rec);
        /* get bounds for canvas genmic browser */
        if (!rec.getReadUnmappedFlag() && rec.getCigar() != null) {
            final long endIndex = convertContigPosToGenomicIndex(rec.getContig(), rec.getUnclippedEnd());
            if (canvasFirstRecordGenomicIndex == null) {
                canvasFirstRecordGenomicIndex = convertContigPosToGenomicIndex(rec.getContig(), rec.getUnclippedStart());
                canvasLastRecordGenomicIndex = endIndex;
            } else if (canvasLastRecordGenomicIndex < endIndex) {
                canvasLastRecordGenomicIndex = endIndex;
            }
        }
        /* FILL pileup */
        if (!rec.getReadUnmappedFlag() && rec.getCigar() != null) {
            int refpos = rec.getUnclippedStart();
            int readpos = 0;
            final byte[] bases = rec.getReadBases();
            final byte[] quals = rec.getOriginalBaseQualities();
            /**
             * function getting the ith base
             */
            final Function<Integer, Character> getBaseAt = new Function<Integer, Character>() {

                @Override
                public Character apply(Integer readPos) {
                    char c;
                    if (bases == null || bases.length <= readPos) {
                        return '?';
                    } else {
                        c = (char) bases[readPos];
                    }
                    c = (rec.getReadNegativeStrandFlag() ? Character.toLowerCase(c) : Character.toUpperCase(c));
                    return c;
                }
            };
            /**
             * function getting the ith base quality
             */
            final Function<Integer, Character> getQualAt = new Function<Integer, Character>() {

                @Override
                public Character apply(Integer readPos) {
                    char c;
                    if (quals == null || quals.length <= readPos) {
                        return '#';
                    } else {
                        c = (char) quals[readPos];
                    }
                    return c;
                }
            };
            for (final CigarElement ce : rec.getCigar()) {
                switch(ce.getOperator()) {
                    case P:
                        break;
                    case N:
                    case D:
                        {
                            refpos += ce.getLength();
                            break;
                        }
                    case H:
                        {
                            for (int i = 0; i < ce.getLength(); ++i) {
                                final Pileup p = getpileup.apply(new ContigPos(rec.getContig(), refpos));
                                p.watch('-', '#', ce.getOperator());
                                ++refpos;
                            }
                            break;
                        }
                    case I:
                        {
                            for (int i = 0; i < ce.getLength(); ++i) {
                                final Pileup p = getpileup.apply(new ContigPos(rec.getContig(), refpos));
                                p.watch('<', getQualAt.apply(readpos), ce.getOperator());
                                readpos++;
                            }
                            break;
                        }
                    case S:
                        for (int i = 0; i < ce.getLength(); ++i) {
                            final Pileup p = getpileup.apply(new ContigPos(rec.getContig(), refpos));
                            p.watch('-', getQualAt.apply(readpos), ce.getOperator());
                            ++readpos;
                            ++refpos;
                        }
                        break;
                    case EQ:
                    case X:
                    case M:
                        for (int i = 0; i < ce.getLength(); ++i) {
                            final Pileup p = getpileup.apply(new ContigPos(rec.getContig(), refpos));
                            p.watch(getBaseAt.apply(readpos), getQualAt.apply(readpos), ce.getOperator());
                            ++readpos;
                            ++refpos;
                        }
                        break;
                }
            }
        }
    }
    if (iter != null)
        iter.close();
    this.canvasScrolVInCoverage.setMin(0);
    final int max_depth = pos2pileup.values().stream().map(P -> P.depth()).max((A, B) -> (A.compareTo(B))).orElse(0);
    this.canvasScrolVInCoverage.setMax(max_depth + 1);
    this.canvasScrolVInCoverage.setValue(0);
    this.recordTable.getItems().setAll(L);
    this.pileupTable.getItems().setAll(pos2pileup.values());
    if (canvasFirstRecordGenomicIndex != null && canvasLastRecordGenomicIndex != null && canvasFirstRecordGenomicIndex.longValue() < canvasLastRecordGenomicIndex.longValue()) {
        this.canvasScrollHGenomicLoc.setMin(canvasFirstRecordGenomicIndex);
        this.canvasScrollHGenomicLoc.setMax(canvasLastRecordGenomicIndex);
        // this.canvasScrollHGenomicLoc.setUnitIncrement(1);
        // this.canvasScrollHGenomicLoc.setBlockIncrement(Math.max(1.0,Math.min( canvasLastRecordGenomicIndex-canvasFirstRecordGenomicIndex, this.canvas.getWidth())));
        this.canvasScrollHGenomicLoc.setValue(canvasFirstRecordGenomicIndex);
    } else {
        this.canvasScrollHGenomicLoc.setMin(0);
        this.canvasScrollHGenomicLoc.setMax(0);
        this.canvasScrollHGenomicLoc.setValue(0);
        this.canvasScrollHGenomicLoc.setBlockIncrement(0);
        this.canvasScrollHGenomicLoc.setUnitIncrement(1);
    }
    if (!this.recordTable.getItems().isEmpty()) {
        final int last_index = this.recordTable.getItems().size() - 1;
        if (!this.recordTable.getItems().get(0).getReadUnmappedFlag() && !this.recordTable.getItems().get(last_index).getReadUnmappedFlag()) {
            super.seqDictionaryCanvas.setItemsInterval(new ContigPos(this.recordTable.getItems().get(0).getContig(), this.recordTable.getItems().get(0).getStart()), new ContigPos(this.recordTable.getItems().get(last_index).getContig(), this.recordTable.getItems().get(last_index).getEnd()));
            if (this.recordTable.getItems().get(0).getContig().equals(this.recordTable.getItems().get(last_index).getContig())) {
                this.gotoField.setText(this.recordTable.getItems().get(0).getContig() + ":" + this.recordTable.getItems().get(0).getStart() + "-" + this.recordTable.getItems().get(last_index).getEnd());
            }
        } else {
            super.seqDictionaryCanvas.setItemsInterval(null, null);
        }
    } else {
        super.seqDictionaryCanvas.setItemsInterval(null, null);
    }
    /* show hide columns for paired end data if no paired data found */
    {
        final boolean contains_paired = this.recordTable.getItems().stream().anyMatch(S -> S.getReadPairedFlag());
        for (final TableColumn<SAMRecord, ?> tc : this.pairedEndColumns) tc.setVisible(contains_paired);
    }
    repaintCanvas();
}

Example 73 with CigarElement

use of htsjdk.samtools.CigarElement in project jvarkit by lindenb.

the class BamStage method repaintCanvas.

/**
 * repaint the canvas area
 */
private void repaintCanvas() {
    final boolean showClip = this.canvasShowClip.isSelected();
    final boolean showReadName = this.canvasShowReadName.isSelected();
    final int baseSize = this.canvasBaseSizeCombo.getValue();
    final double canvaswidth = this.canvas.getCanvas().getWidth();
    final double canvasheight = this.canvas.getCanvas().getHeight();
    final GraphicsContext gc = this.canvas.getCanvas().getGraphicsContext2D();
    gc.setFill(Color.WHITE);
    gc.fillRect(0, 0, canvaswidth, canvasheight);
    double y = baseSize * 2;
    final List<SAMRecord> records = getDisplayableSamRecordStream().collect(Collectors.toList());
    if (records.isEmpty())
        return;
    final long genomicIndex = (long) this.canvasScrollHGenomicLoc.getValue();
    final ContigPos contigStart = super.convertGenomicIndexToContigPos(genomicIndex);
    if (contigStart == null)
        return;
    final int chromStart = contigStart.position;
    final int chromLen = (int) (canvaswidth / baseSize);
    if (chromLen == 0)
        return;
    /* pileup record */
    final List<List<XYRecord>> yxrows = new ArrayList<>();
    for (final SAMRecord rec : records) {
        if (!rec.getReferenceName().equals(contigStart.getContig())) {
            continue;
        }
        final XYRecord newxy = new XYRecord(rec);
        if (rec.getStart() > (chromStart + chromLen))
            break;
        if (rec.getEnd() < (chromStart))
            continue;
        int z = 0;
        for (z = 0; z < yxrows.size(); ++z) {
            final List<XYRecord> row = yxrows.get(z);
            final XYRecord lastyx = row.get(row.size() - 1);
            if (lastyx.getEnd() + 2 < newxy.getStart()) {
                newxy.y = z;
                row.add(newxy);
                break;
            }
        }
        if (z == yxrows.size()) {
            newxy.y = z;
            final List<XYRecord> row = new ArrayList<>();
            row.add(newxy);
            yxrows.add(row);
        }
    }
    final Function<Integer, Double> position2pixel = new Function<Integer, Double>() {

        @Override
        public Double apply(Integer pos) {
            return ((pos - (double) chromStart) / (double) chromLen) * canvaswidth;
        }
    };
    final Hershey hershey = new Hershey();
    // paint contig lines
    hershey.paint(gc, contigStart.getContig(), 1, 0, baseSize * contigStart.getContig().length(), baseSize - 2);
    // paint vertical lines
    for (int x = chromStart; x < chromStart + chromLen; ++x) {
        double px = position2pixel.apply(x);
        gc.setStroke(x % 10 == 0 ? Color.BLACK : Color.GRAY);
        gc.setLineWidth(x % 10 == 0 ? 5 : 0.5);
        gc.strokeLine(px, baseSize, px, canvasheight);
        if (x % 10 == 0) {
            String s = String.valueOf(x);
            gc.setLineWidth(1.0);
            hershey.paint(gc, s, px, baseSize, baseSize * s.length(), baseSize - 1);
        }
    }
    gc.setLineWidth(1);
    for (int yy = (int) this.canvasScrolVInCoverage.getValue(); y < canvasheight && yy < yxrows.size(); ++yy, y += baseSize) {
        final List<XYRecord> row = yxrows.get(yy);
        for (final XYRecord rec : row) {
            /* paint record */
            int baseIndex = 0;
            int refIndex = rec.getStart();
            final byte[] bases = rec.record.getReadBases();
            final Function<Integer, String> getBaseAt = new Function<Integer, String>() {

                @Override
                public String apply(Integer readPos) {
                    char c;
                    if (showReadName) {
                        if (rec.record.getReadNameLength() <= readPos)
                            return "";
                        c = rec.record.getReadName().charAt(readPos);
                    } else if (bases == null || bases.length <= readPos) {
                        return "";
                    } else {
                        c = (char) bases[readPos];
                    }
                    c = (rec.record.getReadNegativeStrandFlag() ? Character.toLowerCase(c) : Character.toUpperCase(c));
                    return String.valueOf(c);
                }
            };
            final Function<Integer, Color> getColorAt = new Function<Integer, Color>() {

                public Color apply(final Integer readPos) {
                    if (bases == null || bases.length <= readPos) {
                        return Color.BLACK;
                    }
                    return JfxNgs.BASE2COLOR.apply((char) bases[readPos]);
                }
            };
            gc.setLineWidth(1.0);
            // arrow end
            {
                double endpos = position2pixel.apply(rec.record.getReadNegativeStrandFlag() ? rec.getStart() : rec.getEnd() + 1);
                double radius = baseSize / 4.0;
                gc.setFill(Color.BLACK);
                gc.fillOval(endpos - radius, y + baseSize / 2.0 - radius, radius * 2, radius * 2);
            }
            final Set<Integer> referenceEvents = new HashSet<>();
            for (final CigarElement ce : rec.record.getCigar()) {
                switch(ce.getOperator()) {
                    case P:
                        break;
                    case I:
                        {
                            baseIndex += ce.getLength();
                            referenceEvents.add(refIndex);
                            break;
                        }
                    case D:
                    case N:
                        {
                            gc.setFill(Color.RED);
                            for (int x = 0; x < ce.getLength(); ++x) {
                                gc.fillRect(position2pixel.apply(refIndex), y, baseSize, baseSize - 1);
                                refIndex++;
                            }
                            break;
                        }
                    case H:
                        {
                            if (showClip) {
                                gc.setFill(Color.YELLOW);
                                for (int x = 0; x < ce.getLength(); ++x) {
                                    gc.fillRect(position2pixel.apply(refIndex), y, baseSize, baseSize - 1);
                                    refIndex++;
                                }
                            } else {
                            // NO refIndex+=ce.getLength();
                            }
                            break;
                        }
                    case S:
                        {
                            if (showClip) {
                                for (int x = 0; x < ce.getLength(); ++x) {
                                    gc.setFill(Color.YELLOW);
                                    gc.fillRect(position2pixel.apply(refIndex), y, baseSize, baseSize - 1);
                                    gc.setStroke(getColorAt.apply(baseIndex));
                                    hershey.paint(gc, getBaseAt.apply(baseIndex), position2pixel.apply(refIndex), y, baseSize - 1, baseSize - 2);
                                    refIndex++;
                                    baseIndex++;
                                }
                            } else {
                                baseIndex += ce.getLength();
                            // NO refIndex+=ce.getLength();
                            }
                            break;
                        }
                    case EQ:
                    case X:
                    case M:
                        {
                            for (int x = 0; x < ce.getLength(); ++x) {
                                gc.setFill(ce.getOperator() == CigarOperator.X ? Color.RED : Color.LIGHTGRAY);
                                gc.fillRect(position2pixel.apply(refIndex), y, baseSize, baseSize - 1);
                                gc.setStroke(getColorAt.apply(baseIndex));
                                hershey.paint(gc, getBaseAt.apply(baseIndex), position2pixel.apply(refIndex), y, baseSize - 1, baseSize - 2);
                                refIndex++;
                                baseIndex++;
                            }
                            break;
                        }
                    default:
                        break;
                }
                if (refIndex > chromStart + chromLen)
                    break;
            }
            gc.setStroke(Color.BLACK);
            gc.strokeRect(position2pixel.apply(rec.getStart()), y, position2pixel.apply(rec.getEnd() + 1) - position2pixel.apply(rec.getStart()), baseSize - 1);
            for (final Integer pos : referenceEvents) {
                double x = position2pixel.apply(pos);
                gc.setStroke(Color.RED);
                gc.setLineWidth(0.5);
                gc.strokeLine(x, y, x, y + baseSize);
            }
        /* end paint record */
        }
    }
    gc.setStroke(Color.BLACK);
    gc.rect(0, 0, canvaswidth - 1, canvasheight - 1);
}

Example 74 with CigarElement

use of htsjdk.samtools.CigarElement in project jvarkit by lindenb.

the class CigarOpPerPositionChartFactory method visit.

@Override
public void visit(final SAMRecord rec) {
    if (rec.getReadUnmappedFlag() || rec.getCigar() == null)
        return;
    int readpos = 0;
    for (final CigarElement ce : rec.getCigar()) {
        switch(ce.getOperator()) {
            case P:
                break;
            case D:
            case N:
                {
                    _visit(ce.getOperator(), readpos);
                    readpos++;
                    break;
                }
            default:
                for (int i = 0; i < ce.getLength(); ++i) {
                    _visit(ce.getOperator(), readpos);
                    readpos++;
                }
                break;
        }
    }
}

Example 75 with CigarElement

use of htsjdk.samtools.CigarElement in project jvarkit by lindenb.

the class VCFAnnoBam method process.

private void process(Rgn rgn, List<SamReader> samReaders) {
    rgn.processed = true;
    int chromStart1 = rgn.interval.getStart();
    int chromEnd1 = rgn.interval.getEnd();
    int[] counts = new int[chromEnd1 - chromStart1 + 1];
    if (counts.length == 0)
        return;
    Arrays.fill(counts, 0);
    for (SamReader samReader : samReaders) {
        /**
         *     start - 1-based, inclusive start of interval of interest. Zero implies start of the reference sequence.
         *	   end - 1-based, inclusive end of interval of interest. Zero implies end of the reference sequence.
         */
        SAMRecordIterator r = samReader.queryOverlapping(rgn.interval.getContig(), chromStart1, chromEnd1);
        while (r.hasNext()) {
            SAMRecord rec = r.next();
            if (rec.getReadUnmappedFlag())
                continue;
            if (this.filter.filterOut(rec))
                continue;
            if (!rec.getReferenceName().equals(rgn.interval.getContig()))
                continue;
            Cigar cigar = rec.getCigar();
            if (cigar == null)
                continue;
            int refpos1 = rec.getAlignmentStart();
            for (CigarElement ce : cigar.getCigarElements()) {
                switch(ce.getOperator()) {
                    case H:
                        break;
                    case S:
                        break;
                    case I:
                        break;
                    case P:
                        break;
                    // reference skip
                    case N:
                    case // deletion in reference
                    D:
                        {
                            refpos1 += ce.getLength();
                            break;
                        }
                    case M:
                    case EQ:
                    case X:
                        {
                            for (int i = 0; i < ce.getLength() && refpos1 <= chromEnd1; ++i) {
                                if (refpos1 >= chromStart1 && refpos1 <= chromEnd1) {
                                    counts[refpos1 - chromStart1]++;
                                }
                                refpos1++;
                            }
                            break;
                        }
                    default:
                        throw new IllegalStateException("Doesn't know how to handle cigar operator:" + ce.getOperator() + " cigar:" + cigar);
                }
            }
        }
        r.close();
    }
    Arrays.sort(counts);
    for (int cov : counts) {
        if (cov <= MIN_COVERAGE)
            rgn.count_no_coverage++;
        rgn.mean += cov;
    }
    rgn.mean /= counts.length;
    rgn.min = counts[0];
    rgn.max = counts[counts.length - 1];
    rgn.percent_covered = (int) (((counts.length - rgn.count_no_coverage) / (double) counts.length) * 100.0);
    rgn.processed = true;
}