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

public int scoreBaseQuality(final SAMRecord barcodedRead) {
	int numBasesBelowQuality=0;
	byte [] qual= barcodedRead.getBaseQualities();
	char [] seq = barcodedRead.getReadString().toUpperCase().toCharArray();

	for (BaseRange b: baseRanges)
		for (int i=b.getStart()-1; i<b.getEnd(); i++) {
			if (i> (seq.length-1))
				throw new TranscriptomeException("Base [" + Integer.toString(i+1) + "] was requested, but the read isn't long enough ["+ barcodedRead.getReadString()+"]");

			byte q = qual[i];
			char s = seq[i];

			if (q < this.baseQualityThrehsold || s=='N')
				numBasesBelowQuality++;
		}

	this.metric.addFailedBase(numBasesBelowQuality);
	return (numBasesBelowQuality);
}
 

public void addQualityScores(final SAMRecord s, final CramCompressionRecord r, final ReferenceTracks t) {
	if (s.getBaseQualities() == SAMRecord.NULL_QUALS) {
		r.qualityScores = SAMRecord.NULL_QUALS;
		r.setForcePreserveQualityScores(false);
		return;
	}

	final byte[] scores = new byte[s.getReadLength()];
	Arrays.fill(scores, (byte) -1);
	for (final PreservationPolicy p : policyList)
		addQS(s, r, scores, t, p);

	if (!r.isForcePreserveQualityScores()) {
		for (int i = 0; i < scores.length; i++) {
			if (scores[i] > -1) {
				if (r.readFeatures == null || r.readFeatures == Collections.EMPTY_LIST)
					r.readFeatures = new LinkedList<ReadFeature>();

				r.readFeatures.add(new BaseQualityScore(i + 1, scores[i]));
			}
		}
		if (r.readFeatures != null)
			Collections.sort(r.readFeatures, readFeaturePositionComparator);
	}
	r.qualityScores = scores;
}
 

@Override
protected void acceptRead(final SAMRecord rec, final ReferenceSequence ref) {
    // Skip unwanted records
    if (PF_READS_ONLY && rec.getReadFailsVendorQualityCheckFlag()) return;
    if (ALIGNED_READS_ONLY && rec.getReadUnmappedFlag()) return;
    if (rec.isSecondaryOrSupplementary()) return;

    final byte[] bases = rec.getReadBases();
    final byte[] quals = rec.getBaseQualities();
    final byte[] oq    = rec.getOriginalBaseQualities();

    final int length = quals.length;

    for (int i=0; i<length; ++i) {
        if (INCLUDE_NO_CALLS || !SequenceUtil.isNoCall(bases[i])) {
            qCounts[quals[i]]++;
            if (oq != null) oqCounts[oq[i]]++;
        }
    }
}
 

private boolean hasHighQualitySoftClipping(SAMRecord read, int start, int length) {
	
	int numHighQualBases = 0;
	int requiredHighQualBases = (int) (softClipFraction * length);
	
	for (int bq = start; bq < start+length; bq++) {
		if (read.getBaseQualities()[bq] >= minBaseQual) {
			numHighQualBases += 1;
			
			if (numHighQualBases >= requiredHighQualBases) {
				return true;
			}
		}
	}
	
	return false;
}
 

private double baseQuality(int startReadIndex, @NotNull final SAMRecord record, int length) {
    int maxIndex = Math.min(startReadIndex + length, record.getBaseQualities().length) - 1;
    int maxLength = maxIndex - startReadIndex + 1;

    double quality = Integer.MAX_VALUE;
    for (int i = 0; i < maxLength; i++) {
        int refPosition = (int) position() + i;
        int readIndex = startReadIndex + i;
        byte rawQuality = record.getBaseQualities()[readIndex];
        byte[] trinucleotideContext = readContext.refTrinucleotideContext(refPosition);
        double recalibratedQuality =
                qualityRecalibrationMap.quality((byte) ref().charAt(i), (byte) alt().charAt(i), trinucleotideContext, rawQuality);
        quality = Math.min(quality, recalibratedQuality);
    }

    return quality;
}
 

