Java SamReader.query方法代码示例

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包/类
private Handler<SAMRecord> filterReads(BamTrack<Read> bamTrack, BamQueryOption options, SamReader reader,
                                       String chromosomeName, boolean coverageOnly, BamTrackEmitter trackEmitter)
        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, trackEmitter); //maxReadCount

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

    return filter;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
/**
 * Given a SamReader, find all reads within a region (ref:start-end) that span a splice junction and tally the number of times
 * each junction is spanned. Optionally tally the average values of "tallyTagName" tags (eg. NM).
 * 
 * @param reader
 * @param ref
 * @param start
 * @param end
 * @param orient
 * @param minOverlap
 * @param tallyTagName
 * @param separateReadCounts
 * @return
 */
public static SortedMap<GenomeSpan, MappedReadCounter> countJunctions(SamReader reader, String ref, int start, int end, Orientation orient, int minOverlap, String tallyTagName) {
    SAMRecordIterator it = reader.query(ref, start, end, true);
    SortedMap<GenomeSpan, MappedReadCounter> counters = new TreeMap<GenomeSpan, MappedReadCounter>();
    
    while (it.hasNext()) {
        SAMRecord read = it.next();
        if (read.isSecondaryOrSupplementary() || read.getDuplicateReadFlag() || read.getNotPrimaryAlignmentFlag() || read.getReadUnmappedFlag() || read.getSupplementaryAlignmentFlag()) {
            // skip all secondary / duplicate / unmapped reads
            continue;
        }
       
        if (isJunctionSpanning(read)) {
            for (GenomeSpan junction: getJunctionsForRead(read, orient, minOverlap)) {
                if (!counters.containsKey(junction)) {
                    counters.put(junction, new MappedReadCounter(tallyTagName));
                }
                
                counters.get(junction).addRead(read);
            }
        }
    }
    it.close();
    return counters;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
void openFile() throws IOException {
  LOG.info("Processing shard " + shard);
  final SamReader reader = BAMIO.openBAM(storageClient, shard.file,
      options.getStringency());
  iterator = null;
  if (reader.hasIndex() && reader.indexing() != null) {
    if (filter == Filter.UNMAPPED_ONLY) {
      LOG.info("Processing unmapped");
      iterator = reader.queryUnmapped();
    } else if (shard.span != null) {
      LOG.info("Processing span for " + shard.contig);
      iterator = reader.indexing().iterator(shard.span);
    } else if (shard.contig.referenceName != null && !shard.contig.referenceName.isEmpty()) {
      LOG.info("Processing all bases for " + shard.contig);
      iterator = reader.query(shard.contig.referenceName, (int) shard.contig.start,
          (int) shard.contig.end, false);
    }
  }
  if (iterator == null) {
    LOG.info("Processing all reads");
    iterator = reader.iterator();
  }
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
/**
 * Count reads in interval using the given filters.
 * @param bam
 * @param gc
 * @param filters List of filters to apply
 * @return
 * @throws MalformedURLException 
 */
public static long countReadsInWindow(String bam, GenomicCoords gc, List<SamRecordFilter> filters) throws MalformedURLException {

	/*  ------------------------------------------------------ */
	/* This chunk prepares SamReader from local bam or URL bam */
	UrlValidator urlValidator = new UrlValidator();
	SamReaderFactory srf=SamReaderFactory.make();
	srf.validationStringency(ValidationStringency.SILENT);
	SamReader samReader;
	if(urlValidator.isValid(bam)){
		samReader = srf.open(SamInputResource.of(new URL(bam)).index(new URL(bam + ".bai")));
	} else {
		samReader= srf.open(new File(bam));
	}
	/*  ------------------------------------------------------ */
	
	long cnt= 0;
	
	Iterator<SAMRecord> sam= samReader.query(gc.getChrom(), gc.getFrom(), gc.getTo(), false);
	AggregateFilter aggregateFilter= new AggregateFilter(filters);
	while(sam.hasNext()){
		SAMRecord rec= sam.next();
		if( !aggregateFilter.filterOut(rec) ){
			cnt++;
		}
	}
	try {
		samReader.close();
	} catch (IOException e) {
		e.printStackTrace();
	}
	return cnt;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
/** 
 * Find reads that overlap a given genomic region
 * @param reader
 * @param pos
 * @param orient
 * @param readLength
 * @param minOverlap
 * @param allowGaps
 * @return
 */
public static List<SAMRecord> findOverlappingReads(SamReader reader, GenomeSpan pos, Orientation orient, int readLength, int minOverlap, boolean allowGaps) {
	List<SAMRecord> out = new ArrayList<SAMRecord>();

	SAMRecordIterator it = reader.query(pos.ref, pos.start - readLength + minOverlap, pos.start + readLength - minOverlap, true);
    while (it.hasNext()) {
        SAMRecord read = it.next();
        if (read.isSecondaryOrSupplementary() || read.getDuplicateReadFlag() || read.getNotPrimaryAlignmentFlag() || read.getReadUnmappedFlag() || read.getSupplementaryAlignmentFlag()) {
            // skip all secondary / duplicate / unmapped reads
            continue;
        }
        
        if (!allowGaps) {
            // skip all reads with gaps
            for (CigarElement el: read.getCigar().getCigarElements()) {
                if (el.getOperator() == CigarOperator.N) {
                    continue;
                }
            }
        }
        
        if (ReadUtils.getFragmentEffectiveStrand(read, orient) != pos.strand) {
        	continue;
        }
        
        for (GenomeSpan region: ReadUtils.getJunctionFlankingRegions(read, orient, minOverlap)) {
        	if (region.start <= (pos.start - minOverlap) && region.end >= (pos.start + minOverlap)) {
        		out.add(read);
        		break;
        	}
        }
    }
    it.close();
	return out;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public static HeaderInfo getHeaderFromBAMFile(Storage.Objects storage, String BAMPath, Iterable<Contig> explicitlyRequestedContigs) throws IOException {
  HeaderInfo result = null;

  // Open and read start of BAM
  LOG.info("Reading header from " + BAMPath);
  final SamReader samReader = BAMIO
      .openBAM(storage, BAMPath, ValidationStringency.DEFAULT_STRINGENCY);
  final SAMFileHeader header = samReader.getFileHeader();
  Contig firstContig = getFirstExplicitContigOrNull(header, explicitlyRequestedContigs);
  if (firstContig == null) {
    final SAMSequenceRecord seqRecord = header.getSequence(0);
    firstContig = new Contig(seqRecord.getSequenceName(), -1, -1);
  }

  LOG.info("Reading first chunk of reads from " + BAMPath);
  final SAMRecordIterator recordIterator = samReader.query(
      firstContig.referenceName, (int)firstContig.start + 1, (int)firstContig.end + 1, false);

  Contig firstShard = null;
  while (recordIterator.hasNext() && result == null) {
    SAMRecord record = recordIterator.next();
    final int alignmentStart = record.getAlignmentStart();
    if (firstShard == null && alignmentStart > firstContig.start &&
        (alignmentStart < firstContig.end || firstContig.end == -1)) {
      firstShard = new Contig(firstContig.referenceName, alignmentStart, alignmentStart);
      LOG.info("Determined first shard to be " + firstShard);
      result = new HeaderInfo(header, firstShard);
    }
  }
  recordIterator.close();
  samReader.close();

  if (result == null) {
    throw new IOException("Did not find reads for the first contig " + firstContig.toString());
  }
  LOG.info("Finished header reading from " + BAMPath);
  return result;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
@Test
public void openBAMTest() throws IOException {
 GCSOptions popts = PipelineOptionsFactory.create().as(GCSOptions.class);
 final Storage.Objects storageClient = Transport.newStorageClient(popts).build().objects();

 SamReader samReader = BAMIO.openBAM(storageClient, TEST_BAM_FNAME, ValidationStringency.DEFAULT_STRINGENCY);
 SAMRecordIterator iterator =  samReader.query("1", 550000, 560000, false);
 int readCount = 0;
 while (iterator.hasNext()) {
     iterator.next();
     readCount++;
 }
 Assert.assertEquals("Unexpected count of unmapped reads",
	  EXPECTED_UNMAPPED_READS_COUNT, readCount);
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
private Optional<Wig> checkAreaForReads(SamReader reader, String sequenceName, int start, int end) {
    try (SAMRecordIterator query = reader.query(sequenceName, start, end, false)) {
        return Optional.ofNullable(query.hasNext() ? new Wig(start, end, 1) : null);
    }
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public boolean importIntoHdfs(String weburl, FileSystem fileSystem,
		String path) throws IOException {

	// check if bai is available else download whole file

	final SamReader reader = SamReaderFactory.makeDefault().validationStringency(htsjdk.samtools.ValidationStringency.SILENT).open(
			SamInputResource.of(new URL(weburl)).index(
					new URL(weburl + ".bai")));

	// path in hdfs
	String[] tiles = weburl.split("/");
	String name = tiles[tiles.length - 1];

	String target = HdfsUtil.path(path, name);

	SAMFileHeader header = reader.getFileHeader();
	SAMSequenceDictionary seqDictionary = header.getSequenceDictionary();

	String referenceName = null;

	for (SAMSequenceRecord record : seqDictionary.getSequences()) {

		if (record.getSequenceLength() == REFERENCE_LENGTH) {
			referenceName = record.getSequenceName();
		}
	}

	if (referenceName == null) {
		reader.close();
		error = "No mitochondrial contig found in " + weburl + ".";
		return false;
	}

	FSDataOutputStream out = fileSystem.create(new Path(target));
	SAMFileWriter writer = new SAMFileWriterFactory().makeBAMWriter(
			reader.getFileHeader(), true, out);

	SAMRecordIterator reads = reader.query(referenceName, 0, 0, false);
	int good = 0;
	int bad = 0;
	int written = 0;
	SAMRecord read = null;
	while (reads.hasNext()) {
		try { // hansi style solution TODO!
			read = reads.next();
			good++;
		} catch (Exception e) {
			// e.printStackTrace(s);
			bad++;
		}
		writer.addAlignment(read);
		written++;
	}

	writer.close();
	reader.close();

	System.out.println("Bad reads: " + bad);
	System.out.println("Good reads: " + good);
	System.out.println("Written reads: " + written);
	
	return true;
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public List<Boolean> filterReads(SamReader samReader, String chrom, int from, int to) throws IOException{

		Iterator<SAMRecord> filterSam= samReader.query(chrom, from, to, false);
		
		AggregateFilter aggregateFilter= new AggregateFilter(this.getFeatureFilter().getSamRecordFilter());

		// This array will contain true/false to indicate whether a record passes the 
		// sam filters AND the awk filter (if given).
		// boolean[] results= new boolean[(int) this.nRecsInWindow];
		List<Boolean> results= new ArrayList<Boolean>();

		List<String> awkDataInput= new ArrayList<String>();
		while(filterSam.hasNext()){ 
			// Record whether a read passes the sam filters. If necessary, we also 
			// store the raw reads for awk.
			SAMRecord rec= filterSam.next();
			boolean passed;
			if(!rec.getReadUnmappedFlag() && 
			        !aggregateFilter.filterOut(rec) &&
			        rec.getAlignmentEnd() >= rec.getAlignmentStart()){
				passed= true;
			} else {
				passed= false;
			}

			// Filter for variant reads: Do it only if there is an intersection between variant interval and current genomic window
			if(this.getFeatureFilter().getVariantChrom().equals(Filter.DEFAULT_VARIANT_CHROM.getValue())){
				// Variant read filter is not set. Nothing to do.
			} else if(passed){ 
				// Memo: Test filter(s) only if a read is passed==true. 
				passed= this.isSNVRead(rec, this.getFeatureFilter().isVariantOnly());
			}
			
			String raw= null;

			if(passed && (! this.getFeatureFilter().getShowRegex().equals(Filter.DEFAULT_SHOW_REGEX.getValue()) || 
					      ! this.getFeatureFilter().getHideRegex().equals(Filter.DEFAULT_HIDE_REGEX.getValue()))){
				// grep
				raw= rec.getSAMString().trim();
				boolean showIt= true;
				if(! this.getFeatureFilter().getShowRegex().equals(Filter.DEFAULT_SHOW_REGEX.getValue())){
					showIt= this.getFeatureFilter().getShowRegex().matcher(raw).find();
				}
				boolean hideIt= false;
				if(! this.getFeatureFilter().getHideRegex().equals(Filter.DEFAULT_HIDE_REGEX.getValue())){
					hideIt= this.getFeatureFilter().getHideRegex().matcher(raw).find();	
				}
				if(!showIt || hideIt){
					passed= false;
				}
			}
			results.add(passed);
			if(passed && this.getAwk() != null && ! this.getAwk().isEmpty()){
				// We pass to awk only records that have been kept so far.
				if(raw == null){
					raw= rec.getSAMString().trim();
				}
				awkDataInput.add(raw);
			}
		} // End loop through reads

		// Apply the awk filter, if given
		if(this.getAwk() != null && ! this.getAwk().equals(Filter.DEFAULT_AWK.getValue())){
			String[] rawLines= new String[awkDataInput.size()];
			rawLines= awkDataInput.toArray(rawLines);
			boolean[] awkResults= Utils.passAwkFilter(rawLines, this.getAwk());
			// Compare the results array with awk filtered. Flip as appropriate the results array
			int awkIdx= 0;
			int i= 0;
			for(boolean isPassed : results){
				if(isPassed){
					if( ! awkResults[awkIdx]){
						results.set(i, false);
					}
					awkIdx++;
				}
				i++;
			}
		}
		return results;
	}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
public void update() throws InvalidGenomicCoordsException, IOException{

		if(this.getyMaxLines() == 0){
			return;
		}
		
		this.userWindowSize= this.getGc().getUserWindowSize();
		
		this.readStack= new ArrayList<List<SamSequenceFragment>>();
		if(this.getGc().getGenomicWindowSize() < this.MAX_REGION_SIZE){

			SamReader samReader= Utils.getSamReader(this.getWorkFilename());
			List<Boolean> passFilter= this.filterReads(samReader, this.getGc().getChrom(), this.getGc().getFrom(), this.getGc().getTo());

			this.nRecsInWindow= 0;
			for(boolean x : passFilter){ // The count of reads in window is the count of reads passing filters
				if(x){
					this.nRecsInWindow++;
				}
			}
			samReader= Utils.getSamReader(this.getWorkFilename());
			Iterator<SAMRecord> sam= samReader.query(this.getGc().getChrom(), this.getGc().getFrom(), this.getGc().getTo(), false);
			
			float max_reads= Float.parseFloat(Config.get(ConfigKey.max_reads_in_stack));
			float probSample= max_reads / this.nRecsInWindow;
			
			// Add this random String to the read name so different screenshot will generate 
			// different samples. 
			String rndOffset= Integer.toString(new Random().nextInt());

			List<TextRead> textReads= new ArrayList<TextRead>();
			ListIterator<Boolean> pass = passFilter.listIterator();
			while(sam.hasNext() && textReads.size() < max_reads){
				SAMRecord rec= sam.next();
				if( pass.next() ){
					String templ_name= Utils.templateNameFromSamReadName(rec.getReadName());
					long v= (templ_name + rndOffset).hashCode(); // Hashing.md5().hashBytes((templ_name + rndOffset).getBytes()).asLong();
					Random rand = new Random(v);
					if(rand.nextFloat() < probSample){ // Downsampler
						TextRead tr= new TextRead(rec, this.getGc());
						textReads.add(tr);
					}
				}
			}
			this.readStack= stackReads(textReads);
		} else {
			this.nRecsInWindow= -1;
		}
	}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
@SuppressWarnings({ "rawtypes", "unchecked" })
@Override
public void update() throws InvalidGenomicCoordsException, IOException{
	
	if(this.getyMaxLines() == 0){
		return;
	}
	if( ! this.useSamtools()){
		String chrom= this.getGc().getChrom();
		
		if(! this.loci.containsKey(chrom)){
			this.loci.put(chrom, new HashMap<Integer, Locus>());
		}
		if(! this.zeroDepthIntervals.containsKey(chrom)){
			this.zeroDepthIntervals.put(chrom, new IntervalTree());
		}

		// Check cache is not growing too much
		if(this.loci.get(chrom).keySet().size() > 500000){
			this.loci.get(chrom).clear();
			this.zeroDepthIntervals.get(chrom).clear();
		}
		
		// Find the positions that we haven't visited before:
		List<Integer> missingPos= new ArrayList<Integer>();
		for(int pos= this.getGc().getFrom(); pos <= this.getGc().getTo(); pos++){
			if( ! this.loci.get(chrom).containsKey(pos) &&  
				! this.zeroDepthIntervals.get(chrom).overlappers(pos, pos).hasNext()){
				missingPos.add(pos);
			}
		}
		for(List<Integer> gap : mergePositionsInIntervals(missingPos)){

			int qryFrom= gap.get(0);
			int qryTo= gap.get(1);
			
			SamReader samReader= Utils.getSamReader(this.getWorkFilename());
			List<Boolean> passFilter= this.filterReads(samReader, chrom, qryFrom, qryTo);
			samReader= Utils.getSamReader(this.getWorkFilename());
			
			Iterator<SAMRecord> sam= samReader.query(chrom, qryFrom, qryTo, false);
	
			ListIterator<Boolean> pass = passFilter.listIterator();
			while(sam.hasNext()){
				SAMRecord rec= sam.next();
				if(pass.next()){
					this.add(rec, qryFrom, qryTo, this.loci.get(chrom));
				}
			}
			
			// Now add the loci that have been collected in this last update
			List<Integer> zeroDepthPos= new ArrayList<Integer>();
			for(int pos= qryFrom; pos <= qryTo; pos++){
				if( ! this.loci.get(chrom).containsKey(pos)){
					zeroDepthPos.add(pos);
				}
			}
			List<List<Integer>> zeroDepthGaps = mergePositionsInIntervals(zeroDepthPos);
			for(List<Integer> z : zeroDepthGaps){
				this.zeroDepthIntervals.get(chrom).put(z.get(0), z.get(1), null);
			}
		}
	}
	List<Float> screenScores= this.prepareScreenScores();
	this.setScreenScores(screenScores);
}
 

import htsjdk.samtools.SamReader; //导入方法依赖的package包/类
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsReadable(REFERENCE_SEQUENCE);
    IOUtil.assertFileIsWritable(OUTPUT);

    final SamReader in = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).open(INPUT);

    final ReferenceSequenceFile reference = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE);
    final SAMSequenceDictionary refDict = reference.getSequenceDictionary();

    if (refDict == null) {
        log.error("No reference sequence dictionary found. Aborting.  You can create a sequence dictionary for the reference fasta using CreateSequenceDictionary.jar.");
        CloserUtil.close(in);
        return 1;
    }

    printDictionary("SAM/BAM file", in.getFileHeader().getSequenceDictionary());
    printDictionary("Reference", refDict);
    final Map<Integer, Integer> newOrder = buildSequenceDictionaryMap(refDict, in.getFileHeader().getSequenceDictionary());

    // has to be after we create the newOrder
    final SAMFileHeader outHeader = in.getFileHeader().clone();
    outHeader.setSequenceDictionary(refDict);

    log.info("Writing reads...");
    if (in.hasIndex()) {
        try( final SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(outHeader, true, OUTPUT)) {

            // write the reads in contig order
            for (final SAMSequenceRecord contig : refDict.getSequences()) {
                final SAMRecordIterator it = in.query(contig.getSequenceName(), 0, 0, false);
                writeReads(out, it, newOrder, contig.getSequenceName());
            }
            // don't forget the unmapped reads
            writeReads(out, in.queryUnmapped(), newOrder, "unmapped");
        }
    } else {
        try (final SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(outHeader, false, OUTPUT)) {
            writeReads(out, in.iterator(), newOrder, "All reads");
        }
    }

    // cleanup
    CloserUtil.close(in);
    return 0;
}