src/api/internal/BamRandomAccessController_p.cpp

   1 // ***************************************************************************
   2 // BamRandomAccessController_p.cpp (c) 2011 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // ---------------------------------------------------------------------------
   5 // Last modified: 5 April 2011(DB)
   6 // ---------------------------------------------------------------------------
   7 // Manages random access operations in a BAM file
   8 // **************************************************************************
   9
  10 #include <api/BamIndex.h>
  11 #include <api/internal/BamRandomAccessController_p.h>
  12 #include <api/internal/BamReader_p.h>
  13 #include <api/internal/BamIndexFactory_p.h>
  14 using namespace BamTools;
  15 using namespace BamTools::Internal;
  16
  17 #include <iostream>
  18 using namespace std;
  19
  20 BamRandomAccessController::BamRandomAccessController(void)
  21     : m_index(0)
  22     , m_indexCacheMode(BamIndex::LimitedIndexCaching)
  23     , m_hasAlignmentsInRegion(true)
  24 { }
  25
  26 BamRandomAccessController::~BamRandomAccessController(void) {
  27     Close();
  28 }
  29
  30 void BamRandomAccessController::AdjustRegion(const int& referenceCount) {
  31
  32     // skip if no index available
  33     if ( m_index == 0 )
  34         return;
  35
  36     // see if any references in region have alignments
  37     m_hasAlignmentsInRegion = false;
  38     int currentId = m_region.LeftRefID;
  39     const int rightBoundRefId = ( m_region.isRightBoundSpecified() ? m_region.RightRefID : referenceCount - 1 );
  40     while ( currentId <= rightBoundRefId ) {
  41         m_hasAlignmentsInRegion = m_index->HasAlignments(currentId);
  42         if ( m_hasAlignmentsInRegion ) break;
  43         ++currentId;
  44     }
  45
  46     // if no data found on any reference in region
  47     if ( !m_hasAlignmentsInRegion )
  48         return;
  49
  50     // if left bound of desired region had no data, use first reference that had data
  51     // otherwise, leave requested region as-is
  52     if ( currentId != m_region.LeftRefID ) {
  53         m_region.LeftRefID = currentId;
  54         m_region.LeftPosition = 0;
  55     }
  56 }
  57
  58 // returns alignments' "RegionState": { Before|Overlaps|After } current region
  59 BamRandomAccessController::RegionState
  60 BamRandomAccessController::AlignmentState(const BamAlignment& alignment) const {
  61
  62     // if region has no left bound at all
  63     if ( !m_region.isLeftBoundSpecified() )
  64         return OverlapsRegion;
  65
  66     // handle unmapped reads - return AFTER region to halt processing
  67     if ( alignment.RefID == -1 )
  68         return AfterRegion;
  69
  70     // if alignment is on any reference before left bound reference
  71     if ( alignment.RefID < m_region.LeftRefID )
  72         return BeforeRegion;
  73
  74     // if alignment is on left bound reference
  75     else if ( alignment.RefID == m_region.LeftRefID ) {
  76
  77         // if alignment starts at or after left bound position
  78         if ( alignment.Position >= m_region.LeftPosition) {
  79
  80             if ( m_region.isRightBoundSpecified() &&            // right bound is specified AND
  81                  m_region.LeftRefID == m_region.RightRefID &&   // left & right bounds on same reference AND
  82                  alignment.Position > m_region.RightPosition )  // alignment starts after right bound position
  83                 return AfterRegion;
  84
  85             // otherwise, alignment overlaps region
  86             else return OverlapsRegion;
  87         }
  88
  89         // alignment starts before left bound position
  90         else {
  91
  92             // if alignment overlaps left bound position
  93             if ( alignment.GetEndPosition() >= m_region.LeftPosition )
  94                 return OverlapsRegion;
  95             else
  96                 return BeforeRegion;
  97         }
  98     }
  99
 100     // otherwise alignment is on a reference after left bound reference
 101     else {
 102
 103         // if region has a right bound
 104         if ( m_region.isRightBoundSpecified() ) {
 105
 106             // alignment is on any reference between boundaries
 107             if ( alignment.RefID < m_region.RightRefID )
 108                 return OverlapsRegion;
 109
 110             // alignment is on any reference after right boundary
 111             else if ( alignment.RefID > m_region.RightRefID )
 112                 return AfterRegion;
 113
 114             // alignment is on right bound reference
 115             else {
 116
 117                 // if alignment starts on or before right bound position
 118                 if ( alignment.Position <= m_region.RightPosition )
 119                     return OverlapsRegion;
 120                 else
 121                     return AfterRegion;
 122             }
 123         }
 124
 125         // otherwise, alignment starts after left bound and there is no right bound
 126         else return OverlapsRegion;
 127     }
 128 }
 129
 130 void BamRandomAccessController::Close(void) {
 131     ClearIndex();
 132     ClearRegion();
 133 }
 134
 135 void BamRandomAccessController::ClearIndex(void) {
 136     delete m_index;
 137     m_index = 0;
 138 }
 139
 140 void BamRandomAccessController::ClearRegion(void) {
 141     m_region.clear();
 142     m_hasAlignmentsInRegion = true;
 143 }
 144
 145 bool BamRandomAccessController::CreateIndex(BamReaderPrivate* reader,
 146                                             const BamIndex::IndexType& type) {
 147
 148     // skip if reader is invalid
 149     if ( reader == 0 )
 150         return false;
 151
 152     // create new index of requested type
 153     BamIndex* newIndex = BamIndexFactory::CreateIndexOfType(type, reader);
 154     if ( newIndex == 0 ) {
 155         cerr << "BamRandomAccessController ERROR: could not create index of type " << type << endl;
 156         return false;
 157     }
 158
 159     // attempt to build index from current BamReader file
 160     if ( !newIndex->Create() ) {
 161         cerr << "BamRandomAccessController ERROR: could not create index for BAM file: "
 162              << reader->Filename() << endl;
 163         return false;
 164     }
 165
 166     // save new index
 167     SetIndex(newIndex);
 168
 169     // set new index's cache mode & return success
 170     newIndex->SetCacheMode(m_indexCacheMode);
 171     return true;
 172 }
 173
 174 bool BamRandomAccessController::HasIndex(void) const {
 175     return ( m_index != 0 );
 176 }
 177
 178 bool BamRandomAccessController::HasRegion(void) const  {
 179     return ( !m_region.isNull() );
 180 }
 181
 182 bool BamRandomAccessController::IndexHasAlignmentsForReference(const int& refId) {
 183     return m_index->HasAlignments(refId);
 184 }
 185
 186 bool BamRandomAccessController::LocateIndex(BamReaderPrivate* reader,
 187                                             const BamIndex::IndexType& preferredType)
 188 {
 189     // look up index filename, deferring to preferredType if possible
 190     const string& indexFilename = BamIndexFactory::FindIndexFilename(reader->Filename(), preferredType);
 191
 192     // if no index file found (of any type)
 193     if ( indexFilename.empty() ) {
 194         cerr << "BamRandomAccessController WARNING: "
 195              << "could not find index file for BAM: "
 196              << reader->Filename() << endl;
 197         return false;
 198     }
 199
 200     // otherwise open & use index file that was found
 201     return OpenIndex(indexFilename, reader);
 202 }
 203
 204 bool BamRandomAccessController::OpenIndex(const string& indexFilename, BamReaderPrivate* reader) {
 205
 206     // attempt create new index of type based on filename
 207     BamIndex* index = BamIndexFactory::CreateIndexFromFilename(indexFilename, reader);
 208     if ( index == 0 ) {
 209         cerr << "BamRandomAccessController ERROR: could not create index for file: " << indexFilename << endl;
 210         return false;
 211     }
 212
 213     // set cache mode
 214     index->SetCacheMode(m_indexCacheMode);
 215
 216     // attempt to load data from index file
 217     if ( !index->Load(indexFilename) ) {
 218         cerr << "BamRandomAccessController ERROR: could not load index data from file: " << indexFilename << endl;
 219         return false;
 220     }
 221
 222     // save new index & return success
 223     SetIndex(index);
 224     return true;
 225 }
 226
 227 bool BamRandomAccessController::RegionHasAlignments(void) const {
 228     return m_hasAlignmentsInRegion;
 229 }
 230
 231 void BamRandomAccessController::SetIndex(BamIndex* index) {
 232     if ( m_index )
 233         ClearIndex();
 234     m_index = index;
 235 }
 236
 237 void BamRandomAccessController::SetIndexCacheMode(const BamIndex::IndexCacheMode& mode) {
 238     m_indexCacheMode = mode;
 239     if ( m_index )
 240         m_index->SetCacheMode(mode);
 241 }
 242
 243 bool BamRandomAccessController::SetRegion(BamReaderPrivate* reader,
 244                                           const BamRegion& region,
 245                                           const int& referenceCount)
 246 {
 247     // store region
 248     m_region = region;
 249
 250     // cannot jump when no index is available
 251     if ( !HasIndex() )
 252         return false;
 253
 254     // adjust region as necessary to reflect where data actually begins
 255     AdjustRegion(referenceCount);
 256
 257     // if no data present, return true
 258     //   * Not an error, but future attempts to access alignments in this region will not return data
 259     //     Returning true is useful in a BamMultiReader setting where some BAM files may
 260     //     lack alignments in regions where other BAMs do have data.
 261     if ( !m_hasAlignmentsInRegion )
 262         return true;
 263
 264     // return success/failure of jump to specified region,
 265     //
 266     //  * Index::Jump() is allowed to modify the m_hasAlignmentsInRegion flag
 267     //    This covers 'corner case' where a region is requested that lies beyond the last
 268     //    alignment on a reference. If this occurs, any subsequent calls to GetNextAlignment[Core]
 269     //    will not return data. BamMultiReader will still be able to successfully pull alignments
 270     //    from a region from multiple files even if one or more have no data.
 271     return m_index->Jump(m_region, &m_hasAlignmentsInRegion);
 272 }