C++ SamFileHeader类代码示例

void testAddHeaderAndTagToFile(const char* inputName, const char* outputName)
{
    SamFile inSam, outSam;
    assert(inSam.OpenForRead(inputName));
    assert(outSam.OpenForWrite(outputName));

    // Read the SAM Header.
    SamFileHeader samHeader;
    assert(inSam.ReadHeader(samHeader));

    // Add a header line.
    assert(samHeader.addHeaderLine("@RG\tID:myID\tSM:mySM") == false);
    assert(samHeader.addHeaderLine("@RG\tID:myID3\tSM:mySM") == true);

    // Write Header
    assert(outSam.WriteHeader(samHeader));

    SamRecord samRecord;
    assert(inSam.ReadRecord(samHeader, samRecord));
    //   validateRead1(samRecord);
    // Add two tags.
    assert(samRecord.addIntTag("XA", 123));
    assert(samRecord.addIntTag("XA", 456));
    assert(samRecord.addTag("RR", 'Z', "myID1"));
    assert(samRecord.addTag("RR", 'Z', "myID2"));

    // Write as Sam.
    assert(outSam.WriteRecord(samHeader, samRecord));

    // TODO, add test to verify it was written correctly.

    // Read a couple of records to make sure it properly can read them even
    // if they are bigger than the original.
    assert(inSam.ReadRecord(samHeader, samRecord));
    assert(inSam.ReadRecord(samHeader, samRecord));

    //  Check the MD tag, which requires the reference.
    GenomeSequence reference("testFiles/chr1_partial.fa");
    assert(SamTags::isMDTagCorrect(samRecord, reference) == false);
    String newMDTag;
    SamTags::createMDTag(newMDTag, samRecord, reference);
    assert(newMDTag == "2T1N0");
    assert(SamTags::updateMDTag(samRecord, reference));
    // Write as Sam.
    assert(outSam.WriteRecord(samHeader, samRecord));
}

// build the read group library map
void Dedup_LowMem::buildReadGroupLibraryMap(SamFileHeader& header) {
    rgidLibMap.clear();
    numLibraries = 0;
    std::map<std::string,uint32_t> libNameMap;

    SamHeaderRecord * headerRecord = header.getNextRGRecord();
    while(headerRecord != NULL) {
        std::string ID = headerRecord->getTagValue("ID");
        std::string LB = headerRecord->getTagValue("LB");

        if ( ID.empty() ) {
            std::string headerRecordString;
            headerRecord->appendString(headerRecordString);
            Logger::gLogger->error("Cannot find readGroup ID information in the header line %s",
                                   headerRecordString.c_str());
        }
        if ( rgidLibMap.find(ID) != rgidLibMap.end() ) {
            Logger::gLogger->error("The readGroup ID %s is not a unique identifier",ID.c_str());
        }

        if ( LB.empty() ) {
            std::string headerRecordString;
            headerRecord->appendString(headerRecordString);
            Logger::gLogger->warning("Cannot find library information in the header line %s. Using empty string for library name",
                                     headerRecordString.c_str());
        }

        if ( libNameMap.find( LB ) != libNameMap.end() ) {
            rgidLibMap[ID] = libNameMap[LB];
        }
        else {
            numLibraries = libNameMap.size()+1;
            libNameMap[LB] = numLibraries;
            rgidLibMap[ID] = numLibraries;
        }
        headerRecord = header.getNextRGRecord();
    }

    if (numLibraries > 0xff) {
        Logger::gLogger->error("More than 255 library names are identified. Dedup_LowMem currently only allows up to 255 library names");
    }
}

void testCopyHeader(SamFileHeader& samHeader)
{
    // Copy the header.
    SamFileHeader samHeader2;
    
    SamHeaderRecord* recPtr = samHeader.getNextHeaderRecord();
    while(recPtr != NULL)
    {
        samHeader2.addRecordCopy(*recPtr);
        recPtr = samHeader.getNextHeaderRecord();
    }
    // Add the comments.
    std::string nextComment = samHeader.getNextComment();
    while(nextComment != SamFileHeader::EMPTY_RETURN)
    {
        samHeader2.addComment(nextComment.c_str());
        nextComment = samHeader.getNextComment();
    }
    // Validate the header.
    validateHeader(samHeader2);
}