private static boolean hasOriginalQualScores( final File bam ) throws IOException {
    try ( final SamReader reader = SamReaderFactory.makeDefault().open(bam) ) {
        for ( SAMRecord read : reader ) {
            final byte[] originalBaseQualities = read.getOriginalBaseQualities();
            Assert.assertNotNull(originalBaseQualities);
            final byte[] baseQualities = read.getBaseQualities();
            Assert.assertEquals(originalBaseQualities.length, baseQualities.length);
            for (int i = 0; i < originalBaseQualities.length; ++i) {
                if (originalBaseQualities[i] != baseQualities[i]) {
                    return false;
                }
            }
        }
    }
    return true;
}
 

int avgCentreQuality(int readIndex, @NotNull final SAMRecord record) {
    int leftOffset = this.readIndex() - readBases.leftCentreIndex();
    int rightOffset = readBases.rightCentreIndex() - this.readIndex();

    int leftIndex = readIndex - leftOffset;
    int rightIndex = readIndex + rightOffset;

    float quality = 0;
    for (int i = leftIndex; i <= rightIndex; i++) {
        quality += record.getBaseQualities()[i];
    }
    return Math.round(quality / (rightIndex - leftIndex + 1));
}
 

@Override
public void handleAlignment(@NotNull final SAMRecord r, @NotNull final CigarElement e, final int startReadIndex, final int refPos) {
    for (int i = 0; i < e.getLength(); i++) {
        int readIndex = startReadIndex + i;
        int position = refPos + i;

        if (position > bounds.end()) {
            return;
        }

        if (position < bounds.start()) {
            continue;
        }

        byte ref = refGenome.base(position);
        byte alt = r.getReadBases()[readIndex];
        byte quality = r.getBaseQualities()[readIndex];
        byte[] trinucleotideContext = refGenome.trinucleotideContext(position);

        if (alt != N && isValid(trinucleotideContext)) {
            final QualityCounterKey key = ImmutableQualityCounterKey.builder()
                    .ref(ref)
                    .alt(alt)
                    .qual(quality)
                    .position(position)
                    .trinucleotideContext(trinucleotideContext)
                    .build();
            qualityMap.computeIfAbsent(key, QualityCounter::new).increment();
        }
    }
}
 

private Object[] getBaseAndQualAtPosition(SAMRecord read, int refPos) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= refPos && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case S:
			case I:
				readPos += elem.getLength();
				break;
			case D:
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				if (refPos < (refPosInRead + elem.getLength())) {
					readPos += refPos - refPosInRead;
					if (readPos < read.getReadLength()) {
						// Found the base.  Return it
						return new Object[] { read.getReadString().charAt(readPos) , (int) read.getBaseQualities()[readPos] };
					}
				} else {
					readPos += elem.getLength();
					refPosInRead += elem.getLength();
				}
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return new Object[] { 'N', (int) 0 };
}
 

