C++ BamReader类代码示例

// old mate-finding code, for position-sorted BAM files, didn't work well
BamAlignment 
lookForMate(BamReader& rdr, BamAlignment& al, RefVector& refs)
{
    if (! rdr.Jump(al.MateRefID, al.MatePosition)) {
        cout << "*** Could not jump to " << al.MateRefID << ":" << al.MatePosition << endl;
        exit(1);
    }
    BamAlignment al_jump;
    while (rdr.GetNextAlignment(al_jump)) {
        if (al_jump.Name == al.Name 
            && al_jump.RefID == al.MateRefID
            && al_jump.Position == al.MatePosition) {
         cout << "MATE FOUND" << endl;
            printAlignmentInfo(al_jump, refs);
            break;
        } else if (al_jump.Position > al.MatePosition) {
         cout << "NO MATE FOUND, beyond MatePosition" << endl;
            break;
        }
    }
    if (! rdr.Jump(al.RefID, al.Position)) {
        cout << "*** Could not return to " << al.RefID << ":" << al.Position << endl;
        exit(1);
    }
    // need to return to the read we were on after the jump, but how?
    rdr.GetNextAlignment(al);
    return(al_jump);
}

bool BamParse::loadIndex(BamReader& reader)
{
	bool bsignal = reader.OpenIndex(this->filename_index);
	if(bsignal==false) 
	{
		cerr << "Index file not found, now create it!!!"<<endl;
		bool bci=reader.CreateIndex();
		if(bci==false)
		{
			cerr << "Index file cannot be created!!!"<<endl;
			return false;
		}
		else
		{
			bool blct=reader.LocateIndex();
			if(blct==false)
			{
				cerr << "Index file cannot be located!!!"<<endl;
				return false;
			}
		}

	}
	return true;
}

/*
Description:
	Load all the bam into memory at one time if no parameters set, otherwise load the needed part of the bam.
	Save the parsed info into vector. 
*/
bool BamParse::parseAlignment(int chrom1, int chrom1_begin, int chrom2, int chrom2_end)
{
	BamReader reader;
    if ( !reader.Open(filename) ) {
        cerr << "Bamtools ERROR: could not open input BAM file: " << filename << endl;
        return false;
    }
		
	//check whether need to set a region.
	if(chrom1>-1 && chrom1_begin>-1 && chrom2>-1 && chrom2_end>-1)
	{
		this->loadIndex(reader);
		BamRegion br(chrom1,chrom1_begin,chrom2,chrom2_end);
		bool is_set=reader.SetRegion(br);
		if(is_set==false)
		{
			return false;//cannot set the region.
		}
	}

	//process input data
    BamAlignment al;   
	while ( reader.GetNextAlignment(al) )
	{
		if(al.Position<0) continue;

		BamAlignmentRecord* bar=new BamAlignmentRecord();
		setAlignmentRecord(al,bar);
		bam_aln_records.push_back(bar);
	}

	reader.Close();
	return true;
}

int main (int argc, char *argv[]) {
    
     string bamfiletopen = string(argv[1]);
     BamReader reader;

     if ( !reader.Open(bamfiletopen) ) {
    	cerr << "Could not open input BAM files." << endl;
    	return 1;
     }

    BamAlignment al;
    while ( reader.GetNextAlignment(al) ) {
	string reconstructedReference = reconstructRef(&al);
	cout<<al.QueryBases<<endl;
	cout<<reconstructedReference<<endl;

	
	pair< string, vector<int> >  reconP = reconstructRefWithPos(&al);
	for(unsigned int i=0;i<reconP.first.size();i++){
	    cout<<reconP.first[i]<<"\t"<<reconP.second[i]<<endl;
	}

    }
    reader.Close();

    return 0;
}

// returns BAM file pointers to beginning of alignment data
bool BamMultiReaderPrivate::RewindReaders()
{

    m_errorString.clear();
    bool errorsEncountered = false;

    // iterate over readers
    std::vector<MergeItem>::iterator readerIter = m_readers.begin();
    std::vector<MergeItem>::iterator readerEnd = m_readers.end();
    for (; readerIter != readerEnd; ++readerIter) {
        MergeItem& item = (*readerIter);
        BamReader* reader = item.Reader;
        if (reader == 0) continue;

        // attempt rewind on BamReader
        if (!reader->Rewind()) {
            m_errorString.append(1, '\t');
            m_errorString.append(reader->GetErrorString());
            m_errorString.append(1, '\n');
            errorsEncountered = true;
        }
    }

    return !errorsEncountered;
}