SamFile::SortedType SamFile::getSortOrderFromHeader(SamFileHeader& header)
{
    const char* tag = header.getSortOrder();
   
    // Default to unsorted since if it is not specified in the header
    // that is the value that should be used.
    SortedType headerSortOrder = UNSORTED;
    if(strcmp(tag, "queryname") == 0)
    {
        headerSortOrder = QUERY_NAME;
    }
    else if(strcmp(tag, "coordinate") == 0)
    {
        headerSortOrder = COORDINATE;
    }
    return(headerSortOrder);
}

// Get the number of unmapped reads in the specified reference id.  
// Returns -1 for out of range refIDs.
int32_t SamFile::getNumUnMappedReadsFromIndex(const char* refName,
                                              SamFileHeader& header)
{
    // The bam index must have already been read.
    if(myBamIndex == NULL)
    {
        myStatus.setStatus(SamStatus::FAIL_ORDER, 
                           "Cannot get num unmapped reads from the index until it has been read.");
        return(false);
    }
    int32_t refID = BamIndex::REF_ID_UNMAPPED;
    if((strcmp(refName, "") != 0) && (strcmp(refName, "*") != 0))
    {
        // Reference name specified, so read just the "-1" entries.
        refID =  header.getReferenceID(refName);
    }
    return(myBamIndex->getNumUnMappedReads(refID));
}

bool SamFile::processNewSection(SamFileHeader &header)
{
    myNewSection = false;

    // If there is no index file, return failure.
    if(myBamIndex == NULL)
    {
        // No bam index has been read.
        myStatus.setStatus(SamStatus::FAIL_ORDER, 
                           "Cannot read section since there is no index file open");
        throw(std::runtime_error("SOFTWARE BUG: trying to read a BAM record by section prior to opening the BAM Index file."));
        return(false);
    }

    // If there is not a BAM file open for reading, return failure.
    if(!myIsBamOpenForRead)
    {
        // There is not a BAM file open for reading.
        myStatus.setStatus(SamStatus::FAIL_ORDER, 
                           "Cannot read section since there is no bam file open");
        throw(std::runtime_error("SOFTWARE BUG: trying to read a BAM record by section without opening a BAM file."));
        return(false);
    }

    if(myHasHeader == false)
    {
        // The header has not yet been read.
        myStatus.setStatus(SamStatus::FAIL_ORDER, 
                           "Cannot read record since the header has not been read.");
        throw(std::runtime_error("SOFTWARE BUG: trying to read a BAM record by section prior to opening the header."));
        return(false);
    }

    // Indexed Bam open for read, so disable read buffering because iftell
    // will be used.
    // Needs to be done here after we already know that the header has been
    // read.
    myFilePtr->disableBuffering();

    myChunksToRead.clear();
    // Reset the end of the current chunk.  We are resetting our read, so
    // we no longer have a "current chunk" that we are reading.
    myCurrentChunkEnd = 0;

    // Check to see if the read section was set based on the reference name
    // but not yet converted to reference id.
    if(!myRefName.empty())
    {
        myRefID = header.getReferenceID(myRefName.c_str());
        // Clear the myRefName length so this code is only executed once.
        myRefName.clear();

        // Check to see if a reference id was found.
        if(myRefID == SamReferenceInfo::NO_REF_ID)
        {
            myStatus = SamStatus::NO_MORE_RECS;
            return(false);
        }
    }

    // Get the chunks associated with this reference region.
    if(myBamIndex->getChunksForRegion(myRefID, myStartPos, myEndPos, 
                                      myChunksToRead) == true)
    {
        myStatus = SamStatus::SUCCESS;
    }
    else
    {
        String errorMsg = "Failed to get the specified region, refID = ";
        errorMsg += myRefID;
        errorMsg += "; startPos = ";
        errorMsg += myStartPos;
        errorMsg += "; endPos = ";
        errorMsg += myEndPos;
        myStatus.setStatus(SamStatus::FAIL_PARSE, 
                           errorMsg);
    }
    return(true);
}

