Examples with SimplePosition com.github.lindenb.jvarkit.samtools.util.SimplePosition used on opensource projects

Example 1 with SimplePosition

use of com.github.lindenb.jvarkit.samtools.util.SimplePosition in project jvarkit by lindenb.

the class Biostar81455 method doWork.

@Override
public int doWork(final List<String> args) {
    BufferedReader r = null;
    String line;
    PrintStream out = null;
    final CharSplitter tab = CharSplitter.TAB;
    try {
        try (final GtfReader gtfReader = new GtfReader(this.gtfPath)) {
            gtfReader.getAllGenes().stream().forEach(G -> this.geneMap.put(new Interval(G), G));
        }
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(r);
    }
    final ContigNameConverter contigNameConverter = ContigNameConverter.fromIntervalTreeMap(this.geneMap);
    try {
        r = super.openBufferedReader(oneFileOrNull(args));
        out = openPathOrStdoutAsPrintStream(this.outputFile);
        while ((line = r.readLine()) != null) {
            if (line.startsWith("#")) {
                out.println(line);
                continue;
            }
            boolean found = false;
            final String[] tokens = tab.split(line);
            final int pos1 = Integer.parseInt(tokens[1]) + (this.one_based ? 0 : 1);
            final String convertCtg = contigNameConverter.apply(tokens[0]);
            if (convertCtg == null) {
                LOG.error("CANNOT FIND contig " + tokens[0]);
                out.println("##UNKNOWN CONTIG " + line);
                continue;
            }
            final SimplePosition position = new SimplePosition(convertCtg, pos1);
            final List<Transcript> transcripts = this.geneMap.getOverlapping(position).stream().flatMap(G -> G.getTranscripts().stream()).filter(T -> T.overlaps(position)).collect(Collectors.toList());
            if (transcripts.isEmpty()) {
                LOG.info("no gene found in chromosome " + tokens[0] + " (check chrom prefix?)");
            } else {
                for (final Transcript kg : transcripts) {
                    Exon bestExon = null;
                    for (final Exon exon : kg.getExons()) {
                        if (bestExon == null || Math.abs(distance(position.getPosition(), exon)) < Math.abs(distance(position.getPosition(), bestExon))) {
                            bestExon = exon;
                        }
                    }
                    if (bestExon != null) {
                        out.print(line);
                        out.print("\t");
                        out.print(bestExon.getTranscript().getId());
                        out.print("\t");
                        out.print(kg.getStart() - 1);
                        out.print("\t");
                        out.print(kg.getEnd());
                        out.print("\t");
                        out.print(kg.getStrand());
                        out.print("\t");
                        out.print(bestExon.getName());
                        out.print("\t");
                        out.print(bestExon.getStart() - 1);
                        out.print("\t");
                        out.print(bestExon.getEnd());
                        out.print("\t");
                        out.print(distance(position.getPosition(), bestExon));
                        out.println();
                        found = true;
                    }
                }
            }
            if (!found) {
                out.println(line + "\tNULL");
            }
        }
        return RETURN_OK;
    } catch (final Exception err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(r);
        CloserUtil.close(out);
    }
}

Example 2 with SimplePosition

use of com.github.lindenb.jvarkit.samtools.util.SimplePosition in project jvarkit by lindenb.

the class FindAVariation method scanRemote.

private void scanRemote(final String url) {
    VCFReader tabix = null;
    try {
        tabix = VCFReaderFactory.makeDefault().open(url);
        final VCFHeader header = tabix.getHeader();
        final SAMSequenceDictionary dict = header.getSequenceDictionary();
        for (final SimplePosition m : convertFromVcfHeader(url, header)) {
            if (dict.getSequence(m.getContig()) == null)
                continue;
            final java.util.Iterator<VariantContext> iter2 = tabix.query(new Interval(m));
            while (iter2.hasNext()) {
                final VariantContext ctx = iter2.next();
                if (this.onlySnp && (ctx.getStart() != m.getPosition() || ctx.getEnd() != m.getPosition()))
                    continue;
                report(url, header, ctx, m);
            }
            CloserUtil.close(iter2);
        }
        tabix.close();
        tabix = null;
    } catch (final htsjdk.tribble.TribbleException.InvalidHeader err) {
        LOG.warn(url + "\t" + err.getMessage());
    } catch (final Throwable err) {
        LOG.severe("cannot read " + url, err);
    } finally {
        CloserUtil.close(tabix);
    }
}

