Java Code Examples for htsjdk.samtools.SamReader#queryOverlapping()

public void slice(@NotNull final SamReader samReader, final List<GenomeRegion> regions, @NotNull final Consumer<SAMRecord> consumer)
{
    mConsumerHalt = false;

    final QueryInterval[] queryIntervals = createIntervals(regions, samReader.getFileHeader());

    try (final SAMRecordIterator iterator = samReader.queryOverlapping(queryIntervals))
    {
        while (!mConsumerHalt && iterator.hasNext())
        {
            final SAMRecord record = iterator.next();

            if (passesFilters(record))
            {
                consumer.accept(record);
            }
        }
    }
}
 

public void slice(@NotNull final SamReader samReader, final QueryInterval[] queryIntervals, @NotNull final Consumer<SAMRecord> consumer)
{
    mConsumerHalt = false;

    try (final SAMRecordIterator iterator = samReader.queryOverlapping(queryIntervals))
    {
        while (!mConsumerHalt && iterator.hasNext())
        {
            final SAMRecord record = iterator.next();

            if (passesFilters(record))
            {
                consumer.accept(record);
            }
        }
    }
}
 

public List<SAMRecord> slice(@NotNull final SamReader samReader, final QueryInterval[] queryIntervals)
{
    final List<SAMRecord> records = Lists.newArrayList();

    try (final SAMRecordIterator iterator = samReader.queryOverlapping(queryIntervals))
    {
        while (iterator.hasNext())
        {
            final SAMRecord record = iterator.next();

            if (passesFilters(record))
            {
                records.add(record);
            }
        }
    }

    return records;
}
 

private static void sliceFromVCF(@NotNull CommandLine cmd) throws IOException {
    String inputPath = cmd.getOptionValue(INPUT);
    String vcfPath = cmd.getOptionValue(VCF);
    int proximity = Integer.parseInt(cmd.getOptionValue(PROXIMITY, "500"));

    SamReaderFactory readerFactory = createFromCommandLine(cmd);
    SamReader reader = readerFactory.open(new File(inputPath));

    QueryInterval[] intervals = getIntervalsFromVCF(vcfPath, reader.getFileHeader(), proximity);
    CloseableIterator<SAMRecord> iterator = reader.queryOverlapping(intervals);
    SAMFileWriter writer = new SAMFileWriterFactory().setCreateIndex(true)
            .makeBAMWriter(reader.getFileHeader(), true, new File(cmd.getOptionValue(OUTPUT)));

    writeToSlice(writer, iterator);

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

@NotNull
private static CloseableIterator<SAMRecord> getIterator(@NotNull SamReader reader, @NotNull QueryInterval[] intervals,
        @NotNull final long[] filePointers) {
    if (reader instanceof SamReader.PrimitiveSamReaderToSamReaderAdapter) {
        SamReader.PrimitiveSamReaderToSamReaderAdapter adapter = (SamReader.PrimitiveSamReaderToSamReaderAdapter) reader;
        if (adapter.underlyingReader() instanceof BAMFileReader) {
            BAMFileReader bamReader = (BAMFileReader) adapter.underlyingReader();
            return bamReader.createIndexIterator(intervals, false, filePointers);
        }
    }
    return reader.queryOverlapping(intervals);
}
 

public void slice(@NotNull final SamReader samReader, @NotNull final Consumer<SAMRecord> consumer) {
    final Set<String> processed = Sets.newHashSet();
    final QueryInterval[] queryIntervals = createIntervals(regions, samReader.getFileHeader());

    try (final SAMRecordIterator iterator = samReader.queryOverlapping(queryIntervals)) {
        while (iterator.hasNext()) {
            final SAMRecord record = iterator.next();
            if (samRecordMeetsQualityRequirements(record)) {
                if (processed.add(record.toString())) {
                    consumer.accept(record);
                }
            }
        }
    }
}
 

public void slice(@NotNull final SamReader samReader, @NotNull final Consumer<SAMRecord> consumer) {
    final QueryInterval[] queryIntervals = createIntervals(regions, samReader.getFileHeader());

    try (final SAMRecordIterator iterator = samReader.queryOverlapping(queryIntervals)) {
        while (iterator.hasNext()) {
            final SAMRecord record = iterator.next();
            if (samRecordMeetsQualityRequirements(record)) {
                consumer.accept(record);
            }
        }
    }
}
 

public ReadLocusIterator(SamReader samReader, Feature region, int maxDepth) {
       
	this.maxDepth = maxDepth;
	
	if (region != null) {
		samIter = new ForwardShiftInsertIterator(samReader.queryOverlapping(region.getSeqname(), (int) region.getStart(), (int) region.getEnd()));
	} else {
		samIter = new ForwardShiftInsertIterator(samReader.iterator());
	}
}
 

public BamOverlapChecker(SamReader bam_file){
    
    SAMFileHeader  header = bam_file.getFileHeader();
    SAMSequenceDictionary dict = header.getSequenceDictionary();
    List<SAMSequenceRecord> sequences = dict.getSequences();
   
    
    booleanMap = new HashMap<String, boolean[]>();
    
    for(SAMSequenceRecord sequence : sequences){
        int sequenceEnd = sequence.getSequenceLength();
        int arrayLength = (int) Math.ceil( (float) sequenceEnd / (float)stepSize );
        boolean[] tempArray;
        tempArray = new boolean[arrayLength];
        
        for(int i=0;i<arrayLength;i++){
            SAMRecordIterator bamQuery = bam_file.queryOverlapping(sequence.getSequenceName(), i*stepSize, (i+1)*stepSize);
            if(bamQuery.hasNext()){
                tempArray[i] = true;
            }else{
                tempArray[i] = false;
            }
            bamQuery.close();
        }
        booleanMap.put(sequence.getSequenceName(),tempArray );
        //System.out.println("Finished checking the bam for chromosome " + sequence.getSequenceName());
    }
    
}
 

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);
	}
 

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.");
}
 

