C++ SamRecord类代码示例

// Soft clip the record from the front and/or the back.
SamFilter::FilterStatus SamFilter::softClip(SamRecord& record, 
                                            int32_t numFrontClips,
                                            int32_t numBackClips)
{
    //////////////////////////////////////////////////////////
    Cigar* cigar = record.getCigarInfo();
    FilterStatus status = NONE;
    int32_t startPos = record.get0BasedPosition();
    CigarRoller updatedCigar;

    status = softClip(*cigar, numFrontClips, numBackClips, 
                      startPos, updatedCigar);

    if(status == FILTERED)
    {
        /////////////////////////////
        // The entire read is clipped, so rather than clipping it,
        // filter it out.
        filterRead(record);
        return(FILTERED);
    }
    else if(status == CLIPPED)
    {
        // Part of the read was clipped, and now that we have
        // an updated cigar, update the read.
        record.setCigar(updatedCigar);

        // Update the starting position.
        record.set0BasedPosition(startPos);
    }
    return(status);
}

bool SamFile::checkRecordInSection(SamRecord& record)
{
    bool recordFound = true;
    if(myRefID == BamIndex::REF_ID_ALL)
    {
        return(true);
    }
    // Check to see if it is in the correct reference/position.
    if(record.getReferenceID() != myRefID)
    {
        // Incorrect reference ID, return no more records.
        myStatus = SamStatus::NO_MORE_RECS;
        return(false);
    }
   
    // Found a record.
    recordFound = true;

    // If start/end position are set, verify that the alignment falls
    // within those.
    // If the alignment start is greater than the end of the region,
    // return NO_MORE_RECS.
    // Since myEndPos is Exclusive 0-based, anything >= myEndPos is outside
    // of the region.
    if((myEndPos != -1) && (record.get0BasedPosition() >= myEndPos))
    {
        myStatus = SamStatus::NO_MORE_RECS;
        return(false);
    }
        
    // We know the start is less than the end position, so the alignment
    // overlaps the region if the alignment end position is greater than the
    // start of the region.
    if((myStartPos != -1) && (record.get0BasedAlignmentEnd() < myStartPos))
    {
        // If it does not overlap the region, so go to the next
        // record...set recordFound back to false.
        recordFound = false;
    }

    if(!myOverlapSection)
    {
        // Needs to be fully contained.  Not fully contained if
        // 1) the record start position is < the region start position.
        // or
        // 2) the end position is specified and the record end position
        //    is greater than or equal to the region end position.
        //    (equal to since the region is exclusive.
        if((record.get0BasedPosition() < myStartPos) ||
           ((myEndPos != -1) && 
            (record.get0BasedAlignmentEnd() >= myEndPos)))
        {
            // This record is not fully contained, so move on to the next
            // record.
            recordFound = false;
        }
    }

    return(recordFound);
}

void Demux::match(){
	BamReader bamReader(bamFilePath);
	SamRecord samRecord;
	while (	bamReader.getNextRecord(samRecord)) {
		string recordName(samRecord.getReadName());
		recordName = cHandler.decrypt(recordName);
		//clean record name
		int len=recordName.find("$");
		recordName=recordName.substr(0,len);
		//clean ended
		if (!isDecoy(recordName)){
			string outputFile = generateFileName(recordName);
			printf("%s\n",outputFile.c_str());
			if (writers.find(outputFile) == writers.end()) //if the BamWriter is not initialized
				writers[outputFile] = new BamWriter(outputFile, bamReader.getHeader());
			BamWriter *writer = writers[outputFile];
			writer->writeRecord(samRecord);
		}
    }
    for(auto it=writers.begin();it!=writers.end();it++){
    	BamWriter *writer = it->second;
    	writer->close();
    	delete writer;
    }
}

void Bam2FastQ::handlePairedRN(SamRecord& samRec)
{
    static SamRecord* prevRec = NULL;
    static std::string prevRN = "";

    if(prevRec == NULL)
    {
        prevRec = &samRec;
    }
    else
    {
        if(strcmp(prevRec->getReadName(), samRec.getReadName()) != 0)
        {
            // Read Name does not match, error, did not find pair.
            std::cerr << "Paired Read, " << prevRec->getReadName()
                      << " but couldn't find mate, so writing as "
                      << "unpaired (single-ended)\n";
            ++myNumMateFailures;
            writeFastQ(*prevRec, myUnpairedFile,
                       myUnpairedFileNameExt);
            // Save this record to check against the next one.
            prevRec = &samRec;
        }
        else
        {
            // Matching ReadNames.
            // Found the mate.
            ++myNumPairs;
            // Check which is the first in the pair.
            if(SamFlag::isFirstFragment(samRec.getFlag()))
            {
                if(SamFlag::isFirstFragment(prevRec->getFlag()))
                {
                    std::cerr << "Both reads of " << samRec.getReadName()
                              << " are first fragment, so "
                              << "splitting one to be in the 2nd fastq.\n";
                }
                writeFastQ(samRec, myFirstFile, myFirstFileNameExt,
                           myFirstRNExt.c_str());
                writeFastQ(*prevRec, mySecondFile, mySecondFileNameExt,
                           mySecondRNExt.c_str());
            }
            else
            {
                if(!SamFlag::isFirstFragment(prevRec->getFlag()))
                {
                    std::cerr << "Neither read of " << samRec.getReadName()
                              << " are first fragment, so "
                              << "splitting one to be in the 2nd fastq.\n";
                }
                writeFastQ(*prevRec, myFirstFile, myFirstFileNameExt, myFirstRNExt.c_str());
                writeFastQ(samRec, mySecondFile, mySecondFileNameExt, mySecondRNExt.c_str());
            }
            // No previous record.
            prevRec = NULL;
        }
    }
}

