Java SAMRecord类代码示例

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
public int processRecordUnpaired(SAMRecord sr){
int totalOp = 0;
String hlagene = HLA.extractHLAGeneName(sr.getReferenceName());
HLAGraph hg = this.hlaName2Graph.get(hlagene);
//hg.traverse();
if(hg != null){
    if(hg.isClassI()){
	boolean qc = this.qcCheck(sr);
	if(!qc)
	    return 0;
    }
    totalOp += hg.addWeight(sr, HLA.readNum);
    HLA.readNum++;
}
return totalOp;
   }
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
private void parseTag(final SAMRecord samRecord, final String tag) {
    Map.Entry<String, Object> entry = null;
    try {
        entry = tagCodec.decode(tag);
    } catch (SAMFormatException e) {
        e.printStackTrace();
    }
    if (entry != null) {
        if (entry.getValue() instanceof TagValueAndUnsignedArrayFlag) {
            final TagValueAndUnsignedArrayFlag valueAndFlag =
                    (TagValueAndUnsignedArrayFlag) entry.getValue();
            if (valueAndFlag.isUnsignedArray) {
                samRecord.setUnsignedArrayAttribute(entry.getKey(),
                        valueAndFlag.value);
            } else {
                samRecord.setAttribute(entry.getKey(), valueAndFlag.value);
            }
        } else {
            samRecord.setAttribute(entry.getKey(), entry.getValue());
        }
    }
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
private static Integer calculateNucleotideFrequencies(final List<SAMRecord> blockReads,
                                                      final Map<String, Integer> atgcMap,
                                                      final NucleotideContent nucleotideContent) {
    for (SAMRecord read : blockReads) {
        char[] nucleotides = read.getReadString().toCharArray();
        if (nucleotides.length != 0) {
            nucleotideContent.incTotalNumberOfNucleotides(nucleotides.length);
            calculateCount(nucleotides, nucleotideContent);
        }
    }
    atgcMap.put("A", nucleotideContent.getACount());
    atgcMap.put("T", nucleotideContent.getTCount());
    atgcMap.put("G", nucleotideContent.getGCount());
    atgcMap.put("C", nucleotideContent.getCCount());
    return nucleotideContent.getTotalNumberOfNucleotides();
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
protected static Cigar reclipCigar(Cigar cigar, SAMRecord read) {
    ArrayList<CigarElement> elements = new ArrayList<CigarElement>();

    int i = 0;
    int n = read.getCigar().numCigarElements();
    while (i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()))
        elements.add(read.getCigar().getCigarElement(i++));

    elements.addAll(cigar.getCigarElements());

    i++;
    while (i < n && !isClipOperator(read.getCigar().getCigarElement(i).getOperator()))
        i++;

    while (i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()))
        elements.add(read.getCigar().getCigarElement(i++));

    return new Cigar(elements);
}
 

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包/类
private Handler<SAMRecord> filterReads(BamTrack<Read> bamTrack, BamQueryOption options, SamReader reader,
                                       String chromosomeName, boolean coverageOnly) throws IOException {
                                                                                        //int maxReadCount
    CloseableIterator<SAMRecord> iterator = reader.query(chromosomeName, bamTrack.getStartIndex(),
                                                               bamTrack.getEndIndex(), false);
    LOG.debug(getMessage(MessagesConstants.DEBUG_GET_ITERATOR_QUERY, iterator.toString()));

    iterator = setIteratorFiltering(iterator, options);

    final Handler<SAMRecord> filter = BamUtil.createSAMRecordHandler(bamTrack, options, referenceManager,
            coverageOnly); //maxReadCount

    while (iterator.hasNext()) {
        final SAMRecord samRecord = iterator.next();
        //if read unmapped
        filter.add(samRecord);
    }

    return filter;
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
private void getReads(final BamFile bamFile, BamTrack<Read> bamTrack, final BamQueryOption options) throws
        IOException {
    Chromosome chromosome = bamTrack.getChromosome();
    try (SamReader reader = makeSamReader(bamFile, Collections.singletonList(chromosome),
            chromosome.getReferenceId())) {
        LOG.debug(getMessage(MessagesConstants.DEBUG_FILE_OPENING, bamFile.getPath()));
        String chromosomeName = options.getChromosomeName();

        if (reader.getFileHeader().getSequence(chromosomeName) == null) {
            chromosomeName = Utils.changeChromosomeName(chromosomeName);
        }

        Handler<SAMRecord> filter = filterReads(bamTrack, options, reader, chromosomeName, options.getMode() ==
                BamTrackMode.COVERAGE);
        //maxReadsCount

        bamTrack.setBlocks(filter.getSifter().getReadListResult());
        bamTrack.setMinPosition(filter.getMinPosition());
        bamTrack.setReferenceBuffer(filter.getReferenceBuff());
        bamTrack.setDownsampleCoverage(filter.getSifter().getDownsampleCoverageResult());
        bamTrack.setBaseCoverage(filter.getBaseCoverage(bamTrack.getScaleFactor()));
        bamTrack.setSpliceJunctions(filter.getSpliceJunctions());
    }
}
 

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包/类
private CloseableIterator<SAMRecord> setIteratorFiltering(final CloseableIterator<SAMRecord> iterator,
                                                          final BamQueryOption options) {
    List<SamRecordFilter> filters = new ArrayList<>();
    if (options.isFilterDuplicate()) {
        filters.add(new DuplicateReadFilter());
    }
    if (options.isFilterNotPrimary()) {
        filters.add(new NotPrimaryAlignmentFilter());
    }
    if (options.isFilterVendorQualityFail()) {
        filters.add(new FailsVendorReadQualityFilter());
    }
    if (options.isFilterSupplementaryAlignment()) {
        filters.add(new SecondaryOrSupplementaryFilter());
    }

    if (!filters.isEmpty()) {
        return new FilteringSamIterator(iterator, new AggregateFilter(filters));
    }

    return iterator;
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
@Nullable
private Read getReadFromBamFile(ReadQuery query, Chromosome chromosome, BamFile bamFile) throws IOException {
    try (SamReader reader = bamHelper.makeSamReader(bamFile, Collections.singletonList(chromosome),
                                                    chromosome.getReferenceId())) {
        String chromosomeName = chromosome.getName();
        if (reader.getFileHeader().getSequence(chromosomeName) == null) {
            chromosomeName = Utils.changeChromosomeName(chromosomeName);
        }

        SAMRecordIterator iterator = reader.query(chromosomeName, query.getStartIndex(), query.getEndIndex(),
                                                  true);
        while (iterator.hasNext()) {
            final SAMRecord samRecord = iterator.next();
            if (samRecord.getReadName().equals(query.getName())) {
                BamQueryOption option = new BamQueryOption();
                option.setRefID(chromosome.getReferenceId());
                option.setChromosomeName(chromosome.getName());
                SAMRecordHandler recordHandler = new SAMRecordHandler(query.getStartIndex(), query.getEndIndex(),
                                                                      referenceManager, null, option);
                List<BasePosition> diffBase = recordHandler.computeDifferentBase(samRecord);
                return BamUtil.createExtendedRead(samRecord, diffBase);
            }
        }
    }
    return null;
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
/**
 * Reads must by added sorted by read's start coordinate. If a read intersects the frame,
 * it will be added to the right frame.
 * @param record representing a read
 * @param start of the read
 * @param end of the read
 * @param differentBase list of read bases that differ from the reference nucleotides
 * @param headStr soft clipped start of the read
 * @param tailStr soft clipped end of the read
 * @throws IOException
 */
@Override
public void add(final SAMRecord record, final int start, final int end,
                final List<BasePosition> differentBase, final String headStr, final String tailStr)
        throws IOException {
    if (statPosition < 0) {
        statPosition = start;
        border = statPosition + frame - 1;
        border = border > endTrack ? endTrack : border;
    }
    if (start > border) {
        flushList();
        refreshForNewFrame(start);
    }
    bufferResult.add(BamUtil.createReadFromRecord(record, start, end, differentBase, headStr, tailStr));
    totalReadsCount++;
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
/**
 * Create Read object from SAMRecord
 * @param samRecord a SAMRecord from BAM file
 * @param start new start of a record, without softClip
 * @param end new end of a record, without softClip
 * @param differentBase list of bases, where variation occurred
 * @param head head sequence of a Record
 * @param tail tail sequence of a Record
 * @return
 */
public static Read createReadFromRecord(final SAMRecord samRecord, final int start, final int end,
                                        final List<BasePosition> differentBase,
                                        final String head, final String tail) {
    final Read read = new Read();
    read.setStartIndex(start);
    read.setEndIndex(end);
    read.setName(samRecord.getReadName());
    read.setCigarString(samRecord.getCigarString());
    read.setStand(!samRecord.getReadNegativeStrandFlag());
    read.setMappingQuality(samRecord.getMappingQuality());
    read.setFlagMask(samRecord.getFlags());
    read.setTLen(samRecord.getInferredInsertSize());
    read.setRNext(samRecord.getMateReferenceName());
    read.setPNext(samRecord.getMateAlignmentStart());
    read.setPairedReadName(samRecord.getPairedReadName());
    read.setRName(samRecord.getReferenceName());
    read.setDifferentBase(differentBase);
    read.setHeadSequence(head);
    read.setTailSequence(tail);
    return read;
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
private String getXSTag(List<SAMRecord.SAMTagAndValue> list) {
    String tagValue = null;
    for (SAMRecord.SAMTagAndValue samTagAndValue : list) {
        if (samTagAndValue.tag.equalsIgnoreCase(XS_TAG)) {
            tagValue = samTagAndValue.value.toString();
        }
    }
    return tagValue;
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
private Map<Sequence, List<SAMRecord>> getReadsFromFile(final Chromosome chromosome, final Track<Sequence> track,
                                                        final BamFile bamFile, final List<Sequence> blocks)
        throws IOException {
    final Map<Sequence, List<SAMRecord>> records = new HashMap<>();
    try (SamReader reader = makeSamReader(bamFile, Collections.singletonList(chromosome),
            chromosome.getReferenceId())) {
        LOG.debug(getMessage(MessagesConstants.DEBUG_FILE_OPENING, bamFile.getPath()));
        String chromosomeName = chromosome.getName();
        final int startIndex = track.getStartIndex();
        final int endIndex = track.getEndIndex();

        if (reader.getFileHeader().getSequence(chromosomeName) == null) {
            chromosomeName = Utils.changeChromosomeName(chromosomeName);
        }
        final SAMRecordIterator iterator = reader.queryOverlapping(chromosomeName, startIndex, endIndex);

        LOG.debug(getMessage(MessagesConstants.DEBUG_GET_ITERATOR_QUERY, iterator.toString()));
        while (iterator.hasNext()) {
            final SAMRecord samRecord = iterator.next();
            //if read unmapped
            if (!samRecord.getSAMFlags().contains(SAMFlag.READ_UNMAPPED) && !samRecord.getCigar().isEmpty() &&
                    samRecord.getEnd() > samRecord.getStart()) {
                addRecordToBlock(blocks, records, samRecord);
            }
        }
    }
    return records;
}
 

import htsjdk.samtools.SAMRecord; //导入依赖的package包/类
/**
 * Judge if a read pass all the filters. It should not be added into a pileup if return false
 *
 * @param read read to check
 * @return if it shouldn't pass all the filter, false is returned
 */
public boolean filter(SAMRecord read) {
    for (Filter f : filters) {
        if (f.filterOut(read))
            return false;
    }
    return true;
}