// get next alignment among all files without parsing character data from alignments
bool BamMultiReader::GetNextAlignmentCore(BamAlignment& nextAlignment) {

    // bail out if we are at EOF in all files, means no more alignments to process
    if (!HasOpenReaders())
        return false;

    // when all alignments have stepped into a new target sequence, update our
    // current reference sequence id
    UpdateReferenceID();

    // our lowest alignment and reader will be at the front of our alignment index
    BamAlignment* alignment = alignments.begin()->second.second;
    BamReader* reader = alignments.begin()->second.first;

    // now that we have the lowest alignment in the set, save it by copy to our argument
    nextAlignment = BamAlignment(*alignment);
    //memcpy(&nextAlignment, alignment, sizeof(BamAlignment));

    // remove this alignment index entry from our alignment index
    alignments.erase(alignments.begin());

    // and add another entry if we can get another alignment from the reader
    if (reader->GetNextAlignmentCore(*alignment)) {
        alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position), 
                                    make_pair(reader, alignment)));
    } else { // do nothing
        //cerr << "reached end of file " << lowestReader->GetFilename() << endl;
    }

    return true;

}

void BamToFastq::SingleFastq() {
    // open the 1st fastq file for writing
    ofstream fq(_fastq1.c_str(), ios::out);
    if ( !fq ) {
        cerr << "Error: The first fastq file (" << _fastq1 << ") could not be opened.  Exiting!" << endl;
        exit (1);
    }
    // open the BAM file
    BamReader reader;
    reader.Open(_bamFile);
    BamAlignment bam;
    while (reader.GetNextAlignment(bam)) {
        // extract the sequence and qualities for the BAM "query"
        string seq  = bam.QueryBases;
        string qual = bam.Qualities;
        if (bam.IsReverseStrand() == true) {
            reverseComplement(seq);
            reverseSequence(qual);
        }
        fq << "@" << bam.Name << endl;
        fq << seq << endl;
        fq << "+" << endl;
        fq << qual << endl;
    }
}

// creates index files for BAM files that don't have them
bool BamMultiReaderPrivate::CreateIndexes(const BamIndex::IndexType& type) {

    bool errorsEncountered = false;
    m_errorString.clear();

    // iterate over readers
    vector<MergeItem>::iterator itemIter = m_readers.begin();
    vector<MergeItem>::iterator itemEnd  = m_readers.end();
    for ( ; itemIter != itemEnd; ++itemIter ) {
        MergeItem& item = (*itemIter);
        BamReader* reader = item.Reader;
        if ( reader == 0 ) continue;

        // if reader doesn't have an index, create one
        if ( !reader->HasIndex() ) {
            if ( !reader->CreateIndex(type) ) {
                m_errorString.append(1, '\t');
                m_errorString.append(reader->GetErrorString());
                m_errorString.append(1, '\n');
                errorsEncountered = true;
            }
        }
    }

    // check for errors encountered before returning success/fail
    if ( errorsEncountered ) {
        const string currentError = m_errorString;
        const string message = string("error while creating index files: ") + "\n" + currentError;
        SetErrorString("BamMultiReader::CreateIndexes", message);
        return false;
    } else
        return true;
}

// close the BAM files
void BamMultiReader::Close(void) {
  
    // close all BAM readers and clean up pointers
    vector<pair<BamReader*, BamAlignment*> >::iterator readerIter = readers.begin();
    vector<pair<BamReader*, BamAlignment*> >::iterator readerEnd  = readers.end();
    for ( ; readerIter != readerEnd; ++readerIter) {
      
        BamReader* reader = (*readerIter).first;
        BamAlignment* alignment = (*readerIter).second;
        
        // close the reader
        if ( reader) reader->Close();  
        
        // delete reader pointer
        delete reader;
        reader = 0;

        // delete alignment pointer
        delete alignment;
        alignment = 0;
    }

    // clear out the container
    readers.clear();
}

// locate (& load) index files for BAM readers that don't already have one loaded
bool BamMultiReaderPrivate::LocateIndexes(const BamIndex::IndexType& preferredType)
{

    bool errorsEncountered = false;
    m_errorString.clear();

    // iterate over readers
    std::vector<MergeItem>::iterator readerIter = m_readers.begin();
    std::vector<MergeItem>::iterator readerEnd = m_readers.end();
    for (; readerIter != readerEnd; ++readerIter) {
        MergeItem& item = (*readerIter);
        BamReader* reader = item.Reader;
        if (reader == 0) continue;

        // if reader has no index, try to locate one
        if (!reader->HasIndex()) {
            if (!reader->LocateIndex(preferredType)) {
                m_errorString.append(1, '\t');
                m_errorString.append(reader->GetErrorString());
                m_errorString.append(1, '\n');
                errorsEncountered = true;
            }
        }
    }

    // check for errors encountered before returning success/fail
    if (errorsEncountered) {
        const std::string currentError = m_errorString;
        const std::string message =
            std::string("error while locating index files: \n") + currentError;
        SetErrorString("BamMultiReader::LocatingIndexes", message);
        return false;
    } else
        return true;
}

