Java SAMRecord.getBaseQualities方法代码示例

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
public static String encodeBQTag(SAMRecord read, byte[] baq) {
    // Offset to base alignment quality (BAQ), of the same length as the read sequence.
    // At the i-th read base, BAQi = Qi - (BQi - 64) where Qi is the i-th base quality.
    // so BQi = Qi - BAQi + 64
    byte[] bqTag = new byte[baq.length];
    for ( int i = 0; i < bqTag.length; i++) {
        final int bq = (int)read.getBaseQualities()[i] + 64;
        final int baq_i = (int)baq[i];
        final int tag = bq - baq_i;
        // problem with the calculation of the correction factor; this is our problem
        if ( tag < 0 )
            throw new ReviewedGATKException("BAQ tag calculation error.  BAQ value above base quality at " + read);
        // the original quality is too high, almost certainly due to using the wrong encoding in the BAM file
        if ( tag > Byte.MAX_VALUE )
            throw new UserException.MisencodedBAM(read, "we encountered an extremely high quality score (" + (int)read.getBaseQualities()[i] + ") with BAQ correction factor of " + baq_i);
        bqTag[i] = (byte)tag;
    }
    return new String(bqTag);
}
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
/**
 * Returns a new qual array for read that includes the BAQ adjustment.  Does not support on-the-fly BAQ calculation
 *
 * @param read the SAMRecord to operate on
 * @param overwriteOriginalQuals If true, we replace the original qualities scores in the read with their BAQ'd version
 * @param useRawQualsIfNoBAQTag If useRawQualsIfNoBAQTag is true, then if there's no BAQ annotation we just use the raw quality scores.  Throws IllegalStateException is false and no BAQ tag is present
 * @return
 */
public static byte[] calcBAQFromTag(SAMRecord read, boolean overwriteOriginalQuals, boolean useRawQualsIfNoBAQTag) {
    byte[] rawQuals = read.getBaseQualities();
    byte[] newQuals = rawQuals;
    byte[] baq = getBAQTag(read);

    if ( baq != null ) {
        // Offset to base alignment quality (BAQ), of the same length as the read sequence.
        // At the i-th read base, BAQi = Qi - (BQi - 64) where Qi is the i-th base quality.
        newQuals = overwriteOriginalQuals ? rawQuals : new byte[rawQuals.length];
        for ( int i = 0; i < rawQuals.length; i++) {
            int rawQual = (int)rawQuals[i];
            int baq_delta = (int)baq[i] - 64;
            int newval =  rawQual - baq_delta;
            if ( newval < 0 )
                throw new UserException.MalformedBAM(read, "BAQ tag error: the BAQ value is larger than the base quality");
            newQuals[i] = (byte)newval;
        }
    } else if ( ! useRawQualsIfNoBAQTag ) {
        throw new IllegalStateException("Required BAQ tag to be present, but none was on read " + read.getReadName());
    }

    return newQuals;
}
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
private int getNumClippedBasesAtStart(final SAMRecord read) {
    // compute total number of clipped bases (soft or hard clipped)
    // check for hard clips (never consider these bases):
    final Cigar c = read.getCigar();
    final CigarElement first = c.getCigarElement(0);

    int numStartClippedBases = 0;
    if (first.getOperator() == CigarOperator.H) {
        numStartClippedBases = first.getLength();
    }
    final byte[] unclippedReadBases = read.getReadBases();
    final byte[] unclippedReadQuals = read.getBaseQualities();

    // Do a stricter base clipping than provided by CIGAR string, since this one may be too conservative,
    // and may leave a string of Q2 bases still hanging off the reads.
    for (int i = numStartClippedBases; i < unclippedReadBases.length; i++) {
        if (unclippedReadQuals[i] < PairHMMIndelErrorModel.BASE_QUAL_THRESHOLD)
            numStartClippedBases++;
        else
            break;

    }

    return numStartClippedBases;
}
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
/**
 * HACK TO CREATE GATKSAMRECORD BASED ONLY A SAMRECORD FOR TESTING PURPOSES ONLY
 * @param read
 */
