Java SAMRecord.getReadName方法代码示例

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包/类
/**
 * Write the record.
 *
 * @param alignment SAMRecord.
 */
public BasicSamRecord transfer(final SAMRecord alignment) {
    String readName = alignment.getReadName();
    int flags = alignment.getFlags();
    int contigId = alignment.getReferenceIndex();
    if (contigId == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
        contigId = FAKE_CONTIG_ID;
    }
    String contigName = alignment.getReferenceName();
    int position = alignment.getAlignmentStart();
    int mapQ = alignment.getMappingQuality();
    String cigar = alignment.getCigarString();
    int mateContigId = alignment.getMateReferenceIndex();
    if (mateContigId == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
        mateContigId = FAKE_CONTIG_ID;
    }
    String mateContigName = alignment.getMateReferenceName();
    int matePosition = alignment.getMateAlignmentStart();
    int inferredSize = alignment.getInferredInsertSize();
    byte[] sequence = alignment.getReadString().getBytes();
    byte[] quality = alignment.getBaseQualityString().getBytes();
    List<SAMRecord.SAMTagAndValue> attributes = alignment.getAttributes();
    List<String> encodedTags = new ArrayList<>(attributes.size());
    for (SAMRecord.SAMTagAndValue attribute : attributes) {
        encodedTags.add(tagCodec.encode(attribute.tag, attribute.value));
    }

    // 在这里统一不压缩，压缩操作是由scala代码中的执行引擎统一调配的。
    return BasicSamRecord.apply(false, readName, flags, contigId, contigName, position, mapQ,
            cigar, mateContigId, mateContigName, matePosition, inferredSize,
            sequence, quality, encodedTags);
}
 

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包/类
/**
 * Determine the appropriate start and stop offsets in the reads for the bases given the cigar string
 * @param read
 * @return
 */
private final Pair<Integer,Integer> calculateQueryRange(SAMRecord read) {
    int queryStart = -1, queryStop = -1;
    int readI = 0;

    // iterate over the cigar elements to determine the start and stop of the read bases for the BAQ calculation
    for ( CigarElement elt : read.getCigar().getCigarElements() ) {
        switch (elt.getOperator()) {
            case N:  return null; // cannot handle these
            case H : case P : case D: break; // ignore pads, hard clips, and deletions
            case I : case S: case M: case EQ: case X:
                int prev = readI;
                readI += elt.getLength();
                if ( includeClippedBases || elt.getOperator() != CigarOperator.S) {
                    if ( queryStart == -1 )
                        queryStart = prev;
                    queryStop = readI;
                }
                // in the else case we aren't including soft clipped bases, so we don't update
                // queryStart or queryStop
                break;
            default: throw new ReviewedGATKException("BUG: Unexpected CIGAR element " + elt + " in read " + read.getReadName());
        }
    }

    if ( queryStop == queryStart ) {
        // this read is completely clipped away, and yet is present in the file for some reason
        // usually they are flagged as non-PF, but it's possible to push them through the BAM
        //System.err.printf("WARNING -- read is completely clipped away: " + read.format());
        return null;
    }

    return new Pair<Integer, Integer>(queryStart, queryStop);
}
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
/**
 * Generic algorithm that takes an iterable over T objects, a getter routine to extract the reads in T,
 * and returns a FragmentCollection that contains the T objects whose underlying reads either overlap (or
 * not) with their mate pairs.
 *
 * @param readContainingObjects An iterator of objects that contain GATKSAMRecords
 * @param nElements the number of elements to be provided by the iterator, which is usually known upfront and
 *                  greatly improves the efficiency of the fragment calculation
 * @param getter a helper function that takes an object of type T and returns is associated GATKSAMRecord
 * @param <T>
 * @return a fragment collection
 */
