From a276e3f3926309a9b566bed52ccb8aba9ee334f3 Mon Sep 17 00:00:00 2001
From: Derek <derekwbarnett@gmail.com>
Date: Sat, 9 Oct 2010 19:28:42 -0400
Subject: [PATCH] Fixed: bug(s) related to empty references and regions.

* NOTE - This fix does introduce a slight modification to the *.bti index format.
  So any existing BTI index files will need to be rebuilt to support the bug fix (apologies).
---
 src/api/BamAux.h      |   10 +-
 src/api/BamIndex.cpp  |  153 ++--
 src/api/BamIndex.h    |   10 +-
 src/api/BamReader.cpp | 1584 +++++++++++++++++++++--------------------
 4 files changed, 930 insertions(+), 827 deletions(-)

diff --git a/src/api/BamAux.h b/src/api/BamAux.h
index 3eff2d7..9d38e7d 100644
--- a/src/api/BamAux.h
+++ b/src/api/BamAux.h
@@ -3,7 +3,7 @@
 // Marth Lab, Department of Biology, Boston College
 // All rights reserved.
 // ---------------------------------------------------------------------------
-// Last modified: 18 September 2010 (DB)
+// Last modified: 9 October 2010 (DB)
 // ---------------------------------------------------------------------------
 // Provides the basic constants, data structures, utilities etc. 
 // used throughout the API for handling BAM files