private static ScoreDiff[] findScoreDiffs(SAMRecord r1, SAMRecord r2, byte[] ref, List<PreservationPolicy> policies) {
	if (!r1.getCigarString().equals(r2.getCigarString()))
		throw new RuntimeException("CIGAR string are different.");

	List<ScoreDiff> diffs = new ArrayList<SamRecordComparision.ScoreDiff>();
	int posInRead = 1, posInRef = r1.getAlignmentStart();

	for (CigarElement ce : r1.getCigar().getCigarElements()) {
		if (ce.getOperator().consumesReadBases()) {
			for (int i = 0; i < ce.getLength(); i++) {
				if (r1.getBaseQualities()[i + posInRead - 1] != r2.getBaseQualities()[i + posInRead - 1]) {
					ScoreDiff d = new ScoreDiff();
					d.base = r1.getReadBases()[i + posInRead - 1];
					d.refBase = ref[i + posInRef - 1];
					ce.getOperator();
					d.cigarOp = CigarOperator.enumToCharacter(ce.getOperator());
					d.oScore = r1.getBaseQualities()[i + posInRead - 1];
					d.score = r2.getBaseQualities()[i + posInRead - 1];
					d.posInRead = i + posInRead;

					if (d.score == Binning.Illumina_binning_matrix[d.oScore])
						d.treatment = QualityScoreTreatmentType.BIN;
					else if (d.score == DEFAULT_SCORE)
						d.treatment = QualityScoreTreatmentType.DROP;

					if (policies == null || !obeysLossyModel(policies, r1, d))
						diffs.add(d);
				}
			}
		}

		posInRead += ce.getOperator().consumesReadBases() ? ce.getLength() : 0;
		posInRef += ce.getOperator().consumesReferenceBases() ? ce.getLength() : 0;
	}
	return diffs.toArray(new ScoreDiff[diffs.size()]);
}
 

public static QENCODING guessEncoding(final SAMRecord read) throws QualityException {
    final byte[] quals = read.getBaseQualities();
    byte max = quals[0];
    for ( int i = 0; i < quals.length; i++ ) {
        if(quals[i] > max) max =quals[i];
    }
    Logger.DEBUG("Max quality: " + max, 3);
    if(max <= REASONABLE_SANGER_THRESHOLD) {
        Logger.DEBUG("SANGER Quality encoding");
        return QENCODING.SANGER;
    } else {
        Logger.DEBUG("ILLUMINA Quality encoding");
        return QENCODING.ILLUMINA;
    }
}
 

/** Get the the number of bases in the given alignment block and record that have base quality greater or equal to the minimum */
public static int getNumBasesPassingMinimumBaseQuality(final SAMRecord record, final AlignmentBlock block, final int minimumBaseQuality) {
    int basesInBlockAtMinimumQuality = 0;
    final byte[] baseQualities = record.getBaseQualities();
    for (int idx = block.getReadStart(); idx <= CoordMath.getEnd(block.getLength(), block.getReadStart()); idx++) { // idx is one-based
        if (minimumBaseQuality <= baseQualities[idx-1]) basesInBlockAtMinimumQuality++;
    }
    return basesInBlockAtMinimumQuality;
}
 

private BaseInfo getBaseAtReferencePosition(SAMRecord read, int refPos) {
	boolean isNearEdge = false;
	boolean isIndel = false;
	int alignmentStart = read.getAlignmentStart();
	Cigar cigar = read.getCigar();
	
	char base = 'N';
	
	int readIdx = 0;
	int currRefPos = alignmentStart;
	
	for (CigarElement element : cigar.getCigarElements()) {
					
		if (element.getOperator() == CigarOperator.M) {
			readIdx += element.getLength();
			currRefPos += element.getLength();
			
			if (currRefPos > refPos) {  // TODO: double check end of read base here...
				
				int offset = currRefPos - refPos;
				
				if ((offset < 3) || (offset+3 >= element.getLength())) {
					// This position is within 3 bases of start/end of alignment or clipping or indel.
					// Tag this base for evaluation downstream.
					isNearEdge = true;
				}
				
				readIdx -= offset;
				if ((readIdx < 0) || (readIdx >= read.getReadBases().length)) {
					System.err.println("Read index out of bounds for read: " + read.getSAMString());
					break;
				}
				
				if (read.getBaseQualities()[readIdx] >= MIN_BASE_QUALITY) {
					base = (char) read.getReadBases()[readIdx];
				}
				break;
			}
		} else if (element.getOperator() == CigarOperator.I) {
			if (currRefPos == refPos) {
				//TODO: Handle insertions
				isIndel = true;
				break;
			}
			readIdx += element.getLength();
		} else if (element.getOperator() == CigarOperator.D) {
			if (refPos >= currRefPos && refPos <= currRefPos+element.getLength()) {
				//TODO: Handle deletions
				isIndel = true;
				break;
			}				
			currRefPos += element.getLength();
		} else if (element.getOperator() == CigarOperator.S) {
			readIdx += element.getLength();
		}			
	}
	
	return new BaseInfo(Character.toUpperCase(base), isNearEdge, isIndel);
}
 

