Java Code Examples for htsjdk.samtools.SAMRecord#getReferenceIndex()

/**
 * Determines if the read represented by a SAM record belongs to the top or bottom strand
 * or if it cannot determine strand position due to one of the reads being unmapped.
 * Top strand is defined as having a read 1 unclipped 5' coordinate
 * less than the read 2 unclipped 5' coordinate.  If a read is unmapped
 * we do not attempt to determine the strand to which the read or its mate belongs.
 * If the mate belongs to a different contig from the read, then the reference
 * index for the read and its mate is used in leu of the unclipped 5' coordinate.
 * @param rec Record to determine top or bottom strand
 * @return Top or bottom strand, unknown if it cannot be determined.
 */
static ReadStrand getStrand(final SAMRecord rec) {
    if (rec.getReadUnmappedFlag() || rec.getMateUnmappedFlag()) {
        return ReadStrand.UNKNOWN;
    }

    // If the read pair are aligned to different contigs we use
    // the reference index to determine relative 5' coordinate ordering.
    // Both the read and its mate should not have their unmapped flag set to true.
    if (!rec.getReferenceIndex().equals(rec.getMateReferenceIndex())) {
        if (rec.getFirstOfPairFlag() == (rec.getReferenceIndex() < rec.getMateReferenceIndex())) {
            return ReadStrand.TOP;
        } else {
            return ReadStrand.BOTTOM;
        }
    }

    final int read5PrimeStart = (rec.getReadNegativeStrandFlag()) ? rec.getUnclippedEnd() : rec.getUnclippedStart();
    final int mate5PrimeStart = (rec.getMateNegativeStrandFlag()) ? SAMUtils.getMateUnclippedEnd(rec) : SAMUtils.getMateUnclippedStart(rec);

    if (rec.getFirstOfPairFlag() == (read5PrimeStart <= mate5PrimeStart)) {
        return ReadStrand.TOP;
    } else {
        return ReadStrand.BOTTOM;
    }
}
 

private SAMRecordWrapper getNext(int idx) {
	SAMRecordWrapper record = null;
	if (iterators[idx].hasNext()) {
		SAMRecord read = iterators[idx].next();
		// If no genomic location is assigned, we've reached the unmapped read pairs.  Do not continue...
		// TODO: Need to include these in final bam files
		if (read.getReferenceIndex() >= 0) {
			record = new SAMRecordWrapper(read, isFiltered(read), shouldAssemble(read), idx);
		}
	}
	
	return record;
}
 

@Override
public boolean addRecord(SAMRecord r) throws IOException {
  // Don't attempt to reorder all the fully-unmapped (i.e. no reference sequence) records
  if (r.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
    if (!mRecordSet.isEmpty()) {
      flush();
    }
    flushRecord(r);
    return true;
  } else {
    return super.addRecord(r);
  }
}
 

/**
   * Compares SAM records on reference index, start position and if paired, mated vs unmated.
   * @param a first record
   * @param b second record
   * @return -1 if first record comes before second record, 1 if the converse. 0 for equal on compared values
   */
  public static int compareSamRecords(SAMRecord a, SAMRecord b) {
    final int thisRef = a.getReferenceIndex() & 0x7FFFFFFF; // makes -1 a bignum
    final int thatRef = b.getReferenceIndex() & 0x7FFFFFFF;
    final int rc = Integer.compare(thisRef, thatRef);
    if (rc != 0) {
      return rc;
    }

    final int ac = Integer.compare(a.getAlignmentStart(), b.getAlignmentStart());
    if (ac != 0) {
      return ac;
    }

    // Do this ... (this one doesn't have same ordering as original)
    //return Integer.compare(~a.getFlags(), ~b.getFlags());

    // or this ...
    // Following compares READ_PAIRED_FLAG, PORPER_PAIRED_FLAG, UNMAPPED_FLAG in that order.
    // The ^3 toggles values to get the ordering we require
    // The reverse makes sure comparison is in the order we require
    return Integer.compare(reverse3bits(a.getFlags() ^ 3), reverse3bits(b.getFlags() ^ 3));

    // or this ...

//    final int rpc = Boolean.compare(b.getReadPairedFlag(), a.getReadPairedFlag());
//    if (rpc != 0) {
//      return rpc;
//    }
//
//    if (a.getReadPairedFlag()) {
//      assert b.getReadPairedFlag();
//      final int mc = Boolean.compare(b.getProperPairFlag(), a.getProperPairFlag());
//      if (mc != 0) {
//        return mc;
//      }
//    }
//
//    return Boolean.compare(a.getReadUnmappedFlag(), b.getReadUnmappedFlag());
  }
 

