Java IntervalList.size方法代码示例

import htsjdk.samtools.util.IntervalList; //导入方法依赖的package包/类
public static OverlapDetector<Interval> makeOverlapDetector(final File samFile, final SAMFileHeader header, final File ribosomalIntervalsFile, final Log log) {

        final OverlapDetector<Interval> ribosomalSequenceOverlapDetector = new OverlapDetector<Interval>(0, 0);
        if (ribosomalIntervalsFile != null) {

            final IntervalList ribosomalIntervals = IntervalList.fromFile(ribosomalIntervalsFile);
            if (ribosomalIntervals.size() == 0) {
                log.warn("The RIBOSOMAL_INTERVALS file, " + ribosomalIntervalsFile.getAbsolutePath() + " does not contain intervals");
            }
            try {
                SequenceUtil.assertSequenceDictionariesEqual(header.getSequenceDictionary(), ribosomalIntervals.getHeader().getSequenceDictionary());
            } catch (SequenceUtil.SequenceListsDifferException e) {
                throw new PicardException("Sequence dictionaries differ in " + samFile.getAbsolutePath() + " and " + ribosomalIntervalsFile.getAbsolutePath(),
                        e);
            }
            final IntervalList uniquedRibosomalIntervals = ribosomalIntervals.uniqued();
            final List<Interval> intervals = uniquedRibosomalIntervals.getIntervals();
            ribosomalSequenceOverlapDetector.addAll(intervals, intervals);
        }
        return ribosomalSequenceOverlapDetector;
    }
 

import htsjdk.samtools.util.IntervalList; //导入方法依赖的package包/类
/** Calculates a few statistics about the bait design that can then be output. */
void calculateStatistics(final IntervalList targets, final IntervalList baits) {
    this.TARGET_TERRITORY = (int) targets.getUniqueBaseCount();
    this.TARGET_COUNT = targets.size();
    this.BAIT_TERRITORY = (int) baits.getUniqueBaseCount();
    this.BAIT_COUNT = baits.size();
    this.DESIGN_EFFICIENCY = this.TARGET_TERRITORY / (double) this.BAIT_TERRITORY;

    // Figure out the intersection between all targets and all baits
    final IntervalList tmp = new IntervalList(targets.getHeader());
    final OverlapDetector<Interval> detector = new OverlapDetector<Interval>(0, 0);
    detector.addAll(baits.getIntervals(), baits.getIntervals());

    for (final Interval target : targets) {
        final Collection<Interval> overlaps = detector.getOverlaps(target);

        if (overlaps.isEmpty()) {
            this.ZERO_BAIT_TARGETS++;
        } else {
            for (final Interval i : overlaps) tmp.add(target.intersect(i));
        }
    }

    tmp.uniqued();
    this.BAIT_TARGET_TERRITORY_INTERSECTION = (int) tmp.getBaseCount();
}
 

import htsjdk.samtools.util.IntervalList; //导入方法依赖的package包/类
/**
 * Returns an {@link AllelicCountCollection} based on the pileup at sites (specified by an interval list)
 * in a sorted BAM file.  Reads and bases below the specified mapping quality and base quality, respectively,
 * are filtered out of the pileup.  The alt count is defined as the total count minus the ref count, and the
 * alt nucleotide is defined as the non-ref base with the highest count, with ties broken by the order of the
 * bases in {@link AllelicCountCollector#BASES}.
 * @param bamFile           sorted BAM file
 * @param siteIntervals     interval list of sites
 * @param minMappingQuality minimum mapping quality required for reads to be included in pileup
 * @param minBaseQuality    minimum base quality required for bases to be included in pileup
 * @return                  AllelicCountCollection of ref/alt counts at sites in BAM file
 */