private int getBaseQuality(SAMRecord read, int index) {
	return (char) read.getBaseQualities()[index];
}
 

/**
 * Unroll both the read bases and quality data for a CG alignment
 * @param rec the alignment
 * @return the unrolled read
 */
public static SamSequence unrollCgRead(SAMRecord rec) {
  final int projectedSplit = rec.getAlignmentStart() * ((rec.getFlags() & SamBamConstants.SAM_MATE_IS_REVERSE) != 0 ? 1 : -1);
  final byte flags = SamSequence.getFlags(rec);
  final byte[] expandedRead;
  final byte[] expandedQual;
  if (rec.getReadUnmappedFlag()) {
    expandedRead = rec.getReadBases();
    expandedQual = rec.getBaseQualities();
  } else {
    final String gc = rec.getStringAttribute(SamUtils.ATTRIBUTE_CG_RAW_READ_INSTRUCTIONS);
    if (gc == null) {
      throw new NoTalkbackSlimException("SAM Record does not contain CGI read reconstruction attribute: " + rec.getSAMString());
    }
    final byte[] gq = FastaUtils.asciiToRawQuality(SamUtils.allowEmpty(rec.getStringAttribute(SamUtils.ATTRIBUTE_CG_OVERLAP_QUALITY)));
    final byte[] gs = SamUtils.allowEmpty(rec.getStringAttribute(SamUtils.ATTRIBUTE_CG_OVERLAP_BASES)).getBytes();
    final boolean legacyLegacy = gq.length == gs.length / 2;
    expandedRead = unrollLegacyRead(rec.getReadBases(), gs, gc);
    if (expandedRead == null) {
      throw new NoTalkbackSlimException("Could not reconstruct read bases for record: " + rec.getSAMString());
    }
    if (rec.getBaseQualities().length == 0) {
      expandedQual = rec.getBaseQualities();
    } else {
      if (!legacyLegacy && gq.length != gs.length) {
        throw new NoTalkbackSlimException("Unexpected length of CG quality information: " + rec.getSAMString());
      }
      final byte[] samQualities = rec.getBaseQualities();
      if (legacyLegacy) {
        expandedQual = unrollLegacyLegacyQualities(samQualities, gq, gc);
      } else {
        final byte[] bytes = unrollLegacyRead(samQualities, gq, gc);
        expandedQual = bytes;
      }
    }
  }

  final byte[] readBytes = CgUtils.unPad(expandedRead, !rec.getReadNegativeStrandFlag());
  final byte[] readQual = CgUtils.unPad(expandedQual, !rec.getReadNegativeStrandFlag());

  return new SamSequence(rec.getReadName(), readBytes, readQual, flags, projectedSplit);
}
 

public void acceptRecord(final SAMRecord rec, final ReferenceSequence ref) {
    if (!this.includeSecondaryAlignments    && rec.getNotPrimaryAlignmentFlag()) return;
    if (!this.includeSupplementalAlignments && rec.getSupplementaryAlignmentFlag()) return;

    final int length = rec.getReadLength();
    metrics.TOTAL_READS++;
    metrics.TOTAL_BASES += length;

    final boolean isPfRead = !rec.getReadFailsVendorQualityCheckFlag();
    if (isPfRead) {
        metrics.PF_READS++;
        metrics.PF_BASES += length;
    }

    final byte[] quals;
    if (this.useOriginalQualities) {
        byte[] tmp = rec.getOriginalBaseQualities();
        if (tmp == null) tmp = rec.getBaseQualities();
        quals = tmp;
    } else {
        quals = rec.getBaseQualities();
    }

    // add up quals, and quals >= 20
    for (final int qual : quals) {
        metrics.Q20_EQUIVALENT_YIELD += qual;

        if (qual >= 30) {
            metrics.Q20_BASES++;
            metrics.Q30_BASES++;
        } else if (qual >= 20) {
            metrics.Q20_BASES++;
        }

        if (isPfRead) {
            metrics.PF_Q20_EQUIVALENT_YIELD += qual;
            if (qual >= 30) {
                metrics.PF_Q20_BASES++;
                metrics.PF_Q30_BASES++;
            } else if (qual >= 20) {
                metrics.PF_Q20_BASES++;
            }
        }
    }
}
 