// Validate that the record is sorted compared to the previously read record
// if there is one, according to the specified sort order.
// If the sort order is UNSORTED, true is returned.
bool SamFile::validateSortOrder(SamRecord& record, SamFileHeader& header)
{
    if(myRefPtr != NULL)
    {
        record.setReference(myRefPtr);
    }
    record.setSequenceTranslation(myReadTranslation);

    bool status = false;
    if(mySortedType == UNSORTED)
    {
        // Unsorted, so nothing to validate, just return true.
        status = true;
    }
    else 
    {
        // Check to see if mySortedType is based on the header.
        if(mySortedType == FLAG)
        {
            // Determine the sorted type from what was read out of the header.
            mySortedType = getSortOrderFromHeader(header);
        }

        if(mySortedType == QUERY_NAME)
        {
            // Validate that it is sorted by query name.
            // Get the query name from the record.
            const char* readName = record.getReadName();

            // Check if it is sorted either in samtools way or picard's way.
            if((myPrevReadName.Compare(readName) > 0) &&
               (strcmp(myPrevReadName.c_str(), readName) > 0))
            {
                // return false.
                String errorMessage = "ERROR: File is not sorted by read name at record ";
                errorMessage += myRecordCount;
                errorMessage += "\n\tPrevious record was ";
                errorMessage += myPrevReadName;
                errorMessage += ", but this record is ";
                errorMessage += readName;
                myStatus.setStatus(SamStatus::INVALID_SORT, 
                                   errorMessage.c_str());
                status = false;
            }
            else
            {
                myPrevReadName = readName;
                status = true;
            }
        }
        else 
        {
            // Validate that it is sorted by COORDINATES.
            // Get the leftmost coordinate and the reference index.
            int32_t refID = record.getReferenceID();
            int32_t coord = record.get0BasedPosition();
            // The unmapped reference id is at the end of a sorted file.
            if(refID == BamIndex::REF_ID_UNMAPPED)
            {
                // A new reference ID that is for the unmapped reads
                // is always valid.
                status = true;
                myPrevRefID = refID;
                myPrevCoord = coord;
            }
            else if(myPrevRefID == BamIndex::REF_ID_UNMAPPED)
            {
                // Previous reference ID was for unmapped reads, but the
                // current one is not, so this is not sorted.
                String errorMessage = "ERROR: File is not coordinate sorted at record ";
                errorMessage += myRecordCount;
                errorMessage += "\n\tPrevious record was unmapped, but this record is ";
                errorMessage += header.getReferenceLabel(refID) + ":" + coord;
                myStatus.setStatus(SamStatus::INVALID_SORT, 
                                   errorMessage.c_str());
                status = false;
            }
            else if(refID < myPrevRefID)
            {
                // Current reference id is less than the previous one, 
                //meaning that it is not sorted.
                String errorMessage = "ERROR: File is not coordinate sorted at record ";
                errorMessage += myRecordCount;
                errorMessage += "\n\tPrevious record was ";
                errorMessage += header.getReferenceLabel(myPrevRefID) + ":" + myPrevCoord;
                errorMessage += ", but this record is ";
                errorMessage += header.getReferenceLabel(refID) + ":" + coord;
                myStatus.setStatus(SamStatus::INVALID_SORT, 
                                   errorMessage.c_str());
                status = false;
            }
            else
            {
                // The reference IDs are in the correct order.
                if(refID > myPrevRefID)
                {
                    // New reference id, so set the previous coordinate to -1
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    {
        // Just remove the extension from the input filename.
        int extStart = inFile.FastFindLastChar('.');
        if(extStart <= 0)
        {
            outBase = inFile;
        }
        else
        {
            outBase = inFile.Left(extStart);
        }
    }

    // Check to see if the first/second/single-ended were specified and
    // if not, set them.
    std::string firstExt = "_1.fastq";
    if(interleave)
    {
        firstExt = "_interleaved.fastq";
    }
    getFileName(firstOut, outBase, firstExt.c_str());
    getFileName(secondOut, outBase, "_2.fastq");
    getFileName(unpairedOut, outBase, ".fastq");

    if(params)
    {
        inputParameters.Status();
    }

    // Open the files for reading/writing.
    // Open prior to opening the output files,
    // so if there is an error, the outputs don't get created.
    SamFile samIn;
    SamFileHeader samHeader;
    samIn.OpenForRead(inFile, &samHeader);

    // Open the output files.
    myUnpairedFile = ifopen(unpairedOut, "w");

    // Only open the first file if it is different than an already opened file.
    if(firstOut != unpairedOut)
    {
        myFirstFile = ifopen(firstOut, "w");
    }
    else
    {
        myFirstFile = myUnpairedFile;
    }

    // If it is interleaved or the 2nd file is not a new name, set it appropriately.
    if(interleave || secondOut == firstOut)
    {
        mySecondFile = myFirstFile;
    }
    else if(secondOut == unpairedOut)
    {
        mySecondFile = myUnpairedFile;
    }
    else
    {
        mySecondFile = ifopen(secondOut, "w");
    }
    
    if(myUnpairedFile == NULL)
    {
        std::cerr << "Failed to open " << unpairedOut

// Read a BAM file's header.
bool BamInterface::readHeader(IFILE filePtr, SamFileHeader& header,
                              SamStatus& status)
{
    if(filePtr == NULL)
    {
        // File is not open, return false.
        status.setStatus(SamStatus::FAIL_ORDER,
                         "Cannot read header since the file pointer is null");
        return(false);
    }
    if(filePtr->isOpen() == false)
    {
        status.setStatus(SamStatus::FAIL_ORDER,
                         "Cannot read header since the file is not open");
        return(false);
    }

    // Clear the passed in header.
    header.resetHeader();

    int32_t headerLength;
    int readSize = ifread(filePtr, &headerLength, sizeof(headerLength));

    if(readSize != sizeof(headerLength))
    {
        String errMsg = "Failed to read the BAM header length, read ";
        errMsg += readSize;
        errMsg += " bytes instead of ";
        errMsg += (unsigned int)sizeof(headerLength);
        status.setStatus(SamStatus::FAIL_IO, errMsg.c_str());
        return(false);
    }

    String headerStr;
    if(headerLength > 0)
    {
        // Read the header.
        readSize =
            ifread(filePtr, headerStr.LockBuffer(headerLength + 1), headerLength);
        headerStr[headerLength] = 0;
        headerStr.UnlockBuffer();
        if(readSize != headerLength)
        {
            // Failed to read the header.
            status.setStatus(SamStatus::FAIL_IO, "Failed to read the BAM header.");
            return(false);
        }
    }

    // Parse the header that was read.
    if(!header.addHeader(headerStr))
    {
        // Status is set in the method on failure.
        status.setStatus(SamStatus::FAIL_PARSE, header.getErrorMessage());
        return(false);
    }

    int referenceCount;
    // Read the number of references sequences.
    ifread(filePtr, &referenceCount, sizeof(int));

    // Get and clear the reference info so it can be set
    // from the bam reference table.
    SamReferenceInfo& refInfo =
        header.getReferenceInfoForBamInterface();
    refInfo.clear();

    CharBuffer refName;

    // Read each reference sequence
    for (int i = 0; i < referenceCount; i++)
    {
        int nameLength;
        int rc;
        // Read the length of the reference name.
        rc = ifread(filePtr, &nameLength, sizeof(int));
        if(rc != sizeof(int))
        {
            status.setStatus(SamStatus::FAIL_IO,
                             "Failed to read the BAM reference dictionary.");
            return(false);
        }

        // Read the name.
        refName.readFromFile(filePtr, nameLength);

        // Read the length of the reference sequence.
        int32_t refLen;
        rc = ifread(filePtr, &refLen, sizeof(int));

        if(rc != sizeof(int)) {
            status.setStatus(SamStatus::FAIL_IO,
                             "Failed to read the BAM reference dictionary.");
            return(false);
        }

        refInfo.add(refName.c_str(), refLen);
    }

//.........这里部分代码省略.........

void parseOutRG(SamFileHeader& header, std::string& noRgPgString, SamFileHeader* newHeader, bool ignorePI)
{
    noRgPgString.clear();
    // strings for comparing if two RGs with same ID are the same.
    static std::string prevString = "";
    static std::string newString = "";

    SamHeaderRecord* rec = header.getNextHeaderRecord();
    while(rec != NULL)
    {
        if(rec->getType() == SamHeaderRecord::RG)
        {
            if(newHeader != NULL)
            {
                // This is an RG line.
                // First check if this RG is already included in the new header.
                SamHeaderRG* prevRG = newHeader->getRG(rec->getTagValue("ID"));
                
                if(prevRG != NULL)
                {
                    // This RG already exists, check that they are the same.
                    // If they are the same, there is nothing to do.
                    bool status = true;
                    prevString.clear();
                    newString.clear();
                    status &= prevRG->appendString(prevString);
                    status &= rec->appendString(newString);

                    if(prevString != newString)
                    {
                        if(!ignorePI)
                        {
                            Logger::gLogger->error("Failed to add readgroup to "
                                                   "header, duplicate, but "
                                                   "non-identical RG ID, %s\n"
                                                   "prev:\t(%s)\nnew:\t(%s)",
                                                   rec->getTagValue("ID"),
                                                   prevString.c_str(),
                                                   newString.c_str());
                        }
                        else
                        {
                            // Check for a PI string.
                            size_t prevPIStart = prevString.find("PI:");
                            size_t newPIStart = newString.find("PI:");

                            // If they are both npos, then PI was not found
                            // so fail.
                            if((prevPIStart == std::string::npos) &&
                               (newPIStart == std::string::npos))
                            {
                                // They are not identical, so report an error.
                                Logger::gLogger->error("Failed to add readgroup"
                                                       " to header, duplicate,"
                                                       " but non-identical RG"
                                                       " ID, %s\n"
                                                       "prev:\t(%s)\nnew:\t(%s)",
                                                       rec->getTagValue("ID"),
                                                       prevString.c_str(),
                                                       newString.c_str());
                            }
                            else
                            {
                                // PI found in one or both strings.
                                size_t prevPIEnd;
                                size_t newPIEnd;
                                if(prevPIStart == std::string::npos)
                                {
                                    // new string has PI, so compare to the start of that.
                                    prevPIStart = newPIStart;
                                    prevPIEnd = newPIStart;
                                }
                                else
                                {
                                    prevPIEnd = prevString.find('\t', prevPIStart) + 1;
                                }
                                if(newPIStart == std::string::npos)
                                {
                                    // new string has PI, so compare to the start of that.
                                    newPIStart = prevPIStart;
                                    newPIEnd = newPIStart;
                                }
                                else
                                {
                                    newPIEnd = newString.find('\t', newPIStart) + 1;
                                }
                                // Compare before PI.
                                if((newString.compare(0, newPIStart, prevString, 0, prevPIStart) != 0) ||
                                   (newString.compare(newPIEnd, std::string::npos, prevString,
                                                      prevPIEnd, std::string::npos) != 0))
                                {
                                    // They are not identical, so report an error.
                                    Logger::gLogger->error("Failed to add readgroup to header, "
                                                           "duplicate, but non-identical RG ID, %s, "
                                                           "even when ignoring PI\n"
                                                           "prev:\t(%s)\nnew:\t(%s)",
                                                           rec->getTagValue("ID"),
                                                           prevString.c_str(),
                                                           newString.c_str());
                                }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        PileupElementBaseQCStats::setPercentStats(false);
    }

    if(baseQCPtr != NULL)
    {
        PileupElementBaseQCStats::setOutputFile(baseQCPtr);
        PileupElementBaseQCStats::printHeader();
    }
    if((baseQCPtr != NULL) || baseSum)
    {
        PileupElementBaseQCStats::setMapQualFilter(minMapQual);
        PileupElementBaseQCStats::setBaseSum(baseSum);
    }

    if(params)
    {
        inputParameters.Status();
    }

    // Open the file for reading.
    SamFile samIn;
    if(!samIn.OpenForRead(inFile))
    {
        fprintf(stderr, "%s\n", samIn.GetStatusMessage());
        return(samIn.GetStatus());
    }

    samIn.SetReadFlags(requiredFlags, excludeFlags);

    // Set whether or not basic statistics should be generated.
    samIn.GenerateStatistics(basic);

    // Read the sam header.
    SamFileHeader samHeader;
    if(!samIn.ReadHeader(samHeader))
    {
        fprintf(stderr, "%s\n", samIn.GetStatusMessage());
        return(samIn.GetStatus());
    }

    // Open the bam index file for reading if we are
    // doing unmapped reads (also set the read section).
    if(useIndex)
    {
        samIn.ReadBamIndex(indexFile);

        if(unmapped)
        {
            samIn.SetReadSection(-1);
        }

        if(!regionList.IsEmpty())
        {
            myRegionList = ifopen(regionList, "r");
        }
    }

    //////////////////////////
    // Read dbsnp if specified and doing baseQC
    if(((baseQCPtr != NULL) || baseSum) && (!dbsnp.IsEmpty()))
    {
        // Read the dbsnp file.
        IFILE fdbSnp;
        fdbSnp = ifopen(dbsnp,"r");
        // Determine how many entries.
        const SamReferenceInfo& refInfo = samHeader.getReferenceInfo();

void testModHeader(SamFileHeader& samHeader)
{
    // Check the header line.
    std::string headerString = "";
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:11\tLN:134452384\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\tLB:library2\[email protected]\tComment 1\[email protected]\tComment 2\n");

    // Remove a tag - by setting it to "".
    assert(samHeader.setRGTag("LB", "", "myID2") == true);


    // Check the header line.
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:11\tLN:134452384\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\[email protected]\tComment 1\[email protected]\tComment 2\n");

    //  Add an HD tag.
    SamHeaderHD* hd = new SamHeaderHD();
    assert(hd->setTag("VN", "1.3") == true);
    assert(samHeader.addHD(hd) == true);
    assert(strcmp(samHeader.getHDTagValue("VN"), "1.3") == 0);
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:11\tLN:134452384\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\[email protected]\tVN:1.3\[email protected]\tComment 1\[email protected]\tComment 2\n");

    // Try adding another HD tag.
    SamHeaderHD* hd2 = new SamHeaderHD();
    assert(hd2->setTag("VN", "1.4") == true);
    assert(samHeader.addHD(hd2) == false);
    assert(strcmp(samHeader.getHDTagValue("VN"), "1.4") != 0);
    assert(strcmp(samHeader.getHDTagValue("VN"), "1.3") == 0);
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:11\tLN:134452384\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\[email protected]\tVN:1.3\[email protected]\tComment 1\[email protected]\tComment 2\n");

    // Remove the entire HD Tag.
    assert(samHeader.removeHD() == true);
    assert(strcmp(samHeader.getHDTagValue("VN"), "") == 0);
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:11\tLN:134452384\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\[email protected]\tComment 1\[email protected]\tComment 2\n");

    // Remove an entire SQ Tag.
    assert(strcmp(samHeader.getSQTagValue("LN", "11"), "134452384") == 0);
    assert(samHeader.removeSQ("11") == true);
    assert(strcmp(samHeader.getSQTagValue("LN", "11"), "") == 0);
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\[email protected]\tComment 1\[email protected]\tComment 2\n");

    // Try adding a null HD tag.
    hd = NULL;
    assert(samHeader.addHD(hd) == false);
    assert(strcmp(samHeader.getHDTagValue("VN"), "") == 0);
    assert(strcmp(samHeader.getHDTagValue("VN"), "1.4") != 0);
    assert(strcmp(samHeader.getHDTagValue("VN"), "1.3") != 0);
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\[email protected]\tComment 1\[email protected]\tComment 2\n");

    // Try adding a null SQ tag.
    SamHeaderSQ* sq = NULL;
    assert(samHeader.addSQ(sq) == false);
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\[email protected]\tComment 1\[email protected]\tComment 2\n");

    // Try adding an HD tag again.
    assert(samHeader.addHD(hd2) == true);
    assert(strcmp(samHeader.getHDTagValue("VN"), "1.4") == 0);
    assert(strcmp(samHeader.getHDTagValue("VN"), "1.3") != 0);
    assert(samHeader.getHeaderString(headerString) == true);
    assert(headerString == "@SQ\tSN:1\tLN:247249719\[email protected]\tSN:2\tLN:242951149\[email protected]\tSN:3\tLN:199501827\[email protected]\tSN:4\tLN:191273063\[email protected]\tSN:5\tLN:180857866\[email protected]\tSN:6\tLN:170899992\[email protected]\tSN:7\tLN:158821424\[email protected]\tSN:8\tLN:146274826\[email protected]\tSN:9\tLN:140273252\[email protected]\tSN:10\tLN:135374737\[email protected]\tSN:12\tLN:132349534\[email protected]\tSN:13\tLN:114142980\[email protected]\tSN:14\tLN:106368585\[email protected]\tSN:15\tLN:100338915\[email protected]\tSN:16\tLN:88827254\[email protected]\tSN:17\tLN:78774742\[email protected]\tSN:18\tLN:76117153\[email protected]\tSN:19\tLN:63811651\[email protected]\tSN:20\tLN:62435964\[email protected]\tSN:21\tLN:46944323\[email protected]\tSN:22\tLN:49691432\[email protected]\tSN:X\tLN:154913754\[email protected]\tID:myID\tLB:library\tSM:sample\[email protected]\tID:myID2\tSM:sample2\[email protected]\tVN:1.4\[email protected]\tComment 1\[email protected]\tComment 2\n");


    // TODO Get the comments.

}

// add readgroup header line to the SamFileHeader
void addReadGroupToHeader(SamFileHeader& header, ReadGroup& rg) {
  if ( !header.addHeaderLine(rg.s_header_line.c_str()) ) {
    Logger::gLogger->error("Failed to add ID = %s, header line %s",rg.s_id.c_str(),rg.s_header_line.c_str());
  }
}

void parseOutRG(SamFileHeader& header, std::string& noRgPgString, SamFileHeader* newHeader)
{
    noRgPgString.clear();
    // strings for comparing if two RGs with same ID are the same.
    static std::string prevString = "";
    static std::string newString = "";

    SamHeaderRecord* rec = header.getNextHeaderRecord();
    while(rec != NULL)
    {
        if(rec->getType() == SamHeaderRecord::RG)
        {
            if(newHeader != NULL)
            {
                // This is an RG line.
                // First check if this RG is already included in the new header.
                SamHeaderRG* prevRG = newHeader->getRG(rec->getTagValue("ID"));
                
                if(prevRG != NULL)
                {
                    // This RG already exists, check that they are the same.
                    // If they are the same, there is nothing to do.
                    bool status = true;
                    prevString.clear();
                    newString.clear();
                    status &= prevRG->appendString(prevString);
                    status &= rec->appendString(newString);
                    if(prevString != newString)
                    {
                        // They are not identical, so report an error.
                        Logger::gLogger->error("Failed to add readgroup to header, "
                                               "duplicate, but non-identical RG ID, %s",
                                               rec->getTagValue("ID"));
                    }
                }
                else
                {
                    // This RG does not exist yet, so add it to the new header.
                    if(!newHeader->addRecordCopy((SamHeaderRG&)(*rec)))
                    {
                        // Failed to add the RG, exit.
                        Logger::gLogger->error("Failed to add readgroup to header, %s",
                                               newHeader->getErrorMessage());
                    }
                }
            }
        }
        else if(rec->getType() == SamHeaderRecord::PG)
        {
            if(newHeader != NULL)
            {
                // This is a PG line.
                // First check if this PG is already included in the new header.
                SamHeaderPG* prevPG = newHeader->getPG(rec->getTagValue("ID"));
                
                if(prevPG != NULL)
                {
                    // This PG already exists, check if they are the same.
                    // If they are the same, there is nothing to do.
                    bool status = true;
                    prevString.clear();
                    newString.clear();
                    status &= prevPG->appendString(prevString);
                    status &= rec->appendString(newString);
                    if(prevString != newString)
                    {
                        // They are not identical, ignore for now.
                        // TODO: change the ID, and add it.
                        Logger::gLogger->warning("Warning: dropping duplicate, "
                                                 "but non-identical PG ID, %s",
                                                 rec->getTagValue("ID"));
                    }
                }
                else
                {
                    // This PG does not exist yet, so add it to the new header.
                    if(!newHeader->addRecordCopy((SamHeaderPG&)(*rec)))
                    {
                        // Failed to add the PG, exit.
                        Logger::gLogger->error("Failed to add PG to header, %s",
                                               newHeader->getErrorMessage());
                    }
                }
            }
        }
        else
        {
            rec->appendString(noRgPgString);
        }
        rec = header.getNextHeaderRecord();
    }

    // Append the comments.
    header.appendCommentLines(noRgPgString);
}

//.........这里部分代码省略.........
      }
      Logger::gLogger->writeLog("Input list file : %s", bamList.c_str());
  }
  Logger::gLogger->writeLog("Output BAM file : %s",s_out.c_str());
  Logger::gLogger->writeLog("Output log file : %s",s_logger.c_str());
  Logger::gLogger->writeLog("Verbose mode    : %s",b_verbose ? "On" : "Off");
  
  vector<ReadGroup> v_readgroups;      // readGroups corresponding to BAM file
  vector<ReadGroup> v_uniq_readgroups; // unique readGroups written to header

  // If the list file is being used instead of the individual bams, parse it.
  if(!s_list.empty())
  {
      // parse the list file and fill the vectors above
      if ( parseListFile(s_list, vs_in_bam_files, v_readgroups, v_uniq_readgroups) == false ) {
          Logger::gLogger->error("Error in parsing the list file %s",s_list.c_str());
      }
      if ( vs_in_bam_files.size() != v_readgroups.size() ) {
          Logger::gLogger->error("parseListFile gave different size for vs_in_bam_files, v_readgroups: %d, %d", vs_in_bam_files.size(), v_readgroups.size());
      }
  }

  // sanity check
  uint32_t n_bams = vs_in_bam_files.size();
  Logger::gLogger->writeLog("Total of %d BAM files are being merged",n_bams);

  if ( n_bams < 2 )
  {
      Logger::gLogger->error("At least two BAM files must be specified for merging");
  }

  // create SamFile and SamFileHeader object for each BAM file
  SamFile *p_in_bams = new SamFile[n_bams];
  SamFileHeader *p_headers = new SamFileHeader[n_bams];

  // read each BAM file and its header, 
  // making sure that the headers are identical

  std::string firstHeaderNoRGPG = "";
  std::string headerNoRGPG = "";
  SamFileHeader newHeader;

  std::string firstHeaderString = "";
  for(uint32_t i=0; i < n_bams; ++i)
  {
      if ( ! p_in_bams[i].OpenForRead(vs_in_bam_files[i].c_str()) )
      {
          Logger::gLogger->error("Cannot open BAM file %s for reading",vs_in_bam_files[i].c_str());
      }
      p_in_bams[i].setSortedValidation(SamFile::COORDINATE);
      
      p_in_bams[i].ReadHeader(p_headers[i]);

      // Extract the RGs from this header.
      if(i == 0)
      {
          // First header, so store it as the first header
          newHeader = p_headers[i];
          // Determine the header without RG.
          parseOutRG(p_headers[i], firstHeaderNoRGPG, NULL);
      }
      else
      {
          parseOutRG(p_headers[i], headerNoRGPG, &newHeader);
          if(firstHeaderNoRGPG != headerNoRGPG)
          {