private void recordPreviousAndAdvance(SAMRecord previousRecord, SAMRecord rec) throws IOException {
  if (previousRecord != null && previousRecord.getReferenceIndex() == mCurrentTemplate) {
    mPreviousAlignmentStart = previousRecord.getAlignmentStart() - 1; // to 0-based
  } else {
    mPreviousAlignmentStart = Integer.MIN_VALUE;
  }
  setBuffered(rec);
  advanceSubIterator();
}
 

/**
 * Calculates the mode for outward-facing pairs, using the first SAMPLE_FOR_MODE
 * outward-facing pairs found in INPUT
 */
private double getOutieMode() {

    int samplePerFile = SAMPLE_FOR_MODE / INPUT.size();

    Histogram<Integer> histo = new Histogram<Integer>();

    for (File f : INPUT) {
        SamReader reader = SamReaderFactory.makeDefault().open(f);
        int sampled = 0;
        for (Iterator<SAMRecord> it = reader.iterator(); it.hasNext() && sampled < samplePerFile; ) {
            SAMRecord sam = it.next();
            if (!sam.getFirstOfPairFlag()) {
                continue;
            }
            // If we get here we've hit the end of the aligned reads
            if (sam.getReadUnmappedFlag() && sam.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
                break;
            } else if (sam.getReadUnmappedFlag() || sam.getMateUnmappedFlag()) {
                continue;
            } else if ((sam.getAttribute(SAMTag.MQ.name()) == null ||
                    sam.getIntegerAttribute(SAMTag.MQ.name()) >= MINIMUM_MAPPING_QUALITY) &&
                    sam.getMappingQuality() >= MINIMUM_MAPPING_QUALITY &&
                    sam.getMateNegativeStrandFlag() != sam.getReadNegativeStrandFlag() &&
                    sam.getMateReferenceIndex().equals(sam.getReferenceIndex()) &&
                    SamPairUtil.getPairOrientation(sam) == PairOrientation.RF) {
                histo.increment(Math.abs(sam.getInferredInsertSize()));
                sampled++;
            }
        }
        CloserUtil.close(reader);
    }

    return histo.size() > 0 ? histo.getMode() : 0;
}
 

@Override
public void acceptRecord(final GcBiasCollectorArgs args) {
    final SAMRecord rec = args.getRec();
    if (logCounter < 100 && rec.getReadBases().length == 0) {
        log.warn("Omitting read " + rec.getReadName() + " with '*' in SEQ field.");
        if (++logCounter == 100) {
            log.warn("There are more than 100 reads with '*' in SEQ field in file.");
        }
        return;
    }
    if (!rec.getReadUnmappedFlag()) {
        if (referenceIndex != rec.getReferenceIndex() || gc == null) {
            final ReferenceSequence ref = args.getRef();
            refBases = ref.getBases();
            StringUtil.toUpperCase(refBases);
            final int refLength = refBases.length;
            final int lastWindowStart = refLength - scanWindowSize;
            gc = GcBiasUtils.calculateAllGcs(refBases, lastWindowStart, scanWindowSize);
            referenceIndex = rec.getReferenceIndex();
        }

        addReadToGcData(rec, this.gcData);
        if (ignoreDuplicates && !rec.getDuplicateReadFlag()) {
            addReadToGcData(rec, this.gcDataNonDups);
        }
    } else {
        updateTotalClusters(rec, this.gcData);
        if (ignoreDuplicates && !rec.getDuplicateReadFlag()) {
            updateTotalClusters(rec, this.gcDataNonDups);
        }
    }
}
 