public GATKSAMRecord(final SAMRecord read) {
    super(read.getHeader());
    setReferenceIndex(read.getReferenceIndex());
    setAlignmentStart(read.getAlignmentStart());
    mReadNameLength = (short)read.getReadNameLength();
    setMappingQuality(read.getMappingQuality());
    mCigarLength = read.getCigarLength();
    setFlags(read.getFlags());
    mReadLength = read.getReadLength();
    setMateReferenceIndex(read.getMateReferenceIndex());
    setMateAlignmentStart(read.getMateAlignmentStart());
    setInferredInsertSize(read.getInferredInsertSize());
    mRestOfBinaryData = null;

    super.setReadName(read.getReadName());
    super.setCigarString(read.getCigarString());
    super.setReadBases(read.getReadBases());
    super.setBaseQualities(read.getBaseQualities());

    SAMReadGroupRecord samRG = read.getReadGroup();
    clearAttributes();
    if (samRG != null) {
        GATKSAMReadGroupRecord rg = new GATKSAMReadGroupRecord(samRG);
        setReadGroup(rg);
    }
    List<SAMTagAndValue> attributes = read.getAttributes();
    for(SAMTagAndValue tagAndValue : attributes) {
        setAttribute(tagAndValue.tag, tagAndValue.value);
    }
}
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
/**
 * Modifies read in place so that the base quality scores are capped by the BAQ calculation.  Uses the BAQ
 * tag if present already and alwaysRecalculate is false, otherwise fires up the HMM and does the BAQ on the fly
 * using the refReader to obtain the reference bases as needed.
 *
 * @param read
 * @param refContentProvider
 * @param calculationType
 * @return BQ qualities for use, in case qmode is DONT_MODIFY
 */
public byte[] baqRead(SAMRecord read, RefContentProvider refContentProvider, CalculationMode calculationType, QualityMode qmode ) {
    if ( DEBUG ) System.out.printf("BAQ %s read %s%n", calculationType, read.getReadName());

    byte[] BAQQuals = read.getBaseQualities();      // in general we are overwriting quals, so just get a pointer to them
    if ( calculationType == CalculationMode.OFF) { // we don't want to do anything
        ; // just fall though
    } else if ( excludeReadFromBAQ(read) ) {
        ; // just fall through
    } else {
        final boolean readHasBAQTag = hasBAQTag(read);

        if ( calculationType == CalculationMode.RECALCULATE || ! readHasBAQTag ) {
            if ( DEBUG ) System.out.printf("  Calculating BAQ on the fly%n");
            BAQCalculationResult hmmResult = calcBAQFromHMM(read, refContentProvider);
            if ( hmmResult != null ) {
                switch ( qmode ) {
                    case ADD_TAG:         addBAQTag(read, hmmResult.bq); break;
                    case OVERWRITE_QUALS: System.arraycopy(hmmResult.bq, 0, read.getBaseQualities(), 0, hmmResult.bq.length); break;
                    case DONT_MODIFY:     BAQQuals = hmmResult.bq; break;
                    default:              throw new ReviewedGATKException("BUG: unexpected qmode " + qmode);
                }
            } else if ( readHasBAQTag ) {
                // remove the BAQ tag if it's there because we cannot trust it
                read.setAttribute(BAQ_TAG, null);
            }
        } else if ( qmode == QualityMode.OVERWRITE_QUALS ) { // only makes sense if we are overwriting quals
            if ( DEBUG ) System.out.printf("  Taking BAQ from tag%n");
            // this overwrites the original qualities
            calcBAQFromTag(read, true, false);
        }
    }

    return BAQQuals;
}
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
public BAQCalculationResult(SAMRecord read, byte[] ref) {
    this(read.getBaseQualities(), read.getReadBases(), ref);
}