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

@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]]++;
        }
    }
}
 

void addRecord(final SAMRecord rec) {
    final byte[] quals = (useOriginalQualities ? rec.getOriginalBaseQualities() : rec.getBaseQualities());
    if (quals == null) return;

    final int length = quals.length;
    final boolean rc = rec.getReadNegativeStrandFlag();
    ensureArraysBigEnough(length+1);

    for (int i=0; i<length; ++i) {
        final int cycle = rc ? length-i : i+1;

        if (rec.getReadPairedFlag() && rec.getSecondOfPairFlag()) {
            secondReadTotalsByCycle[cycle] += quals[i];
            secondReadCountsByCycle[cycle] += 1;
        }
        else {
            firstReadTotalsByCycle[cycle] += quals[i];
            firstReadCountsByCycle[cycle] += 1;
        }
    }
}
 

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

/**
 * Takes an individual SAMRecord and applies the set of changes/reversions to it that
 * have been requested by program level options.
 */
public void revertSamRecord(final SAMRecord rec) {
    if (RESTORE_ORIGINAL_QUALITIES) {
        final byte[] oq = rec.getOriginalBaseQualities();
        if (oq != null) {
            rec.setBaseQualities(oq);
            rec.setOriginalBaseQualities(null);
        }
    }

    if (REMOVE_DUPLICATE_INFORMATION) {
        rec.setDuplicateReadFlag(false);
    }

    if (REMOVE_ALIGNMENT_INFORMATION) {
        if (rec.getReadNegativeStrandFlag()) {
            rec.reverseComplement(true);
            rec.setReadNegativeStrandFlag(false);
        }

        // Remove all alignment based information about the read itself
        rec.setReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
        rec.setAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
        rec.setCigarString(SAMRecord.NO_ALIGNMENT_CIGAR);
        rec.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);

        rec.setInferredInsertSize(0);
        rec.setNotPrimaryAlignmentFlag(false);
        rec.setProperPairFlag(false);
        rec.setReadUnmappedFlag(true);

        // Then remove any mate flags and info related to alignment
        rec.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
        rec.setMateNegativeStrandFlag(false);
        rec.setMateReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
        rec.setMateUnmappedFlag(rec.getReadPairedFlag());

        if (RESTORE_HARDCLIPS) {
            String hardClippedBases = rec.getStringAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG);
            String hardClippedQualities = rec.getStringAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG);
            if (hardClippedBases != null && hardClippedQualities != null) {
                // Record has already been reverse complemented if this was on the negative strand
                rec.setReadString(rec.getReadString() + hardClippedBases);
                rec.setBaseQualities(SAMUtils.fastqToPhred(SAMUtils.phredToFastq(rec.getBaseQualities()) + hardClippedQualities));

                // Remove hard clipping storage tags
                rec.setAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG, null);
                rec.setAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG, null);
            }
        }

        // And then remove any tags that are calculated from the alignment
        ATTRIBUTE_TO_CLEAR.forEach(tag -> rec.setAttribute(tag, null));
    }
}
 

/**
 * Checks if we can skip the SAM/BAM file when reverting origin base qualities and adding mate cigars.
 *
 * @param inputFile                   the SAM/BAM input file
 * @param maxRecordsToExamine         the maximum number of records to examine before quitting
 * @param revertOriginalBaseQualities true if we are to revert original base qualities, false otherwise
 * @return whether we can skip or not, and the explanation why.
 */
public static CanSkipSamFile canSkipSAMFile(final File inputFile, final int maxRecordsToExamine, boolean revertOriginalBaseQualities,
                                            final File referenceFasta) {
    final SamReader in = SamReaderFactory.makeDefault().referenceSequence(referenceFasta).enable(SamReaderFactory.Option.EAGERLY_DECODE).open(inputFile);
    final Iterator<SAMRecord> iterator = in.iterator();
    int numRecordsExamined = 0;
    CanSkipSamFile returnType = CanSkipSamFile.FOUND_NO_EVIDENCE;

    while (iterator.hasNext() && numRecordsExamined < maxRecordsToExamine) {
        final SAMRecord record = iterator.next();

        if (revertOriginalBaseQualities && null != record.getOriginalBaseQualities()) {
            // has OQ, break and return case #2
            returnType = CanSkipSamFile.CANNOT_SKIP_FOUND_OQ;
            break;
        }

        // check if mate pair and its mate is mapped
        if (record.getReadPairedFlag() && !record.getMateUnmappedFlag()) {
            if (null == SAMUtils.getMateCigar(record)) {
                // has no MC, break and return case #2
                returnType = CanSkipSamFile.CANNOT_SKIP_FOUND_NO_MC;
                break;
            } else {
                // has MC, previously checked that it does not have OQ, break and return case #1
                returnType = CanSkipSamFile.CAN_SKIP;
                break;
            }
        }

        numRecordsExamined++;
    }

    // no more records anyhow, so we can skip
    if (!iterator.hasNext() && CanSkipSamFile.FOUND_NO_EVIDENCE == returnType) {
        returnType = CanSkipSamFile.CAN_SKIP;
    }

    CloserUtil.close(in);

    return returnType;
}
 

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

/**
 *
 */
private Counts computeAlleleFraction(final SamLocusIterator.LocusInfo info, final byte refBase) {
    final Counts counts = new Counts();
    final byte altBase = (refBase == 'C') ? (byte) 'A' : (byte) 'T';

    for (final SamLocusIterator.RecordAndOffset rec : info.getRecordAndOffsets()) {
        final byte qual;
        final SAMRecord samrec = rec.getRecord();

        if (USE_OQ) {
            final byte[] oqs = samrec.getOriginalBaseQualities();
            if (oqs != null) qual = oqs[rec.getOffset()];
            else qual = rec.getBaseQuality();
        } else {
            qual = rec.getBaseQuality();
        }

        // Skip if below qual, or if library isn't a match
        if (qual < MINIMUM_QUALITY_SCORE) continue;
        if (!this.library.equals(Optional.ofNullable(samrec.getReadGroup().getLibrary()).orElse(UNKNOWN_LIBRARY))) continue;

        // Get the read base, and get it in "as read" orientation
        final byte base = rec.getReadBase();
        final byte baseAsRead = samrec.getReadNegativeStrandFlag() ? SequenceUtil.complement(base) : base;
        final int read = samrec.getReadPairedFlag() && samrec.getSecondOfPairFlag() ? 2 : 1;

        // Figure out how to count the alternative allele. If the damage is caused by oxidation of G
        // during shearing (in non-rnaseq data), then we know that:
        //     G>T observation is always in read 1
        //     C>A observation is always in read 2
        // But if the substitution is from other causes the distribution of A/T across R1/R2 will be
        // random.
        if (base == refBase) {
            if (baseAsRead == 'G' && read == 1) ++counts.oxidatedC;
            else if (baseAsRead == 'G' && read == 2) ++counts.controlC;
            else if (baseAsRead == 'C' && read == 1) ++counts.controlC;
            else if (baseAsRead == 'C' && read == 2) ++counts.oxidatedC;
        } else if (base == altBase) {
            if (baseAsRead == 'T' && read == 1) ++counts.oxidatedA;
            else if (baseAsRead == 'T' && read == 2) ++counts.controlA;
            else if (baseAsRead == 'A' && read == 1) ++counts.controlA;
            else if (baseAsRead == 'A' && read == 2) ++counts.oxidatedA;
        }
    }

    return counts;
}
 

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