import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private static void sliceFromS3(@NotNull final CommandLine cmd) throws IOException, EmptyFileException {
    final OkHttpClient httpClient = SbpS3Client.create(Integer.parseInt(cmd.getOptionValue(MAX_CONCURRENT_REQUESTS)));
    final URL bamUrl = SbpS3UrlGenerator.generateUrl(cmd.getOptionValue(BUCKET), cmd.getOptionValue(INPUT));
    final URL indexUrl = SbpS3UrlGenerator.generateUrl(cmd.getOptionValue(BUCKET), cmd.getOptionValue(INDEX));
    final String outputPath = cmd.getOptionValue(OUTPUT);
    final String bedPath = cmd.getOptionValue(BED);
    final int maxBufferSize = readMaxBufferSize(cmd);
    final File indexFile = downloadIndex(indexUrl);
    final SamReader reader = SamReaderFactory.makeDefault().open(SamInputResource.of(bamUrl).index(indexFile));
    LOGGER.info("Generating query intervals from BED file: {}", bedPath);
    final QueryInterval[] intervals = getIntervalsFromBED(bedPath, reader.getFileHeader());
    final BAMFileSpan span = bamSpanForIntervals(indexFile, reader.getFileHeader(), intervals);
    final List<Chunk> expandedChunks = expandChunks(span.getChunks());
    LOGGER.info("Generated {} query intervals which map to {} bam chunks", intervals.length, expandedChunks.size());
    final SamInputResource bamResource =
            SamInputResource.of(new CachingSeekableHTTPStream(httpClient, bamUrl, expandedChunks, maxBufferSize)).index(indexFile);
    final SamReader cachingReader = SamReaderFactory.makeDefault().open(bamResource);

    LOGGER.info("Slicing bam...");
    final CloseableIterator<SAMRecord> iterator = getIterator(cachingReader, intervals, span.toCoordinateArray());
    writeToSlice(outputPath, cachingReader.getFileHeader(), iterator);
    cachingReader.close();
    reader.close();
    indexFile.deleteOnExit();
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private static File queryNameSortedBAM(final SamReader reader, final QueryInterval[] intervals, final String name) throws IOException {

        final SAMFileHeader header = reader.getFileHeader().clone();
        header.setSortOrder(SAMFileHeader.SortOrder.queryname);

        final File file = File.createTempFile(name, ".bam");
        final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(header, false, file);

        final SAMRecordIterator iterator = reader.queryOverlapping(intervals);
        while (iterator.hasNext()) {
            writer.addAlignment(iterator.next());
        }

        iterator.close();
        writer.close();

        return file;
    }
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
/**
 * Converts a List of SimpleIntervals into the format required by the SamReader query API
 * @param rawIntervals SimpleIntervals to be converted
 * @return A sorted, merged list of QueryIntervals suitable for passing to the SamReader query API
 */
static QueryInterval[] prepareQueryIntervals( final List<Interval>
		rawIntervals, final SAMSequenceDictionary sequenceDictionary ) {
	if ( rawIntervals == null || rawIntervals.isEmpty() ) {
		return null;
	}

	// Convert each SimpleInterval to a QueryInterval
	final QueryInterval[] convertedIntervals =
			rawIntervals.stream()
					.map(rawInterval -> convertSimpleIntervalToQueryInterval(rawInterval, sequenceDictionary))
					.toArray(QueryInterval[]::new);

	// Intervals must be optimized (sorted and merged) in order to use the htsjdk query API
	return QueryInterval.optimizeIntervals(convertedIntervals);
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
/**
 * Converts an interval in SimpleInterval format into an htsjdk QueryInterval.
 *
 * In doing so, a header lookup is performed to convert from contig name to index
 *
 * @param interval interval to convert
 * @param sequenceDictionary sequence dictionary used to perform the conversion
 * @return an equivalent interval in QueryInterval format
 */
private static QueryInterval convertSimpleIntervalToQueryInterval( final Interval interval,	final SAMSequenceDictionary sequenceDictionary ) {
	if (interval == null) {
		throw new IllegalArgumentException("interval may not be null");
	}
	if (sequenceDictionary == null) {
		throw new IllegalArgumentException("sequence dictionary may not be null");
	}

	final int contigIndex = sequenceDictionary.getSequenceIndex(interval.getContig());
	if ( contigIndex == -1 ) {
		throw new IllegalArgumentException("Contig " + interval.getContig() + " not present in reads sequence " +
				"dictionary");
	}

	return new QueryInterval(contigIndex, interval.getStart(), interval.getEnd());
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
/**
 * Converts a List of SimpleIntervals into the format required by the SamReader query API
 * @param rawIntervals SimpleIntervals to be converted
 * @return A sorted, merged list of QueryIntervals suitable for passing to the SamReader query API
 */
private QueryInterval[] prepareQueryIntervals( final List<SimpleInterval> rawIntervals ) {
    if ( rawIntervals == null || rawIntervals.isEmpty() ) {
        return null;
    }

    // This might take a while with large interval lists, so log a status message
    logger.debug("Preparing intervals for traversal");

    // Convert each SimpleInterval to a QueryInterval
    final QueryInterval[] convertedIntervals =
            rawIntervals.stream()
            .map(rawInterval -> IntervalUtils.convertSimpleIntervalToQueryInterval(rawInterval, reader.getFileHeader().getSequenceDictionary()))
            .toArray(QueryInterval[]::new);

    // Intervals must be optimized (sorted and merged) in order to use the htsjdk query API
    return QueryInterval.optimizeIntervals(convertedIntervals);
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
/** Parses the BAM file into SAMRecords. Will be distributed onto at least 'numPartitions' partitions. **/
public JavaRDD<SAMRecord> getReads(JavaSparkContext ctx, int numPartitions) {
    try {
        Path bamPath = IOUtils.getPath(bam);
        ChannelAsSeekableStream bamOverNIO = new ChannelAsSeekableStream(Files.newByteChannel(bamPath), bamPath.toString());
        final byte[] index = getIndex();
        SeekableStream indexInMemory = new ByteArraySeekableStream(index);

        SamReader bam3 = SamReaderFactory.makeDefault()
                .validationStringency(ValidationStringency.LENIENT)
                .enable(SamReaderFactory.Option.CACHE_FILE_BASED_INDEXES)
                .open(SamInputResource.of(bamOverNIO).index(indexInMemory));
        List<QueryInterval> chunks = getAllChunksBalanced(bam3, numPartitions);

        // Ideally we'd get exactly the number of chunks the user is asking for, but until then...
        logger.debug("We got: " + chunks.size() + " chunks.");

        return ctx.parallelize(chunks, chunks.size()).flatMap(qi -> new ReadsIterable(bam, index, qi).iterator());
    }
    catch ( IOException e ) {
        throw new GATKException("I/O error loading reads", e);
    }
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
public ReadsIterator() throws IOException {
    Path fpath = IOUtils.getPath(path);
    byte[] indexData = index;
    SeekableStream indexInMemory = new ByteArraySeekableStream(indexData);
    // set high-level retries to deal with servers that might be temporarily overloaded
    // while we're reading a very long file from them.
    SeekableByteChannelPrefetcher chan = new SeekableByteChannelPrefetcher(
        Files.newByteChannel(fpath), BUFSIZE);
    ChannelAsSeekableStream bamOverNIO = new ChannelAsSeekableStream(chan, path);
    bam = SamReaderFactory.makeDefault()
            .validationStringency(ValidationStringency.LENIENT)
            .enable(SamReaderFactory.Option.CACHE_FILE_BASED_INDEXES)
            .open(SamInputResource.of(bamOverNIO).index(indexInMemory));

    QueryInterval[] array = new QueryInterval[1];
    array[0] = interval;
    query = bam.query(array, false);
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
@Override
public SAMRecordIterator query(QueryInterval[] intervals, boolean contained) {
  final GA4GHQueryInterval[] myIntervals = new GA4GHQueryInterval[intervals.length];
  final GA4GHQueryInterval.ReadPositionConstraint constraint = contained ?
    GA4GHQueryInterval.ReadPositionConstraint.CONTAINED : 
    GA4GHQueryInterval.ReadPositionConstraint.OVERLAPPING;
  final SAMFileHeader header = getFileHeader();
  for (int i = 0; i < intervals.length; i++) {
    final QueryInterval interval = intervals[i];
    final String sequence = header.getSequence(
        interval.referenceIndex).getSequenceName();
    myIntervals[i] = new GA4GHQueryInterval(sequence, interval.start, 
        interval.end, constraint);
  }
  return query(myIntervals);
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private List<SAMRecord> extractSupportingReads(
	final VariantContext ctx,
	final String sample,
	final SamReader reader,
	QueryInterval intervals[]
	) {
intervals = QueryInterval.optimizeIntervals(intervals);
final List<SAMRecord> L = new ArrayList<>();
final CloseableIterator<SAMRecord> iter= reader.query(intervals,false);
while(iter.hasNext())
	{
	final SAMRecord rec=iter.next();
	if(rec.getReadUnmappedFlag()) continue;
	if(rec.getDuplicateReadFlag()) continue;
	if(rec.getCigar()==null || rec.getCigar().isEmpty()) continue;
	final SAMReadGroupRecord rg = rec.getReadGroup();
	if(rg==null) continue;
	if(!sample.equals(rg.getSample())) continue;
	L.add(rec);
	}
iter.close();
return L;
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private QueryInterval[] getRegionsAsQueryIntervalArray(final SamReader r) 
{
final List<QueryInterval> queryIntervals=new ArrayList<>();
for(final String intervalStr:this.intervalStrList)
	{
	if(StringUtil.isBlank(intervalStr)) {
		throw new IllegalStateException("empty interval string in "+this.intervalStrList);
	}
	final SAMSequenceDictionary dict =r.getFileHeader().getSequenceDictionary();
	if(dict==null) throw new JvarkitException.BamDictionaryMissing(r.getResourceDescription());
	final Interval i= new IntervalParser(dict).
		setContigNameIsWholeContig(true).
		parse(intervalStr);
	if(i==null) throw new IllegalArgumentException("Cannot parse interval "+intervalStr);
	final int referenceIndex = dict.getSequenceIndex(i.getContig());
	if(referenceIndex<0) throw new IllegalArgumentException("tid<0 ??! for "+i);
	queryIntervals.add(new QueryInterval(referenceIndex, i.getStart(), i.getEnd()));
	}
Collections.sort(queryIntervals);

return queryIntervals.toArray(QueryInterval.optimizeIntervals(new QueryInterval[queryIntervals.size()]));
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private static void sliceFromVCF(@NotNull final CommandLine cmd) throws IOException {
    final String inputPath = cmd.getOptionValue(INPUT);
    final String outputPath = cmd.getOptionValue(OUTPUT);
    final String vcfPath = cmd.getOptionValue(VCF);
    final int proximity = Integer.parseInt(cmd.getOptionValue(PROXIMITY, "500"));
    final SamReader reader = SamReaderFactory.makeDefault().open(new File(inputPath));
    final QueryInterval[] intervals = getIntervalsFromVCF(vcfPath, reader.getFileHeader(), proximity);
    final CloseableIterator<SAMRecord> iterator = getIterator(reader, intervals);
    writeToSlice(outputPath, reader.getFileHeader(), iterator);
    reader.close();
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
@NotNull
private static QueryInterval[] getIntervalsFromBED(@NotNull final String bedPath, @NotNull final SAMFileHeader header)
        throws IOException, EmptyFileException {
    final Slicer bedSlicer = SlicerFactory.fromBedFile(bedPath);
    final List<QueryInterval> queryIntervals = Lists.newArrayList();
    for (final GenomeRegion region : bedSlicer.regions()) {
        queryIntervals.add(new QueryInterval(header.getSequenceIndex(region.chromosome()), (int) region.start(), (int) region.end()));
    }
    return QueryInterval.optimizeIntervals(queryIntervals.toArray(new QueryInterval[queryIntervals.size()]));
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
@NotNull
private static CloseableIterator<SAMRecord> getIterator(@NotNull final SamReader reader, @NotNull final QueryInterval[] intervals,
        @NotNull final long[] filePointers) {
    if (reader instanceof SamReader.PrimitiveSamReaderToSamReaderAdapter) {
        final SamReader.PrimitiveSamReaderToSamReaderAdapter adapter = (SamReader.PrimitiveSamReaderToSamReaderAdapter) reader;
        if (adapter.underlyingReader() instanceof BAMFileReader) {
            final BAMFileReader bamReader = (BAMFileReader) adapter.underlyingReader();
            return bamReader.createIndexIterator(intervals, false, filePointers);
        }
    }
    return reader.queryOverlapping(intervals);
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private static void writeToSlice(final String path, final SamReader reader, final QueryInterval[] intervals) {
    final File outputBAM = new File(path);
    final SAMFileWriter writer = new SAMFileWriterFactory().makeBAMWriter(reader.getFileHeader(), true, outputBAM);
    final SAMRecordIterator iterator = reader.queryOverlapping(intervals);
    while (iterator.hasNext()) {
        writer.addAlignment(iterator.next());
    }
    iterator.close();
    writer.close();
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
@NotNull
private SAMRecordIterator queryBam(@NotNull final PotentialMNVRegion potentialMnvRegion) {
    final int referenceIndex = getReferenceIndex(potentialMnvRegion);
    final QueryInterval[] queryIntervals =
            new QueryInterval[] { new QueryInterval(referenceIndex, potentialMnvRegion.start(), potentialMnvRegion.end() - 1) };
    return tumorReader().queryOverlapping(queryIntervals);
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private SortedMap<QueryInterval, List<Allele>> getQueryIntervalsMap(final File vcf) {

        final Map<String, Integer> contigIndexMap = new HashMap<>();
        final VCFFileReader vcfReader = new VCFFileReader(vcf, false);

        // We want to look at unfiltered SNP sites for which the sample is genotyped as a het
        // with high quality.
        final CompoundFilter compoundFilter = new CompoundFilter(true);
        compoundFilter.add(new SnpFilter());
        compoundFilter.add(new PassingVariantFilter());
        compoundFilter.add(new GenotypeQualityFilter(MINIMUM_GQ, SAMPLE));
        compoundFilter.add(new HeterozygosityFilter(true, SAMPLE));

        final Iterator<VariantContext> hetIterator = new FilteringVariantContextIterator(vcfReader.iterator(), compoundFilter);

        for (final VCFContigHeaderLine vcfContig : vcfReader.getFileHeader().getContigLines()) {
            contigIndexMap.put(vcfContig.getID(), vcfContig.getContigIndex());
        }

        // return a TreeMap since the keys are comparable, and this will use their order in the iteration
        final SortedMap<QueryInterval, List<Allele>> map = new TreeMap<>();

        while (hetIterator.hasNext()) {
            final VariantContext vc = hetIterator.next();
            map.put(new QueryInterval(contigIndexMap.get(vc.getContig()), vc.getStart(), vc.getEnd()), vc.getGenotype(SAMPLE).getAlleles());
        }

        vcfReader.close();

        return map;
    }
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
/**
 * Converts an interval in SimpleInterval format into an htsjdk QueryInterval.
 *
 * In doing so, a header lookup is performed to convert from contig name to index
 *
 * @param interval interval to convert
 * @param sequenceDictionary sequence dictionary used to perform the conversion
 * @return an equivalent interval in QueryInterval format
 */
public static QueryInterval convertSimpleIntervalToQueryInterval( final SimpleInterval interval, final SAMSequenceDictionary sequenceDictionary ) {
    Utils.nonNull(interval);
    Utils.nonNull(sequenceDictionary);

    final int contigIndex = sequenceDictionary.getSequenceIndex(interval.getContig());
    if ( contigIndex == -1 ) {
        throw new UserException("Contig " + interval.getContig() + " not present in reads sequence dictionary");
    }

    return new QueryInterval(contigIndex, interval.getStart(), interval.getEnd());
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private static List<QueryInterval> getAllChunksBalanced(SamReader bam, int countPerContig) {
    List<QueryInterval> ret = new ArrayList<>();
    SAMFileHeader header = bam.getFileHeader();
    for (SAMSequenceRecord s : header.getSequenceDictionary().getSequences()) {
        ret.addAll(getChunksBalanced(bam, s.getSequenceIndex(), countPerContig));
    }
    return ret;
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
private static List<QueryInterval> getChunksBalanced(SamReader bam, int sequenceIndex, int retCount) {
    List<QueryInterval> ret = new ArrayList<>();
    BAMIndex index = bam.indexing().getIndex();
    SAMFileHeader header = bam.getFileHeader();
    SAMSequenceRecord s = header.getSequence(sequenceIndex);
    long totalLength = chunksLength(getChunks(index, sequenceIndex, 1, s.getSequenceLength() + 1));
    if (totalLength == 0) {
        return ret;
    }
    int sofar = 0;
    long targetLength = totalLength / retCount;
    int end = s.getSequenceLength();
    int step = s.getSequenceLength() / (100 * retCount);
    if (step < 1) step = 1;
    int start = 1;
    for (int j = step; j < end; j += step) {
        if (j > end) j = end;
        List<Chunk> candidate = getChunks(index, sequenceIndex, start, j);
        long size = chunksLength(candidate);
        if (size < targetLength) {
            // not big enough yet
            continue;
        }
        if (size > targetLength * 2) {
            // too large, search for a good separation point
            // TODO
        }
        // good, emit.
        ret.add(new QueryInterval(sequenceIndex, start, j + 1));
        start = j;
        sofar += size;
        if (ret.size() < retCount) {
            targetLength = (totalLength - sofar) / (retCount - ret.size());
        } else {
            targetLength = totalLength / retCount;
        }

    }
    return ret;
}
 

import htsjdk.samtools.QueryInterval; //导入依赖的package包/类
@Test
public void testConvertSimpleIntervalToQueryInterval() {
    final SAMSequenceRecord contigRecord = new SAMSequenceRecord("1", 100);
    final SAMSequenceDictionary dictionary = new SAMSequenceDictionary(Arrays.asList(contigRecord));
    final SimpleInterval originalInterval = new SimpleInterval("1", 5, 10);

    final QueryInterval convertedInterval = IntervalUtils.convertSimpleIntervalToQueryInterval(originalInterval, dictionary);
    Assert.assertEquals(convertedInterval.referenceIndex, 0);
    Assert.assertEquals(convertedInterval.start, 5);
    Assert.assertEquals(convertedInterval.end, 10);
}