// make a 64-bit genomic coordinate [24bit-chr][32bit-pos][8bit-orientation]
uint64_t getGenomicCoordinate(SamRecord& r) {
  // 64bit string consisting of
  // 24bit refID, 32bit pos, 8bit orientation
  if ( ( r.getReferenceID() < 0 ) || ( r.get0BasedPosition() < 0 ) ) {
    return UNMAPPED_GENOMIC_COORDINATE;
  }
  else {
    return ( ( static_cast<uint64_t>(r.getReferenceID()) << 40 ) | ( static_cast<uint64_t>(r.get0BasedPosition()) << 8 ) | static_cast<uint64_t>( r.getFlag() & 0x0010 ) );
  }
}

// Returns the number of bases in the passed in read that overlap the
// region that is currently set.
uint32_t SamFile::GetNumOverlaps(SamRecord& samRecord)
{
    if(myRefPtr != NULL)
    {
        samRecord.setReference(myRefPtr);
    }
    samRecord.setSequenceTranslation(myReadTranslation);

    // Get the overlaps in the sam record for the region currently set
    // for this file.
    return(samRecord.getNumOverlaps(myStartPos, myEndPos));
}

void BamInterface::readRecord(IFILE filePtr, SamFileHeader& header,
                              SamRecord& record,
                              SamStatus& samStatus)
{
    // TODO - need to validate there are @SQ lines in both sam/bam - MAYBE!

    // SetBufferFromFile will reset the record prior to reading a new one.
    if(record.setBufferFromFile(filePtr, header) != SamStatus::SUCCESS)
    {
        // Failed, so add the error message.
        samStatus.addError(record.getStatus());
    }
}

void ClipOverlap::cleanupMateMap(MateMapByCoord& mateMap,
                                 SamCoordOutput* outputBufferPtr,
                                 int32_t chrom, int32_t position)
{
    // Cleanup any reads in the mateMap whose mates are prior to the position
    // currently being processed in the file.  It means the mate was not found 
    // as expected.  Stop cleaning up once one is found that is not passed.
    uint64_t chromPos = 0;
    if((chrom != -1) && (position != -1))

    {
        chromPos = SamHelper::combineChromPos(chrom, position);
    }
    else
    {
        chrom = -1;
    }
    
    // Stop after the first read is found whose mate has not yet been reached.
    SamRecord* firstRec = mateMap.first();
    while(firstRec != NULL)
    {
        uint64_t firstMateChromPos = 
            SamHelper::combineChromPos(firstRec->getMateReferenceID(),
                                       firstRec->get0BasedMatePosition());
        if((firstMateChromPos < chromPos) || (chrom == -1))
        {
            // Already past the mate's position, so note this read and
            // write it.
            ++myNumMateFailures;
            if((outputBufferPtr != NULL) && !myOverlapsOnly)
            {
                outputBufferPtr->add(firstRec);
            }
            else
            {
                myPool.releaseRecord(firstRec);
            }
            // Remove this record.
            mateMap.popFirst();
            // Get the next record to check.
            firstRec = mateMap.first();
        }
        else
        {
            // The first record's mate position has not yet been passed, so
            // stop cleaning up the buffer.
            break;
        }
    }
}

// Finds the total base quality of a read
int Dedup_LowMem::getBaseQuality(SamRecord & record) {
    const char* baseQualities = record.getQuality();
    int readLength = record.getReadLength();
    int quality = 0.;
    if(strcmp(baseQualities, "*") == 0)
    {
        return(0);
    }
    for(int i=0; i < readLength; ++i) {
        int q = static_cast<int>(baseQualities[i])-33;
        if ( q >= myMinQual ) quality += q;
    }
    return quality;
}

