use of htsjdk.samtools.AlignmentBlock in project gatk by broadinstitute.
the class CollectSequencingArtifactMetrics method acceptRead.
@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 UserException("Record contains library that is missing from header: " + library);
}
// 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 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;
final int contextFullLength = 2 * CONTEXT_SIZE + 1;
if (contextStartIndex < 0 || contextStartIndex + contextFullLength > ref.length())
continue;
// skip contexts with N bases
final String context = StringUtil.bytesToString(ref.getBases(), contextStartIndex, contextFullLength).toUpperCase();
if (context.contains("N"))
continue;
// skip low BQ sites
if (failsBaseQualityCutoff(readPos, rec))
continue;
// skip N bases in read
final char readBase = Character.toUpperCase((char) rec.getReadBases()[readPos - 1]);
if (readBase == 'N')
continue;
// count the base!
artifactCounters.get(library).countRecord(context, readBase, rec);
}
}
}
Aggregations