@@ -111,6 +111,14 @@ struct BamRegion {
         , RightPosition(rightPos)
     { }
     
+    // copy constructor
+    BamRegion(const BamRegion& other)
+	: LeftRefID(other.LeftRefID)
+	, LeftPosition(other.LeftPosition)
+	, RightRefID(other.RightRefID)
+	, RightPosition(other.RightPosition)
+    { }
+    
     // member functions
     void clear(void) { LeftRefID = -1; LeftPosition = -1; RightRefID = -1; RightPosition = -1; }
     bool isLeftBoundSpecified(void) const { return ( LeftRefID != -1 && LeftPosition != -1 ); }
diff --git a/src/api/BamIndex.cpp b/src/api/BamIndex.cpp
index f79b2aa..51d7850 100644
--- a/src/api/BamIndex.cpp
+++ b/src/api/BamIndex.cpp
@@ -3,7 +3,7 @@
 // Marth Lab, Department of Biology, Boston College
 // All rights reserved.
 // ---------------------------------------------------------------------------
-// Last modified: 18 September 2010 (DB)
+// Last modified: 9 October 2010 (DB)
 // ---------------------------------------------------------------------------
 // Provides index functionality - both for the default (standardized) BAM 
 // index format (.bai) as well as a BamTools-specific (nonstandard) index 
@@ -33,17 +33,6 @@ BamIndex::BamIndex(BamTools::BgzfData* bgzf, BamTools::BamReader* reader, bool i
         m_references = m_reader->GetReferenceData();
 }
 
-bool BamIndex::HasAlignments(const int& referenceID) {
-    
-    // return false if invalid ID
-    if ( (referenceID < 0) || (referenceID >= (int)m_references.size()) ) 
-        return false;
-    
-    // else return status of reference (has alignments?)
-    else
-        return m_references.at(referenceID).RefHasAlignments;
-}
-
 // #########################################################################################
 // #########################################################################################
 
@@ -120,9 +109,11 @@ struct BamStandardIndex::BamStandardIndexPrivate {
     
     // creates index data (in-memory) from current reader data
     bool Build(void);
+    // returns whether reference has alignments or no
+    bool HasAlignments(const int& referenceID) const;
     // attempts to use index to jump to region; returns success/fail
     bool Jump(const BamTools::BamRegion& region, bool* hasAlignmentsInRegion);
-      // loads existing data from file into memory
+    // loads existing data from file into memory
     bool Load(const std::string& filename);
     // writes in-memory index data out to file 
     // N.B. - (this is the original BAM filename, method will modify it to use applicable extension)
@@ -134,7 +125,7 @@ struct BamStandardIndex::BamStandardIndexPrivate {
     // calculate bins that overlap region
     int BinsFromRegion(const BamRegion& region, const bool isRightBoundSpecified, uint16_t bins[BamTools::MAX_BIN]);
     // calculates offset(s) for a given region
-    bool GetOffsets(const BamRegion& region, const bool isRightBoundSpecified, vector<int64_t>& offsets);
+    bool GetOffsets(const BamRegion& region, const bool isRightBoundSpecified, vector<int64_t>& offsets, bool* hasAlignmentsInRegion);
     // saves BAM bin entry for index
     void InsertBinEntry(BamBinMap& binMap, const uint32_t& saveBin, const uint64_t& saveOffset, const uint64_t& lastOffset);
     // saves linear offset entry for index
@@ -175,13 +166,13 @@ int BamStandardIndex::BamStandardIndexPrivate::BinsFromRegion(const BamRegion& r
     return i;
 }
 
+// creates index data (in-memory) from current reader data
 bool BamStandardIndex::BamStandardIndexPrivate::Build(void) { 
   
     // localize parent data
     if ( m_parent == 0 ) return false;
-    BamReader* reader = m_parent->m_reader;
-    BgzfData*  mBGZF  = m_parent->m_BGZF;
-    RefVector& references = m_parent->m_references;
+    BamReader* reader     = m_parent->m_reader;
+    BgzfData*  mBGZF      = m_parent->m_BGZF;
   
     // be sure reader & BGZF file are valid & open for reading
     if ( reader == 0 || mBGZF == 0 || !mBGZF->IsOpen ) 
@@ -191,7 +182,7 @@ bool BamStandardIndex::BamStandardIndexPrivate::Build(void) {
     reader->Rewind();
 
     // get reference count, reserve index space
-    int numReferences = (int)references.size();
+    const int numReferences = (int)m_parent->m_references.size();
     for ( int i = 0; i < numReferences; ++i ) 
         m_indexData.push_back(ReferenceIndex());
     
@@ -233,7 +224,7 @@ bool BamStandardIndex::BamStandardIndexPrivate::Build(void) {
         if ( (bAlignment.RefID >= 0) && (bAlignment.Bin < 4681) ) {
 
             // save linear offset entry (matched to BAM entry refID)
-            ReferenceIndex& refIndex = m_indexData.at(bAlignment.RefID);
+            ReferenceIndex& refIndex    = m_indexData.at(bAlignment.RefID);
             LinearOffsetVector& offsets = refIndex.Offsets;
             InsertLinearOffset(offsets, bAlignment, lastOffset);
         }
@@ -289,7 +280,6 @@ bool BamStandardIndex::BamStandardIndexPrivate::Build(void) {
     MergeChunks();
     
     // iterate through references in index
-    // store whether reference has data &
     // sort offsets in linear offset vector
     BamStandardIndexData::iterator indexIter = m_indexData.begin();
     BamStandardIndexData::iterator indexEnd  = m_indexData.end();
@@ -297,12 +287,8 @@ bool BamStandardIndex::BamStandardIndexPrivate::Build(void) {
 
         // get reference index data
         ReferenceIndex& refIndex = (*indexIter);
-        BamBinMap& binMap = refIndex.Bins;
         LinearOffsetVector& offsets = refIndex.Offsets;
 
-        // store whether reference has alignments or no
-        references[i].RefHasAlignments = ( binMap.size() > 0 );
-
         // sort linear offsets
         sort(offsets.begin(), offsets.end());
     }
@@ -311,7 +297,8 @@ bool BamStandardIndex::BamStandardIndexPrivate::Build(void) {
     return reader->Rewind();
 }
 
-bool BamStandardIndex::BamStandardIndexPrivate::GetOffsets(const BamRegion& region, const bool isRightBoundSpecified, vector<int64_t>& offsets) { 
+// calculates offset(s) for a given region
+bool BamStandardIndex::BamStandardIndexPrivate::GetOffsets(const BamRegion& region, const bool isRightBoundSpecified, vector<int64_t>& offsets, bool* hasAlignmentsInRegion) { 
   
     // calculate which bins overlap this region
     uint16_t* bins = (uint16_t*)calloc(MAX_BIN, 2);
@@ -352,8 +339,21 @@ bool BamStandardIndex::BamStandardIndexPrivate::GetOffsets(const BamRegion& regi
     // sort the offsets before returning
     sort(offsets.begin(), offsets.end());
     
-    // return whether any offsets were found
-    return ( offsets.size() != 0 );
+    // set flag & return success
+    *hasAlignmentsInRegion = (offsets.size() != 0 );
+    return true;
+}
+
+// returns whether reference has alignments or no
+bool BamStandardIndex::BamStandardIndexPrivate::HasAlignments(const int& referenceID) const {
+    
+    // ID not found
+    if ( referenceID < 0 || referenceID >= (int)m_indexData.size() ) 
+	return false;
+	
+    // return whether reference contains data
+    const ReferenceIndex& refIndex = m_indexData.at(referenceID);
+    return !( refIndex.Bins.empty() );
 }
 
 // saves BAM bin entry for index
@@ -404,6 +404,7 @@ void BamStandardIndex::BamStandardIndexPrivate::InsertLinearOffset(LinearOffsetV
     }
 }      
 
+// attempts to use index to jump to region; returns success/fail
 bool BamStandardIndex::BamStandardIndexPrivate::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) {
   
     // localize parent data
@@ -416,17 +417,14 @@ bool BamStandardIndex::BamStandardIndexPrivate::Jump(const BamRegion& region, bo
     if ( reader == 0 || mBGZF == 0 || !mBGZF->IsOpen ) 
         return false;
     
-    // see if left-bound reference of region has alignments
-    if ( !m_parent->HasAlignments(region.LeftRefID) ) return false; 
-    
     // make sure left-bound position is valid
     if ( region.LeftPosition > references.at(region.LeftRefID).RefLength ) return false;
         
     // calculate offsets for this region
     // if failed, print message, set flag, and return failure
     vector<int64_t> offsets;
-    if ( !GetOffsets(region, region.isRightBoundSpecified(), offsets) ) {
-        fprintf(stderr, "ERROR: Could not jump: unable to calculate offset for specified region.\n");
+    if ( !GetOffsets(region, region.isRightBoundSpecified(), offsets, hasAlignmentsInRegion) ) {
+        fprintf(stderr, "ERROR: Could not jump: unable to calculate offset(s) for specified region.\n");
         *hasAlignmentsInRegion = false;
         return false;
     }
@@ -451,7 +449,7 @@ bool BamStandardIndex::BamStandardIndexPrivate::Jump(const BamRegion& region, bo
     
     // if error in jumping, print message & set flag
     if ( !result ) {
-        fprintf(stderr, "ERROR: Could not jump: unable to calculate offset for specified region.\n");
+        fprintf(stderr, "ERROR: Could not jump: unable to determine correct offset for specified region.\n");
         *hasAlignmentsInRegion = false;
     }
     
@@ -459,12 +457,12 @@ bool BamStandardIndex::BamStandardIndexPrivate::Jump(const BamRegion& region, bo
     return result;
 }
 
+// loads existing data from file into memory
 bool BamStandardIndex::BamStandardIndexPrivate::Load(const string& filename)  { 
     
     // localize parent data
     if ( m_parent == 0 ) return false;
     const bool isBigEndian = m_parent->m_isBigEndian;
-    RefVector& references  = m_parent->m_references;
   
     // open index file, abort on error
     FILE* indexStream = fopen(filename.c_str(), "rb");
@@ -501,11 +499,6 @@ bool BamStandardIndex::BamStandardIndexPrivate::Load(const string& filename)  {
         elementsRead = fread(&numBins, 4, 1, indexStream);
         if ( isBigEndian ) SwapEndian_32(numBins);
 
-        if ( numBins > 0 ) {
-            RefData& refEntry = references[i];
-            refEntry.RefHasAlignments = true;
-        }
-
         // intialize BinVector
         BamBinMap binMap;
 
@@ -652,6 +645,7 @@ bool BamStandardIndex::BamStandardIndexPrivate::Write(const std::string& bamFile
     if ( m_parent == 0 ) return false;
     const bool isBigEndian = m_parent->m_isBigEndian;
   
+    // open index file
     string indexFilename = bamFilename + ".bai";
     FILE* indexStream = fopen(indexFilename.c_str(), "wb");
     if ( indexStream == 0 ) {
@@ -761,6 +755,9 @@ BamStandardIndex::~BamStandardIndex(void) {
 // creates index data (in-memory) from current reader data
 bool BamStandardIndex::Build(void) { return d->Build(); }
 
+// returns whether reference has alignments or no
+bool BamStandardIndex::HasAlignments(const int& referenceID) const { return d->HasAlignments(referenceID); }
+
 // attempts to use index to jump to region; returns success/fail
 bool BamStandardIndex::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) { return d->Jump(region, hasAlignmentsInRegion); }
 
@@ -809,7 +806,7 @@ struct BamToolsIndex::BamToolsIndexPrivate {
     
     // keep a list of any supported versions here 
     // (might be useful later to handle any 'legacy' versions if the format changes)
-    // listed for example like: BTI_1_1 = 2, BTI_1_2 = 3, BTI_2_0 = 4, and so on
+    // listed for example like: BTI_1_0 = 1, BTI_1_1 = 2, BTI_1_2 = 3, BTI_2_0 = 4, and so on
     //
     // so a change introduced in (hypothetical) BTI_1_2 would be handled from then on by: 
     //
@@ -817,12 +814,15 @@ struct BamToolsIndex::BamToolsIndexPrivate {
     //   do something new 
     // else 
     //   do the old thing
+    
     enum Version { BTI_1_0 = 1 
                  , BTI_1_1
+                 , BTI_1_2
                  };  
   
     // -------------------------
     // data members
+    
     BamToolsIndexData m_indexData;
     BamToolsIndex*    m_parent;
     int32_t           m_blockSize;
@@ -834,7 +834,7 @@ struct BamToolsIndex::BamToolsIndexPrivate {
     BamToolsIndexPrivate(BamToolsIndex* parent) 
         : m_parent(parent)
         , m_blockSize(1000)
-        , m_version(BTI_1_1) // latest version - used for writing new index files
+        , m_version(BTI_1_2) // latest version - used for writing new index files
     { }
     
     ~BamToolsIndexPrivate(void) { }
@@ -845,10 +845,12 @@ struct BamToolsIndex::BamToolsIndexPrivate {
     // creates index data (in-memory) from current reader data
     bool Build(void);
     // returns supported file extension
-    const std::string Extension(void) const { return std::string(".bti"); }
+    const string Extension(void) const { return string(".bti"); }
+    // returns whether reference has alignments or no
+    bool HasAlignments(const int& referenceID) const;
     // attempts to use index to jump to region; returns success/fail
     bool Jump(const BamTools::BamRegion& region, bool* hasAlignmentsInRegion);
-      // loads existing data from file into memory
+    // loads existing data from file into memory
     bool Load(const std::string& filename);
     // writes in-memory index data out to file 
     // N.B. - (this is the original BAM filename, method will modify it to use applicable extension)
@@ -861,13 +863,13 @@ struct BamToolsIndex::BamToolsIndexPrivate {
     bool GetOffset(const BamRegion& region, int64_t& offset, bool* hasAlignmentsInRegion);
 };
 
+// creates index data (in-memory) from current reader data
 bool BamToolsIndex::BamToolsIndexPrivate::Build(void) { 
   
     // localize parent data
     if ( m_parent == 0 ) return false;
-    BamReader* reader = m_parent->m_reader;
-    BgzfData*  mBGZF  = m_parent->m_BGZF;
-    RefVector& references = m_parent->m_references;
+    BamReader* reader     = m_parent->m_reader;
+    BgzfData*  mBGZF      = m_parent->m_BGZF;
   
     // be sure reader & BGZF file are valid & open for reading
     if ( reader == 0 || mBGZF == 0 || !mBGZF->IsOpen ) 
@@ -878,6 +880,11 @@ bool BamToolsIndex::BamToolsIndexPrivate::Build(void) {
     if ( !reader->Rewind() ) 
         return false;
     
+    // initialize index data structure with space for all references
+    const int numReferences = (int)m_parent->m_references.size();
+    for ( int i = 0; i < numReferences; ++i )
+	m_indexData[i].clear();    
+    
     // set up counters and markers
     int32_t currentBlockCount      = 0;
     int64_t currentAlignmentOffset = mBGZF->Tell();
@@ -893,13 +900,10 @@ bool BamToolsIndex::BamToolsIndexPrivate::Build(void) {
         // if block contains data (not the first time through) AND alignment is on a new reference
         if ( currentBlockCount > 0 && al.RefID != blockRefId ) {
           
-            // make sure reference flag is set
-            references[blockRefId].RefHasAlignments = true;
-          
-            // store previous data
-            m_indexData[blockRefId].push_back( BamToolsIndexEntry(blockMaxEndPosition, blockStartOffset, blockStartPosition) );
-            
-            // intialize new block
+	    // store previous data
+	    m_indexData[blockRefId].push_back( BamToolsIndexEntry(blockMaxEndPosition, blockStartOffset, blockStartPosition) );
+
+            // intialize new block for current alignment's reference
             currentBlockCount   = 0;
             blockMaxEndPosition = al.GetEndPosition();
             blockStartOffset    = currentAlignmentOffset;
@@ -913,7 +917,7 @@ bool BamToolsIndex::BamToolsIndexPrivate::Build(void) {
       
         // increment block counter
         ++currentBlockCount;
-        
+
         // check end position
         int32_t alignmentEndPosition = al.GetEndPosition();
         if ( alignmentEndPosition > blockMaxEndPosition ) 
@@ -921,10 +925,6 @@ bool BamToolsIndex::BamToolsIndexPrivate::Build(void) {
         
         // if block is full, get offset for next block, reset currentBlockCount
         if ( currentBlockCount == m_blockSize ) {
-          
-            // make sure reference flag is set
-            references[blockRefId].RefHasAlignments = true;
-          
             m_indexData[blockRefId].push_back( BamToolsIndexEntry(blockMaxEndPosition, blockStartOffset, blockStartPosition) );
             blockStartOffset  = mBGZF->Tell();
             currentBlockCount = 0;
@@ -935,7 +935,7 @@ bool BamToolsIndex::BamToolsIndexPrivate::Build(void) {
         currentAlignmentOffset = mBGZF->Tell();  
     }
     
-    // store final block
+    // store final block with data
     m_indexData[blockRefId].push_back( BamToolsIndexEntry(blockMaxEndPosition, blockStartOffset, blockStartPosition) );
     
     // attempt to rewind back to begnning of alignments
@@ -973,6 +973,15 @@ bool BamToolsIndex::BamToolsIndexPrivate::GetOffset(const BamRegion& region, int
     return true; 
 }
 
+// returns whether reference has alignments or no
+bool BamToolsIndex::BamToolsIndexPrivate::HasAlignments(const int& referenceID) const {  
+    if ( m_indexData.find(referenceID) == m_indexData.end() ) 
+	return false;
+    else 
+      return ( !m_indexData.at(referenceID).empty() );
+}
+
+// attempts to use index to jump to region; returns success/fail
 bool BamToolsIndex::BamToolsIndexPrivate::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) {
   
     // clear flag
@@ -990,9 +999,6 @@ bool BamToolsIndex::BamToolsIndexPrivate::Jump(const BamRegion& region, bool* ha
         return false;
     }
   
-    // see if left-bound reference of region has alignments
-    if ( !m_parent->HasAlignments(region.LeftRefID) ) return false; 
-    
     // make sure left-bound position is valid
     if ( region.LeftPosition > references.at(region.LeftRefID).RefLength ) return false;
   
@@ -1002,17 +1008,17 @@ bool BamToolsIndex::BamToolsIndexPrivate::Jump(const BamRegion& region, bool* ha
         fprintf(stderr, "ERROR: Could not jump - unable to calculate offset for specified region.\n");
         return false;
     }
-  
+    
     // attempt seek in file, return success/fail
     return mBGZF->Seek(offset);    
 }
 
+// loads existing data from file into memory
 bool BamToolsIndex::BamToolsIndexPrivate::Load(const string& filename) { 
   
     // localize parent data
     if ( m_parent == 0 ) return false;
     const bool isBigEndian = m_parent->m_isBigEndian;
-    RefVector& references  = m_parent->m_references;
   
     // open index file, abort on error
     FILE* indexStream = fopen(filename.c_str(), "rb");
@@ -1032,7 +1038,7 @@ bool BamToolsIndex::BamToolsIndexPrivate::Load(const string& filename) {
         fclose(indexStream);
         return false;
     }
-
+    
     // read in version
     int32_t indexVersion;
     elementsRead = fread(&indexVersion, sizeof(indexVersion), 1, indexStream);
@@ -1043,12 +1049,19 @@ bool BamToolsIndex::BamToolsIndexPrivate::Load(const string& filename) {
         return false;
     }
 
-    if ( (Version)indexVersion < BTI_1_1 ) {
+    if ( (Version)indexVersion == BTI_1_0 ) {
         fprintf(stderr, "\nProblem with index file - this version of the index contains a bug related to accessing data near reference ends.\n");
         fprintf(stderr, "\nPlease run \'bamtools index -bti -in yourData.bam\' to generate an up-to-date BamToolsIndex.\n\n");
         fclose(indexStream);
         exit(1);
     }
+    
+    if ( (Version)indexVersion == BTI_1_1 ) {
+        fprintf(stderr, "\nProblem with index file - this version of the index contains a bug related to handling empty references.\n");
+        fprintf(stderr, "\nPlease run \'bamtools index -bti -in yourData.bam\' to generate an up-to-date BamToolsIndex.\n\n");
+        fclose(indexStream);
+        exit(1);
+    }
 
     // read in block size
     elementsRead = fread(&m_blockSize, sizeof(m_blockSize), 1, indexStream);
@@ -1058,7 +1071,7 @@ bool BamToolsIndex::BamToolsIndexPrivate::Load(const string& filename) {
     int32_t numReferences;
     elementsRead = fread(&numReferences, sizeof(numReferences), 1, indexStream);
     if ( isBigEndian ) SwapEndian_32(numReferences);
-    
+
     // iterate over reference entries
     for ( int i = 0; i < numReferences; ++i ) {
       
@@ -1097,9 +1110,6 @@ bool BamToolsIndex::BamToolsIndexPrivate::Load(const string& filename) {
         
         // save refID, offsetVector entry into data structure
         m_indexData.insert( make_pair(i, offsets) );
-        
-        // set ref.HasAlignments flag
-        references[i].RefHasAlignments = (numOffsets != 0);
     }
 
     // close index file and return
@@ -1203,6 +1213,9 @@ BamToolsIndex::~BamToolsIndex(void) {
 // creates index data (in-memory) from current reader data
 bool BamToolsIndex::Build(void) { return d->Build(); }
 
+// returns whether reference has alignments or no
+bool BamToolsIndex::HasAlignments(const int& referenceID) const { return d->HasAlignments(referenceID); }
+
 // attempts to use index to jump to region; returns success/fail
 // a "successful" jump indicates no error, but not whether this region has data
 //   * thus, the method sets a flag to indicate whether there are alignments 
diff --git a/src/api/BamIndex.h b/src/api/BamIndex.h
index 710893a..9da858f 100644
--- a/src/api/BamIndex.h
+++ b/src/api/BamIndex.h
@@ -3,7 +3,7 @@
 // Marth Lab, Department of Biology, Boston College
 // All rights reserved.
 // ---------------------------------------------------------------------------
-// Last modified: 18 September 2010 (DB)
+// Last modified: 8 October 2010 (DB)
 // ---------------------------------------------------------------------------
 // Provides index functionality - both for the standardized BAM index format
 // (".bai") as well as a BamTools-specific (nonstandard) index format (".bti").
@@ -43,7 +43,7 @@ class BamIndex {
         //     available after the jump position
         virtual bool Jump(const BamTools::BamRegion& region, bool* hasAlignmentsInRegion) =0;
         // returns whether reference has alignments or no
-        virtual bool HasAlignments(const int& referenceID); 
+        virtual bool HasAlignments(const int& referenceID) const =0;
         // loads existing data from file into memory
         virtual bool Load(const std::string& filename)  =0;
         // writes in-memory index data out to file 
@@ -58,7 +58,7 @@ class BamIndex {
         // (if not found, attmempts to load other type(s), returns 0 if NONE found)
         //
         // ** default preferred type is BamToolsIndex ** use this anytime it exists
-         enum PreferredIndexType { BAMTOOLS = 0, STANDARD };
+        enum PreferredIndexType { BAMTOOLS = 0, STANDARD };
         static BamIndex* FromBamFilename(const std::string& bamFilename, 
                                          BamTools::BgzfData* bgzf, 
                                          BamTools::BamReader* reader, 
@@ -99,6 +99,8 @@ class BamStandardIndex : public BamIndex {
         bool Build(void);
         // returns supported file extension
         const std::string Extension(void) const { return std::string(".bai"); }
+        // returns whether reference has alignments or no
+        bool HasAlignments(const int& referenceID) const;
         // attempts to use index to jump to region; returns success/fail
         // a "successful" jump indicates no error, but not whether this region has data
         //   * thus, the method sets a flag to indicate whether there are alignments 
@@ -135,6 +137,8 @@ class BamToolsIndex : public BamIndex {
         bool Build(void);
         // returns supported file extension
         const std::string Extension(void) const { return std::string(".bti"); }
+        // returns whether reference has alignments or no
+        bool HasAlignments(const int& referenceID) const;
         // attempts to use index to jump to region; returns success/fail
         // a "successful" jump indicates no error, but not whether this region has data
         //   * thus, the method sets a flag to indicate whether there are alignments 
diff --git a/src/api/BamReader.cpp b/src/api/BamReader.cpp
index f16d98d..ce21fab 100644
--- a/src/api/BamReader.cpp
+++ b/src/api/BamReader.cpp
@@ -1,753 +1,831 @@
-// ***************************************************************************
-// BamReader.cpp (c) 2009 Derek Barnett, Michael Str�mberg
-// Marth Lab, Department of Biology, Boston College
-// All rights reserved.
-// ---------------------------------------------------------------------------
-// Last modified: 17 September 2010 (DB)
-// ---------------------------------------------------------------------------
-// Uses BGZF routines were adapted from the bgzf.c code developed at the Broad
-// Institute.
-// ---------------------------------------------------------------------------
-// Provides the basic functionality for reading BAM files
-// ***************************************************************************
-
-// C++ includes
-#include <algorithm>
-#include <iterator>
-#include <string>
-#include <vector>
-#include <iostream>
-#include "BGZF.h"
-#include "BamReader.h"
-#include "BamIndex.h"
-using namespace BamTools;
-using namespace std;
-
-struct BamReader::BamReaderPrivate {
-
-    // -------------------------------
-    // structs, enums, typedefs
-    // -------------------------------
-    enum RegionState { BEFORE_REGION = 0
-                     , WITHIN_REGION
-                     , AFTER_REGION
-                     };
-
-    // -------------------------------
-    // data members
-    // -------------------------------
-
-    // general file data
-    BgzfData  mBGZF;
-    string    HeaderText;
-    BamIndex* Index;
-    RefVector References;
-    bool      IsIndexLoaded;
-    int64_t   AlignmentsBeginOffset;
-    string    Filename;
-    string    IndexFilename;
-    
-    // system data
-    bool IsBigEndian;
-
-    // user-specified region values
-    BamRegion Region;
-    bool HasAlignmentsInRegion;
-
-    // parent BamReader
-    BamReader* Parent;
-    
-    // BAM character constants
-    const char* DNA_LOOKUP;
-    const char* CIGAR_LOOKUP;
-
-    // -------------------------------
-    // constructor & destructor
-    // -------------------------------
-    BamReaderPrivate(BamReader* parent);
-    ~BamReaderPrivate(void);
-
-    // -------------------------------
-    // "public" interface
-    // -------------------------------
-
-    // file operations
-    void Close(void);
-    bool Open(const std::string& filename, 
-              const std::string& indexFilename, 
-              const bool lookForIndex, 
-              const bool preferStandardIndex);
-    bool Rewind(void);
-    bool SetRegion(const BamRegion& region);
-
-    // access alignment data
-    bool GetNextAlignment(BamAlignment& bAlignment);
-    bool GetNextAlignmentCore(BamAlignment& bAlignment);
-
-    // access auxiliary data
-    int GetReferenceID(const string& refName) const;
-
-    // index operations
-    bool CreateIndex(bool useStandardIndex);
-
-    // -------------------------------
-    // internal methods
-    // -------------------------------
-
-    // *** reading alignments and auxiliary data *** //
-
-    // fills out character data for BamAlignment data
-    bool BuildCharData(BamAlignment& bAlignment);
-    // checks to see if alignment overlaps current region
-    RegionState IsOverlap(BamAlignment& bAlignment);
-    // retrieves header text from BAM file
-    void LoadHeaderData(void);
-    // retrieves BAM alignment under file pointer
-    bool LoadNextAlignment(BamAlignment& bAlignment);
-    // builds reference data structure from BAM file
-    void LoadReferenceData(void);
-
-    // *** index file handling *** //
-
-    // clear out inernal index data structure
-    void ClearIndex(void);
-    // loads index from BAM index file
-    bool LoadIndex(const bool lookForIndex, const bool preferStandardIndex);
-};
-
-// -----------------------------------------------------
-// BamReader implementation (wrapper around BRPrivate)
-// -----------------------------------------------------
-// constructor
-BamReader::BamReader(void) {
-    d = new BamReaderPrivate(this);
-}
-
-// destructor
-BamReader::~BamReader(void) {
-    delete d;
-    d = 0;
-}
-
-// file operations
-void BamReader::Close(void) { d->Close(); }
-bool BamReader::IsIndexLoaded(void) const { return d->IsIndexLoaded; }
-bool BamReader::IsOpen(void) const { return d->mBGZF.IsOpen; }
-bool BamReader::Jump(int refID, int position)  { return d->SetRegion( BamRegion(refID, position) ); }
-bool BamReader::Open(const std::string& filename, 
-                     const std::string& indexFilename, 
-                     const bool lookForIndex, 
-                     const bool preferStandardIndex) 
-{ 
-    return d->Open(filename, indexFilename, lookForIndex, preferStandardIndex); 
-}
-bool BamReader::Rewind(void) { return d->Rewind(); }
-bool BamReader::SetRegion(const BamRegion& region) { return d->SetRegion(region); }
-bool BamReader::SetRegion(const int& leftRefID, const int& leftBound, const int& rightRefID, const int& rightBound) {
-    return d->SetRegion( BamRegion(leftRefID, leftBound, rightRefID, rightBound) );
-}
-
-// access alignment data
-bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }
-bool BamReader::GetNextAlignmentCore(BamAlignment& bAlignment) { return d->GetNextAlignmentCore(bAlignment); }
-
-// access auxiliary data
-const string BamReader::GetHeaderText(void) const { return d->HeaderText; }
-int BamReader::GetReferenceCount(void) const { return d->References.size(); }
-const RefVector& BamReader::GetReferenceData(void) const { return d->References; }
-int BamReader::GetReferenceID(const string& refName) const { return d->GetReferenceID(refName); }
-const std::string BamReader::GetFilename(void) const { return d->Filename; }
-
-// index operations
-bool BamReader::CreateIndex(bool useStandardIndex) { return d->CreateIndex(useStandardIndex); }
-
-// -----------------------------------------------------
-// BamReaderPrivate implementation
-// -----------------------------------------------------
-
-// constructor
-BamReader::BamReaderPrivate::BamReaderPrivate(BamReader* parent)
-    : Index(0)
-    , IsIndexLoaded(false)
-    , AlignmentsBeginOffset(0)
-    , Parent(parent)
-    , DNA_LOOKUP("=ACMGRSVTWYHKDBN")
-    , CIGAR_LOOKUP("MIDNSHP")
-{ 
-    IsBigEndian = SystemIsBigEndian();
-}
-
-// destructor
-BamReader::BamReaderPrivate::~BamReaderPrivate(void) {
-    Close();
-}
-
-bool BamReader::BamReaderPrivate::BuildCharData(BamAlignment& bAlignment) {
-  
-    // calculate character lengths/offsets
-    const unsigned int dataLength      = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
-    const unsigned int seqDataOffset   = bAlignment.SupportData.QueryNameLength + (bAlignment.SupportData.NumCigarOperations * 4);
-    const unsigned int qualDataOffset  = seqDataOffset + (bAlignment.SupportData.QuerySequenceLength+1)/2;
-    const unsigned int tagDataOffset   = qualDataOffset + bAlignment.SupportData.QuerySequenceLength;
-    const unsigned int tagDataLength   = dataLength - tagDataOffset;
-      
-    // set up char buffers
-    const char*     allCharData = bAlignment.SupportData.AllCharData.data();
-    const char*     seqData     = ((const char*)allCharData) + seqDataOffset;
-    const char*     qualData    = ((const char*)allCharData) + qualDataOffset;
-          char*     tagData     = ((char*)allCharData) + tagDataOffset;
-  
-    // store alignment name (relies on null char in name as terminator)
-    bAlignment.Name.assign((const char*)(allCharData));    
-
-    // save query sequence
-    bAlignment.QueryBases.clear();
-    bAlignment.QueryBases.reserve(bAlignment.SupportData.QuerySequenceLength);
-    for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
-        char singleBase = DNA_LOOKUP[ ( ( seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
-        bAlignment.QueryBases.append(1, singleBase);
-    }
-  
-    // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
-    bAlignment.Qualities.clear();
-    bAlignment.Qualities.reserve(bAlignment.SupportData.QuerySequenceLength);
-    for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
-        char singleQuality = (char)(qualData[i]+33);
-        bAlignment.Qualities.append(1, singleQuality);
-    }
-    
-    // if QueryBases is empty (and this is a allowed case)
-    if ( bAlignment.QueryBases.empty() ) 
-        bAlignment.AlignedBases = bAlignment.QueryBases;
-    
-    // if QueryBases contains data, then build AlignedBases using CIGAR data
-    else {
-    
-        // resize AlignedBases
-        bAlignment.AlignedBases.clear();
-        bAlignment.AlignedBases.reserve(bAlignment.SupportData.QuerySequenceLength);
-      
-        // iterate over CigarOps
-        int k = 0;
-        vector<CigarOp>::const_iterator cigarIter = bAlignment.CigarData.begin();
-        vector<CigarOp>::const_iterator cigarEnd  = bAlignment.CigarData.end();
-        for ( ; cigarIter != cigarEnd; ++cigarIter ) {
-            
-            const CigarOp& op = (*cigarIter);
-            switch(op.Type) {
-              
-                case ('M') :
-                case ('I') :
-                    bAlignment.AlignedBases.append(bAlignment.QueryBases.substr(k, op.Length)); // for 'M', 'I' - write bases
-                    // fall through
-                
-                case ('S') :
-                    k += op.Length;                                     // for 'S' - soft clip, skip over query bases
-                    break;
-                    
-                case ('D') :
-                    bAlignment.AlignedBases.append(op.Length, '-');     // for 'D' - write gap character
-                    break;
-                    
-                case ('P') :
-                    bAlignment.AlignedBases.append( op.Length, '*' );   // for 'P' - write padding character
-                    break;
-                    
-                case ('N') :
-                    bAlignment.AlignedBases.append( op.Length, 'N' );  // for 'N' - write N's, skip bases in original query sequence
-                    break;
-                    
-                case ('H') :
-                    break;  // for 'H' - hard clip, do nothing to AlignedBases, move to next op
-                    
-                default:
-                    fprintf(stderr, "ERROR: Invalid Cigar op type\n"); // shouldn't get here
-                    exit(1);
-            }
-        }
-    }
- 
-    // -----------------------
-    // Added: 3-25-2010 DB
-    // Fixed: endian-correctness for tag data
-    // -----------------------
-    if ( IsBigEndian ) {
-        int i = 0;
-        while ( (unsigned int)i < tagDataLength ) {
-          
-            i += 2; // skip tag type (e.g. "RG", "NM", etc)
-            uint8_t type = toupper(tagData[i]);     // lower & upper case letters have same meaning 
-            ++i;                                    // skip value type
-    
-            switch (type) {
-                
-                case('A') :
-                case('C') : 
-                    ++i;
-                    break;
-
-                case('S') : 
-                    SwapEndian_16p(&tagData[i]); 
-                    i += sizeof(uint16_t);
-                    break;
-                    
-                case('F') :
-                case('I') : 
-                    SwapEndian_32p(&tagData[i]);
-                    i += sizeof(uint32_t);
-                    break;
-                
-                case('D') : 
-                    SwapEndian_64p(&tagData[i]);
-                    i += sizeof(uint64_t);
-                    break;
-                
-                case('H') :
-                case('Z') : 
-                    while (tagData[i]) { ++i; }
-                    ++i; // increment one more for null terminator
-                    break;
-                
-                default : 
-                    fprintf(stderr, "ERROR: Invalid tag value type\n"); // shouldn't get here
-                    exit(1);
-            }
-        }
-    }
-    
-    // store TagData
-    bAlignment.TagData.clear();
-    bAlignment.TagData.resize(tagDataLength);
-    memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLength);
-    
-    // clear the core-only flag
-    bAlignment.SupportData.HasCoreOnly = false;
-    
-    // return success
-    return true;
-}
-
-// clear index data structure
-void BamReader::BamReaderPrivate::ClearIndex(void) {
-    delete Index;
-    Index = 0;
-    IsIndexLoaded = false;
-}
-
-// closes the BAM file
-void BamReader::BamReaderPrivate::Close(void) {
-    
-    // close BGZF file stream
-    mBGZF.Close();
-    
-    // clear out index data
-    ClearIndex();
-    
-    // clear out header data
-    HeaderText.clear();
-    
-    // clear out region flags
-    Region.clear();
-}
-
-// creates index for BAM file, saves to file
-// default behavior is to create the BAM standard index (".bai")
-// set flag to false to create the BamTools-specific index (".bti")
-bool BamReader::BamReaderPrivate::CreateIndex(bool useStandardIndex) {
-
-    // clear out prior index data
-    ClearIndex();
-    
-    // create index based on type requested
-    if ( useStandardIndex ) 
-        Index = new BamStandardIndex(&mBGZF, Parent, IsBigEndian);
-    // create BamTools 'custom' index
-    else
-        Index = new BamToolsIndex(&mBGZF, Parent, IsBigEndian);
-    
-    // build new index
-    bool ok = true;
-    ok &= Index->Build();
-    IsIndexLoaded = ok;
-    
-    // attempt to save index data to file
-    ok &= Index->Write(Filename); 
-    
-    // return success/fail of both building & writing index
-    return ok;
-}
-
-// get next alignment (from specified region, if given)
-bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {
-
-    // if valid alignment found, attempt to parse char data, and return success/failure
-    if ( GetNextAlignmentCore(bAlignment) )
-        return BuildCharData(bAlignment);
-    
-    // no valid alignment found
-    else
-        return false;
-}
-
-// retrieves next available alignment core data (returns success/fail)
-// ** DOES NOT parse any character data (read name, bases, qualities, tag data)
-//    these can be accessed, if necessary, from the supportData 
-// useful for operations requiring ONLY positional or other alignment-related information
-bool BamReader::BamReaderPrivate::GetNextAlignmentCore(BamAlignment& bAlignment) {
-
-    // if region is set but has no alignments
-    if ( !Region.isNull() && !HasAlignmentsInRegion ) 
-        return false;
-  
-    // if valid alignment available
-    if ( LoadNextAlignment(bAlignment) ) {
-
-        // set core-only flag
-        bAlignment.SupportData.HasCoreOnly = true;
-      
-        // if region not specified, return success
-        if ( !Region.isLeftBoundSpecified() ) return true;
-
-        // determine region state (before, within, after)
-        BamReader::BamReaderPrivate::RegionState state = IsOverlap(bAlignment);
-      
-        // if alignment lies after region, return false
-        if ( state == AFTER_REGION ) return false;
-
-        while ( state != WITHIN_REGION ) {
-            // if no valid alignment available (likely EOF) return failure
-            if ( !LoadNextAlignment(bAlignment) ) return false;
-            // if alignment lies after region, return false (no available read within region)
-            state = IsOverlap(bAlignment);
-            if ( state == AFTER_REGION ) return false;
-        }
-
-        // return success (alignment found that overlaps region)
-        return true;
-    }
-
-    // no valid alignment
-    else
-        return false;
-}
-
-// returns RefID for given RefName (returns References.size() if not found)
-int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {
-
-    // retrieve names from reference data
-    vector<string> refNames;
-    RefVector::const_iterator refIter = References.begin();
-    RefVector::const_iterator refEnd  = References.end();
-    for ( ; refIter != refEnd; ++refIter) 
-        refNames.push_back( (*refIter).RefName );
-
-    // return 'index-of' refName ( if not found, returns refNames.size() )
-    return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));
-}
-
-// returns region state - whether alignment ends before, overlaps, or starts after currently specified region
-// this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true
-BamReader::BamReaderPrivate::RegionState BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {
-    
-    // --------------------------------------------------
-    // check alignment start against right bound cutoff
-    
-    // if full region of interest was given
-    if ( Region.isRightBoundSpecified() ) {
-      
-        // read starts on right bound reference, but AFTER right bound position
-        if ( bAlignment.RefID == Region.RightRefID && bAlignment.Position > Region.RightPosition )
-            return AFTER_REGION;
-      
-        // if read starts on reference AFTER right bound, return false
-        if ( bAlignment.RefID > Region.RightRefID ) 
-            return AFTER_REGION;
-    }
-  
-    // --------------------------------------------------------
-    // no right bound given OR read starts before right bound
-    // so, check if it overlaps left bound 
-  
-    // if read starts on left bound reference AND after left boundary, return success
-    if ( bAlignment.RefID == Region.LeftRefID && bAlignment.Position >= Region.LeftPosition)
-        return WITHIN_REGION;
-  
-    // if read is on any reference sequence before left bound, return false
-    if ( bAlignment.RefID < Region.LeftRefID )
-        return BEFORE_REGION;
-
-    // --------------------------------------------------------
-    // read is on left bound reference, but starts before left bound position
-
-    // if it overlaps, return WITHIN_REGION
-    if ( bAlignment.GetEndPosition() >= Region.LeftPosition )
-        return WITHIN_REGION;
-    // else begins before left bound position
-    else
-        return BEFORE_REGION;
-}
-
-// load BAM header data
-void BamReader::BamReaderPrivate::LoadHeaderData(void) {
-
-    // check to see if proper BAM header
-    char buffer[4];
-    if (mBGZF.Read(buffer, 4) != 4) {
-        fprintf(stderr, "Could not read header type\n");
-        exit(1);
-    }
-
-    if (strncmp(buffer, "BAM\001", 4)) {
-        fprintf(stderr, "wrong header type!\n");
-        exit(1);
-    }
-
-    // get BAM header text length
-    mBGZF.Read(buffer, 4);
-    unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);
-    if ( IsBigEndian ) SwapEndian_32(headerTextLength); 
-    
-    // get BAM header text
-    char* headerText = (char*)calloc(headerTextLength + 1, 1);
-    mBGZF.Read(headerText, headerTextLength);
-    HeaderText = (string)((const char*)headerText);
-
-    // clean up calloc-ed temp variable
-    free(headerText);
-}
-
-// load existing index data from BAM index file (".bti" OR ".bai"), return success/fail
-bool BamReader::BamReaderPrivate::LoadIndex(const bool lookForIndex, const bool preferStandardIndex) {
-
-    // clear out any existing index data
-    ClearIndex();
-
-    // if no index filename provided, so we need to look for available index files
-    if ( IndexFilename.empty() ) {
-      
-        // attempt to load BamIndex based on current Filename provided & preferStandardIndex flag
-        const BamIndex::PreferredIndexType type = (preferStandardIndex ? BamIndex::STANDARD : BamIndex::BAMTOOLS);
-        Index = BamIndex::FromBamFilename(Filename, &mBGZF, Parent, IsBigEndian, type);
-        
-        // if null, return failure
-        if ( Index == 0 ) return false;
-        
-        // generate proper IndexFilename based on type of index created
-        IndexFilename = Filename + Index->Extension();
-    }
-    
-    else {
-        // attempt to load BamIndex based on IndexFilename provided by client
-        Index = BamIndex::FromIndexFilename(IndexFilename, &mBGZF, Parent, IsBigEndian);
-        
-        // if null, return failure
-        if ( Index == 0 ) return false; 
-    }
-    
-    // an index file was found
-    // return success of loading the index data from file
-    IsIndexLoaded = Index->Load(IndexFilename);
-    return IsIndexLoaded;
-}
-
-// populates BamAlignment with alignment data under file pointer, returns success/fail
-bool BamReader::BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {
-
-    // read in the 'block length' value, make sure it's not zero
-    char buffer[4];
-    mBGZF.Read(buffer, 4);
-    bAlignment.SupportData.BlockLength = BgzfData::UnpackUnsignedInt(buffer);
-    if ( IsBigEndian ) { SwapEndian_32(bAlignment.SupportData.BlockLength); }
-    if ( bAlignment.SupportData.BlockLength == 0 ) return false;
-
-    // read in core alignment data, make sure the right size of data was read
-    char x[BAM_CORE_SIZE];
-    if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) return false; 
-
-    if ( IsBigEndian ) {
-        for ( int i = 0; i < BAM_CORE_SIZE; i+=sizeof(uint32_t) ) 
-            SwapEndian_32p(&x[i]); 
-    }
-    
-    // set BamAlignment 'core' and 'support' data
-    bAlignment.RefID    = BgzfData::UnpackSignedInt(&x[0]);  
-    bAlignment.Position = BgzfData::UnpackSignedInt(&x[4]);
-    
-    unsigned int tempValue = BgzfData::UnpackUnsignedInt(&x[8]);
-    bAlignment.Bin        = tempValue >> 16;
-    bAlignment.MapQuality = tempValue >> 8 & 0xff;
-    bAlignment.SupportData.QueryNameLength = tempValue & 0xff;
-
-    tempValue = BgzfData::UnpackUnsignedInt(&x[12]);
-    bAlignment.AlignmentFlag = tempValue >> 16;
-    bAlignment.SupportData.NumCigarOperations = tempValue & 0xffff;
-
-    bAlignment.SupportData.QuerySequenceLength = BgzfData::UnpackUnsignedInt(&x[16]);
-    bAlignment.MateRefID    = BgzfData::UnpackSignedInt(&x[20]);
-    bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[24]);
-    bAlignment.InsertSize   = BgzfData::UnpackSignedInt(&x[28]);
-    
-    // set BamAlignment length
-    bAlignment.Length = bAlignment.SupportData.QuerySequenceLength;
-    
-    // read in character data - make sure proper data size was read
-    bool readCharDataOK = false;
-    const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
-    char* allCharData = (char*)calloc(sizeof(char), dataLength);
-    
-    if ( mBGZF.Read(allCharData, dataLength) == (signed int)dataLength) { 
-      
-        // store 'allCharData' in supportData structure
-        bAlignment.SupportData.AllCharData.assign((const char*)allCharData, dataLength);
-        
-        // set success flag
-        readCharDataOK = true;
-        
-        // save CIGAR ops 
-        // need to calculate this here so that  BamAlignment::GetEndPosition() performs correctly, 
-        // even when BamReader::GetNextAlignmentCore() is called 
-        const unsigned int cigarDataOffset = bAlignment.SupportData.QueryNameLength;
-        uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);
-        CigarOp op;
-        bAlignment.CigarData.clear();
-        bAlignment.CigarData.reserve(bAlignment.SupportData.NumCigarOperations);
-        for (unsigned int i = 0; i < bAlignment.SupportData.NumCigarOperations; ++i) {
-
-            // swap if necessary
-            if ( IsBigEndian ) SwapEndian_32(cigarData[i]);
-          
-            // build CigarOp structure
-            op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);
-            op.Type   = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];
-
-            // save CigarOp
-            bAlignment.CigarData.push_back(op);
-        }
-    }
-
-    free(allCharData);
-    return readCharDataOK;
-}
-
-// loads reference data from BAM file
-void BamReader::BamReaderPrivate::LoadReferenceData(void) {
-
-    // get number of reference sequences
-    char buffer[4];
-    mBGZF.Read(buffer, 4);
-    unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);
-    if ( IsBigEndian ) SwapEndian_32(numberRefSeqs);
-    if ( numberRefSeqs == 0 ) return;
-    References.reserve((int)numberRefSeqs);
-
-    // iterate over all references in header
-    for (unsigned int i = 0; i != numberRefSeqs; ++i) {
-
-        // get length of reference name
-        mBGZF.Read(buffer, 4);
-        unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);
-        if ( IsBigEndian ) SwapEndian_32(refNameLength);
-        char* refName = (char*)calloc(refNameLength, 1);
-
-        // get reference name and reference sequence length
-        mBGZF.Read(refName, refNameLength);
-        mBGZF.Read(buffer, 4);
-        int refLength = BgzfData::UnpackSignedInt(buffer);
-        if ( IsBigEndian ) SwapEndian_32(refLength); 
-
-        // store data for reference
-        RefData aReference;
-        aReference.RefName   = (string)((const char*)refName);
-        aReference.RefLength = refLength;
-        References.push_back(aReference);
-
-        // clean up calloc-ed temp variable
-        free(refName);
-    }
-}
-
-// opens BAM file (and index)
-bool BamReader::BamReaderPrivate::Open(const string& filename, const string& indexFilename, const bool lookForIndex, const bool preferStandardIndex) {
-
-    // store filenames
-    Filename = filename;
-    IndexFilename = indexFilename;
-
-    // open the BGZF file for reading, return false on failure
-    if ( !mBGZF.Open(filename, "rb") ) return false; 
-    
-    // retrieve header text & reference data
-    LoadHeaderData();
-    LoadReferenceData();
-
-    // store file offset of first alignment
-    AlignmentsBeginOffset = mBGZF.Tell();
-
-    // if no index filename provided 
-    if ( IndexFilename.empty() ) {
-     
-        // client did not specify that index SHOULD be found
-        // useful for cases where sequential access is all that is required
-        if ( !lookForIndex ) return true; 
-          
-        // otherwise, look for index file, return success/fail
-        return LoadIndex(lookForIndex, preferStandardIndex) ;
-    }
-    
-    // client supplied an index filename
-    // attempt to load index data, return success/fail
-    return LoadIndex(lookForIndex, preferStandardIndex);    
-}
-
-// returns BAM file pointer to beginning of alignment data
-bool BamReader::BamReaderPrivate::Rewind(void) {
-   
-    // rewind to first alignment, return false if unable to seek
-    if ( !mBGZF.Seek(AlignmentsBeginOffset) ) return false;
-  
-    // retrieve first alignment data, return false if unable to read
-    BamAlignment al;
-    if ( !LoadNextAlignment(al) ) return false;
-      
-    // reset default region info using first alignment in file
-    Region.clear();
-    Region.LeftRefID    = al.RefID;
-    Region.LeftPosition = al.Position;
-    HasAlignmentsInRegion = true;
-
-    // rewind back to beginning of first alignment
-    // return success/fail of seek
-    return mBGZF.Seek(AlignmentsBeginOffset);
-}
-
-// asks Index to attempt a Jump() to specified region
-// returns success/failure
-bool BamReader::BamReaderPrivate::SetRegion(const BamRegion& region) {
-      
-    // clear out any prior BamReader region data
-    //
-    // N.B. - this is cleared so that BamIndex now has free reign to call
-    // GetNextAlignmentCore() and do overlap checking without worrying about BamReader 
-    // performing any overlap checking of its own and moving on to the next read... Calls 
-    // to GetNextAlignmentCore() with no Region set, simply return the next alignment.
-    // This ensures that the Index is able to do just that. (All without exposing 
-    // LoadNextAlignment() to the public API, and potentially confusing clients with the nomenclature)
-    Region.clear();
-  
-    // check for existing index 
-    if ( !IsIndexLoaded || Index == 0 ) return false; 
-    
-    // attempt jump to user-specified region return false if jump could not be performed at all
-    // (invalid index, unknown reference, etc)
-    //
-    // Index::Jump() is allowed to modify the HasAlignmentsInRegion flag
-    //  * This covers case where a region is requested that lies beyond the last alignment on a reference
-    //    If this occurs, any subsequent calls to GetNexAlignment[Core] simply return false
-    //    BamMultiReader is then able to successfully pull alignments from a region from multiple files
-    //    even if one or more have no data.
-    if ( !Index->Jump(region, &HasAlignmentsInRegion) ) return false;
-      
-    // if jump successful, save region data & return success
-    Region = region;
-    return true;
-}
+// ***************************************************************************
+// BamReader.cpp (c) 2009 Derek Barnett, Michael Str�mberg
+// Marth Lab, Department of Biology, Boston College
+// All rights reserved.
+// ---------------------------------------------------------------------------
+// Last modified: 9 October 2010 (DB)
+// ---------------------------------------------------------------------------
+// Uses BGZF routines were adapted from the bgzf.c code developed at the Broad
+// Institute.
+// ---------------------------------------------------------------------------
+// Provides the basic functionality for reading BAM files
+// ***************************************************************************
+
+// C++ includes
+#include <algorithm>
+#include <iterator>
+#include <string>
+#include <vector>
+#include <iostream>
+#include "BGZF.h"
+#include "BamReader.h"
+#include "BamIndex.h"
+using namespace BamTools;
+using namespace std;
+
+struct BamReader::BamReaderPrivate {
+
+    // -------------------------------
+    // structs, enums, typedefs
+    // -------------------------------
+    enum RegionState { BEFORE_REGION = 0
+                     , WITHIN_REGION
+                     , AFTER_REGION
+                     };
+
+    // -------------------------------
+    // data members
+    // -------------------------------
+
+    // general file data
+    BgzfData  mBGZF;
+    string    HeaderText;
+    BamIndex* Index;
+    RefVector References;
+    bool      IsIndexLoaded;
+    int64_t   AlignmentsBeginOffset;
+    string    Filename;
+    string    IndexFilename;
+    
+    // system data
+    bool IsBigEndian;
+
+    // user-specified region values
+    BamRegion Region;
+    bool HasAlignmentsInRegion;
+
+    // parent BamReader
+    BamReader* Parent;
+    
+    // BAM character constants
+    const char* DNA_LOOKUP;
+    const char* CIGAR_LOOKUP;
+
+    // -------------------------------
+    // constructor & destructor
+    // -------------------------------
+    BamReaderPrivate(BamReader* parent);
+    ~BamReaderPrivate(void);
+
+    // -------------------------------
+    // "public" interface
+    // -------------------------------
+
+    // file operations
+    void Close(void);
+    bool Open(const std::string& filename, 
+              const std::string& indexFilename, 
+              const bool lookForIndex, 
+              const bool preferStandardIndex);
+    bool Rewind(void);
+    bool SetRegion(const BamRegion& region);
+
+    // access alignment data
+    bool GetNextAlignment(BamAlignment& bAlignment);
+    bool GetNextAlignmentCore(BamAlignment& bAlignment);
+
+    // access auxiliary data
+    int GetReferenceID(const string& refName) const;
+
+    // index operations
+    bool CreateIndex(bool useStandardIndex);
+
+    // -------------------------------
+    // internal methods
+    // -------------------------------
+
+    // *** reading alignments and auxiliary data *** //
+
+    // adjusts requested region if necessary (depending on where data actually begins)
+    void AdjustRegion(BamRegion& region);
+    // fills out character data for BamAlignment data
+    bool BuildCharData(BamAlignment& bAlignment);
+    // checks to see if alignment overlaps current region
+    RegionState IsOverlap(BamAlignment& bAlignment);
+    // retrieves header text from BAM file
+    void LoadHeaderData(void);
+    // retrieves BAM alignment under file pointer
+    bool LoadNextAlignment(BamAlignment& bAlignment);
+    // builds reference data structure from BAM file
+    void LoadReferenceData(void);
+    // mark references with 'HasAlignments' status
+    void MarkReferences(void);
+
+    // *** index file handling *** //
+
+    // clear out inernal index data structure
+    void ClearIndex(void);
+    // loads index from BAM index file
+    bool LoadIndex(const bool lookForIndex, const bool preferStandardIndex);
+};
+
+// -----------------------------------------------------
+// BamReader implementation (wrapper around BRPrivate)
+// -----------------------------------------------------
+// constructor
+BamReader::BamReader(void) {
+    d = new BamReaderPrivate(this);
+}
+
+// destructor
+BamReader::~BamReader(void) {
+    delete d;
+    d = 0;
+}
+
+// file operations
+void BamReader::Close(void) { d->Close(); }
+bool BamReader::IsIndexLoaded(void) const { return d->IsIndexLoaded; }
+bool BamReader::IsOpen(void) const { return d->mBGZF.IsOpen; }
+bool BamReader::Jump(int refID, int position)  { return d->SetRegion( BamRegion(refID, position) ); }
+bool BamReader::Open(const std::string& filename, 
+                     const std::string& indexFilename, 
+                     const bool lookForIndex, 
+                     const bool preferStandardIndex) 
+{ 
+    return d->Open(filename, indexFilename, lookForIndex, preferStandardIndex); 
+}
+bool BamReader::Rewind(void) { return d->Rewind(); }
+bool BamReader::SetRegion(const BamRegion& region) { return d->SetRegion(region); }
+bool BamReader::SetRegion(const int& leftRefID, const int& leftBound, const int& rightRefID, const int& rightBound) {
+    return d->SetRegion( BamRegion(leftRefID, leftBound, rightRefID, rightBound) );
+}
+
+// access alignment data
+bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }
+bool BamReader::GetNextAlignmentCore(BamAlignment& bAlignment) { return d->GetNextAlignmentCore(bAlignment); }
+
+// access auxiliary data
+const string BamReader::GetHeaderText(void) const { return d->HeaderText; }
+int BamReader::GetReferenceCount(void) const { return d->References.size(); }
+const RefVector& BamReader::GetReferenceData(void) const { return d->References; }
+int BamReader::GetReferenceID(const string& refName) const { return d->GetReferenceID(refName); }
+const std::string BamReader::GetFilename(void) const { return d->Filename; }
+
+// index operations
+bool BamReader::CreateIndex(bool useStandardIndex) { return d->CreateIndex(useStandardIndex); }
+
+// -----------------------------------------------------
+// BamReaderPrivate implementation
+// -----------------------------------------------------
+
+// constructor
+BamReader::BamReaderPrivate::BamReaderPrivate(BamReader* parent)
+    : Index(0)
+    , IsIndexLoaded(false)
+    , AlignmentsBeginOffset(0)
+    , HasAlignmentsInRegion(true)
+    , Parent(parent)
+    , DNA_LOOKUP("=ACMGRSVTWYHKDBN")
+    , CIGAR_LOOKUP("MIDNSHP")
+{ 
+    IsBigEndian = SystemIsBigEndian();
+}
+
+// destructor
+BamReader::BamReaderPrivate::~BamReaderPrivate(void) {
+    Close();
+}
+
+// adjusts requested region if necessary (depending on where data actually begins)
+void BamReader::BamReaderPrivate::AdjustRegion(BamRegion& region) {
+  
+    // check for valid index first
+    if ( Index == 0 ) return;
+  
+    // see if any references in region have alignments
+    HasAlignmentsInRegion = false;
+    int currentId = region.LeftRefID;
+    while ( currentId <= region.RightRefID ) { 
+	HasAlignmentsInRegion = Index->HasAlignments(currentId);
+	if ( HasAlignmentsInRegion ) break;
+	++currentId;
+    }
+    
+    // if no data found on any reference in region
+    if ( !HasAlignmentsInRegion ) return;
+
+    // if left bound of desired region had no data, use first reference that had data
+    // otherwise, leave requested region as-is
+    if ( currentId != region.LeftRefID ) {
+	region.LeftRefID = currentId;
+	region.LeftPosition = 0;
+    }
+}
+
+// fills out character data for BamAlignment data
+bool BamReader::BamReaderPrivate::BuildCharData(BamAlignment& bAlignment) {
+  
+    // calculate character lengths/offsets
+    const unsigned int dataLength      = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
+    const unsigned int seqDataOffset   = bAlignment.SupportData.QueryNameLength + (bAlignment.SupportData.NumCigarOperations * 4);
+    const unsigned int qualDataOffset  = seqDataOffset + (bAlignment.SupportData.QuerySequenceLength+1)/2;
+    const unsigned int tagDataOffset   = qualDataOffset + bAlignment.SupportData.QuerySequenceLength;
+    const unsigned int tagDataLength   = dataLength - tagDataOffset;
+      
+    // set up char buffers
+    const char*     allCharData = bAlignment.SupportData.AllCharData.data();
+    const char*     seqData     = ((const char*)allCharData) + seqDataOffset;
+    const char*     qualData    = ((const char*)allCharData) + qualDataOffset;
+          char*     tagData     = ((char*)allCharData) + tagDataOffset;
+  
+    // store alignment name (relies on null char in name as terminator)
+    bAlignment.Name.assign((const char*)(allCharData));    
+
+    // save query sequence
+    bAlignment.QueryBases.clear();
+    bAlignment.QueryBases.reserve(bAlignment.SupportData.QuerySequenceLength);
+    for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
+        char singleBase = DNA_LOOKUP[ ( ( seqData[(i/2)] >> (4*(1-(i%2)))) & 0xf ) ];
+        bAlignment.QueryBases.append(1, singleBase);
+    }
+  
+    // save qualities, converting from numeric QV to 'FASTQ-style' ASCII character
+    bAlignment.Qualities.clear();
+    bAlignment.Qualities.reserve(bAlignment.SupportData.QuerySequenceLength);
+    for (unsigned int i = 0; i < bAlignment.SupportData.QuerySequenceLength; ++i) {
+        char singleQuality = (char)(qualData[i]+33);
+        bAlignment.Qualities.append(1, singleQuality);
+    }
+    
+    // if QueryBases is empty (and this is a allowed case)
+    if ( bAlignment.QueryBases.empty() ) 
+        bAlignment.AlignedBases = bAlignment.QueryBases;
+    
+    // if QueryBases contains data, then build AlignedBases using CIGAR data
+    else {
+    
+        // resize AlignedBases
+        bAlignment.AlignedBases.clear();
+        bAlignment.AlignedBases.reserve(bAlignment.SupportData.QuerySequenceLength);
+      
+        // iterate over CigarOps
+        int k = 0;
+        vector<CigarOp>::const_iterator cigarIter = bAlignment.CigarData.begin();
+        vector<CigarOp>::const_iterator cigarEnd  = bAlignment.CigarData.end();
+        for ( ; cigarIter != cigarEnd; ++cigarIter ) {
+            
+            const CigarOp& op = (*cigarIter);
+            switch(op.Type) {
+              
+                case ('M') :
+                case ('I') :
+                    bAlignment.AlignedBases.append(bAlignment.QueryBases.substr(k, op.Length)); // for 'M', 'I' - write bases
+                    // fall through
+                
+                case ('S') :
+                    k += op.Length;                                     // for 'S' - soft clip, skip over query bases
+                    break;
+                    
+                case ('D') :
+                    bAlignment.AlignedBases.append(op.Length, '-');     // for 'D' - write gap character
+                    break;
+                    
+                case ('P') :
+                    bAlignment.AlignedBases.append( op.Length, '*' );   // for 'P' - write padding character
+                    break;
+                    
+                case ('N') :
+                    bAlignment.AlignedBases.append( op.Length, 'N' );  // for 'N' - write N's, skip bases in original query sequence
+                    break;
+                    
+                case ('H') :
+                    break;  // for 'H' - hard clip, do nothing to AlignedBases, move to next op
+                    
+                default:
+                    fprintf(stderr, "ERROR: Invalid Cigar op type\n"); // shouldn't get here
+                    exit(1);
+            }
+        }
+    }
+ 
+    // -----------------------
+    // Added: 3-25-2010 DB
+    // Fixed: endian-correctness for tag data
+    // -----------------------
+    if ( IsBigEndian ) {
+        int i = 0;
+        while ( (unsigned int)i < tagDataLength ) {
+          
+            i += 2; // skip tag type (e.g. "RG", "NM", etc)
+            uint8_t type = toupper(tagData[i]);     // lower & upper case letters have same meaning 
+            ++i;                                    // skip value type
+    
+            switch (type) {
+                
+                case('A') :
+                case('C') : 
+                    ++i;
+                    break;
+
+                case('S') : 
+                    SwapEndian_16p(&tagData[i]); 
+                    i += sizeof(uint16_t);
+                    break;
+                    
+                case('F') :
+                case('I') : 
+                    SwapEndian_32p(&tagData[i]);
+                    i += sizeof(uint32_t);
+                    break;
+                
+                case('D') : 
+                    SwapEndian_64p(&tagData[i]);
+                    i += sizeof(uint64_t);
+                    break;
+                
+                case('H') :
+                case('Z') : 
+                    while (tagData[i]) { ++i; }
+                    ++i; // increment one more for null terminator
+                    break;
+                
+                default : 
+                    fprintf(stderr, "ERROR: Invalid tag value type\n"); // shouldn't get here
+                    exit(1);
+            }
+        }
+    }
+    
+    // store TagData
+    bAlignment.TagData.clear();
+    bAlignment.TagData.resize(tagDataLength);
+    memcpy((char*)bAlignment.TagData.data(), tagData, tagDataLength);
+    
+    // clear the core-only flag
+    bAlignment.SupportData.HasCoreOnly = false;
+    
+    // return success
+    return true;
+}
+
+// clear index data structure
+void BamReader::BamReaderPrivate::ClearIndex(void) {
+    delete Index;
+    Index = 0;
+    IsIndexLoaded = false;
+}
+
+// closes the BAM file
+void BamReader::BamReaderPrivate::Close(void) {
+    
+    // close BGZF file stream
+    mBGZF.Close();
+    
+    // clear out index data
+    ClearIndex();
+    
+    // clear out header data
+    HeaderText.clear();
+    
+    // clear out region flags
+    Region.clear();
+}
+
+// creates index for BAM file, saves to file
+// default behavior is to create the BAM standard index (".bai")
+// set flag to false to create the BamTools-specific index (".bti")
+bool BamReader::BamReaderPrivate::CreateIndex(bool useStandardIndex) {
+
+    // clear out prior index data
+    ClearIndex();
+    
+    // create index based on type requested
+    if ( useStandardIndex ) 
+        Index = new BamStandardIndex(&mBGZF, Parent, IsBigEndian);
+    // create BamTools 'custom' index
+    else
+        Index = new BamToolsIndex(&mBGZF, Parent, IsBigEndian);
+    
+    // build new index
+    bool ok = true;
+    ok &= Index->Build();
+    IsIndexLoaded = ok;
+    
+    // mark empty references
+    MarkReferences();
+    
+    // attempt to save index data to file
+    ok &= Index->Write(Filename); 
+    
+    // return success/fail of both building & writing index
+    return ok;
+}
+
+// get next alignment (from specified region, if given)
+bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {
+
+    // if valid alignment found, attempt to parse char data, and return success/failure
+    if ( GetNextAlignmentCore(bAlignment) )
+        return BuildCharData(bAlignment);
+    
+    // no valid alignment found
+    else return false;
+}
+
+// retrieves next available alignment core data (returns success/fail)
+// ** DOES NOT parse any character data (read name, bases, qualities, tag data)
+//    these can be accessed, if necessary, from the supportData 
+// useful for operations requiring ONLY positional or other alignment-related information
+bool BamReader::BamReaderPrivate::GetNextAlignmentCore(BamAlignment& bAlignment) {
+
+    // if region is set but has no alignments
+    if ( !Region.isNull() && !HasAlignmentsInRegion )
+        return false;
+  
+    // if valid alignment available
+    if ( LoadNextAlignment(bAlignment) ) {
+
+        // set core-only flag
+        bAlignment.SupportData.HasCoreOnly = true;
+      
+	// if region not specified with at least a left boundary, return success
+	if ( !Region.isLeftBoundSpecified() ) return true;
+	
+        // determine region state (before, within, after)
+        BamReader::BamReaderPrivate::RegionState state = IsOverlap(bAlignment);
+
+        // if alignment lies after region, return false
+        if ( state == AFTER_REGION ) return false;
+
+        while ( state != WITHIN_REGION ) {
+            // if no valid alignment available (likely EOF) return failure
+            if ( !LoadNextAlignment(bAlignment) ) return false;
+            // if alignment lies after region, return false (no available read within region)
+            state = IsOverlap(bAlignment);
+            if ( state == AFTER_REGION ) return false;
+        }
+
+        // return success (alignment found that overlaps region)
+        return true;
+    }
+
+    // no valid alignment
+    else return false;
+}
+
+// returns RefID for given RefName (returns References.size() if not found)
+int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {
+
+    // retrieve names from reference data
+    vector<string> refNames;
+    RefVector::const_iterator refIter = References.begin();
+    RefVector::const_iterator refEnd  = References.end();
+    for ( ; refIter != refEnd; ++refIter) 
+        refNames.push_back( (*refIter).RefName );
+
+    // return 'index-of' refName ( if not found, returns refNames.size() )
+    return distance(refNames.begin(), find(refNames.begin(), refNames.end(), refName));
+}
+
+// returns region state - whether alignment ends before, overlaps, or starts after currently specified region
+// this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true
+BamReader::BamReaderPrivate::RegionState BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {
+    
+    // if alignment is on any reference sequence before left bound
+    if ( bAlignment.RefID < Region.LeftRefID ) return BEFORE_REGION;
+  
+    // if alignment starts on left bound reference 
+    else if ( bAlignment.RefID == Region.LeftRefID ) {
+      
+	// if alignment starts at or after left boundary
+	if ( bAlignment.Position >= Region.LeftPosition) {
+	  
+	    // if right boundary is specified AND 
+	    // left/right boundaries are on same reference AND 
+	    // alignment starts past right boundary
+	    if ( Region.isRightBoundSpecified() && 
+	         Region.LeftRefID == Region.RightRefID && 
+	         bAlignment.Position > Region.RightPosition )
+		return AFTER_REGION;
+	    
+	    // otherwise, alignment is within region
+	    return WITHIN_REGION;
+	}
+	
+	// alignment starts before left boundary
+	else {
+	    // check if alignment overlaps left boundary
+	    if ( bAlignment.GetEndPosition() >= Region.LeftPosition ) return WITHIN_REGION;
+	    else return BEFORE_REGION;
+	}
+    }
+  
+    // alignment starts on a reference after the left bound
+    else {
+      
+	// if region has a right boundary
+	if ( Region.isRightBoundSpecified() ) {
+	  
+	    // alignment is on reference between boundaries
+	    if ( bAlignment.RefID < Region.RightRefID ) return WITHIN_REGION;
+	  
+	    // alignment is on reference after right boundary
+	    else if ( bAlignment.RefID > Region.RightRefID ) return AFTER_REGION;
+	  
+	    // alignment is on right bound reference
+	    else {
+		// check if alignment starts before or at right boundary
+		if ( bAlignment.Position <= Region.RightPosition ) return WITHIN_REGION;		
+		else return AFTER_REGION;
+	    }
+	}
+      
+	// otherwise, alignment is after left bound reference, but there is no right boundary
+	else return WITHIN_REGION;
+    }
+}
+
+// load BAM header data
+void BamReader::BamReaderPrivate::LoadHeaderData(void) {
+
+    // check to see if proper BAM header
+    char buffer[4];
+    if (mBGZF.Read(buffer, 4) != 4) {
+        fprintf(stderr, "Could not read header type\n");
+        exit(1);
+    }
+
+    if (strncmp(buffer, "BAM\001", 4)) {
+        fprintf(stderr, "wrong header type!\n");
+        exit(1);
+    }
+
+    // get BAM header text length
+    mBGZF.Read(buffer, 4);
+    unsigned int headerTextLength = BgzfData::UnpackUnsignedInt(buffer);
+    if ( IsBigEndian ) SwapEndian_32(headerTextLength); 
+    
+    // get BAM header text
+    char* headerText = (char*)calloc(headerTextLength + 1, 1);
+    mBGZF.Read(headerText, headerTextLength);
+    HeaderText = (string)((const char*)headerText);
+
+    // clean up calloc-ed temp variable
+    free(headerText);
+}
+
+// load existing index data from BAM index file (".bti" OR ".bai"), return success/fail
+bool BamReader::BamReaderPrivate::LoadIndex(const bool lookForIndex, const bool preferStandardIndex) {
+
+    // clear out any existing index data
+    ClearIndex();
+
+    // if no index filename provided, so we need to look for available index files
+    if ( IndexFilename.empty() ) {
+      
+        // attempt to load BamIndex based on current Filename provided & preferStandardIndex flag
+        const BamIndex::PreferredIndexType type = (preferStandardIndex ? BamIndex::STANDARD : BamIndex::BAMTOOLS);
+        Index = BamIndex::FromBamFilename(Filename, &mBGZF, Parent, IsBigEndian, type);
+        
+        // if null, return failure
+        if ( Index == 0 ) return false;
+        
+        // generate proper IndexFilename based on type of index created
+        IndexFilename = Filename + Index->Extension();
+    }
+    
+    else {
+      
+        // attempt to load BamIndex based on IndexFilename provided by client
+        Index = BamIndex::FromIndexFilename(IndexFilename, &mBGZF, Parent, IsBigEndian);
+        
+        // if null, return failure
+        if ( Index == 0 ) return false;
+    }
+    
+    // an index file was found
+    // return success of loading the index data from file
+    IsIndexLoaded = Index->Load(IndexFilename);
+    
+    // mark empty references
+    MarkReferences();
+    
+    // return index status
+    return IsIndexLoaded;
+}
+
+// populates BamAlignment with alignment data under file pointer, returns success/fail
+bool BamReader::BamReaderPrivate::LoadNextAlignment(BamAlignment& bAlignment) {
+
+    // read in the 'block length' value, make sure it's not zero
+    char buffer[4];
+    mBGZF.Read(buffer, 4);
+    bAlignment.SupportData.BlockLength = BgzfData::UnpackUnsignedInt(buffer);
+    if ( IsBigEndian ) { SwapEndian_32(bAlignment.SupportData.BlockLength); }
+    if ( bAlignment.SupportData.BlockLength == 0 ) return false;
+
+    // read in core alignment data, make sure the right size of data was read
+    char x[BAM_CORE_SIZE];
+    if ( mBGZF.Read(x, BAM_CORE_SIZE) != BAM_CORE_SIZE ) return false; 
+
+    if ( IsBigEndian ) {
+        for ( int i = 0; i < BAM_CORE_SIZE; i+=sizeof(uint32_t) ) 
+            SwapEndian_32p(&x[i]); 
+    }
+    
+    // set BamAlignment 'core' and 'support' data
+    bAlignment.RefID    = BgzfData::UnpackSignedInt(&x[0]);  
+    bAlignment.Position = BgzfData::UnpackSignedInt(&x[4]);
+    
+    unsigned int tempValue = BgzfData::UnpackUnsignedInt(&x[8]);
+    bAlignment.Bin        = tempValue >> 16;
+    bAlignment.MapQuality = tempValue >> 8 & 0xff;
+    bAlignment.SupportData.QueryNameLength = tempValue & 0xff;
+
+    tempValue = BgzfData::UnpackUnsignedInt(&x[12]);
+    bAlignment.AlignmentFlag = tempValue >> 16;
+    bAlignment.SupportData.NumCigarOperations = tempValue & 0xffff;
+
+    bAlignment.SupportData.QuerySequenceLength = BgzfData::UnpackUnsignedInt(&x[16]);
+    bAlignment.MateRefID    = BgzfData::UnpackSignedInt(&x[20]);
+    bAlignment.MatePosition = BgzfData::UnpackSignedInt(&x[24]);
+    bAlignment.InsertSize   = BgzfData::UnpackSignedInt(&x[28]);
+    
+    // set BamAlignment length
+    bAlignment.Length = bAlignment.SupportData.QuerySequenceLength;
+    
+    // read in character data - make sure proper data size was read
+    bool readCharDataOK = false;
+    const unsigned int dataLength = bAlignment.SupportData.BlockLength - BAM_CORE_SIZE;
+    char* allCharData = (char*)calloc(sizeof(char), dataLength);
+    
+    if ( mBGZF.Read(allCharData, dataLength) == (signed int)dataLength) { 
+      
+        // store 'allCharData' in supportData structure
+        bAlignment.SupportData.AllCharData.assign((const char*)allCharData, dataLength);
+        
+        // set success flag
+        readCharDataOK = true;
+        
+        // save CIGAR ops 
+        // need to calculate this here so that  BamAlignment::GetEndPosition() performs correctly, 
+        // even when BamReader::GetNextAlignmentCore() is called 
+        const unsigned int cigarDataOffset = bAlignment.SupportData.QueryNameLength;
+        uint32_t* cigarData = (uint32_t*)(allCharData + cigarDataOffset);
+        CigarOp op;
+        bAlignment.CigarData.clear();
+        bAlignment.CigarData.reserve(bAlignment.SupportData.NumCigarOperations);
+        for (unsigned int i = 0; i < bAlignment.SupportData.NumCigarOperations; ++i) {
+
+            // swap if necessary
+            if ( IsBigEndian ) SwapEndian_32(cigarData[i]);
+          
+            // build CigarOp structure
+            op.Length = (cigarData[i] >> BAM_CIGAR_SHIFT);
+            op.Type   = CIGAR_LOOKUP[ (cigarData[i] & BAM_CIGAR_MASK) ];
+
+            // save CigarOp
+            bAlignment.CigarData.push_back(op);
+        }
+    }
+
+    free(allCharData);
+    return readCharDataOK;
+}
+
+// loads reference data from BAM file
+void BamReader::BamReaderPrivate::LoadReferenceData(void) {
+
+    // get number of reference sequences
+    char buffer[4];
+    mBGZF.Read(buffer, 4);
+    unsigned int numberRefSeqs = BgzfData::UnpackUnsignedInt(buffer);
+    if ( IsBigEndian ) SwapEndian_32(numberRefSeqs);
+    if ( numberRefSeqs == 0 ) return;
+    References.reserve((int)numberRefSeqs);
+
+    // iterate over all references in header
+    for (unsigned int i = 0; i != numberRefSeqs; ++i) {
+
+        // get length of reference name
+        mBGZF.Read(buffer, 4);
+        unsigned int refNameLength = BgzfData::UnpackUnsignedInt(buffer);
+        if ( IsBigEndian ) SwapEndian_32(refNameLength);
+        char* refName = (char*)calloc(refNameLength, 1);
+
+        // get reference name and reference sequence length
+        mBGZF.Read(refName, refNameLength);
+        mBGZF.Read(buffer, 4);
+        int refLength = BgzfData::UnpackSignedInt(buffer);
+        if ( IsBigEndian ) SwapEndian_32(refLength); 
+
+        // store data for reference
+        RefData aReference;
+        aReference.RefName   = (string)((const char*)refName);
+        aReference.RefLength = refLength;
+        References.push_back(aReference);
+
+        // clean up calloc-ed temp variable
+        free(refName);
+    }
+}
+
+// mark references with no alignment data
+void BamReader::BamReaderPrivate::MarkReferences(void) {
+    
+    // ensure index is available
+    if ( Index == 0 || !IsIndexLoaded ) return;
+    
+    // mark empty references
+    for ( int i = 0; i < (int)References.size(); ++i ) 
+	References.at(i).RefHasAlignments = Index->HasAlignments(i);
+}
+
+// opens BAM file (and index)
+bool BamReader::BamReaderPrivate::Open(const string& filename, const string& indexFilename, const bool lookForIndex, const bool preferStandardIndex) {
+
+    // store filenames
+    Filename = filename;
+    IndexFilename = indexFilename;
+
+    // open the BGZF file for reading, return false on failure
+    if ( !mBGZF.Open(filename, "rb") ) return false; 
+    
+    // retrieve header text & reference data
+    LoadHeaderData();
+    LoadReferenceData();
+
+    // store file offset of first alignment
+    AlignmentsBeginOffset = mBGZF.Tell();
+
+    // if no index filename provided 
+    if ( IndexFilename.empty() ) {
+     
+        // client did not specify that index SHOULD be found
+        // useful for cases where sequential access is all that is required
+        if ( !lookForIndex ) return true; 
+          
+        // otherwise, look for index file, return success/fail
+        return LoadIndex(lookForIndex, preferStandardIndex) ;
+    }
+    
+    // client supplied an index filename
+    // attempt to load index data, return success/fail
+    return LoadIndex(lookForIndex, preferStandardIndex);    
+}
+
+// returns BAM file pointer to beginning of alignment data
+bool BamReader::BamReaderPrivate::Rewind(void) {
+   
+    // rewind to first alignment, return false if unable to seek
+    if ( !mBGZF.Seek(AlignmentsBeginOffset) ) return false;
+  
+    // retrieve first alignment data, return false if unable to read
+    BamAlignment al;
+    if ( !LoadNextAlignment(al) ) return false;
+      
+    // reset default region info using first alignment in file
+    Region.clear();
+    HasAlignmentsInRegion = true;
+
+    // rewind back to beginning of first alignment
+    // return success/fail of seek
+    return mBGZF.Seek(AlignmentsBeginOffset);
+}
+
+// asks Index to attempt a Jump() to specified region
+// returns success/failure
+bool BamReader::BamReaderPrivate::SetRegion(const BamRegion& region) {
+      
+    // clear out any prior BamReader region data
+    //
+    // N.B. - this is cleared so that BamIndex now has free reign to call
+    // GetNextAlignmentCore() and do overlap checking without worrying about BamReader 
+    // performing any overlap checking of its own and moving on to the next read... Calls 
+    // to GetNextAlignmentCore() with no Region set, simply return the next alignment.
+    // This ensures that the Index is able to do just that. (All without exposing 
+    // LoadNextAlignment() to the public API, and potentially confusing clients with the nomenclature)
+    Region.clear();
+  
+    // check for existing index 
+    if ( !IsIndexLoaded || Index == 0 ) return false; 
+    
+    // adjust region if necessary to reflect where data actually begins
+    BamRegion adjustedRegion(region);
+    AdjustRegion(adjustedRegion);
+    
+    // if no data present, return true
+    // not an error, but BamReader knows that no data is there for future alignment access
+    // (this is useful in a MultiBamReader setting where some BAM files may lack data in regions
+    // that other BAMs have data)
+    if ( !HasAlignmentsInRegion ) { 
+	Region = adjustedRegion;
+	return true;
+    }
+    
+    // attempt jump to user-specified region return false if jump could not be performed at all
+    // (invalid index, unknown reference, etc)
+    //
+    // Index::Jump() is allowed to modify the HasAlignmentsInRegion flag
+    //  * This covers case where a region is requested that lies beyond the last alignment on a reference
+    //    If this occurs, any subsequent calls to GetNexAlignment[Core] simply return false
+    //    BamMultiReader is then able to successfully pull alignments from a region from multiple files
+    //    even if one or more have no data.
+    if ( !Index->Jump(adjustedRegion, &HasAlignmentsInRegion) )
+	return false;
+    
+    // save region and return success
+    Region = adjustedRegion;
+    return true;
+}
-- 
2.39.5