Example 3 with SimplePosition

use of com.github.lindenb.jvarkit.samtools.util.SimplePosition in project jvarkit by lindenb.

the class FindAVariation method scanPath.

private void scanPath(final String vcfPathString) {
    final Path vcfPath = Paths.get(vcfPathString);
    if (Files.isDirectory(vcfPath))
        return;
    if (!Files.isReadable(vcfPath))
        return;
    VCFIterator iter = null;
    VCFReader r = null;
    try {
        if (vcfPathString.endsWith(FileExtensions.BCF) && BcfToolsUtils.isBcfToolsRequired(Paths.get(vcfPathString))) {
            if (!this.use_bcf)
                return;
        }
        if (isIndexed(vcfPath)) {
            r = VCFReaderFactory.makeDefault().open(vcfPath, true);
            final VCFHeader header = r.getHeader();
            for (final SimplePosition m : convertFromVcfHeader(vcfPath.toString(), header)) {
                try (CloseableIterator<VariantContext> iter2 = r.query(m)) {
                    while (iter2.hasNext()) {
                        final VariantContext ctx = iter2.next();
                        if (this.onlySnp && (ctx.getStart() != m.getPosition() || ctx.getEnd() != m.getPosition()))
                            continue;
                        report(vcfPathString, header, ctx, m);
                    }
                }
            }
            r.close();
            r = null;
        } else if (!this.indexedOnly) {
            iter = VCFUtils.createVCFIteratorFromPath(vcfPath);
            final VCFHeader header = iter.getHeader();
            final Set<SimplePosition> mutlist = convertFromVcfHeader(vcfPath.toString(), iter.getHeader());
            while (iter.hasNext()) {
                final VariantContext ctx = iter.next();
                if (this.onlySnp && ctx.getLengthOnReference() != 1)
                    continue;
                for (final SimplePosition m2 : mutlist) {
                    if (!m2.overlaps(ctx))
                        continue;
                    if (this.onlySnp && (ctx.getStart() != m2.getPosition() || ctx.getEnd() != m2.getPosition()))
                        continue;
                    report(vcfPathString, header, ctx, m2);
                }
            }
            iter.close();
            iter = null;
        }
    } catch (final htsjdk.tribble.TribbleException.InvalidHeader err) {
        LOG.warn(vcfPathString + "\t" + err.getMessage());
    } catch (final Throwable err) {
        LOG.severe("cannot read " + vcfPathString, err);
    } finally {
        CloserUtil.close(r);
        CloserUtil.close(iter);
    }
}

Example 4 with SimplePosition

use of com.github.lindenb.jvarkit.samtools.util.SimplePosition in project jvarkit by lindenb.

the class ValidateCnv method doWork.