/**
 * 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;
}
 

public static String get_allele_specific_overlap_at_snp(GeneticVariant this_variant,
                                                        int sample_index, 
                                                        String chromosome,
                                                        String position, 
                                                        SamReader bam_file){
    
    int pos_int = Integer.parseInt(position);
    
    Alleles all_variants = this_variant.getVariantAlleles();
    Character ref_allele_char = all_variants.getAllelesAsChars()[0];
    String ref_allele = ref_allele_char.toString(); 
    //System.out.println("ref_allele: " + ref_allele);
    Character alt_allele_char = all_variants.getAllelesAsChars()[1];
    String alt_allele = alt_allele_char.toString(); 
    //System.out.println("alt_allele: " + alt_allele);
   
    
    
    int ref_overlap = 0;
    int alt_overlap = 0;
    int no__overlap = 0;
    
    

   
    // now determine per individual the sample variants.
    // I'm assuming the ordering is the same as the individual names created 
    // by the  getSampleNames() method. 
    // Otherwise the data will be nicely permuted, and I will have to convert some stuff.
    
    int position_of_snp = Integer.parseInt(position);                       
    
    //Check to make sure the variant position is not 0.
    if(position_of_snp <= 0){
        System.out.println("A SNP was read with a position lower than 1. This is illegal");
        System.out.println("Please adapt your genotype files by removing SNPs with these illegal positions");
        System.out.println("\tchr: " + chromosome + " pos: " + position);
        throw new IllegalDataException("Variant Position was less than 1");
    }
    
    SAMRecordIterator all_reads_in_region;
    try{
        all_reads_in_region = bam_file.queryOverlapping(chromosome, position_of_snp, position_of_snp);
    } catch(IllegalArgumentException e){
        System.out.println("Found an error when trying the following input:");
        System.out.println("chr:\t"+chromosome);
        System.out.println("pos:\t"+ position);
        System.out.println("If these values look correct, please make sure your bam file is sorted AND indexed by samtools.");
        System.out.println("If the problem persists, perhaps the chromosome (or sequence) are not the same in the genotype or bam file");
        all_reads_in_region = null;
        System.exit(0);
        
    }
    
    String bases = "";
    
    while(all_reads_in_region.hasNext()){
        
        SAMRecord read_in_region = all_reads_in_region.next();
        
        Character base_in_read = get_base_at_position(read_in_region, pos_int);
        if(GlobalVariables.verbosity >= 100){
            System.out.println("base_in_read: " + base_in_read);
        }
        
        if( base_in_read == ref_allele.charAt(0) ){
            ref_overlap++;            
        } else if(base_in_read == alt_allele.charAt(0) ){
            alt_overlap++;
        }else if(base_in_read == '!'  || base_in_read == 'N'){
            continue;
        }else{
            no__overlap++;
        }
    }
    
    //This line below cost me a day to figure out the error.
    all_reads_in_region.close();
    
    String string_for_output;
    string_for_output = Integer.toString(ref_overlap) + "\t" + 
                        Integer.toString(alt_overlap) + "\t" + 
                        Integer.toString(no__overlap);
    
    return(string_for_output);
}