private static void updateTracks(final List<SAMRecord> samRecords, final ReferenceTracks tracks) {
	for (final SAMRecord samRecord : samRecords) {
		if (samRecord.getAlignmentStart() != SAMRecord.NO_ALIGNMENT_START
				&& samRecord.getReferenceIndex() == tracks.getSequenceId()) {
			int refPos = samRecord.getAlignmentStart();
			int readPos = 0;
			for (final CigarElement cigarElement : samRecord.getCigar().getCigarElements()) {
				if (cigarElement.getOperator().consumesReferenceBases()) {
					for (int elementIndex = 0; elementIndex < cigarElement.getLength(); elementIndex++)
						tracks.addCoverage(refPos + elementIndex, 1);
				}
				switch (cigarElement.getOperator()) {
				case M:
				case X:
				case EQ:
					for (int pos = readPos; pos < cigarElement.getLength(); pos++) {
						final byte readBase = samRecord.getReadBases()[readPos + pos];
						final byte refBase = tracks.baseAt(refPos + pos);
						if (readBase != refBase)
							tracks.addMismatches(refPos + pos, 1);
					}
					break;

				default:
					break;
				}

				readPos += cigarElement.getOperator().consumesReadBases() ? cigarElement.getLength() : 0;
				refPos += cigarElement.getOperator().consumesReferenceBases() ? cigarElement.getLength() : 0;
			}
		}
	}
}
 

/** Note: this is the only getKey function that handles unmapped reads
 * specially!
 */
public static long getKey(final SAMRecord rec) {
	final int refIdx = rec.getReferenceIndex();
	final int start  = rec.getAlignmentStart();

	if (!(rec.getReadUnmappedFlag() || refIdx < 0 || start < 0))
		return getKey(refIdx, start);

	// Put unmapped reads at the end, but don't give them all the exact same
	// key so that they can be distributed to different reducers.
	//
	// A random number would probably be best, but to ensure that the same
	// record always gets the same key we use a fast hash instead.
	//
	// We avoid using hashCode(), because it's not guaranteed to have the
	// same value across different processes.

	int hash = 0;
	byte[] var;
	if ((var = rec.getVariableBinaryRepresentation()) != null) {
		// Undecoded BAM record: just hash its raw data.
		hash = (int)MurmurHash3.murmurhash3(var, hash);
	} else {
		// Decoded BAM record or any SAM record: hash a few representative
		// fields together.
		hash = (int)MurmurHash3.murmurhash3(rec.getReadName(), hash);
		hash = (int)MurmurHash3.murmurhash3(rec.getReadBases(), hash);
		hash = (int)MurmurHash3.murmurhash3(rec.getBaseQualities(), hash);
		hash = (int)MurmurHash3.murmurhash3(rec.getCigarString(), hash);
	}
	return getKey0(Integer.MAX_VALUE, hash);
}
 

public int writePairedSAMRecordToContext(SAMRecord sam, boolean useCompact) throws IOException, InterruptedException {
    int count = 0;
    int read1Ref = sam.getReferenceIndex();
    int read2Ref = sam.getMateReferenceIndex();
    if (!sam.getReadUnmappedFlag() && 
            (read1Ref == read2Ref || keepChrSplitPairs) && 
            (read1Ref >= 0 || read2Ref >= 0)) {
        context.getCounter(HalvadeCounters.OUT_BWA_READS).increment(1);
        writableRecord.set(sam);
        int beginpos = sam.getAlignmentStart();
        HashSet<Integer> keys = null;
        if (!mergeBam)
             keys = splitter.getRegions(sam, read1Ref, read2Ref);
        else {
            keys = new HashSet<>();
            keys.add(0);
        }
        for(Integer key : keys) {
            if(useCompact) {
                writeableCompactRegion.setRegion(key, beginpos);
                context.write(writeableCompactRegion, stub);
            } else {
                writableRegion.setChromosomeRegion(read1Ref, beginpos, key);
                context.write(writableRegion, writableRecord);
            }
            count++;
        }
    } else {
        if(sam.getReadUnmappedFlag()) 
            context.getCounter(HalvadeCounters.OUT_UNMAPPED_READS).increment(1);
        else
            context.getCounter(HalvadeCounters.OUT_DIFF_CHR_READS).increment(1);
    }
    return count;
}
 