@Override
protected void acceptRead(final SAMRecord rec, final ReferenceSequence ref) {
    // see if the whole read should be skipped
    if (recordFilter.filterOut(rec)) return;

    // check read group + library
    final String library = (rec.getReadGroup() == null) ? UNKNOWN_LIBRARY : getOrElse(rec.getReadGroup().getLibrary(), UNKNOWN_LIBRARY);
    if (!libraries.contains(library)) {
        // should never happen if SAM is valid
        throw new PicardException("Record contains library that is missing from header: " + library);
    }

    // set up some constants that don't change in the loop below
    final int contextFullLength = 2 * CONTEXT_SIZE + 1;
    final ArtifactCounter counter = artifactCounters.get(library);
    final byte[] readBases = rec.getReadBases();
    final byte[] readQuals;
    if (USE_OQ) {
        final byte[] tmp = rec.getOriginalBaseQualities();
        readQuals = tmp == null ? rec.getBaseQualities() : tmp;
    } else {
        readQuals = rec.getBaseQualities();
    }

    // iterate over aligned positions
    for (final AlignmentBlock block : rec.getAlignmentBlocks()) {
        for (int offset = 0; offset < block.getLength(); offset++) {
            // remember, these are 1-based!
            final int readPos = block.getReadStart() + offset;
            final int refPos = block.getReferenceStart() + offset;

            // skip low BQ sites
            final byte qual = readQuals[readPos - 1];
            if (qual < MINIMUM_QUALITY_SCORE) continue;

            // skip N bases in read
            final char readBase = Character.toUpperCase((char)readBases[readPos - 1]);
            if (readBase == 'N') continue;

            /**
             * Skip regions outside of intervals.
             *
             * NB: IntervalListReferenceSequenceMask.get() has side-effects which assume
             * that successive ReferenceSequence's passed to this method will be in-order
             * (e.g. it will break if you call acceptRead() with chr1, then chr2, then chr1
             * again). So this only works if the underlying iteration is coordinate-sorted.
             */
            if (intervalMask != null && !intervalMask.get(ref.getContigIndex(), refPos)) continue;

            // skip dbSNP sites
            if (dbSnpMask != null && dbSnpMask.isDbSnpSite(ref.getName(), refPos)) continue;

            // skip the ends of the reference
            final int contextStartIndex = refPos - CONTEXT_SIZE - 1;
            if (contextStartIndex < 0 || contextStartIndex + contextFullLength > ref.length()) continue;

            // skip contexts with N bases
            final String context = getRefContext(ref, contextStartIndex, contextFullLength);
            if (context.contains("N")) continue;

            // skip non-ACGT bases
            if (!SequenceUtil.isUpperACGTN((byte)readBase))
                continue;

            // count the base!
            counter.countRecord(context, readBase, rec);
        }
    }
}
 

