Java SamReader.queryOverlapping方法代码示例

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
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.SamReader; //导入方法依赖的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.SamReader; //导入方法依赖的package包/类
public ReadLocusIterator(SamReader samReader, Feature region) {
       	  
	if (region != null) {
		samIter = new ForwardShiftInsertIterator(samReader.queryOverlapping(region.getSeqname(), (int) region.getStart(), (int) region.getEnd()));
	} else {
		samIter = new ForwardShiftInsertIterator(samReader.iterator());
	}
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
private Counts getCounts(SamReader reader, LocusInfo locus) {
		
		int depth = 0;
		int altCount = 0;
		int refCount = 0;
		
		CloseableIterator<SAMRecord> iter = reader.queryOverlapping(locus.chromosome, locus.posStart, locus.posStop);
		while (iter.hasNext()) {
			SAMRecord read = iter.next();
			if (!read.getDuplicateReadFlag()) {
				depth += 1;
				
				Object[] baseAndQual = getBaseAndQualAtPosition(read, locus.posStart);
				Character base = (Character) baseAndQual[0];
				int baseQual = (Integer) baseAndQual[1];
				Character refBase = c2r.getSequence(locus.chromosome, locus.posStart, 1).charAt(0);
				
				// Override input with actual reference
//				locus.ref = new String(new char[] { refBase });
				locus.actualRef = new String(new char[] { refBase });
				
				if (locus.isIndel()) {
					if (hasIndel(read, locus)) {
						altCount += 1;
					} else if (!base.equals('N') && base.equals(refBase)) {
						refCount += 1;
					}
				} else if (baseQual >= minBaseQual) {
					
					if (!base.equals('N') && base.equals(refBase)) {
						refCount += 1;
					} else if (base.equals(locus.alt.charAt(0))) {
						altCount += 1;
					}
				}
			}
		}
		
		iter.close();
		
		return new Counts(refCount, altCount, depth);
	}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public void runDiff() throws Exception {
  SamReaderFactory readerFactory =  SamReaderFactory
      .makeDefault()
      .validationStringency(ValidationStringency.SILENT)
      .enable(SamReaderFactory.Option.CACHE_FILE_BASED_INDEXES);
  LOG.info("Opening file 1 for diff: " + BAMFile1);
  SamReader reader1 = readerFactory.open(new File(BAMFile1));
  LOG.info("Opening file 2 for diff: " + BAMFile2);
  SamReader reader2 = readerFactory.open(new File(BAMFile2));


  try {
    Iterator<Contig> contigsToProcess = null;
    if (options.contigsToProcess != null && !options.contigsToProcess.isEmpty()) {
      Iterable<Contig> parsedContigs = Contig.parseContigsFromCommandLine(options.contigsToProcess);
      referencesToProcess = Sets.newHashSet();
      for (Contig c : parsedContigs) {
        referencesToProcess.add(c.referenceName);
      }
      contigsToProcess = parsedContigs.iterator();
      if (!contigsToProcess.hasNext()) {
        return;
      }
    }
    LOG.info("Comparing headers");
    if (!compareHeaders(reader1.getFileHeader(), reader2.getFileHeader())) {
      error("Headers are not equal");
      return;
    }
    LOG.info("Headers are equal");
    do {
      SAMRecordIterator it1;
      SAMRecordIterator it2;
      if (contigsToProcess == null) {
        LOG.info("Checking all the reads");
        it1 = reader1.iterator();
        it2 = reader2.iterator();
      } else {
        Contig contig = contigsToProcess.next();
        LOG.info("Checking contig " + contig.toString());
        processedContigs++;
        it1 = reader1.queryOverlapping(contig.referenceName, (int)contig.start,  (int)contig.end);
        it2 = reader2.queryOverlapping(contig.referenceName, (int)contig.start,  (int)contig.end);
      }

      if (!compareRecords(it1, it2)) {
        break;
      }

      it1.close();
      it2.close();
    } while (contigsToProcess != null && contigsToProcess.hasNext());
  } catch (Exception ex) {
    throw ex;
  } finally {
    reader1.close();
    reader2.close();
  }
  LOG.info("Processed " + processedContigs + " contigs, " +
      processedLoci + " loci, " + processedReads + " reads.");
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
/**
 * To compare how sharded reading works vs. plain HTSJDK sequential iteration,
 * this method implements such iteration.
 * This makes it easier to discover errors such as reads that are somehow
 * skipped by a sharded approach.
 */
public static Iterable<Read> readSequentiallyForTesting(Objects storageClient,
    String storagePath, Contig contig,
    ReaderOptions options) throws IOException {
  Stopwatch timer = Stopwatch.createStarted();
  SamReader samReader = BAMIO.openBAM(storageClient, storagePath, options.getStringency());
  SAMRecordIterator iterator =  samReader.queryOverlapping(contig.referenceName,
      (int) contig.start + 1,
      (int) contig.end);
  List<Read> reads = new ArrayList<Read>();

  int recordsBeforeStart = 0;
  int recordsAfterEnd = 0;
  int mismatchedSequence = 0;
  int recordsProcessed = 0;
  Filter filter = setupFilter(options, contig.referenceName);
  while (iterator.hasNext()) {
    SAMRecord record = iterator.next();
    final boolean passesFilter = passesFilter(record, filter, contig.referenceName);

    if (!passesFilter) {
      mismatchedSequence++;
      continue;
    }
    if (record.getAlignmentStart() < contig.start) {
      recordsBeforeStart++;
      continue;
    }
    if (record.getAlignmentStart() > contig.end) {
      recordsAfterEnd++;
      continue;
    }
    reads.add(ReadUtils.makeReadGrpc(record));
    recordsProcessed++;
  }
  timer.stop();
  LOG.info("NON SHARDED: Processed " + recordsProcessed +
      " in " + timer +
      ". Speed: " + (recordsProcessed*1000)/timer.elapsed(TimeUnit.MILLISECONDS) + " reads/sec"
      + ", skipped other sequences " + mismatchedSequence
      + ", skippedBefore " + recordsBeforeStart
      + ", skipped after " + recordsAfterEnd);
  return reads;
}