public int writeSAMRecordToContext(SAMRecord sam, boolean useCompact) throws IOException, InterruptedException {
    int count = 0;
    if (!sam.getReadUnmappedFlag()){
        int read1Ref = sam.getReferenceIndex();
        context.getCounter(HalvadeCounters.OUT_BWA_READS).increment(1);
        writableRecord.set(sam);
        int beginpos = sam.getAlignmentStart();
        HashSet<Integer> keys = null;
        if (!mergeBam)
             keys = splitter.getRegions(sam, read1Ref);
        else {
            keys = new HashSet<>();
            keys.add(0);
        }
        for(Integer key : keys) {
            if(useCompact) {
                writeableCompactRegion.setRegion(key, beginpos);
                context.write(writeableCompactRegion, stub);
            } else {
                writableRegion.setChromosomeRegion(read1Ref, beginpos, key);
                context.write(writableRegion, writableRecord);
            }
            count++;
        }
    } else {
        context.getCounter(HalvadeCounters.OUT_UNMAPPED_READS).increment(1);
    }
    return count;
}
 

public static Map<Integer, AlignmentSpan> getSpans(List<SAMRecord> samRecords) {
	Map<Integer, AlignmentSpan> spans = new HashMap<Integer, AlignmentSpan>();
	int unmapped = 0;
	for (final SAMRecord r : samRecords) {

		int refId = r.getReferenceIndex();
		if (refId == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
			unmapped++;
			continue;
		}

		int start = r.getAlignmentStart();
		int end = r.getAlignmentEnd();

		if (spans.containsKey(refId)) {
			spans.get(refId).add(start, end - start, 1);
		} else {
			spans.put(refId, new AlignmentSpan(start, end - start));
		}
	}
	if (unmapped > 0) {
		AlignmentSpan span = new AlignmentSpan(AlignmentSpan.UNMAPPED_SPAN.getStart(),
				AlignmentSpan.UNMAPPED_SPAN.getSpan(), unmapped);
		spans.put(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX, span);
	}
	return spans;
}
 

private void populateNext() {
  mNextRecord = null;
  while (mIterator.hasNext()) {
    final SAMRecord rec = mIterator.next();
    final int refId = rec.getReferenceIndex();

    if (refId != mTemplate) { // Get current sequence ranges

      if (refId > mEndTemplate) {  // Past the last template in requested list so we can bail out
        //System.out.println("Aborting iterator for record at " + rec.getReferenceName() + ":" + rec.getAlignmentStart() + " as it is after last template");
        break;
      }

      mTemplate = refId;
      final RangeList<String> seqList = mRegions.get(mTemplate);
      mSeqRegions = (seqList == null || seqList.getRangeList().size() == 0) ? null : seqList.getRangeList();
      mCurrentRegionIdx = 0;
      mCurrentRegion = mSeqRegions == null ? null : mSeqRegions.get(0);
    }

    if (mSeqRegions == null) { // No regions on this template
      continue;
    }

    final int alignmentStart = rec.getAlignmentStart() - 1; // to 0-based
    int alignmentEnd = rec.getAlignmentEnd(); // to 0-based exclusive = noop
    if (alignmentEnd == 0) {
      // Use the read length to get an end point for unmapped reads
      alignmentEnd = rec.getAlignmentStart() + rec.getReadLength();
    }

    // Advance active region if need be
    while (mCurrentRegion.getEnd() <= alignmentStart && mCurrentRegionIdx < mSeqRegions.size() - 1) {
      ++mCurrentRegionIdx;
      mCurrentRegion = mSeqRegions.get(mCurrentRegionIdx);
    }

    if (alignmentEnd <= mCurrentRegion.getStart()) {  // before region
      continue;
    }

    if (alignmentStart >= mCurrentRegion.getEnd()) {  // past region
      if (mTemplate == mEndTemplate) {
        //System.out.println("Aborting iterator for record at " + rec.getReferenceName() + ":" + rec.getAlignmentStart() + " as it is after last region");
        break;
      } else {
        continue;
      }
    }

    mNextRecord = rec;
    break;
  }
}
 