private static <T> FragmentCollection<T> create(final Iterable<T> readContainingObjects, final int nElements, final ReadGetter<T> getter) {
    Collection<T> singletons = null;
    Collection<List<T>> overlapping = null;
    Map<String, T> nameMap = null;

    int lastStart = -1;

    // build an initial map, grabbing all of the multi-read fragments
    for ( final T p : readContainingObjects ) {
        final SAMRecord read = getter.get(p);

        if ( read.getAlignmentStart() < lastStart ) {
            throw new IllegalArgumentException(String.format(
                    "FragmentUtils.create assumes that the incoming objects are ordered by " +
                            "SAMRecord alignment start, but saw a read %s with alignment start " +
                            "%d before the previous start %d", read.getSAMString(), read.getAlignmentStart(), lastStart));
        }
        lastStart = read.getAlignmentStart();

        final int mateStart = read.getMateAlignmentStart();
        if ( mateStart == 0 || mateStart > read.getAlignmentEnd() ) {
            // if we know that this read won't overlap its mate, or doesn't have one, jump out early
            if ( singletons == null ) singletons = new ArrayList<T>(nElements); // lazy init
            singletons.add(p);
        } else {
            // the read might overlap it's mate, or is the rightmost read of a pair
            final String readName = read.getReadName();
            final T pe1 = nameMap == null ? null : nameMap.get(readName);
            if ( pe1 != null ) {
                // assumes we have at most 2 reads per fragment
                if ( overlapping == null ) overlapping = new ArrayList<List<T>>(); // lazy init
                overlapping.add(Arrays.asList(pe1, p));
                nameMap.remove(readName);
            } else {
                if ( nameMap == null ) nameMap = new HashMap<String, T>(nElements); // lazy init
                nameMap.put(readName, p);
            }
        }
    }

    // add all of the reads that are potentially overlapping but whose mate never showed
    // up to the oneReadPile
    if ( nameMap != null && ! nameMap.isEmpty() ) {
        if ( singletons == null )
            singletons = nameMap.values();
        else
            singletons.addAll(nameMap.values());
    }

    return new FragmentCollection<T>(singletons, overlapping);
}
 

import htsjdk.samtools.SAMRecord; //导入方法依赖的package包/类
public BAQCalculationResult calcBAQFromHMM(SAMRecord read, byte[] ref, int refOffset) {
    // todo -- need to handle the case where the cigar sum of lengths doesn't cover the whole read
    Pair<Integer, Integer> queryRange = calculateQueryRange(read);
    if ( queryRange == null ) return null; // read has Ns, or is completely clipped away

    int queryStart = queryRange.getFirst();
    int queryEnd = queryRange.getSecond();

    BAQCalculationResult baqResult = calcBAQFromHMM(ref, read.getReadBases(), read.getBaseQualities(), queryStart, queryEnd);

    // cap quals
    int readI = 0, refI = 0;
    for ( CigarElement elt : read.getCigar().getCigarElements() ) {
        int l = elt.getLength();
        switch (elt.getOperator()) {
            case N: // cannot handle these
                return null;
            case H : case P : // ignore pads and hard clips
                break;
            case S : refI += l; // move the reference too, in addition to I
            case I :
                // todo -- is it really the case that we want to treat I and S the same?
                for ( int i = readI; i < readI + l; i++ ) baqResult.bq[i] = baqResult.rawQuals[i];
                readI += l;
                break;
            case D : refI += l; break;
            case M :
                for (int i = readI; i < readI + l; i++) {
                    int expectedPos = refI - refOffset + (i - readI);
                    baqResult.bq[i] = capBaseByBAQ( baqResult.rawQuals[i], baqResult.bq[i], baqResult.state[i], expectedPos );
                }
                readI += l; refI += l;
                break;
            default:
                throw new ReviewedGATKException("BUG: Unexpected CIGAR element " + elt + " in read " + read.getReadName());
        }
    }
    if ( readI != read.getReadLength() ) // odd cigar string
        System.arraycopy(baqResult.rawQuals, 0, baqResult.bq, 0, baqResult.bq.length);

    return baqResult;
}