@Override
public int doWork(final List<String> args) {
    if (this.extendFactor <= 0) {
        LOG.error("bad extend factor " + this.extendFactor);
        return -1;
    }
    if (this.treshold < 0 || this.treshold >= 0.25) {
        LOG.error("Bad treshold 0 < " + this.treshold + " >=0.25 ");
        return -1;
    }
    final Map<String, BamInfo> sample2bam = new HashMap<>();
    VariantContextWriter out = null;
    Iterator<? extends Locatable> iterIn = null;
    try {
        final SAMSequenceDictionary dict = SequenceDictionaryUtils.extractRequired(this.rererencePath);
        final CRAMReferenceSource cramReferenceSource = new ReferenceSource(this.rererencePath);
        final List<Path> bamPaths = IOUtils.unrollPaths(this.bamFiles);
        final String input = oneFileOrNull(args);
        if (input == null) {
            iterIn = IntervalListProvider.empty().dictionary(dict).skipUnknownContigs().fromInputStream(stdin(), "bed").iterator();
        } else {
            final IntervalListProvider ilp = IntervalListProvider.from(input).setVariantPredicate(CTX -> {
                if (CTX.isSNP())
                    return false;
                final String svType = CTX.getAttributeAsString(VCFConstants.SVTYPE, "");
                if (svType != null && (svType.equals("INV") || svType.equals("BND")))
                    return false;
                return true;
            }).dictionary(dict).skipUnknownContigs();
            iterIn = ilp.stream().iterator();
        }
        /* register each bam */
        for (final Path p2 : bamPaths) {
            final BamInfo bi = new BamInfo(p2, cramReferenceSource);
            if (sample2bam.containsKey(bi.sampleName)) {
                LOG.error("sample " + bi.sampleName + " specified twice.");
                bi.close();
                return -1;
            }
            sample2bam.put(bi.sampleName, bi);
        }
        if (sample2bam.isEmpty()) {
            LOG.error("no bam was defined");
            return -1;
        }
        final Set<VCFHeaderLine> metadata = new HashSet<>();
        final VCFInfoHeaderLine infoSVSamples = new VCFInfoHeaderLine("N_SAMPLES", 1, VCFHeaderLineType.Integer, "Number of Samples that could carry a SV");
        metadata.add(infoSVSamples);
        final VCFInfoHeaderLine infoSvLen = new VCFInfoHeaderLine("SVLEN", 1, VCFHeaderLineType.Integer, "SV length");
        metadata.add(infoSvLen);
        final BiFunction<String, String, VCFFormatHeaderLine> makeFmt = (TAG, DESC) -> new VCFFormatHeaderLine(TAG, 1, VCFHeaderLineType.Integer, DESC);
        final VCFFormatHeaderLine formatCN = new VCFFormatHeaderLine("CN", 1, VCFHeaderLineType.Float, "normalized copy-number. Treshold was " + this.treshold);
        metadata.add(formatCN);
        final VCFFormatHeaderLine nReadsSupportingSv = makeFmt.apply("RSD", "number of split reads supporting SV.");
        metadata.add(nReadsSupportingSv);
        final VCFFilterHeaderLine filterAllDel = new VCFFilterHeaderLine("ALL_DEL", "number of samples greater than 1 and all are deletions");
        metadata.add(filterAllDel);
        final VCFFilterHeaderLine filterAllDup = new VCFFilterHeaderLine("ALL_DUP", "number of samples  greater than  1 and all are duplication");
        metadata.add(filterAllDup);
        final VCFFilterHeaderLine filterNoSV = new VCFFilterHeaderLine("NO_SV", "There is no DUP or DEL in this variant");
        metadata.add(filterNoSV);
        final VCFFilterHeaderLine filterHomDel = new VCFFilterHeaderLine("HOM_DEL", "There is one Homozygous deletion.");
        metadata.add(filterHomDel);
        final VCFFilterHeaderLine filterHomDup = new VCFFilterHeaderLine("HOM_DUP", "There is one Homozygous duplication.");
        metadata.add(filterHomDup);
        VCFStandardHeaderLines.addStandardFormatLines(metadata, true, VCFConstants.DEPTH_KEY, VCFConstants.GENOTYPE_KEY, VCFConstants.GENOTYPE_FILTER_KEY, VCFConstants.GENOTYPE_QUALITY_KEY);
        VCFStandardHeaderLines.addStandardInfoLines(metadata, true, VCFConstants.DEPTH_KEY, VCFConstants.END_KEY, VCFConstants.ALLELE_COUNT_KEY, VCFConstants.ALLELE_FREQUENCY_KEY, VCFConstants.ALLELE_NUMBER_KEY);
        final VCFHeader header = new VCFHeader(metadata, sample2bam.keySet());
        if (dict != null)
            header.setSequenceDictionary(dict);
        JVarkitVersion.getInstance().addMetaData(this, header);
        final ProgressFactory.Watcher<VariantContext> progress = ProgressFactory.newInstance().dictionary(dict).logger(LOG).build();
        out = this.writingVariantsDelegate.dictionary(dict).open(this.outputFile);
        out.writeHeader(header);
        final Allele DUP_ALLELE = Allele.create("<DUP>", false);
        final Allele DEL_ALLELE = Allele.create("<DEL>", false);
        final Allele REF_ALLELE = Allele.create("N", true);
        while (iterIn.hasNext()) {
            final Locatable ctx = iterIn.next();
            if (ctx == null)
                continue;
            final SAMSequenceRecord ssr = dict.getSequence(ctx.getContig());
            if (ssr == null || ctx.getStart() >= ssr.getSequenceLength())
                continue;
            final int svLen = ctx.getLengthOnReference();
            if (svLen < this.min_abs_sv_size)
                continue;
            if (svLen > this.max_abs_sv_size)
                continue;
            int n_samples_with_cnv = 0;
            final SimplePosition breakPointLeft = new SimplePosition(ctx.getContig(), ctx.getStart());
            final SimplePosition breakPointRight = new SimplePosition(ctx.getContig(), ctx.getEnd());
            final int extend = 1 + (int) (svLen * this.extendFactor);
            final int leftPos = Math.max(1, breakPointLeft.getPosition() - extend);
            final int array_mid_start = breakPointLeft.getPosition() - leftPos;
            final int array_mid_end = breakPointRight.getPosition() - leftPos;
            final int rightPos = Math.min(breakPointRight.getPosition() + extend, ssr.getSequenceLength());
            final VariantContextBuilder vcb = new VariantContextBuilder();
            vcb.chr(ctx.getContig());
            vcb.start(ctx.getStart());
            vcb.stop(ctx.getEnd());
            vcb.attribute(VCFConstants.END_KEY, ctx.getEnd());
            final Set<Allele> alleles = new HashSet<>();
            alleles.add(REF_ALLELE);
            int count_dup = 0;
            int count_del = 0;
            int an = 0;
            final Counter<Allele> countAlleles = new Counter<>();
            final List<Genotype> genotypes = new ArrayList<>(sample2bam.size());
            Double badestGQ = null;
            final double[] raw_coverage = new double[CoordMath.getLength(leftPos, rightPos)];
            for (final String sampleName : sample2bam.keySet()) {
                final BamInfo bi = sample2bam.get(sampleName);
                Arrays.fill(raw_coverage, 0.0);
                int n_reads_supporting_sv = 0;
                try (CloseableIterator<SAMRecord> iter2 = bi.samReader.queryOverlapping(ctx.getContig(), leftPos, rightPos)) {
                    while (iter2.hasNext()) {
                        final SAMRecord rec = iter2.next();
                        if (!SAMRecordDefaultFilter.accept(rec, this.min_mapq))
                            continue;
                        final Cigar cigar = rec.getCigar();
                        if (cigar == null || cigar.isEmpty())
                            continue;
                        // any clip supporting deletion ?
                        boolean read_supports_cnv = false;
                        final int breakpoint_distance = 10;
                        // any clip on left ?
                        if (cigar.isLeftClipped() && rec.getUnclippedStart() < rec.getAlignmentStart() && new SimpleInterval(ctx.getContig(), rec.getUnclippedStart(), rec.getAlignmentStart()).withinDistanceOf(breakPointLeft, breakpoint_distance)) {
                            read_supports_cnv = true;
                        }
                        // any clip on right ?
                        if (!read_supports_cnv && cigar.isRightClipped() && rec.getAlignmentEnd() < rec.getUnclippedEnd() && new SimpleInterval(ctx.getContig(), rec.getAlignmentEnd(), rec.getUnclippedEnd()).withinDistanceOf(breakPointRight, breakpoint_distance)) {
                            read_supports_cnv = true;
                        }
                        if (read_supports_cnv) {
                            n_reads_supporting_sv++;
                        }
                        int ref = rec.getStart();
                        for (final CigarElement ce : cigar) {
                            final CigarOperator op = ce.getOperator();
                            if (op.consumesReferenceBases()) {
                                if (op.consumesReadBases()) {
                                    for (int x = 0; x < ce.getLength() && ref + x - leftPos < raw_coverage.length; ++x) {
                                        final int p = ref + x - leftPos;
                                        if (p < 0 || p >= raw_coverage.length)
                                            continue;
                                        raw_coverage[p]++;
                                    }
                                }
                                ref += ce.getLength();
                            }
                        }
                    }
                // end while iter record
                }
                // end try query for iterator
                // test for great difference between DP left and DP right
                final OptionalDouble medianDepthLeft = Percentile.median().evaluate(raw_coverage, 0, array_mid_start);
                final OptionalDouble medianDepthRight = Percentile.median().evaluate(raw_coverage, array_mid_end, raw_coverage.length - array_mid_end);
                // any is just too low
                if (!medianDepthLeft.isPresent() || medianDepthLeft.getAsDouble() < this.min_depth || !medianDepthRight.isPresent() || medianDepthRight.getAsDouble() < this.min_depth) {
                    final Genotype gt2 = new GenotypeBuilder(sampleName, Arrays.asList(Allele.NO_CALL, Allele.NO_CALL)).filter("LowDp").make();
                    genotypes.add(gt2);
                    continue;
                }
                final double difference_factor = 2.0;
                // even if a value is divided , it remains greater than the other size
                if (medianDepthLeft.getAsDouble() / difference_factor > medianDepthRight.getAsDouble() || medianDepthRight.getAsDouble() / difference_factor > medianDepthLeft.getAsDouble()) {
                    final Genotype gt2 = new GenotypeBuilder(sampleName, Arrays.asList(Allele.NO_CALL, Allele.NO_CALL)).filter("DiffLR").make();
                    genotypes.add(gt2);
                    continue;
                }
                // run median to smooth spline
                final double[] smoothed_cov = new RunMedian(RunMedian.getTurlachSize(raw_coverage.length)).apply(raw_coverage);
                final double[] bounds_cov = IntStream.concat(IntStream.range(0, array_mid_start), IntStream.range(array_mid_end, smoothed_cov.length)).mapToDouble(IDX -> raw_coverage[IDX]).toArray();
                final OptionalDouble optMedianBound = Percentile.median().evaluate(bounds_cov);
                if (!optMedianBound.isPresent() || optMedianBound.getAsDouble() == 0) {
                    final Genotype gt2 = new GenotypeBuilder(sampleName, Arrays.asList(Allele.NO_CALL, Allele.NO_CALL)).filter("MedZero").make();
                    genotypes.add(gt2);
                    continue;
                }
                final double medianBound = optMedianBound.getAsDouble();
                // divide coverage per medianBound
                final double[] normalized_mid_coverage = new double[array_mid_end - array_mid_start];
                for (int i = 0; i < normalized_mid_coverage.length; ++i) {
                    normalized_mid_coverage[i] = smoothed_cov[array_mid_start + i] / medianBound;
                }
                final double normDepth = Percentile.median().evaluate(normalized_mid_coverage).getAsDouble();
                final boolean is_sv;
                final boolean is_hom_deletion = Math.abs(normDepth - 0.0) <= this.treshold;
                final boolean is_het_deletion = Math.abs(normDepth - 0.5) <= this.treshold || (!is_hom_deletion && normDepth <= 0.5);
                final boolean is_hom_dup = Math.abs(normDepth - 2.0) <= this.treshold || normDepth > 2.0;
                final boolean is_het_dup = Math.abs(normDepth - 1.5) <= this.treshold || (!is_hom_dup && normDepth >= 1.5);
                final boolean is_ref = Math.abs(normDepth - 1.0) <= this.treshold;
                final double theoritical_depth;
                final GenotypeBuilder gb;
                if (is_ref) {
                    gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, REF_ALLELE));
                    is_sv = false;
                    theoritical_depth = 1.0;
                    an += 2;
                } else if (is_het_deletion) {
                    gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, DEL_ALLELE));
                    alleles.add(DEL_ALLELE);
                    is_sv = true;
                    theoritical_depth = 0.5;
                    count_del++;
                    an += 2;
                    countAlleles.incr(DEL_ALLELE);
                } else if (is_hom_deletion) {
                    gb = new GenotypeBuilder(sampleName, Arrays.asList(DEL_ALLELE, DEL_ALLELE));
                    alleles.add(DEL_ALLELE);
                    vcb.filter(filterHomDel.getID());
                    is_sv = true;
                    theoritical_depth = 0.0;
                    count_del++;
                    an += 2;
                    countAlleles.incr(DEL_ALLELE, 2);
                } else if (is_het_dup) {
                    gb = new GenotypeBuilder(sampleName, Arrays.asList(REF_ALLELE, DUP_ALLELE));
                    alleles.add(DUP_ALLELE);
                    is_sv = true;
                    theoritical_depth = 1.5;
                    count_dup++;
                    an += 2;
                    countAlleles.incr(DUP_ALLELE);
                } else if (is_hom_dup) {
                    gb = new GenotypeBuilder(sampleName, Arrays.asList(DUP_ALLELE, DUP_ALLELE));
                    alleles.add(DUP_ALLELE);
                    vcb.filter(filterHomDup.getID());
                    is_sv = true;
                    theoritical_depth = 2.0;
                    count_dup++;
                    an += 2;
                    countAlleles.incr(DUP_ALLELE, 2);
                } else {
                    gb = new GenotypeBuilder(sampleName, Arrays.asList(Allele.NO_CALL, Allele.NO_CALL)).filter("Ambigous");
                    is_sv = false;
                    theoritical_depth = 1.0;
                }
                if (is_sv) {
                    n_samples_with_cnv++;
                }
                double gq = Math.abs(theoritical_depth - normDepth);
                gq = Math.min(0.5, gq);
                gq = gq * gq;
                gq = gq / 0.25;
                gq = 99 * (1.0 - gq);
                gb.GQ((int) gq);
                if (badestGQ == null || badestGQ.compareTo(gq) > 0) {
                    badestGQ = gq;
                }
                gb.attribute(formatCN.getID(), normDepth);
                gb.attribute(nReadsSupportingSv.getID(), n_reads_supporting_sv);
                genotypes.add(gb.make());
            }
            vcb.attribute(VCFConstants.ALLELE_NUMBER_KEY, an);
            final List<Allele> orderedAlleles = new ArrayList<>(alleles);
            Collections.sort(orderedAlleles);
            if (orderedAlleles.size() > 1) {
                final List<Integer> acL = new ArrayList<>();
                final List<Double> afL = new ArrayList<>();
                for (int i = 1; i < orderedAlleles.size(); i++) {
                    final Allele a = orderedAlleles.get(i);
                    final int c = (int) countAlleles.count(a);
                    acL.add(c);
                    if (an > 0)
                        afL.add(c / (double) an);
                }
                vcb.attribute(VCFConstants.ALLELE_COUNT_KEY, acL);
                if (an > 0)
                    vcb.attribute(VCFConstants.ALLELE_FREQUENCY_KEY, afL);
            }
            // if(alleles.size()<=1) continue;
            vcb.alleles(orderedAlleles);
            vcb.noID();
            vcb.genotypes(genotypes);
            vcb.attribute(infoSVSamples.getID(), n_samples_with_cnv);
            vcb.attribute(infoSvLen.getID(), svLen);
            if (count_dup == sample2bam.size() && sample2bam.size() != 1) {
                vcb.filter(filterAllDup.getID());
            }
            if (count_del == sample2bam.size() && sample2bam.size() != 1) {
                vcb.filter(filterAllDel.getID());
            }
            if (n_samples_with_cnv == 0) {
                vcb.filter(filterNoSV.getID());
            }
            if (badestGQ != null) {
                vcb.log10PError(badestGQ / -10.0);
            }
            out.add(vcb.make());
        }
        progress.close();
        out.close();
        return 0;
    } catch (final Throwable err) {
        LOG.error(err);
        return -1;
    } finally {
        CloserUtil.close(iterIn);
        CloserUtil.close(out);
        sample2bam.values().forEach(F -> CloserUtil.close(F));
    }
}