/**
 * Copied from net.sf.picard.sam.SamPairUtil. This is a more permissive
 * version of the method, which does not reset alignment start and reference
 * for unmapped reads.
 * 
 * @param rec1
 * @param rec2
 * @param header
 */
public static void setLooseMateInfo(final SAMRecord rec1, final SAMRecord rec2, final SAMFileHeader header) {
	if (rec1.getReferenceName() != SAMRecord.NO_ALIGNMENT_REFERENCE_NAME
			&& rec1.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
		rec1.setReferenceIndex(header.getSequenceIndex(rec1.getReferenceName()));
	if (rec2.getReferenceName() != SAMRecord.NO_ALIGNMENT_REFERENCE_NAME
			&& rec2.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
		rec2.setReferenceIndex(header.getSequenceIndex(rec2.getReferenceName()));

	// If neither read is unmapped just set their mate info
	if (!rec1.getReadUnmappedFlag() && !rec2.getReadUnmappedFlag()) {

		rec1.setMateReferenceIndex(rec2.getReferenceIndex());
		rec1.setMateAlignmentStart(rec2.getAlignmentStart());
		rec1.setMateNegativeStrandFlag(rec2.getReadNegativeStrandFlag());
		rec1.setMateUnmappedFlag(false);
		rec1.setAttribute(SAMTag.MQ.name(), rec2.getMappingQuality());

		rec2.setMateReferenceIndex(rec1.getReferenceIndex());
		rec2.setMateAlignmentStart(rec1.getAlignmentStart());
		rec2.setMateNegativeStrandFlag(rec1.getReadNegativeStrandFlag());
		rec2.setMateUnmappedFlag(false);
		rec2.setAttribute(SAMTag.MQ.name(), rec1.getMappingQuality());
	}
	// Else if they're both unmapped set that straight
	else if (rec1.getReadUnmappedFlag() && rec2.getReadUnmappedFlag()) {
		rec1.setMateNegativeStrandFlag(rec2.getReadNegativeStrandFlag());
		rec1.setMateUnmappedFlag(true);
		rec1.setAttribute(SAMTag.MQ.name(), null);
		rec1.setInferredInsertSize(0);

		rec2.setMateNegativeStrandFlag(rec1.getReadNegativeStrandFlag());
		rec2.setMateUnmappedFlag(true);
		rec2.setAttribute(SAMTag.MQ.name(), null);
		rec2.setInferredInsertSize(0);
	}
	// And if only one is mapped copy it's coordinate information to the
	// mate
	else {
		final SAMRecord mapped = rec1.getReadUnmappedFlag() ? rec2 : rec1;
		final SAMRecord unmapped = rec1.getReadUnmappedFlag() ? rec1 : rec2;

		mapped.setMateReferenceIndex(unmapped.getReferenceIndex());
		mapped.setMateAlignmentStart(unmapped.getAlignmentStart());
		mapped.setMateNegativeStrandFlag(unmapped.getReadNegativeStrandFlag());
		mapped.setMateUnmappedFlag(true);
		mapped.setInferredInsertSize(0);

		unmapped.setMateReferenceIndex(mapped.getReferenceIndex());
		unmapped.setMateAlignmentStart(mapped.getAlignmentStart());
		unmapped.setMateNegativeStrandFlag(mapped.getReadNegativeStrandFlag());
		unmapped.setMateUnmappedFlag(false);
		unmapped.setInferredInsertSize(0);
	}

	boolean firstIsFirst = rec1.getAlignmentStart() < rec2.getAlignmentStart();
	int insertSize = firstIsFirst ? SamPairUtil.computeInsertSize(rec1, rec2) : SamPairUtil.computeInsertSize(rec2,
			rec1);

	rec1.setInferredInsertSize(firstIsFirst ? insertSize : -insertSize);
	rec2.setInferredInsertSize(firstIsFirst ? -insertSize : insertSize);

}
 

SAMRecord advance() {
	while (true) {
		// Pull next record from stream
		if (!wrappedIterator.hasNext())
			return null;

		final SAMRecord record = wrappedIterator.next();
		// If beyond the end of this reference sequence, end iteration
		final int referenceIndex = record.getReferenceIndex();
		if (referenceIndex != mReferenceIndex) {
			if (referenceIndex < 0 || referenceIndex > mReferenceIndex) {
				return null;
			}
			// If before this reference sequence, continue
			continue;
		}
		if (mRegionStart == 0 && mRegionEnd == Integer.MAX_VALUE) {
			// Quick exit to avoid expensive alignment end calculation
			return record;
		}
		final int alignmentStart = record.getAlignmentStart();
		// If read is unmapped but has a coordinate, return it if the
		// coordinate is within
		// the query region, regardless of whether the mapped mate will
		// be returned.
		final int alignmentEnd;
		if (mQueryType == QueryType.STARTING_AT) {
			alignmentEnd = -1;
		} else {
			alignmentEnd = (record.getAlignmentEnd() != SAMRecord.NO_ALIGNMENT_START ? record.getAlignmentEnd()
					: alignmentStart);
		}

		if (alignmentStart > mRegionEnd) {
			// If scanned beyond target region, end iteration
			return null;
		}
		// Filter for overlap with region
		if (mQueryType == QueryType.CONTAINED) {
			if (alignmentStart >= mRegionStart && alignmentEnd <= mRegionEnd) {
				return record;
			}
		} else if (mQueryType == QueryType.OVERLAPPING) {
			if (alignmentEnd >= mRegionStart && alignmentStart <= mRegionEnd) {
				return record;
			}
		} else {
			if (alignmentStart == mRegionStart) {
				return record;
			}
		}
	}
}
 

private static QueryInterval queryIntervalFromSamRecord(final SAMRecord samRecord) {
    return new QueryInterval(samRecord.getReferenceIndex(), samRecord.getStart(), samRecord.getEnd());
}
 

/**
 * Helper function that writes reads from iterator it into writer out, updating each SAMRecord along the way
 * according to the newOrder mapping from dictionary index -> index.  Name is used for printing only.
 */
private void writeReads(final SAMFileWriter out,
                        final SAMRecordIterator it,
                        final Map<Integer, Integer> newOrder,
                        final String name) {
    long counter = 0;
    log.info("  Processing " + name);

    while (it.hasNext()) {
        counter++;
        final SAMRecord read = it.next();
        final int oldRefIndex = read.getReferenceIndex();
        final int oldMateIndex = read.getMateReferenceIndex();
        final int newRefIndex = newOrderIndex(read, oldRefIndex, newOrder);

        read.setHeader(out.getFileHeader());
        read.setReferenceIndex(newRefIndex);

        // read becoming unmapped
        if (oldRefIndex != NO_ALIGNMENT_REFERENCE_INDEX &&
                newRefIndex == NO_ALIGNMENT_REFERENCE_INDEX) {
            read.setAlignmentStart(NO_ALIGNMENT_START);
            read.setReadUnmappedFlag(true);
            read.setCigarString(SAMRecord.NO_ALIGNMENT_CIGAR);
            read.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);
        }

        final int newMateIndex = newOrderIndex(read, oldMateIndex, newOrder);
        if (oldMateIndex != NO_ALIGNMENT_REFERENCE_INDEX &&
                newMateIndex == NO_ALIGNMENT_REFERENCE_INDEX) { // mate becoming unmapped
            read.setMateAlignmentStart(NO_ALIGNMENT_START);
            read.setMateUnmappedFlag(true);
            read.setAttribute(SAMTag.MC.name(), null);      // Set the Mate Cigar String to null
        }
        read.setMateReferenceIndex(newMateIndex);

        out.addAlignment(read);
    }

    it.close();
    log.info("Wrote " + counter + " reads");
}
 

private void populateNext(boolean force) throws IOException {
  final SAMRecord previousRecord = mNextRecord;
  mNextRecord = null;
  if (force) {
    advanceSubIterator();
  }
  while (mCurrentIt != null) {
    if (!mCurrentIt.hasNext()) {  // Only happens when stream is exhausted, so effectively just closes things out.
      advanceSubIterator();
    } else {
      final SAMRecord rec = nextOrBuffered();
      final int refId = rec.getReferenceIndex();

      if (refId > mCurrentTemplate) { // current offset has exceeded region and block overlapped next template
        recordPreviousAndAdvance(previousRecord, rec);
      } else {

        if (refId < mCurrentTemplate) { // Current block may occasionally return records from the previous template if the block overlaps
          continue;
        }

        final int alignmentStart = rec.getAlignmentStart() - 1; // to 0-based
        int alignmentEnd = rec.getAlignmentEnd(); // to 0-based exclusive = noop
        if (alignmentEnd == 0) {
          // Use the read length to get an end point for unmapped reads
          alignmentEnd = rec.getAlignmentStart() + rec.getReadLength();
        }

        if (alignmentEnd <= mCurrentRegion.getStart()) {  // before region
          continue;
        }

        if (alignmentStart <= mPreviousAlignmentStart) { // this record would have been already returned by an earlier region
          //Diagnostic.developerLog("Ignoring record from earlier block at " + rec.getReferenceName() + ":" + rec.getAlignmentStart());
          ++mDoubleFetched;
          if (mDoubleFetched % 100000 == 0) {
            Diagnostic.developerLog("Many double-fetched records for source " + mLabel + " noticed at " + rec.getReferenceName() + ":" + rec.getAlignmentStart() + " in region " + mCurrentRegion + " (skipping through to " + mPreviousAlignmentStart + ")");
          }
          continue;
        }

        if (alignmentStart >= mCurrentRegion.getEnd()) {  // past current region, advance the iterator and record the furtherest we got
          recordPreviousAndAdvance(previousRecord, rec);
          continue;
        }

        mNextRecord = rec;
        break;
      }
    }
  }
}
 

/**
 * Convenience method. See
 * {@link net.sf.cram.AlignmentsTags#calculateMdAndNmTags(SAMRecord, byte[], int, boolean, boolean)}
 * for details.
 * 
 * @param samRecord
 * @param referenceSource
 * @param samSequenceDictionary
 * @param calculateMdTag
 * @param calculateNmTag
 * @throws IOException
 */
public static void calculateMdAndNmTags(SAMRecord samRecord, ReferenceSource referenceSource,
		SAMSequenceDictionary samSequenceDictionary, boolean calculateMdTag, boolean calculateNmTag)
		throws IOException {
	if (!samRecord.getReadUnmappedFlag() && samRecord.getReferenceIndex() != SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
		SAMSequenceRecord sequence = samSequenceDictionary.getSequence(samRecord.getReferenceIndex());
		ReferenceRegion region = referenceSource.getRegion(sequence, samRecord.getAlignmentStart(),
				samRecord.getAlignmentEnd());
		calculateMdAndNmTags(samRecord, region.array, samRecord.getAlignmentStart() - 1, calculateMdTag,
				calculateNmTag);
	}
}