public void acceptRecord(final SAMRecordAndReference args) {
	mappedRecordCount++;

	final SAMRecord samRecord = args.getSamRecord();
	final ReferenceSequence referenceSequence = args.getReferenceSequence();

	final byte[] readBases = samRecord.getReadBases();
	final byte[] readQualities = samRecord.getBaseQualities();
	final byte[] refBases = referenceSequence.getBases();

	if (samRecord.getReadLength() < minReadLength) {
		smallReadCount++;
		return;
	} else if (SequenceUtil.countMismatches(samRecord, refBases, true) > Math.round(samRecord.getReadLength() * maxMismatchRate)) {
		mismatchCount++;
		return;
	}

	// We only record non-CpG C sites if there was at least one CpG in the read, keep track of
	// the values for this record and then apply to the global totals if valid
	int recordCpgs = 0;

	for (final AlignmentBlock alignmentBlock : samRecord.getAlignmentBlocks()) {
		final int blockLength = alignmentBlock.getLength();
		final int refFragmentStart = alignmentBlock.getReferenceStart() - 1;
		final int readFragmentStart = alignmentBlock.getReadStart() - 1;

		final byte[] refFragment = getFragment(refBases, refFragmentStart, blockLength);
		final byte[] readFragment = getFragment(readBases, readFragmentStart, blockLength);
		final byte[] readQualityFragment = getFragment(readQualities, readFragmentStart, blockLength);

		if (samRecord.getReadNegativeStrandFlag()) {
			// In the case of a negative strand, reverse (and complement for base arrays) the reference,
			// reads & qualities so that it can be treated as a positive strand for the rest of the process
			SequenceUtil.reverseComplement(refFragment);
			SequenceUtil.reverseComplement(readFragment);
			SequenceUtil.reverseQualities(readQualityFragment);
		}

		for (int i=0; i < blockLength-1; i++) {
			final int curRefIndex = getCurRefIndex(refFragmentStart, blockLength, i, samRecord.getReadNegativeStrandFlag());

			// Look at a 2-base window to see if we're on a CpG site, and if so check for conversion
			// (CG -> TG). We do not consider ourselves to be on a CpG site if we're on the last base of a read
			if ((SequenceUtil.basesEqual(refFragment[i], SequenceUtil.C)) &&
					(SequenceUtil.basesEqual(refFragment[i+1], SequenceUtil.G))) {
				// We want to catch the case where there's a CpG in the reference, even if it is not valid
				// to prevent the C showing up as a non-CpG C down below. Otherwise this could have been all
				// in one if statement
				if (isValidCpg(refFragment, readFragment, readQualityFragment, i)) {
					recordCpgs++;
					final CpgLocation curLocation = new CpgLocation(samRecord.getReferenceName(), curRefIndex);
					cpgTotal.increment(curLocation);
					if (SequenceUtil.isBisulfiteConverted(readFragment[i], refFragment[i])) {
						cpgConverted.increment(curLocation);
					}
				}
				i++;
			} else if (isC(refFragment[i], readFragment[i]) && isAboveCytoQcThreshold(readQualities, i) &&
					SequenceUtil.bisulfiteBasesEqual(false, readFragment[i+1], refFragment[i+1])) {
				// C base in the reference that's not associated with a CpG
				nCytoTotal++;
				if (SequenceUtil.isBisulfiteConverted(readFragment[i], refFragment[i])) {
					nCytoConverted++;
				}
			}
		}
	}

	if (recordCpgs == 0) {
		noCpgCount++;
	}
}
 