bool BamMultiReaderPrivate::PopNextCachedAlignment(BamAlignment& al, const bool needCharData)
{

    // skip if no alignments available
    if (m_alignmentCache == 0 || m_alignmentCache->IsEmpty()) return false;

    // pop next merge item entry from cache
    MergeItem item = m_alignmentCache->TakeFirst();
    BamReader* reader = item.Reader;
    BamAlignment* alignment = item.Alignment;
    if (reader == 0 || alignment == 0) return false;

    // set char data if requested
    if (needCharData) {
        alignment->BuildCharData();
        alignment->Filename = reader->GetFilename();
    }

    // store cached alignment into destination parameter (by copy)
    al = *alignment;

    // load next alignment from reader & store in cache
    SaveNextAlignment(reader, alignment);
    return true;
}

// opens BAM files
bool BamMultiReaderPrivate::Open(const vector<string>& filenames) {

    m_errorString.clear();

    // put all current readers back at beginning (refreshes alignment cache)
    if ( !Rewind() ) {
        const string currentError = m_errorString;
        const string message = string("unable to rewind existing readers: \n\t") + currentError;
        SetErrorString("BamMultiReader::Open", message);
        return false;
    }

    // iterate over filenames
    bool errorsEncountered = false;
    vector<string>::const_iterator filenameIter = filenames.begin();
    vector<string>::const_iterator filenameEnd  = filenames.end();
    for ( ; filenameIter != filenameEnd; ++filenameIter ) {
        const string& filename = (*filenameIter);
        if ( filename.empty() ) continue;

        // attempt to open BamReader
        BamReader* reader = new BamReader;
        const bool readerOpened = reader->Open(filename);

        // if opened OK, store it
        if ( readerOpened )
            m_readers.push_back( MergeItem(reader, new BamAlignment) );

        // otherwise store error & clean up invalid reader
        else {
            m_errorString.append(1, '\t');
            m_errorString += string("unable to open file: ") + filename;
            m_errorString.append(1, '\n');
            errorsEncountered = true;

            delete reader;
            reader = 0;
        }
    }

    // check for errors while opening
    if ( errorsEncountered ) {
        const string currentError = m_errorString;
        const string message = string("unable to open all files: \t\n") + currentError;
        SetErrorString("BamMultiReader::Open", message);
        return false;
    }

    // check for BAM file consistency
    if ( !ValidateReaders() ) {
        const string currentError = m_errorString;
        const string message = string("unable to open inconsistent files: \t\n") + currentError;
        SetErrorString("BamMultiReader::Open", message);
        return false;
    }

    // update alignment cache
    return UpdateAlignmentCache();
}

// saves index data to BAM index files (".bai"/".bti") where necessary, returns success/fail
bool BamMultiReader::CreateIndexes(bool useStandardIndex) {
    bool result = true;
    for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
        BamReader* reader = it->first;
        result &= reader->CreateIndex(useStandardIndex);
    }
    return result;
}

// returns BAM file pointers to beginning of alignment data
bool BamMultiReader::Rewind(void) { 
    bool result = true;
    for (vector<pair<BamReader*, BamAlignment*> >::iterator it = readers.begin(); it != readers.end(); ++it) {
        BamReader* reader = it->first;
        result &= reader->Rewind();
    }
    return result;
}

/**
 * Main work method.  Reads the BAM file once and collects sorted information about
 * the 5' ends of both ends of each read (or just one end in the case of pairs).
 * Then makes a pass through those determining duplicates before re-reading the
 * input file and writing it out with duplication flags set correctly.
 */
int MarkDuplicates::runInternal() {
    
    ogeNameThread("am_MarkDuplicates");

    if(verbose)
        cerr << "Reading input file and constructing read end information." << endl;
    buildSortedReadEndLists();
    
    generateDuplicateIndexes();
    
    if(verbose)
        cerr << "Marking " << numDuplicateIndices << " records as duplicates." << endl;
    
    BamReader in;

    in.Open(getBufferFileName());

    // Now copy over the file while marking all the necessary indexes as duplicates
    long recordInFileIndex = 0;
    
    long written = 0;
    while (true) {
        BamAlignment * prec = in.GetNextAlignment();
        if(!prec)
            break;
        
        if (prec->IsPrimaryAlignment()) {
            if (duplicateIndexes.count(recordInFileIndex) == 1)
                prec->SetIsDuplicate(true);
            else
                prec->SetIsDuplicate(false);
        }
        recordInFileIndex++;
        
        if (removeDuplicates && prec->IsDuplicate()) {
            // do nothing
        }
        else {
            putOutputAlignment(prec);
            if (verbose && read_count && ++written % 100000 == 0) {
                cerr << "\rWritten " << written << " records (" << written * 100 / read_count <<"%)." << std::flush;
            }
        }
    }

    if (verbose && read_count)
        cerr << "\rWritten " << written << " records (" << written * 100 / read_count <<"%)." << endl;

    in.Close();

    remove(getBufferFileName().c_str());
    
    return 0;
}