void SamFilter::filterRead(SamRecord& record)
{
    // Filter the read by marking it as unmapped.
    uint16_t flag = record.getFlag(); 
    SamFlag::setUnmapped(flag);
    // Clear N/A flags.
    flag &= ~SamFlag::PROPER_PAIR;
    flag &= ~SamFlag::SECONDARY_ALIGNMENT;
    flag &= ~SamFlag::SUPPLEMENTARY_ALIGNMENT;
    record.setFlag(flag);
    // Clear Cigar
    record.setCigar("*");
    // Clear mapping quality
    record.setMapQuality(0);
}

// determine whether the record's position is different from the previous record
bool Dedup_LowMem::hasPositionChanged(SamRecord& record)
{
    if (lastReference != record.getReferenceID() ||
            lastCoordinate < record.get0BasedPosition())
    {
        if (lastReference != record.getReferenceID())
        {
            lastReference = record.getReferenceID();
            Logger::gLogger->writeLog("Reading ReferenceID %d\n", lastReference);
        }
        lastCoordinate = record.get0BasedPosition();
        return true;
    }
    return false;
}

void Bam2FastQ::writeFastQ(SamRecord& samRec, IFILE filePtr,
                             const char* readNameExt)
{
    static int16_t flag;
    static std::string sequence;
    static String quality;

    if(filePtr == NULL)
    {
        return;
    }

    flag = samRec.getFlag();
    const char* readName = samRec.getReadName();
    sequence = samRec.getSequence();
    quality = samRec.getQuality();
    
    if(SamFlag::isReverse(flag) && myReverseComp)
    {
        // It is reverse, so reverse compliment the sequence
        BaseUtilities::reverseComplement(sequence);
        // Reverse the quality.
        quality.Reverse();
    }
    else
    {
        // Ensure it is all capitalized.
        int seqLen = sequence.size();
        for (int i = 0; i < seqLen; i++)
        {
            sequence[i] = (char)toupper(sequence[i]);
        }
    }
    
    if(myRNPlus)
    {

        ifprintf(filePtr, "@%s%s\n%s\n+%s%s\n%s\n", readName, readNameExt,
                 sequence.c_str(), readName, readNameExt, quality.c_str());
    }
    else
    {
        ifprintf(filePtr, "@%s%s\n%s\n+\n%s\n", readName, readNameExt,
                 sequence.c_str(), quality.c_str());
    }
    // Release the record.
    myPool.releaseRecord(&samRec);
}

void testRead(SamFile &inSam)
{
    // Read the SAM Header.
    SamFileHeader samHeader;
    assert(inSam.ReadHeader(samHeader));

    validateHeader(samHeader);

    testCopyHeader(samHeader);    

    testModHeader(samHeader);

    SamRecord samRecord;
    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead1(samRecord);

    // Set a new quality and get the buffer.
    samRecord.setQuality("ABCDE");
    validateRead1ModQuality(samRecord);
    //   void* buffer = samRecord.getRecordBuffer();

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead2(samRecord);

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead3(samRecord);

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead4(samRecord);

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead5(samRecord);

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead6(samRecord);

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead7(samRecord);

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead8(samRecord);

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead9(samRecord);

    assert(inSam.ReadRecord(samHeader, samRecord) == true);
    validateRead10(samRecord);
}

// NOTE: Only positions where the reference and read both have bases that
//       are different and not 'N' are considered mismatches.
uint32_t SamFilter::sumMismatchQuality(SamRecord& record, 
                                       GenomeSequence& refSequence,
                                       uint8_t defaultQualityInt)
{
    // Track the mismatch info.
    int mismatchQual = 0;
    int numMismatch = 0;

    SamQuerySeqWithRefIter sequenceIter(record, refSequence);

    SamSingleBaseMatchInfo baseMatchInfo;
    while(sequenceIter.getNextMatchMismatch(baseMatchInfo))
    {
        if(baseMatchInfo.getType() == SamSingleBaseMatchInfo::MISMATCH)
        {
            // Got a mismatch, get the associated quality.
            char readQualityChar = 
                record.getQuality(baseMatchInfo.getQueryIndex());
            uint8_t readQualityInt = 
                BaseUtilities::getPhredBaseQuality(readQualityChar);
            
            if(readQualityInt == BaseUtilities::UNKNOWN_QUALITY_INT)
            {
                // Quality was not specified, so use the configured setting.
                readQualityInt = defaultQualityInt;
            }
            mismatchQual += readQualityInt;
            ++numMismatch;
        }
    }

    return(mismatchQual);
}

SamStatus::Status BamInterface::writeRecord(IFILE filePtr,
        SamFileHeader& header,
        SamRecord& record,
        SamRecord::SequenceTranslation translation)
{
    // Write the file, returning the status.
    return(record.writeRecordBuffer(filePtr, translation));
}