@NotNull
private List<AltRead> processAlignment(@NotNull final SAMRecord record, int readBasesStartIndex, int refPositionStart,
        int alignmentLength, final IndexedBases refBases, int numberOfEvents) {
    final List<AltRead> result = Lists.newArrayList();

    int refIndex = refBases.index(refPositionStart);

    for (int i = 0; i < alignmentLength; i++) {
        int refPosition = refPositionStart + i;
        int readBaseIndex = readBasesStartIndex + i;
        int refBaseIndex = refIndex + i;

        if (!inBounds(refPosition)) {
            continue;
        }

        final byte refByte = refBases.bases()[refBaseIndex];
        final String ref = String.valueOf((char) refByte);
        final byte readByte = record.getReadBases()[readBaseIndex];
        boolean findReadContext = findReadContext(readBaseIndex, record);

        final RefContext refContext = candidates.refContext(record.getContig(), refPosition);
        if (!refContext.reachedLimit()) {
            int baseQuality = record.getBaseQualities()[readBaseIndex];
            if (readByte != refByte) {
                final String alt = String.valueOf((char) readByte);
                final ReadContext readContext =
                        findReadContext ? readContextFactory.createSNVContext(refPosition, readBaseIndex, record, refBases) : null;

                result.add(new AltRead(refContext, ref, alt, baseQuality, numberOfEvents, readContext));

                if (config.mnvEnabled()) {
                    int mnvMaxLength = mnvLength(readBaseIndex, refBaseIndex, record.getReadBases(), refBases.bases());
                    for (int mnvLength = 2; mnvLength <= mnvMaxLength; mnvLength++) {

                        final String mnvRef = new String(refBases.bases(), refBaseIndex, mnvLength);
                        final String mnvAlt = new String(record.getReadBases(), readBaseIndex, mnvLength);

                        // Only check last base because some subsets may not be valid,
                        // ie CA > TA is not a valid subset of CAC > TAT
                        if (mnvRef.charAt(mnvLength - 1) != mnvAlt.charAt(mnvLength - 1)) {
                            final ReadContext mnvReadContext = findReadContext ? readContextFactory.createMNVContext(refPosition,
                                    readBaseIndex,
                                    mnvLength,
                                    record,
                                    refBases) : null;

                            result.add(new AltRead(refContext,
                                    mnvRef,
                                    mnvAlt,
                                    baseQuality,
                                    NumberEvents.numberOfEventsWithMNV(numberOfEvents, mnvRef, mnvAlt),
                                    mnvReadContext));
                        }
                    }
                }
            } else {
                refContext.refRead();
            }
        }
    }

    return result;
}
 

/**
 * For a read on this pileup, get the base and quality of the base that is at the same
 * position as the SNP.  If the read does not overlap the interval, then return an empty array.
 * @param r The read
 * @return A length 2 byte array containing the base and quality.  Empty if the read does not overlap the interval.
 */
public byte [] getBaseAndQualityOverlappingInterval (final SAMRecord r) {
	byte [] result = new byte [2];

	List<CigarElement> elements = r.getCigar().getCigarElements();
	Iterator<AlignmentBlock> blocks = r.getAlignmentBlocks().iterator();

	int finalSNPLocalPosition=-1;
	int lengthTraversed=0;

	for (CigarElement ce: elements) {
		CigarOperator co = ce.getOperator();
		// you're in an alignment block
		if (co==CigarOperator.M) {
			// get the next alignment block
			AlignmentBlock b = blocks.next();
			int refStart = b.getReferenceStart();
			int snpLocalPos=this.getPosition() - refStart +1;
			int blockLength=b.getLength();

			// is the local position inside this alignment block?
			// if not, move onto the next block.
			if (snpLocalPos >0 && snpLocalPos<=blockLength) {
				// found it!  Done.
				finalSNPLocalPosition=snpLocalPos+lengthTraversed;
				break;
			}
		}
		// consume the bases if necessary and move on to the next element.
		if (co.consumesReadBases())
			lengthTraversed+=ce.getLength();
	}

	// if the position is assigned, then add to the pileup.
	if (finalSNPLocalPosition!=-1) {
		// arrays 0 based.
		byte [] quals = r.getBaseQualities();
		byte [] bases = r.getReadBases();
		byte base = bases[finalSNPLocalPosition-1];
		// char baseChar = StringUtil.byteToChar(base);
		byte qual = quals[finalSNPLocalPosition-1];
		result[0]=base;
		result[1]=qual;
		return (result);
	}
	return (ArrayUtils.EMPTY_BYTE_ARRAY);
}