Java SamReader类代码示例

import htsjdk.samtools.SamReader; //导入依赖的package包/类
@Nullable
private Read getReadFromBamFile(ReadQuery query, Chromosome chromosome, BamFile bamFile) throws IOException {
    try (SamReader reader = bamHelper.makeSamReader(bamFile, Collections.singletonList(chromosome),
                                                    chromosome.getReferenceId())) {
        String chromosomeName = chromosome.getName();
        if (reader.getFileHeader().getSequence(chromosomeName) == null) {
            chromosomeName = Utils.changeChromosomeName(chromosomeName);
        }

        SAMRecordIterator iterator = reader.query(chromosomeName, query.getStartIndex(), query.getEndIndex(),
                                                  true);
        while (iterator.hasNext()) {
            final SAMRecord samRecord = iterator.next();
            if (samRecord.getReadName().equals(query.getName())) {
                BamQueryOption option = new BamQueryOption();
                option.setRefID(chromosome.getReferenceId());
                option.setChromosomeName(chromosome.getName());
                SAMRecordHandler recordHandler = new SAMRecordHandler(query.getStartIndex(), query.getEndIndex(),
                                                                      referenceManager, null, option);
                List<BasePosition> diffBase = recordHandler.computeDifferentBase(samRecord);
                return BamUtil.createExtendedRead(samRecord, diffBase);
            }
        }
    }
    return null;
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
private void getReads(final BamFile bamFile, BamTrack<Read> bamTrack, final BamQueryOption options) throws
        IOException {
    Chromosome chromosome = bamTrack.getChromosome();
    try (SamReader reader = makeSamReader(bamFile, Collections.singletonList(chromosome),
            chromosome.getReferenceId())) {
        LOG.debug(getMessage(MessagesConstants.DEBUG_FILE_OPENING, bamFile.getPath()));
        String chromosomeName = options.getChromosomeName();

        if (reader.getFileHeader().getSequence(chromosomeName) == null) {
            chromosomeName = Utils.changeChromosomeName(chromosomeName);
        }

        Handler<SAMRecord> filter = filterReads(bamTrack, options, reader, chromosomeName, options.getMode() ==
                BamTrackMode.COVERAGE);
        //maxReadsCount

        bamTrack.setBlocks(filter.getSifter().getReadListResult());
        bamTrack.setMinPosition(filter.getMinPosition());
        bamTrack.setReferenceBuffer(filter.getReferenceBuff());
        bamTrack.setDownsampleCoverage(filter.getSifter().getDownsampleCoverageResult());
        bamTrack.setBaseCoverage(filter.getBaseCoverage(bamTrack.getScaleFactor()));
        bamTrack.setSpliceJunctions(filter.getSpliceJunctions());
    }
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
private Handler<SAMRecord> filterReads(BamTrack<Read> bamTrack, BamQueryOption options, SamReader reader,
                                       String chromosomeName, boolean coverageOnly) throws IOException {
                                                                                        //int maxReadCount
    CloseableIterator<SAMRecord> iterator = reader.query(chromosomeName, bamTrack.getStartIndex(),
                                                               bamTrack.getEndIndex(), false);
    LOG.debug(getMessage(MessagesConstants.DEBUG_GET_ITERATOR_QUERY, iterator.toString()));

    iterator = setIteratorFiltering(iterator, options);

    final Handler<SAMRecord> filter = BamUtil.createSAMRecordHandler(bamTrack, options, referenceManager,
            coverageOnly); //maxReadCount

    while (iterator.hasNext()) {
        final SAMRecord samRecord = iterator.next();
        //if read unmapped
        filter.add(samRecord);
    }

    return filter;
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
private void parseBam(final SamReader reader) {
    Assert.notNull(reader, getMessage(RESOURCE_NOT_FOUND));
    final SAMFileHeader samFileHeader = reader.getFileHeader();

    //check we can read this Bam-file
    //get list of chromosome
    final List<SAMSequenceRecord> list = samFileHeader.getSequenceDictionary().getSequences();
    Assert.notEmpty(list, getMessage(MessagesConstants.WRONG_HEADER_BAM_FILE_EMPTY_FILE));

    //get first chromosome and make a request to the file with this chromosome
    final SAMSequenceRecord samSequenceRecord = list.get(0);
    SAMRecordIterator iterator = reader.query(samSequenceRecord.getSequenceName(),
            Constants.BAM_START_INDEX_TEST, Math.min(Constants.MAX_BAM_END_INDEX_TEST,
                    samSequenceRecord.getSequenceLength()), false);
    Assert.notNull(iterator);
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
/**
 * A method that creates a SamReader object that's passed on to the HTSJDK.
 * Each time someone tries to open a SamReader on an URL,
 * HTSJDK checks if there's a custom reader factory and if it's there, this method is called.
 *
 * @param url target file URL
 * @return A SamReader object on a specified file URL
 */
@Override
public SamReader open(URL url) {
    PerformanceMonitor.start();
    AmazonS3URI amazonURI = new AmazonS3URI(url.toString());
    S3Client client = new S3Client();
    S3InputStreamFactory streamFactory = new S3InputStreamFactory(client);

    //download index file if is possible, and then start download .bam file
    final Optional<SeekableStream> indexStream;
    try {
        IndexLoader loader = new IndexLoader(client);
        indexStream = loader.loadIndex(amazonURI);
    } catch (IOException e) {
        throw new RuntimeIOException(e.getMessage() + " failed to download index", e);
    }

    SeekableStream stream = new S3SeekableStream(amazonURI, client, streamFactory);
    SamReaderFactory factory = SamReaderFactory.makeDefault();
    SamInputResource inputResource = SamInputResource.of(stream);

    indexStream.ifPresent(inputResource::index);

    return factory.open(inputResource);
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
public static void main(String args[]) throws IOException, ParseException {
      Options options = new Options();
      options.addOption("u","uniquehits",false,"only output hits with a single mapping");
      options.addOption("s","nosuboptimal",false,"do not include hits whose score is not equal to the best score for the read");
      options.addOption("p","pairedend",false,"output paired-end hits");
      options.addOption("j","junctions",false,"output junction mapping reads (reads with gaps)");
      CommandLineParser parser = new GnuParser();
      CommandLine cl = parser.parse( options, args, false );            
  	uniqueOnly = cl.hasOption("uniquehits");
  	filterSubOpt = cl.hasOption("nosuboptimal");
  	inclPairedEnd = cl.hasOption("pairedend");
  	inclJunction = cl.hasOption("junctions");
  	SamReaderFactory factory =
          SamReaderFactory.makeDefault()
              .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS, SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS)
              .validationStringency(ValidationStringency.SILENT);
SamReader reader = factory.open(SamInputResource.of(System.in));
      CloseableIterator<SAMRecord> iter = reader.iterator();
      while (iter.hasNext()) {
          SAMRecord record = iter.next();
          if (record.getReadUnmappedFlag()) {continue; }
          processRecord(record);
      }
      iter.close();
      reader.close();
  }
 

import htsjdk.samtools.SamReader; //导入依赖的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.SamReader; //导入依赖的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.SamReader; //导入依赖的package包/类
@NotNull
public Multimap<Chromosome, CobaltCount> fromBam(@NotNull final String referenceBam, @NotNull final String tumorBam)
        throws IOException, ExecutionException, InterruptedException {
    final File tumorFile = new File(tumorBam);
    final File referenceFile = new File(referenceBam);

    final SamReaderFactory readerFactory = SamReaderFactory.make();
    final String chromosomeLengthFileName = ChromosomeLengthFile.generateFilename(outputDirectory, tumor);
    final List<ChromosomeLength> lengths;
    try (SamReader reader = readerFactory.open(tumorFile)) {
        lengths = ChromosomeLengthFactory.create(reader.getFileHeader());
    }
    ChromosomeLengthFile.write(chromosomeLengthFileName, lengths);

    LOGGER.info("Calculating Read Count from {}", tumorFile.toString());
    final List<Future<ChromosomeReadCount>> tumorFutures = createFutures(readerFactory, tumorFile, lengths);

    LOGGER.info("Calculating Read Count from {}", referenceFile.toString());
    final List<Future<ChromosomeReadCount>> referenceFutures = createFutures(readerFactory, referenceFile, lengths);

    final Multimap<String, ReadCount> tumorCounts = fromFutures(tumorFutures);
    final Multimap<String, ReadCount> referenceCounts = fromFutures(referenceFutures);

    LOGGER.info("Read Count Complete");
    return CobaltCountFactory.merge(referenceCounts, tumorCounts);
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
private Map<String, Sequence> parseContigFromBam() {
	// TODO Auto-generated method stub
	
	final Map<String, Sequence> contig_list = new HashMap<String, Sequence>();
	final SamReader in1 = factory.open(new File(this.bam_list[0]));
	try {
		List<SAMSequenceRecord> seqs = 
				in1.getFileHeader().getSequenceDictionary().getSequences();
		for(SAMSequenceRecord seq : seqs) 
			contig_list.put( seq.getSequenceName(), 
					new Sequence(seq.getSequenceIndex(), seq.getSequenceLength()) );
		in1.close();
	} catch (IOException e) {
		// TODO Auto-generated catch block
		e.printStackTrace();
	}
	
	return contig_list;
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
public void countSplices(String input) {
	SamReader reader = SAMRecordUtils.getSamReader(input);
	
	for (SAMRecord read : reader) {
		if (read.getCigarString().contains("N")) {
			for (SpliceJunction junc : getJunctions(read)) {
				incrementCount(junc, read);
			}
		}
	}
	
	List<SpliceJunction> junctions = new ArrayList<SpliceJunction>(uniqueReads.keySet());
	Collections.sort(junctions, new SpliceJunctionComparator(reader.getFileHeader()));
	
	for (SpliceJunction junction : junctions) {
		String rec = String.format("%s\t%d\t%d\t.\t.\t.\t%d\t%d\t.", junction.chrom, junction.start, junction.stop,
				uniqueReads.get(junction), multiMapReads.get(junction));
		
		System.out.println(rec);
	}
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
public void call(String normal, String tumor, String vcf, String reference, int minBaseQual, double minMaf) throws IOException {
	loadLoci(vcf);
	c2r = new CompareToReference2();
	c2r.init(reference);
	this.minBaseQual = minBaseQual;
	this.minMaf = minMaf;
	
	System.err.println("Processing positions");

	SamReader normalReader = SAMRecordUtils.getSamReader(normal);
	
	SamReader tumorReader = SAMRecordUtils.getSamReader(tumor);
       
	for (LocusInfo locus : loci) {
		locus.normalCounts = getCounts(normalReader, locus);
		locus.tumorCounts = getCounts(tumorReader, locus);
	}
	
	normalReader.close();
	tumorReader.close();
	
	System.err.println("Writing results");
	
	outputResults();
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
public MultiSamReader(String[] inputBams, int minMapqForAssembly, boolean isPairedEnd, Feature region) {
	
	//TODO: Assert all SAM Headers have same sequence dict
	readers = new SamReader[inputBams.length];
	nextRecord = new SAMRecordWrapper[inputBams.length];
	iterators = new Iterator[inputBams.length];
	this.minMapqForAssembly = minMapqForAssembly;
	this.isPairedEnd = isPairedEnd;
	
	int idx = 0;
	for (String bamFileName : inputBams) {
		SamReader reader = SAMRecordUtils.getSamReader(bamFileName);
		
		readers[idx] = reader;
		
		// TODO: Pad by region size?
		iterators[idx] = readers[idx].queryOverlapping(region.getSeqname(), (int) region.getStart(), (int) region.getEnd());
		
		// cache next record
		cacheNextRecord(idx);
		
		idx += 1;
	}
	
	clientIterator = new MultiSamReaderIterator(this);
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
public static void main(String[] args) throws IOException {
	File bamFile = new File(
			"./data/miniCaviar_IDT_NEB.runA.NA12878.bwa.bam");
	File newBamFile = new File(
			"./data/miniCaviar_IDT_NEB.runA.NA12878.bwa.chrom1.bam");
	SAMFileWriter outputSam = null;
	final SamReader reader = SamReaderFactory.makeDefault()
			.validationStringency(ValidationStringency.SILENT)
			.open(bamFile);

	outputSam = new SAMFileWriterFactory().makeBAMWriter(
			reader.getFileHeader(), true, newBamFile);

	int currentReads = 0;
	for (final SAMRecord samRecord : reader) {
		if (samRecord.getReferenceIndex() == 0)
			writeBam(samRecord, outputSam);

	}

	reader.close();
	outputSam.close();
}
 

import htsjdk.samtools.SamReader; //导入依赖的package包/类
public static void main(String[] args) {
	String bam = "./data/HG00115.chrom11.ILLUMINA.bwa.GBR.exome.20130415.bam";
	SamReader reader = SamReaderFactory.makeDefault()
			.validationStringency(ValidationStringency.SILENT)
			.open(new File(bam));
	int length = 0;
	for (final SAMRecord record : reader) {
		if (record.getReadUnmappedFlag())
			continue;
		length += record.getReadLength();
	}

	ReferenceSequenceFile rsf = ReferenceSequenceFileFactory
			.getReferenceSequenceFile(new File("data/hs37d5.fa"));
	ReferenceSequence rs = rsf.getSequence("11");

	System.out.println(length * 1.0 / rs.length());
}