Example 5 with SimplePosition

use of com.github.lindenb.jvarkit.samtools.util.SimplePosition in project jvarkit by lindenb.

the class StructuralVariantComparator method testBNDDistance.

private boolean testBNDDistance(final Locatable a, final Locatable b) {
    final Function<Locatable, Locatable> bnd2interval = LOC -> {
        if (!(LOC instanceof VariantContext)) {
            return new SimplePosition(LOC.getContig(), LOC.getStart());
        }
        final VariantContext CTX = VariantContext.class.cast(LOC);
        // can be negative
        int extend_5 = 0;
        // can be negative
        int extend_3 = 0;
        if (CTX.hasAttribute("CIPOS")) {
            try {
                // can be negative
                final List<Integer> xList = CTX.getAttributeAsIntList("CIPOS", 0);
                if (xList.size() == 2) {
                    extend_5 = xList.get(0).intValue();
                    extend_3 = xList.get(1).intValue();
                    if (extend_5 > 0)
                        extend_5 = 0;
                    if (extend_3 < 0)
                        extend_3 = 0;
                }
            } catch (final Throwable err) {
                extend_5 = 0;
                extend_3 = 0;
            }
        }
        return new SimpleInterval(CTX.getContig(), CTX.getStart() + extend_5, /* can be negative */
        CTX.getEnd() + extend_3);
    };
    // we use SimplePosition because getEnd() can be large when the BND is on the same chromosome
    return this.withinDistanceOf(bnd2interval.apply(a), bnd2interval.apply(b), this.bnd_max_distance);
}