Search in sources :

Example 1 with AlignmentBlock

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);
        }
    }
}
Also used : UserException(org.broadinstitute.hellbender.exceptions.UserException) AlignmentBlock(htsjdk.samtools.AlignmentBlock)

Aggregations

AlignmentBlock (htsjdk.samtools.AlignmentBlock)1 UserException (org.broadinstitute.hellbender.exceptions.UserException)1