Added concept of a fully specified region of interest to the BamReader API. Added...

diff --git a/BamAux.h b/BamAux.h
index 6777dbc07809c64b47ff0f9af9f47806a71fffda..c079e1dacf4d20975428be3daee638934114d5e7 100644 (file)
--- a/BamAux.h
+++ b/BamAux.h
@@ -3,7 +3,7 @@
 // Marth Lab, Department of Biology, Boston College\r
 // All rights reserved.\r
 // ---------------------------------------------------------------------------\r
-// Last modified: 8 June 2010 (DB)\r
+// Last modified: 16 June 2010 (DB)\r
 // ---------------------------------------------------------------------------\r
 // Provides the basic constants, data structures, etc. for using BAM files\r
 // ***************************************************************************\r
@@ -208,6 +208,26 @@ struct RefData {
 typedef std::vector<RefData>      RefVector;\r
 typedef std::vector<BamAlignment> BamAlignmentVector;\r
 \r
+struct BamRegion {\r
+  \r
+    // data members\r
+    int LeftRefID;\r
+    int LeftPosition;\r
+    int RightRefID;\r
+    int RightPosition;\r
+    \r
+    // constructor\r
+    BamRegion(const int& leftID   = -1, \r
+              const int& leftPos  = -1,\r
+              const int& rightID  = -1,\r
+              const int& rightPos = -1)\r
+        : LeftRefID(leftID)\r
+        , LeftPosition(leftPos)\r
+        , RightRefID(rightID)\r
+        , RightPosition(rightPos)\r
+    { }\r
+};\r
+\r
 // ----------------------------------------------------------------\r
 // Indexing structs & typedefs\r
 \r

diff --git a/BamReader.cpp b/BamReader.cpp
index e5265e0ca0446318ee9e3955a4db2defbe3aefa4..e3aca5c09bd038371e8c90f6489d61af978fb700 100644 (file)
--- a/BamReader.cpp
+++ b/BamReader.cpp
@@ -3,7 +3,7 @@
 // Marth Lab, Department of Biology, Boston College\r
 // All rights reserved.\r
 // ---------------------------------------------------------------------------\r
-// Last modified: 8 June 2010 (DB)\r
+// Last modified: 16 June 2010 (DB)\r
 // ---------------------------------------------------------------------------\r
 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad\r
 // Institute.\r
@@ -43,6 +43,10 @@ struct BamReader::BamReaderPrivate {
     bool IsBigEndian;\r
 \r
     // user-specified region values\r
+    BamRegion Region;\r
+    bool IsLeftBoundSpecified;\r
+    bool IsRightBoundSpecified;\r
+    \r
     bool IsRegionSpecified;\r
     int  CurrentRefID;\r
     int  CurrentLeft;\r
@@ -66,6 +70,7 @@ struct BamReader::BamReaderPrivate {
     bool Jump(int refID, int position = 0);\r
     void Open(const string& filename, const string& indexFilename = "");\r
     bool Rewind(void);\r
+    bool SetRegion(const BamRegion& region);\r
 \r
     // access alignment data\r
     bool GetNextAlignment(BamAlignment& bAlignment);\r
@@ -133,9 +138,19 @@ BamReader::~BamReader(void) {
 \r
 // file operations\r
 void BamReader::Close(void) { d->Close(); }\r
-bool BamReader::Jump(int refID, int position) { return d->Jump(refID, position); }\r
+bool BamReader::Jump(int refID, int position) { \r
+    d->Region.LeftRefID = refID;\r
+    d->Region.LeftPosition = position;\r
+    d->IsLeftBoundSpecified = true;\r
+    d->IsRightBoundSpecified = false;\r
+    return d->Jump(refID, position); \r
+}\r
 void BamReader::Open(const string& filename, const string& indexFilename) { d->Open(filename, indexFilename); }\r
 bool BamReader::Rewind(void) { return d->Rewind(); }\r
+bool BamReader::SetRegion(const BamRegion& region) { return d->SetRegion(region); }\r
+bool BamReader::SetRegion(const int& leftRefID, const int& leftBound, const int& rightRefID, const int& rightBound) {\r
+    return d->SetRegion( BamRegion(leftRefID, leftBound, rightRefID, rightBound) );\r
+}\r
 \r
 // access alignment data\r
 bool BamReader::GetNextAlignment(BamAlignment& bAlignment) { return d->GetNextAlignment(bAlignment); }\r
@@ -159,6 +174,8 @@ bool BamReader::CreateIndex(void) { return d->CreateIndex(); }
 BamReader::BamReaderPrivate::BamReaderPrivate(void)\r
     : IsIndexLoaded(false)\r
     , AlignmentsBeginOffset(0)\r
+    , IsLeftBoundSpecified(false)\r
+    , IsRightBoundSpecified(false)\r
     , IsRegionSpecified(false)\r
     , CurrentRefID(0)\r
     , CurrentLeft(0)\r
@@ -471,6 +488,8 @@ void BamReader::BamReaderPrivate::Close(void) {
     mBGZF.Close();\r
     ClearIndex();\r
     HeaderText.clear();\r
+    IsLeftBoundSpecified = false;\r
+    IsRightBoundSpecified = false;\r
     IsRegionSpecified = false;\r
 }\r
 \r
@@ -507,28 +526,12 @@ int BamReader::BamReaderPrivate::GetReferenceID(const string& refName) const {
 // get next alignment (from specified region, if given)\r
 bool BamReader::BamReaderPrivate::GetNextAlignment(BamAlignment& bAlignment) {\r
 \r
-    // if valid alignment available\r
-    if ( LoadNextAlignment(bAlignment) ) {\r
-\r
-        // if region not specified, return success\r
-        if ( !IsRegionSpecified ) { \r
-          bool ok = BuildCharData(bAlignment);\r
-          return ok; \r
-        }\r
-\r
-        // load next alignment until region overlap is found\r
-        while ( !IsOverlap(bAlignment) ) {\r
-            // if no valid alignment available (likely EOF) return failure\r
-            if ( !LoadNextAlignment(bAlignment) ) return false;\r
-        }\r
-\r
-        // return success (alignment found that overlaps region)\r
-        bool ok = BuildCharData(bAlignment);\r
-        return ok;\r
-    }\r
-\r
-    // no valid alignment\r
-    else \r
+    // if valid alignment found, attempt to parse char data, and return success/failure\r
+    if ( GetNextAlignmentCore(bAlignment) )\r
+        return BuildCharData(bAlignment);\r
+    \r
+    // no valid alignment found\r
+    else\r
         return false;\r
 }\r
 \r
@@ -542,7 +545,7 @@ bool BamReader::BamReaderPrivate::GetNextAlignmentCore(BamAlignment& bAlignment)
     if ( LoadNextAlignment(bAlignment) ) {\r
 \r
         // if region not specified, return success\r
-        if ( !IsRegionSpecified ) return true;\r
+        if ( !IsLeftBoundSpecified ) return true;\r
 \r
         // load next alignment until region overlap is found\r
         while ( !IsOverlap(bAlignment) ) {\r
@@ -652,17 +655,40 @@ void BamReader::BamReaderPrivate::InsertLinearOffset(LinearOffsetVector& offsets
     }\r
 }\r
 \r
-// returns whether alignment overlaps currently specified region (refID, leftBound)\r
+// returns whether alignment overlaps currently specified region\r
+// this *internal* method should ONLY called when (at least) IsLeftBoundSpecified == true\r
 bool BamReader::BamReaderPrivate::IsOverlap(BamAlignment& bAlignment) {\r
+    \r
+    // --------------------------------------------------\r
+    // check alignment start against right bound cutoff\r
+    \r
+    // if full region of interest was given\r
+    if ( IsRightBoundSpecified ) {\r
+      \r
+        // read starts on right bound reference, but AFTER right bound position\r
+        if ( bAlignment.RefID == Region.RightRefID && bAlignment.Position > Region.RightPosition )\r
+            return false;\r
+      \r
+        // if read starts on reference AFTER right bound, return false\r
+        if ( bAlignment.RefID > Region.RightRefID ) \r
+            return false;\r
+    }\r
+  \r
+    // --------------------------------------------------------\r
+    // no right bound given OR read starts before right bound\r
+    // so, check if it overlaps left bound \r
+  \r
+    // if read starts on left bound reference AND after left boundary, return success\r
+    if ( bAlignment.RefID == Region.LeftRefID && bAlignment.Position >= Region.LeftPosition)\r
+        return true;\r
+  \r
+    // if read is on any reference sequence before left bound, return false\r
+    if ( bAlignment.RefID < Region.LeftRefID )\r
+        return false;\r
 \r
-    // if on different reference sequence, quit\r
-    if ( bAlignment.RefID != CurrentRefID ) { return false; }\r
-\r
-    // read starts after left boundary\r
-    if ( bAlignment.Position >= CurrentLeft) { return true; }\r
-\r
-    // return whether alignment end overlaps left boundary\r
-    return ( bAlignment.GetEndPosition() >= CurrentLeft );\r
+    // read is on left bound reference, but starts before left bound position\r
+    // calculate the alignment end to see if it overlaps\r
+    return ( bAlignment.GetEndPosition() >= Region.LeftPosition );\r
 }\r
 \r
 // jumps to specified region(refID, leftBound) in BAM file, returns success/fail\r
@@ -671,19 +697,15 @@ bool BamReader::BamReaderPrivate::Jump(int refID, int position) {
     // if data exists for this reference and position is valid    \r
     if ( References.at(refID).RefHasAlignments && (position <= References.at(refID).RefLength) ) {\r
 \r
-        // set current region\r
-        CurrentRefID = refID;\r
-        CurrentLeft  = position;\r
-        IsRegionSpecified = true;\r
-\r
         // calculate offset\r
-        int64_t offset = GetOffset(CurrentRefID, CurrentLeft);\r
+        int64_t offset = GetOffset(Region.LeftRefID, Region.LeftPosition);\r
 \r
         // if in valid offset, return failure\r
-        if ( offset == -1 ) { return false; }\r
-\r
         // otherwise return success of seek operation\r
-        else { return mBGZF.Seek(offset); }\r
+        if ( offset == -1 ) \r
+            return false;\r
+        else \r
+            return mBGZF.Seek(offset);\r
     }\r
 \r
     // invalid jump request parameters, return failure\r
@@ -1048,15 +1070,36 @@ bool BamReader::BamReaderPrivate::Rewind(void) {
         if ( References.at(refID).RefHasAlignments ) { break; }\r
     }\r
 \r
-    // store default bounds for first alignment\r
-    CurrentRefID = refID;\r
-    CurrentLeft = 0;\r
-    IsRegionSpecified = false;\r
+    // reset default region info\r
+    Region.LeftRefID = refID;\r
+    Region.LeftPosition = 0;\r
+    Region.RightRefID = -1;\r
+    Region.RightPosition = -1;\r
+    IsLeftBoundSpecified = false;\r
+    IsRightBoundSpecified = false;\r
 \r
     // return success/failure of seek\r
     return mBGZF.Seek(AlignmentsBeginOffset);\r
 }\r
 \r
+// sets a region of interest (with left & right bound reference/position)\r
+// attempts a Jump() to left bound as well\r
+// returns success/failure of Jump()\r
+bool BamReader::BamReaderPrivate::SetRegion(const BamRegion& region) {\r
+    \r
+    // save region of interest\r
+    Region = region;\r
+    \r
+    // set flags\r
+    if ( region.LeftRefID >= 0 && region.LeftPosition >= 0 ) \r
+        IsLeftBoundSpecified = true;\r
+    if ( region.RightRefID >= 0 && region.RightPosition >= 0 ) \r
+        IsRightBoundSpecified = true;\r
+    \r
+    // attempt jump to beginning of region, return success/fail of Jump()\r
+    return Jump( Region.LeftRefID, Region.LeftPosition );\r
+}\r
+\r
 // saves index data to BAM index file (".bai"), returns success/fail\r
 bool BamReader::BamReaderPrivate::WriteIndex(void) {\r
 \r

diff --git a/BamReader.h b/BamReader.h
index 0b13786b1b112a48cccae944ac73f14c4f7e6e53..fc9a00354b5fa031e685ac78343ceac6df373b2c 100644 (file)
--- a/BamReader.h
+++ b/BamReader.h
@@ -3,7 +3,7 @@
 // Marth Lab, Department of Biology, Boston College\r
 // All rights reserved.\r
 // ---------------------------------------------------------------------------\r
-// Last modified: 8 June 2010 (DB)\r
+// Last modified: 16 June 2010 (DB)\r
 // ---------------------------------------------------------------------------\r
 // Uses BGZF routines were adapted from the bgzf.c code developed at the Broad\r
 // Institute.\r
@@ -21,7 +21,7 @@
 #include "BamAux.h"\r
 \r
 namespace BamTools {\r
-\r
+  \r
 class BamReader {\r
 \r
     // constructor / destructor\r
@@ -44,6 +44,11 @@ class BamReader {
         void Open(const std::string& filename, const std::string& indexFilename = "");\r
         // returns file pointer to beginning of alignments\r
         bool Rewind(void);\r
+        // sets a region of interest (with left & right bound reference/position)\r
+        // attempts a Jump() to left bound as well\r
+        // returns success/failure of Jump()\r
+        bool SetRegion(const BamRegion& region);\r
+        bool SetRegion(const int& leftRefID, const int& leftBound, const int& rightRefID, const int& rightBound);\r
 \r
         // ----------------------\r
         // access alignment data\r