public AllelicCountCollection collect(final File bamFile,
                                      final IntervalList siteIntervals,
                                      final int minMappingQuality,
                                      final int minBaseQuality) {
    try (final SamReader reader = readerFactory.open(bamFile)) {
        ParamUtils.isPositiveOrZero(minMappingQuality, "Minimum mapping quality must be nonnegative.");
        ParamUtils.isPositiveOrZero(minBaseQuality, "Minimum base quality must be nonnegative.");
        if (reader.getFileHeader().getSortOrder() != SAMFileHeader.SortOrder.coordinate) {
            throw new UserException.BadInput("BAM file " + bamFile.toString() + " must be coordinate sorted.");
        }

        final int numberOfSites = siteIntervals.size();
        final boolean useIndex = numberOfSites < MAX_INTERVALS_FOR_INDEX;
        final SamLocusIterator locusIterator = new SamLocusIterator(reader, siteIntervals, useIndex);

        //set read and locus filters [note: read counts match IGV, but off by a few from pysam.mpileup]
        final List<SamRecordFilter> samFilters = Arrays.asList(new NotPrimaryAlignmentFilter(), new DuplicateReadFilter());
        locusIterator.setSamFilters(samFilters);
        locusIterator.setEmitUncoveredLoci(true);
        locusIterator.setIncludeNonPfReads(false);
        locusIterator.setMappingQualityScoreCutoff(minMappingQuality);
        locusIterator.setQualityScoreCutoff(minBaseQuality);

        logger.info("Examining " + numberOfSites + " sites in total...");
        int locusCount = 0;
        final AllelicCountCollection counts = new AllelicCountCollection();
        for (final SamLocusIterator.LocusInfo locus : locusIterator) {
            if (locusCount % NUMBER_OF_SITES_PER_LOGGED_STATUS_UPDATE == 0) {
                logger.info("Examined " + locusCount + " sites.");
            }
            locusCount++;

            final Nucleotide refBase = Nucleotide.valueOf(referenceWalker.get(locus.getSequenceIndex()).getBases()[locus.getPosition() - 1]);
            if (!BASES.contains(refBase)) {
                logger.warn(String.format("The reference position at %d has an unknown base call (value: %s). Skipping...",
                        locus.getPosition(), refBase.toString()));
                continue;
            }

            final Nucleotide.Counter baseCounts = getPileupBaseCounts(locus);
            final int totalBaseCount = BASES.stream().mapToInt(b -> (int) baseCounts.get(b)).sum(); //only include total ACGT counts in binomial test (exclude N, etc.)
            final int refReadCount = (int) baseCounts.get(refBase);
            final int altReadCount = totalBaseCount - refReadCount;                                 //we take alt = total - ref instead of the actual alt count
            final Nucleotide altBase = inferAltFromPileupBaseCounts(baseCounts, refBase);

            counts.add(new AllelicCount(
                    new SimpleInterval(locus.getSequenceName(), locus.getPosition(), locus.getPosition()),
                    refReadCount, altReadCount, refBase, altBase));
        }
        logger.info(locusCount + " sites out of " + numberOfSites + " total sites were examined.");
        return counts;
    } catch (final IOException | SAMFormatException e) {
        throw new UserException("Unable to collect allelic counts from " + bamFile);
    }
}
 

import htsjdk.samtools.util.IntervalList; //导入方法依赖的package包/类
/**
 * Writes out fasta files for each pool and also agilent format files if requested.
 *
 * @param dir      the directory to output files into
 * @param basename the basename of each file
 * @param baits    the set of baits to write out
 */
void writePoolFiles(final File dir, final String basename, final IntervalList baits) {
    final int copies;
    if (FILL_POOLS && baits.size() < POOL_SIZE) copies = (int) Math.floor(POOL_SIZE / (double) baits.size());
    else copies = 1;

    int written = 0;
    int nextPool = 0;
    BufferedWriter out = null;
    BufferedWriter agilentOut = null;
    final String prefix = DESIGN_NAME.substring(0, Math.min(DESIGN_NAME.length(), 8)) + "_"; // prefix for 15 digit bait id
    final NumberFormat fmt = new DecimalFormat("000000");

    try {
        for (int i = 0; i < copies; ++i) {
            final boolean rc = i % 2 == 1;

            int baitId = 1;
            for (final Interval interval : baits) {
                final Bait bait = (Bait) interval;

                if (written++ % POOL_SIZE == 0) {
                    if (out != null) out.close();
                    if (agilentOut != null) agilentOut.close();

                    final String filename = basename + ".pool" + nextPool++ + ".design.";
                    out = IOUtil.openFileForBufferedWriting(new File(dir, filename + "fasta"));
                    if (OUTPUT_AGILENT_FILES) {
                        agilentOut = IOUtil.openFileForBufferedWriting(new File(dir, filename + "txt"));
                    }
                }

                writeBaitFasta(out, interval, rc);
                if (OUTPUT_AGILENT_FILES) {
                    agilentOut.append(prefix).append(fmt.format(baitId++));
                    agilentOut.append("\t");
                    agilentOut.append(getBaitSequence(bait, rc).toUpperCase());
                    agilentOut.newLine();
                }
            }
        }

        CloserUtil.close(out);
        CloserUtil.close(agilentOut);
    } catch (Exception e) {
        throw new PicardException("Error while writing pool files.", e);
    }
}