C++ BamReader::HasIndex方法代码示例

// 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;
}

// 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;
}

void getRefAltCounts(BamReader& reader, const BamRegion &r, int& nReadsRef, int& nReadsAlt, vcfRecord record)
{
	// Return the number of reads aligning reference and alternate allele
	if (!reader.IsOpen())
	{
		cout << "Can't access BAM file" << endl;
		return;
	}
	if (!reader.HasIndex())
	{
		cout << "Can't access BAM file index" << endl;
		return;
	}
	
	
	BamAlignment al;
	assert(reader.SetRegion(r));
		
	while ( reader.GetNextAlignment(al) )
	{
		// int baseAlignment = al.Position;
		
		bool found = false;
		int index = 0;
		int indexPosition = al.Position;
		for ( size_t i = 0; i < al.CigarData.size(); i++ )
		{
			switch ( al.CigarData.at(i).Type )
			{
				case 'H' : break; // Hard clip
				case 'S' : break; // Soft clip
				case 'M' : // Match or mismatch
				case 'D' :
				case 'N' :
				case 'P' : 
					if ( indexPosition + al.CigarData.at(i).Length > (record.pos-1) )
					{
						index += (record.pos-1) - indexPosition;
						found = true;
					}
					else
					{
						index += al.CigarData.at(i).Length;
						indexPosition += al.CigarData.at(i).Length;
					}
					break;
				case 'I' :
					index += al.CigarData.at(i).Length;
					break;
			}
			if ( found )
				break;
		}

		if ( !found || index >= al.AlignedBases.size())
		{
			cerr << "Problem with finding variant in read:\n";
			for ( unsigned int i = 0; i < al.AlignedBases.size(); i++ )
				cout << al.AlignedBases.at(i) << " ";
			cout << endl;
			for ( unsigned int i = 0; i < al.CigarData.size(); i++ )
				cout << al.CigarData.at(i).Type << al.CigarData.at(i).Length << " ";
			cout << endl;
			cout << "Read position: " << al.RefID+1 << " " << al.Position+1 << endl;
			cout << "Variant position: " << record.chr << " " << record.pos << endl;
			exit(1);
		}
		
		if ( al.AlignedBases.at(index) == record.ref.at(0) )
			nReadsRef++;
		else if ( al.AlignedBases.at(index) == record.alt.at(0) )
			nReadsAlt++;
	}
}