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