|
6
|
1 // ***************************************************************************
|
|
|
2 // BamStandardIndex.cpp (c) 2010 Derek Barnett
|
|
|
3 // Marth Lab, Department of Biology, Boston College
|
|
|
4 // All rights reserved.
|
|
|
5 // ---------------------------------------------------------------------------
|
|
|
6 // Last modified: 22 November 2010 (DB)
|
|
|
7 // ---------------------------------------------------------------------------
|
|
|
8 // Provides index operations for the standardized BAM index format (".bai")
|
|
|
9 // ***************************************************************************
|
|
|
10
|
|
|
11 #include <BamAlignment.h>
|
|
|
12 #include <BamReader.h>
|
|
|
13 #include <BGZF.h>
|
|
|
14 #include <BamStandardIndex_p.h>
|
|
|
15 using namespace BamTools;
|
|
|
16 using namespace BamTools::Internal;
|
|
|
17
|
|
|
18 #include <cstdio>
|
|
|
19 #include <cstdlib>
|
|
|
20 #include <algorithm>
|
|
|
21 #include <iostream>
|
|
|
22 #include <map>
|
|
|
23 using namespace std;
|
|
|
24
|
|
|
25 BamStandardIndex::BamStandardIndex(BgzfData* bgzf, BamReader* reader)
|
|
|
26 : BamIndex(bgzf, reader)
|
|
|
27 , m_dataBeginOffset(0)
|
|
|
28 , m_hasFullDataCache(false)
|
|
|
29 {
|
|
|
30 m_isBigEndian = BamTools::SystemIsBigEndian();
|
|
|
31 }
|
|
|
32
|
|
|
33 BamStandardIndex::~BamStandardIndex(void) {
|
|
|
34 ClearAllData();
|
|
|
35 }
|
|
|
36
|
|
|
37 // calculate bins that overlap region
|
|
|
38 int BamStandardIndex::BinsFromRegion(const BamRegion& region,
|
|
|
39 const bool isRightBoundSpecified,
|
|
|
40 uint16_t bins[MAX_BIN])
|
|
|
41 {
|
|
|
42 // get region boundaries
|
|
|
43 uint32_t begin = (unsigned int)region.LeftPosition;
|
|
|
44 uint32_t end;
|
|
|
45
|
|
|
46 // if right bound specified AND left&right bounds are on same reference
|
|
|
47 // OK to use right bound position
|
|
|
48 if ( isRightBoundSpecified && ( region.LeftRefID == region.RightRefID ) )
|
|
|
49 end = (unsigned int)region.RightPosition;
|
|
|
50
|
|
|
51 // otherwise, use end of left bound reference as cutoff
|
|
|
52 else
|
|
|
53 end = (unsigned int)m_references.at(region.LeftRefID).RefLength - 1;
|
|
|
54
|
|
|
55 // initialize list, bin '0' always a valid bin
|
|
|
56 int i = 0;
|
|
|
57 bins[i++] = 0;
|
|
|
58
|
|
|
59 // get rest of bins that contain this region
|
|
|
60 unsigned int k;
|
|
|
61 for (k = 1 + (begin>>26); k <= 1 + (end>>26); ++k) { bins[i++] = k; }
|
|
|
62 for (k = 9 + (begin>>23); k <= 9 + (end>>23); ++k) { bins[i++] = k; }
|
|
|
63 for (k = 73 + (begin>>20); k <= 73 + (end>>20); ++k) { bins[i++] = k; }
|
|
|
64 for (k = 585 + (begin>>17); k <= 585 + (end>>17); ++k) { bins[i++] = k; }
|
|
|
65 for (k = 4681 + (begin>>14); k <= 4681 + (end>>14); ++k) { bins[i++] = k; }
|
|
|
66
|
|
|
67 // return number of bins stored
|
|
|
68 return i;
|
|
|
69 }
|
|
|
70
|
|
|
71 // creates index data (in-memory) from current reader data
|
|
|
72 bool BamStandardIndex::Build(void) {
|
|
|
73
|
|
|
74 // be sure reader & BGZF file are valid & open for reading
|
|
|
75 if ( m_reader == 0 || m_BGZF == 0 || !m_BGZF->IsOpen )
|
|
|
76 return false;
|
|
|
77
|
|
|
78 // move file pointer to beginning of alignments
|
|
|
79 m_reader->Rewind();
|
|
|
80
|
|
|
81 // get reference count, reserve index space
|
|
|
82 const int numReferences = (int)m_references.size();
|
|
|
83 m_indexData.clear();
|
|
|
84 m_hasFullDataCache = false;
|
|
|
85 SetReferenceCount(numReferences);
|
|
|
86
|
|
|
87 // sets default constant for bin, ID, offset, coordinate variables
|
|
|
88 const uint32_t defaultValue = 0xffffffffu;
|
|
|
89
|
|
|
90 // bin data
|
|
|
91 uint32_t saveBin(defaultValue);
|
|
|
92 uint32_t lastBin(defaultValue);
|
|
|
93
|
|
|
94 // reference ID data
|
|
|
95 int32_t saveRefID(defaultValue);
|
|
|
96 int32_t lastRefID(defaultValue);
|
|
|
97
|
|
|
98 // offset data
|
|
|
99 uint64_t saveOffset = m_BGZF->Tell();
|
|
|
100 uint64_t lastOffset = saveOffset;
|
|
|
101
|
|
|
102 // coordinate data
|
|
|
103 int32_t lastCoordinate = defaultValue;
|
|
|
104
|
|
|
105 BamAlignment bAlignment;
|
|
|
106 while ( m_reader->GetNextAlignmentCore(bAlignment) ) {
|
|
|
107
|
|
|
108 // change of chromosome, save ID, reset bin
|
|
|
109 if ( lastRefID != bAlignment.RefID ) {
|
|
|
110 lastRefID = bAlignment.RefID;
|
|
|
111 lastBin = defaultValue;
|
|
|
112 }
|
|
|
113
|
|
|
114 // if lastCoordinate greater than BAM position - file not sorted properly
|
|
|
115 else if ( lastCoordinate > bAlignment.Position ) {
|
|
|
116 fprintf(stderr, "BAM file not properly sorted:\n");
|
|
|
117 fprintf(stderr, "Alignment %s : %d > %d on reference (id = %d)", bAlignment.Name.c_str(),
|
|
|
118 lastCoordinate, bAlignment.Position, bAlignment.RefID);
|
|
|
119 exit(1);
|
|
|
120 }
|
|
|
121
|
|
|
122 // if valid reference && BAM bin spans some minimum cutoff (smaller bin ids span larger regions)
|
|
|
123 if ( (bAlignment.RefID >= 0) && (bAlignment.Bin < 4681) ) {
|
|
|
124
|
|
|
125 // save linear offset entry (matched to BAM entry refID)
|
|
|
126 BamStandardIndexData::iterator indexIter = m_indexData.find(bAlignment.RefID);
|
|
|
127 if ( indexIter == m_indexData.end() ) return false; // error
|
|
|
128 ReferenceIndex& refIndex = (*indexIter).second;
|
|
|
129 LinearOffsetVector& offsets = refIndex.Offsets;
|
|
|
130 SaveLinearOffset(offsets, bAlignment, lastOffset);
|
|
|
131 }
|
|
|
132
|
|
|
133 // if current BamAlignment bin != lastBin, "then possibly write the binning index"
|
|
|
134 if ( bAlignment.Bin != lastBin ) {
|
|
|
135
|
|
|
136 // if not first time through
|
|
|
137 if ( saveBin != defaultValue ) {
|
|
|
138
|
|
|
139 // save Bam bin entry
|
|
|
140 BamStandardIndexData::iterator indexIter = m_indexData.find(saveRefID);
|
|
|
141 if ( indexIter == m_indexData.end() ) return false; // error
|
|
|
142 ReferenceIndex& refIndex = (*indexIter).second;
|
|
|
143 BamBinMap& binMap = refIndex.Bins;
|
|
|
144 SaveBinEntry(binMap, saveBin, saveOffset, lastOffset);
|
|
|
145 }
|
|
|
146
|
|
|
147 // update saveOffset
|
|
|
148 saveOffset = lastOffset;
|
|
|
149
|
|
|
150 // update bin values
|
|
|
151 saveBin = bAlignment.Bin;
|
|
|
152 lastBin = bAlignment.Bin;
|
|
|
153
|
|
|
154 // update saveRefID
|
|
|
155 saveRefID = bAlignment.RefID;
|
|
|
156
|
|
|
157 // if invalid RefID, break out
|
|
|
158 if ( saveRefID < 0 ) break;
|
|
|
159 }
|
|
|
160
|
|
|
161 // make sure that current file pointer is beyond lastOffset
|
|
|
162 if ( m_BGZF->Tell() <= (int64_t)lastOffset ) {
|
|
|
163 fprintf(stderr, "Error in BGZF offsets.\n");
|
|
|
164 exit(1);
|
|
|
165 }
|
|
|
166
|
|
|
167 // update lastOffset
|
|
|
168 lastOffset = m_BGZF->Tell();
|
|
|
169
|
|
|
170 // update lastCoordinate
|
|
|
171 lastCoordinate = bAlignment.Position;
|
|
|
172 }
|
|
|
173
|
|
|
174 // save any leftover BAM data (as long as refID is valid)
|
|
|
175 if ( saveRefID >= 0 ) {
|
|
|
176 // save Bam bin entry
|
|
|
177 BamStandardIndexData::iterator indexIter = m_indexData.find(saveRefID);
|
|
|
178 if ( indexIter == m_indexData.end() ) return false; // error
|
|
|
179 ReferenceIndex& refIndex = (*indexIter).second;
|
|
|
180 BamBinMap& binMap = refIndex.Bins;
|
|
|
181 SaveBinEntry(binMap, saveBin, saveOffset, lastOffset);
|
|
|
182 }
|
|
|
183
|
|
|
184 // simplify index by merging chunks
|
|
|
185 MergeChunks();
|
|
|
186
|
|
|
187 // iterate through references in index
|
|
|
188 // sort offsets in linear offset vector
|
|
|
189 BamStandardIndexData::iterator indexIter = m_indexData.begin();
|
|
|
190 BamStandardIndexData::iterator indexEnd = m_indexData.end();
|
|
|
191 for ( int i = 0; indexIter != indexEnd; ++indexIter, ++i ) {
|
|
|
192
|
|
|
193 // get reference index data
|
|
|
194 ReferenceIndex& refIndex = (*indexIter).second;
|
|
|
195 LinearOffsetVector& offsets = refIndex.Offsets;
|
|
|
196
|
|
|
197 // sort linear offsets
|
|
|
198 sort(offsets.begin(), offsets.end());
|
|
|
199 }
|
|
|
200
|
|
|
201 // rewind file pointer to beginning of alignments, return success/fail
|
|
|
202 return m_reader->Rewind();
|
|
|
203 }
|
|
|
204
|
|
|
205 // check index file magic number, return true if OK
|
|
|
206 bool BamStandardIndex::CheckMagicNumber(void) {
|
|
|
207
|
|
|
208 // read in magic number
|
|
|
209 char magic[4];
|
|
|
210 size_t elementsRead = fread(magic, sizeof(char), 4, m_indexStream);
|
|
|
211
|
|
|
212 // compare to expected value
|
|
|
213 if ( strncmp(magic, "BAI\1", 4) != 0 ) {
|
|
|
214 fprintf(stderr, "Problem with index file - invalid format.\n");
|
|
|
215 fclose(m_indexStream);
|
|
|
216 return false;
|
|
|
217 }
|
|
|
218
|
|
|
219 // return success/failure of load
|
|
|
220 return (elementsRead == 4);
|
|
|
221 }
|
|
|
222
|
|
|
223 // clear all current index offset data in memory
|
|
|
224 void BamStandardIndex::ClearAllData(void) {
|
|
|
225 BamStandardIndexData::const_iterator indexIter = m_indexData.begin();
|
|
|
226 BamStandardIndexData::const_iterator indexEnd = m_indexData.end();
|
|
|
227 for ( ; indexIter != indexEnd; ++indexIter ) {
|
|
|
228 const int& refId = (*indexIter).first;
|
|
|
229 ClearReferenceOffsets(refId);
|
|
|
230 }
|
|
|
231 }
|
|
|
232
|
|
|
233 // clear all index offset data for desired reference
|
|
|
234 void BamStandardIndex::ClearReferenceOffsets(const int& refId) {
|
|
|
235
|
|
|
236 // look up refId, skip if not found
|
|
|
237 BamStandardIndexData::iterator indexIter = m_indexData.find(refId);
|
|
|
238 if ( indexIter == m_indexData.end() ) return ;
|
|
|
239
|
|
|
240 // clear reference data
|
|
|
241 ReferenceIndex& refEntry = (*indexIter).second;
|
|
|
242 refEntry.Bins.clear();
|
|
|
243 refEntry.Offsets.clear();
|
|
|
244
|
|
|
245 // set flag
|
|
|
246 m_hasFullDataCache = false;
|
|
|
247 }
|
|
|
248
|
|
|
249 // return file position after header metadata
|
|
|
250 const off_t BamStandardIndex::DataBeginOffset(void) const {
|
|
|
251 return m_dataBeginOffset;
|
|
|
252 }
|
|
|
253
|
|
|
254 // calculates offset(s) for a given region
|
|
|
255 bool BamStandardIndex::GetOffsets(const BamRegion& region,
|
|
|
256 const bool isRightBoundSpecified,
|
|
|
257 vector<int64_t>& offsets,
|
|
|
258 bool* hasAlignmentsInRegion)
|
|
|
259 {
|
|
|
260 // return false if leftBound refID is not found in index data
|
|
|
261 if ( m_indexData.find(region.LeftRefID) == m_indexData.end() )
|
|
|
262 return false;
|
|
|
263
|
|
|
264 // load index data for region if not already cached
|
|
|
265 if ( !IsDataLoaded(region.LeftRefID) ) {
|
|
|
266 bool loadedOk = true;
|
|
|
267 loadedOk &= SkipToReference(region.LeftRefID);
|
|
|
268 loadedOk &= LoadReference(region.LeftRefID);
|
|
|
269 if ( !loadedOk ) return false;
|
|
|
270 }
|
|
|
271
|
|
|
272 // calculate which bins overlap this region
|
|
|
273 uint16_t* bins = (uint16_t*)calloc(MAX_BIN, 2);
|
|
|
274 int numBins = BinsFromRegion(region, isRightBoundSpecified, bins);
|
|
|
275
|
|
|
276 // get bins for this reference
|
|
|
277 BamStandardIndexData::const_iterator indexIter = m_indexData.find(region.LeftRefID);
|
|
|
278 if ( indexIter == m_indexData.end() ) return false; // error
|
|
|
279 const ReferenceIndex& refIndex = (*indexIter).second;
|
|
|
280 const BamBinMap& binMap = refIndex.Bins;
|
|
|
281
|
|
|
282 // get minimum offset to consider
|
|
|
283 const LinearOffsetVector& linearOffsets = refIndex.Offsets;
|
|
|
284 const uint64_t minOffset = ( (unsigned int)(region.LeftPosition>>BAM_LIDX_SHIFT) >= linearOffsets.size() )
|
|
|
285 ? 0 : linearOffsets.at(region.LeftPosition>>BAM_LIDX_SHIFT);
|
|
|
286
|
|
|
287 // store all alignment 'chunk' starts (file offsets) for bins in this region
|
|
|
288 for ( int i = 0; i < numBins; ++i ) {
|
|
|
289
|
|
|
290 const uint16_t binKey = bins[i];
|
|
|
291 map<uint32_t, ChunkVector>::const_iterator binIter = binMap.find(binKey);
|
|
|
292 if ( (binIter != binMap.end()) && ((*binIter).first == binKey) ) {
|
|
|
293
|
|
|
294 // iterate over chunks
|
|
|
295 const ChunkVector& chunks = (*binIter).second;
|
|
|
296 std::vector<Chunk>::const_iterator chunksIter = chunks.begin();
|
|
|
297 std::vector<Chunk>::const_iterator chunksEnd = chunks.end();
|
|
|
298 for ( ; chunksIter != chunksEnd; ++chunksIter) {
|
|
|
299
|
|
|
300 // if valid chunk found, store its file offset
|
|
|
301 const Chunk& chunk = (*chunksIter);
|
|
|
302 if ( chunk.Stop > minOffset )
|
|
|
303 offsets.push_back( chunk.Start );
|
|
|
304 }
|
|
|
305 }
|
|
|
306 }
|
|
|
307
|
|
|
308 // clean up memory
|
|
|
309 free(bins);
|
|
|
310
|
|
|
311 // sort the offsets before returning
|
|
|
312 sort(offsets.begin(), offsets.end());
|
|
|
313
|
|
|
314 // set flag & return success
|
|
|
315 *hasAlignmentsInRegion = (offsets.size() != 0 );
|
|
|
316
|
|
|
317 // if cache mode set to none, dump the data we just loaded
|
|
|
318 if (m_cacheMode == BamIndex::NoIndexCaching )
|
|
|
319 ClearReferenceOffsets(region.LeftRefID);
|
|
|
320
|
|
|
321 // return succes
|
|
|
322 return true;
|
|
|
323 }
|
|
|
324
|
|
|
325 // returns whether reference has alignments or no
|
|
|
326 bool BamStandardIndex::HasAlignments(const int& refId) const {
|
|
|
327 BamStandardIndexData::const_iterator indexIter = m_indexData.find(refId);
|
|
|
328 if ( indexIter == m_indexData.end() ) return false; // error
|
|
|
329 const ReferenceIndex& refEntry = (*indexIter).second;
|
|
|
330 return refEntry.HasAlignments;
|
|
|
331 }
|
|
|
332
|
|
|
333 // return true if all index data is cached
|
|
|
334 bool BamStandardIndex::HasFullDataCache(void) const {
|
|
|
335 return m_hasFullDataCache;
|
|
|
336 }
|
|
|
337
|
|
|
338 // returns true if index cache has data for desired reference
|
|
|
339 bool BamStandardIndex::IsDataLoaded(const int& refId) const {
|
|
|
340
|
|
|
341 // look up refId, return false if not found
|
|
|
342 BamStandardIndexData::const_iterator indexIter = m_indexData.find(refId);
|
|
|
343 if ( indexIter == m_indexData.end() ) return false;
|
|
|
344
|
|
|
345 // see if reference has alignments
|
|
|
346 // if not, it's not a problem to have no offset data
|
|
|
347 const ReferenceIndex& refEntry = (*indexIter).second;
|
|
|
348 if ( !refEntry.HasAlignments ) return true;
|
|
|
349
|
|
|
350 // return whether bin map contains data
|
|
|
351 return ( !refEntry.Bins.empty() );
|
|
|
352 }
|
|
|
353
|
|
|
354 // attempts to use index to jump to region; returns success/fail
|
|
|
355 bool BamStandardIndex::Jump(const BamRegion& region, bool* hasAlignmentsInRegion) {
|
|
|
356
|
|
|
357 // be sure reader & BGZF file are valid & open for reading
|
|
|
358 if ( m_reader == 0 || m_BGZF == 0 || !m_BGZF->IsOpen )
|
|
|
359 return false;
|
|
|
360
|
|
|
361 // make sure left-bound position is valid
|
|
|
362 if ( region.LeftPosition > m_references.at(region.LeftRefID).RefLength )
|
|
|
363 return false;
|
|
|
364
|
|
|
365 // calculate offsets for this region
|
|
|
366 // if failed, print message, set flag, and return failure
|
|
|
367 vector<int64_t> offsets;
|
|
|
368 if ( !GetOffsets(region, region.isRightBoundSpecified(), offsets, hasAlignmentsInRegion) ) {
|
|
|
369 fprintf(stderr, "ERROR: Could not jump: unable to calculate offset(s) for specified region.\n");
|
|
|
370 *hasAlignmentsInRegion = false;
|
|
|
371 return false;
|
|
|
372 }
|
|
|
373
|
|
|
374 // iterate through offsets
|
|
|
375 BamAlignment bAlignment;
|
|
|
376 bool result = true;
|
|
|
377 for ( vector<int64_t>::const_iterator o = offsets.begin(); o != offsets.end(); ++o) {
|
|
|
378
|
|
|
379 // attempt seek & load first available alignment
|
|
|
380 // set flag to true if data exists
|
|
|
381 result &= m_BGZF->Seek(*o);
|
|
|
382 *hasAlignmentsInRegion = m_reader->GetNextAlignmentCore(bAlignment);
|
|
|
383
|
|
|
384 // if this alignment corresponds to desired position
|
|
|
385 // return success of seeking back to the offset before the 'current offset' (to cover overlaps)
|
|
|
386 if ( ((bAlignment.RefID == region.LeftRefID) &&
|
|
|
387 ((bAlignment.Position + bAlignment.Length) > region.LeftPosition)) ||
|
|
|
388 (bAlignment.RefID > region.LeftRefID) )
|
|
|
389 {
|
|
|
390 if ( o != offsets.begin() ) --o;
|
|
|
391 return m_BGZF->Seek(*o);
|
|
|
392 }
|
|
|
393 }
|
|
|
394
|
|
|
395 // if error in jumping, print message & set flag
|
|
|
396 if ( !result ) {
|
|
|
397 fprintf(stderr, "ERROR: Could not jump: unable to determine correct offset for specified region.\n");
|
|
|
398 *hasAlignmentsInRegion = false;
|
|
|
399 }
|
|
|
400
|
|
|
401 // return success/failure
|
|
|
402 return result;
|
|
|
403 }
|
|
|
404
|
|
|
405 // clears index data from all references except the first
|
|
|
406 void BamStandardIndex::KeepOnlyFirstReferenceOffsets(void) {
|
|
|
407 BamStandardIndexData::const_iterator indexBegin = m_indexData.begin();
|
|
|
408 KeepOnlyReferenceOffsets((*indexBegin).first);
|
|
|
409 }
|
|
|
410
|
|
|
411 // clears index data from all references except the one specified
|
|
|
412 void BamStandardIndex::KeepOnlyReferenceOffsets(const int& refId) {
|
|
|
413 BamStandardIndexData::iterator mapIter = m_indexData.begin();
|
|
|
414 BamStandardIndexData::iterator mapEnd = m_indexData.end();
|
|
|
415 for ( ; mapIter != mapEnd; ++mapIter ) {
|
|
|
416 const int entryRefId = (*mapIter).first;
|
|
|
417 if ( entryRefId != refId )
|
|
|
418 ClearReferenceOffsets(entryRefId);
|
|
|
419 }
|
|
|
420 }
|
|
|
421
|
|
|
422 bool BamStandardIndex::LoadAllReferences(bool saveData) {
|
|
|
423
|
|
|
424 // skip if data already loaded
|
|
|
425 if ( m_hasFullDataCache ) return true;
|
|
|
426
|
|
|
427 // get number of reference sequences
|
|
|
428 uint32_t numReferences;
|
|
|
429 if ( !LoadReferenceCount((int&)numReferences) )
|
|
|
430 return false;
|
|
|
431
|
|
|
432 // iterate over reference entries
|
|
|
433 bool loadedOk = true;
|
|
|
434 for ( int i = 0; i < (int)numReferences; ++i )
|
|
|
435 loadedOk &= LoadReference(i, saveData);
|
|
|
436
|
|
|
437 // set flag
|
|
|
438 if ( loadedOk && saveData )
|
|
|
439 m_hasFullDataCache = true;
|
|
|
440
|
|
|
441 // return success/failure of loading references
|
|
|
442 return loadedOk;
|
|
|
443 }
|
|
|
444
|
|
|
445 // load header data from index file, return true if loaded OK
|
|
|
446 bool BamStandardIndex::LoadHeader(void) {
|
|
|
447
|
|
|
448 bool loadedOk = CheckMagicNumber();
|
|
|
449
|
|
|
450 // store offset of beginning of data
|
|
|
451 m_dataBeginOffset = ftell64(m_indexStream);
|
|
|
452
|
|
|
453 // return success/failure of load
|
|
|
454 return loadedOk;
|
|
|
455 }
|
|
|
456
|
|
|
457 // load a single index bin entry from file, return true if loaded OK
|
|
|
458 // @saveData - save data in memory if true, just read & discard if false
|
|
|
459 bool BamStandardIndex::LoadBin(ReferenceIndex& refEntry, bool saveData) {
|
|
|
460
|
|
|
461 size_t elementsRead = 0;
|
|
|
462
|
|
|
463 // get bin ID
|
|
|
464 uint32_t binId;
|
|
|
465 elementsRead += fread(&binId, sizeof(binId), 1, m_indexStream);
|
|
|
466 if ( m_isBigEndian ) SwapEndian_32(binId);
|
|
|
467
|
|
|
468 // load alignment chunks for this bin
|
|
|
469 ChunkVector chunks;
|
|
|
470 bool chunksOk = LoadChunks(chunks, saveData);
|
|
|
471
|
|
|
472 // store bin entry
|
|
|
473 if ( chunksOk && saveData )
|
|
|
474 refEntry.Bins.insert(pair<uint32_t, ChunkVector>(binId, chunks));
|
|
|
475
|
|
|
476 // return success/failure of load
|
|
|
477 return ( (elementsRead == 1) && chunksOk );
|
|
|
478 }
|
|
|
479
|
|
|
480 bool BamStandardIndex::LoadBins(ReferenceIndex& refEntry, bool saveData) {
|
|
|
481
|
|
|
482 size_t elementsRead = 0;
|
|
|
483
|
|
|
484 // get number of bins
|
|
|
485 int32_t numBins;
|
|
|
486 elementsRead += fread(&numBins, sizeof(numBins), 1, m_indexStream);
|
|
|
487 if ( m_isBigEndian ) SwapEndian_32(numBins);
|
|
|
488
|
|
|
489 // set flag
|
|
|
490 refEntry.HasAlignments = ( numBins != 0 );
|
|
|
491
|
|
|
492 // iterate over bins
|
|
|
493 bool binsOk = true;
|
|
|
494 for ( int i = 0; i < numBins; ++i )
|
|
|
495 binsOk &= LoadBin(refEntry, saveData);
|
|
|
496
|
|
|
497 // return success/failure of load
|
|
|
498 return ( (elementsRead == 1) && binsOk );
|
|
|
499 }
|
|
|
500
|
|
|
501 // load a single index bin entry from file, return true if loaded OK
|
|
|
502 // @saveData - save data in memory if true, just read & discard if false
|
|
|
503 bool BamStandardIndex::LoadChunk(ChunkVector& chunks, bool saveData) {
|
|
|
504
|
|
|
505 size_t elementsRead = 0;
|
|
|
506
|
|
|
507 // read in chunk data
|
|
|
508 uint64_t start;
|
|
|
509 uint64_t stop;
|
|
|
510 elementsRead += fread(&start, sizeof(start), 1, m_indexStream);
|
|
|
511 elementsRead += fread(&stop, sizeof(stop), 1, m_indexStream);
|
|
|
512
|
|
|
513 // swap endian-ness if necessary
|
|
|
514 if ( m_isBigEndian ) {
|
|
|
515 SwapEndian_64(start);
|
|
|
516 SwapEndian_64(stop);
|
|
|
517 }
|
|
|
518
|
|
|
519 // save data if requested
|
|
|
520 if ( saveData ) chunks.push_back( Chunk(start, stop) );
|
|
|
521
|
|
|
522 // return success/failure of load
|
|
|
523 return ( elementsRead == 2 );
|
|
|
524 }
|
|
|
525
|
|
|
526 bool BamStandardIndex::LoadChunks(ChunkVector& chunks, bool saveData) {
|
|
|
527
|
|
|
528 size_t elementsRead = 0;
|
|
|
529
|
|
|
530 // read in number of chunks
|
|
|
531 uint32_t numChunks;
|
|
|
532 elementsRead += fread(&numChunks, sizeof(numChunks), 1, m_indexStream);
|
|
|
533 if ( m_isBigEndian ) SwapEndian_32(numChunks);
|
|
|
534
|
|
|
535 // initialize space for chunks if we're storing this data
|
|
|
536 if ( saveData ) chunks.reserve(numChunks);
|
|
|
537
|
|
|
538 // iterate over chunks
|
|
|
539 bool chunksOk = true;
|
|
|
540 for ( int i = 0; i < (int)numChunks; ++i )
|
|
|
541 chunksOk &= LoadChunk(chunks, saveData);
|
|
|
542
|
|
|
543 // sort chunk vector
|
|
|
544 sort( chunks.begin(), chunks.end(), ChunkLessThan );
|
|
|
545
|
|
|
546 // return success/failure of load
|
|
|
547 return ( (elementsRead == 1) && chunksOk );
|
|
|
548 }
|
|
|
549
|
|
|
550 // load a single index linear offset entry from file, return true if loaded OK
|
|
|
551 // @saveData - save data in memory if true, just read & discard if false
|
|
|
552 bool BamStandardIndex::LoadLinearOffsets(ReferenceIndex& refEntry, bool saveData) {
|
|
|
553
|
|
|
554 size_t elementsRead = 0;
|
|
|
555
|
|
|
556 // read in number of linear offsets
|
|
|
557 int32_t numLinearOffsets;
|
|
|
558 elementsRead += fread(&numLinearOffsets, sizeof(numLinearOffsets), 1, m_indexStream);
|
|
|
559 if ( m_isBigEndian ) SwapEndian_32(numLinearOffsets);
|
|
|
560
|
|
|
561 // set up destination vector (if we're saving the data)
|
|
|
562 LinearOffsetVector linearOffsets;
|
|
|
563 if ( saveData ) linearOffsets.reserve(numLinearOffsets);
|
|
|
564
|
|
|
565 // iterate over linear offsets
|
|
|
566 uint64_t linearOffset;
|
|
|
567 for ( int i = 0; i < numLinearOffsets; ++i ) {
|
|
|
568 elementsRead += fread(&linearOffset, sizeof(linearOffset), 1, m_indexStream);
|
|
|
569 if ( m_isBigEndian ) SwapEndian_64(linearOffset);
|
|
|
570 if ( saveData ) linearOffsets.push_back(linearOffset);
|
|
|
571 }
|
|
|
572
|
|
|
573 // sort linear offsets
|
|
|
574 sort ( linearOffsets.begin(), linearOffsets.end() );
|
|
|
575
|
|
|
576 // save in reference index entry if desired
|
|
|
577 if ( saveData ) refEntry.Offsets = linearOffsets;
|
|
|
578
|
|
|
579 // return success/failure of load
|
|
|
580 return ( elementsRead == (size_t)(numLinearOffsets + 1) );
|
|
|
581 }
|
|
|
582
|
|
|
583 bool BamStandardIndex::LoadFirstReference(bool saveData) {
|
|
|
584 BamStandardIndexData::const_iterator indexBegin = m_indexData.begin();
|
|
|
585 return LoadReference((*indexBegin).first, saveData);
|
|
|
586 }
|
|
|
587
|
|
|
588 // load a single reference from file, return true if loaded OK
|
|
|
589 // @saveData - save data in memory if true, just read & discard if false
|
|
|
590 bool BamStandardIndex::LoadReference(const int& refId, bool saveData) {
|
|
|
591
|
|
|
592 // look up refId
|
|
|
593 BamStandardIndexData::iterator indexIter = m_indexData.find(refId);
|
|
|
594
|
|
|
595 // if reference not previously loaded, create new entry
|
|
|
596 if ( indexIter == m_indexData.end() ) {
|
|
|
597 ReferenceIndex newEntry;
|
|
|
598 newEntry.HasAlignments = false;
|
|
|
599 m_indexData.insert( pair<int32_t, ReferenceIndex>(refId, newEntry) );
|
|
|
600 }
|
|
|
601
|
|
|
602 // load reference data
|
|
|
603 indexIter = m_indexData.find(refId);
|
|
|
604 ReferenceIndex& entry = (*indexIter).second;
|
|
|
605 bool loadedOk = true;
|
|
|
606 loadedOk &= LoadBins(entry, saveData);
|
|
|
607 loadedOk &= LoadLinearOffsets(entry, saveData);
|
|
|
608 return loadedOk;
|
|
|
609 }
|
|
|
610
|
|
|
611 // loads number of references, return true if loaded OK
|
|
|
612 bool BamStandardIndex::LoadReferenceCount(int& numReferences) {
|
|
|
613
|
|
|
614 size_t elementsRead = 0;
|
|
|
615
|
|
|
616 // read reference count
|
|
|
617 elementsRead += fread(&numReferences, sizeof(numReferences), 1, m_indexStream);
|
|
|
618 if ( m_isBigEndian ) SwapEndian_32(numReferences);
|
|
|
619
|
|
|
620 // return success/failure of load
|
|
|
621 return ( elementsRead == 1 );
|
|
|
622 }
|
|
|
623
|
|
|
624 // merges 'alignment chunks' in BAM bin (used for index building)
|
|
|
625 void BamStandardIndex::MergeChunks(void) {
|
|
|
626
|
|
|
627 // iterate over reference enties
|
|
|
628 BamStandardIndexData::iterator indexIter = m_indexData.begin();
|
|
|
629 BamStandardIndexData::iterator indexEnd = m_indexData.end();
|
|
|
630 for ( ; indexIter != indexEnd; ++indexIter ) {
|
|
|
631
|
|
|
632 // get BAM bin map for this reference
|
|
|
633 ReferenceIndex& refIndex = (*indexIter).second;
|
|
|
634 BamBinMap& bamBinMap = refIndex.Bins;
|
|
|
635
|
|
|
636 // iterate over BAM bins
|
|
|
637 BamBinMap::iterator binIter = bamBinMap.begin();
|
|
|
638 BamBinMap::iterator binEnd = bamBinMap.end();
|
|
|
639 for ( ; binIter != binEnd; ++binIter ) {
|
|
|
640
|
|
|
641 // get chunk vector for this bin
|
|
|
642 ChunkVector& binChunks = (*binIter).second;
|
|
|
643 if ( binChunks.size() == 0 ) continue;
|
|
|
644
|
|
|
645 ChunkVector mergedChunks;
|
|
|
646 mergedChunks.push_back( binChunks[0] );
|
|
|
647
|
|
|
648 // iterate over chunks
|
|
|
649 int i = 0;
|
|
|
650 ChunkVector::iterator chunkIter = binChunks.begin();
|
|
|
651 ChunkVector::iterator chunkEnd = binChunks.end();
|
|
|
652 for ( ++chunkIter; chunkIter != chunkEnd; ++chunkIter) {
|
|
|
653
|
|
|
654 // get 'currentChunk' based on numeric index
|
|
|
655 Chunk& currentChunk = mergedChunks[i];
|
|
|
656
|
|
|
657 // get iteratorChunk based on vector iterator
|
|
|
658 Chunk& iteratorChunk = (*chunkIter);
|
|
|
659
|
|
|
660 // if chunk ends where (iterator) chunk starts, then merge
|
|
|
661 if ( currentChunk.Stop>>16 == iteratorChunk.Start>>16 )
|
|
|
662 currentChunk.Stop = iteratorChunk.Stop;
|
|
|
663
|
|
|
664 // otherwise
|
|
|
665 else {
|
|
|
666 // set currentChunk + 1 to iteratorChunk
|
|
|
667 mergedChunks.push_back(iteratorChunk);
|
|
|
668 ++i;
|
|
|
669 }
|
|
|
670 }
|
|
|
671
|
|
|
672 // saved merged chunk vector
|
|
|
673 (*binIter).second = mergedChunks;
|
|
|
674 }
|
|
|
675 }
|
|
|
676 }
|
|
|
677
|
|
|
678 // saves BAM bin entry for index
|
|
|
679 void BamStandardIndex::SaveBinEntry(BamBinMap& binMap,
|
|
|
680 const uint32_t& saveBin,
|
|
|
681 const uint64_t& saveOffset,
|
|
|
682 const uint64_t& lastOffset)
|
|
|
683 {
|
|
|
684 // look up saveBin
|
|
|
685 BamBinMap::iterator binIter = binMap.find(saveBin);
|
|
|
686
|
|
|
687 // create new chunk
|
|
|
688 Chunk newChunk(saveOffset, lastOffset);
|
|
|
689
|
|
|
690 // if entry doesn't exist
|
|
|
691 if ( binIter == binMap.end() ) {
|
|
|
692 ChunkVector newChunks;
|
|
|
693 newChunks.push_back(newChunk);
|
|
|
694 binMap.insert( pair<uint32_t, ChunkVector>(saveBin, newChunks));
|
|
|
695 }
|
|
|
696
|
|
|
697 // otherwise
|
|
|
698 else {
|
|
|
699 ChunkVector& binChunks = (*binIter).second;
|
|
|
700 binChunks.push_back( newChunk );
|
|
|
701 }
|
|
|
702 }
|
|
|
703
|
|
|
704 // saves linear offset entry for index
|
|
|
705 void BamStandardIndex::SaveLinearOffset(LinearOffsetVector& offsets,
|
|
|
706 const BamAlignment& bAlignment,
|
|
|
707 const uint64_t& lastOffset)
|
|
|
708 {
|
|
|
709 // get converted offsets
|
|
|
710 int beginOffset = bAlignment.Position >> BAM_LIDX_SHIFT;
|
|
|
711 int endOffset = (bAlignment.GetEndPosition() - 1) >> BAM_LIDX_SHIFT;
|
|
|
712
|
|
|
713 // resize vector if necessary
|
|
|
714 int oldSize = offsets.size();
|
|
|
715 int newSize = endOffset + 1;
|
|
|
716 if ( oldSize < newSize )
|
|
|
717 offsets.resize(newSize, 0);
|
|
|
718
|
|
|
719 // store offset
|
|
|
720 for( int i = beginOffset + 1; i <= endOffset; ++i ) {
|
|
|
721 if ( offsets[i] == 0 )
|
|
|
722 offsets[i] = lastOffset;
|
|
|
723 }
|
|
|
724 }
|
|
|
725
|
|
|
726 // initializes index data structure to hold @count references
|
|
|
727 void BamStandardIndex::SetReferenceCount(const int& count) {
|
|
|
728 for ( int i = 0; i < count; ++i )
|
|
|
729 m_indexData[i].HasAlignments = false;
|
|
|
730 }
|
|
|
731
|
|
|
732 bool BamStandardIndex::SkipToFirstReference(void) {
|
|
|
733 BamStandardIndexData::const_iterator indexBegin = m_indexData.begin();
|
|
|
734 return SkipToReference( (*indexBegin).first );
|
|
|
735 }
|
|
|
736
|
|
|
737 // position file pointer to desired reference begin, return true if skipped OK
|
|
|
738 bool BamStandardIndex::SkipToReference(const int& refId) {
|
|
|
739
|
|
|
740 // attempt rewind
|
|
|
741 if ( !Rewind() ) return false;
|
|
|
742
|
|
|
743 // read in number of references
|
|
|
744 uint32_t numReferences;
|
|
|
745 size_t elementsRead = fread(&numReferences, sizeof(numReferences), 1, m_indexStream);
|
|
|
746 if ( elementsRead != 1 ) return false;
|
|
|
747 if ( m_isBigEndian ) SwapEndian_32(numReferences);
|
|
|
748
|
|
|
749 // iterate over reference entries
|
|
|
750 bool skippedOk = true;
|
|
|
751 int currentRefId = 0;
|
|
|
752 while (currentRefId != refId) {
|
|
|
753 skippedOk &= LoadReference(currentRefId, false);
|
|
|
754 ++currentRefId;
|
|
|
755 }
|
|
|
756
|
|
|
757 // return success
|
|
|
758 return skippedOk;
|
|
|
759 }
|
|
|
760
|
|
|
761 // write header to new index file
|
|
|
762 bool BamStandardIndex::WriteHeader(void) {
|
|
|
763
|
|
|
764 size_t elementsWritten = 0;
|
|
|
765
|
|
|
766 // write magic number
|
|
|
767 elementsWritten += fwrite("BAI\1", sizeof(char), 4, m_indexStream);
|
|
|
768
|
|
|
769 // store offset of beginning of data
|
|
|
770 m_dataBeginOffset = ftell64(m_indexStream);
|
|
|
771
|
|
|
772 // return success/failure of write
|
|
|
773 return (elementsWritten == 4);
|
|
|
774 }
|
|
|
775
|
|
|
776 // write index data for all references to new index file
|
|
|
777 bool BamStandardIndex::WriteAllReferences(void) {
|
|
|
778
|
|
|
779 size_t elementsWritten = 0;
|
|
|
780
|
|
|
781 // write number of reference sequences
|
|
|
782 int32_t numReferenceSeqs = m_indexData.size();
|
|
|
783 if ( m_isBigEndian ) SwapEndian_32(numReferenceSeqs);
|
|
|
784 elementsWritten += fwrite(&numReferenceSeqs, sizeof(numReferenceSeqs), 1, m_indexStream);
|
|
|
785
|
|
|
786 // iterate over reference sequences
|
|
|
787 bool refsOk = true;
|
|
|
788 BamStandardIndexData::const_iterator indexIter = m_indexData.begin();
|
|
|
789 BamStandardIndexData::const_iterator indexEnd = m_indexData.end();
|
|
|
790 for ( ; indexIter != indexEnd; ++ indexIter )
|
|
|
791 refsOk &= WriteReference( (*indexIter).second );
|
|
|
792
|
|
|
793 // return success/failure of write
|
|
|
794 return ( (elementsWritten == 1) && refsOk );
|
|
|
795 }
|
|
|
796
|
|
|
797 // write index data for bin to new index file
|
|
|
798 bool BamStandardIndex::WriteBin(const uint32_t& binId, const ChunkVector& chunks) {
|
|
|
799
|
|
|
800 size_t elementsWritten = 0;
|
|
|
801
|
|
|
802 // write BAM bin ID
|
|
|
803 uint32_t binKey = binId;
|
|
|
804 if ( m_isBigEndian ) SwapEndian_32(binKey);
|
|
|
805 elementsWritten += fwrite(&binKey, sizeof(binKey), 1, m_indexStream);
|
|
|
806
|
|
|
807 // write chunks
|
|
|
808 bool chunksOk = WriteChunks(chunks);
|
|
|
809
|
|
|
810 // return success/failure of write
|
|
|
811 return ( (elementsWritten == 1) && chunksOk );
|
|
|
812 }
|
|
|
813
|
|
|
814 // write index data for bins to new index file
|
|
|
815 bool BamStandardIndex::WriteBins(const BamBinMap& bins) {
|
|
|
816
|
|
|
817 size_t elementsWritten = 0;
|
|
|
818
|
|
|
819 // write number of bins
|
|
|
820 int32_t binCount = bins.size();
|
|
|
821 if ( m_isBigEndian ) SwapEndian_32(binCount);
|
|
|
822 elementsWritten += fwrite(&binCount, sizeof(binCount), 1, m_indexStream);
|
|
|
823
|
|
|
824 // iterate over bins
|
|
|
825 bool binsOk = true;
|
|
|
826 BamBinMap::const_iterator binIter = bins.begin();
|
|
|
827 BamBinMap::const_iterator binEnd = bins.end();
|
|
|
828 for ( ; binIter != binEnd; ++binIter )
|
|
|
829 binsOk &= WriteBin( (*binIter).first, (*binIter).second );
|
|
|
830
|
|
|
831 // return success/failure of write
|
|
|
832 return ( (elementsWritten == 1) && binsOk );
|
|
|
833 }
|
|
|
834
|
|
|
835 // write index data for chunk entry to new index file
|
|
|
836 bool BamStandardIndex::WriteChunk(const Chunk& chunk) {
|
|
|
837
|
|
|
838 size_t elementsWritten = 0;
|
|
|
839
|
|
|
840 // localize alignment chunk offsets
|
|
|
841 uint64_t start = chunk.Start;
|
|
|
842 uint64_t stop = chunk.Stop;
|
|
|
843
|
|
|
844 // swap endian-ness if necessary
|
|
|
845 if ( m_isBigEndian ) {
|
|
|
846 SwapEndian_64(start);
|
|
|
847 SwapEndian_64(stop);
|
|
|
848 }
|
|
|
849
|
|
|
850 // write to index file
|
|
|
851 elementsWritten += fwrite(&start, sizeof(start), 1, m_indexStream);
|
|
|
852 elementsWritten += fwrite(&stop, sizeof(stop), 1, m_indexStream);
|
|
|
853
|
|
|
854 // return success/failure of write
|
|
|
855 return ( elementsWritten == 2 );
|
|
|
856 }
|
|
|
857
|
|
|
858 // write index data for chunk entry to new index file
|
|
|
859 bool BamStandardIndex::WriteChunks(const ChunkVector& chunks) {
|
|
|
860
|
|
|
861 size_t elementsWritten = 0;
|
|
|
862
|
|
|
863 // write chunks
|
|
|
864 int32_t chunkCount = chunks.size();
|
|
|
865 if ( m_isBigEndian ) SwapEndian_32(chunkCount);
|
|
|
866 elementsWritten += fwrite(&chunkCount, sizeof(chunkCount), 1, m_indexStream);
|
|
|
867
|
|
|
868 // iterate over chunks
|
|
|
869 bool chunksOk = true;
|
|
|
870 ChunkVector::const_iterator chunkIter = chunks.begin();
|
|
|
871 ChunkVector::const_iterator chunkEnd = chunks.end();
|
|
|
872 for ( ; chunkIter != chunkEnd; ++chunkIter )
|
|
|
873 chunksOk &= WriteChunk( (*chunkIter) );
|
|
|
874
|
|
|
875 // return success/failure of write
|
|
|
876 return ( (elementsWritten == 1) && chunksOk );
|
|
|
877 }
|
|
|
878
|
|
|
879 // write index data for linear offsets entry to new index file
|
|
|
880 bool BamStandardIndex::WriteLinearOffsets(const LinearOffsetVector& offsets) {
|
|
|
881
|
|
|
882 size_t elementsWritten = 0;
|
|
|
883
|
|
|
884 // write number of linear offsets
|
|
|
885 int32_t offsetCount = offsets.size();
|
|
|
886 if ( m_isBigEndian ) SwapEndian_32(offsetCount);
|
|
|
887 elementsWritten += fwrite(&offsetCount, sizeof(offsetCount), 1, m_indexStream);
|
|
|
888
|
|
|
889 // iterate over linear offsets
|
|
|
890 LinearOffsetVector::const_iterator offsetIter = offsets.begin();
|
|
|
891 LinearOffsetVector::const_iterator offsetEnd = offsets.end();
|
|
|
892 for ( ; offsetIter != offsetEnd; ++offsetIter ) {
|
|
|
893
|
|
|
894 // write linear offset
|
|
|
895 uint64_t linearOffset = (*offsetIter);
|
|
|
896 if ( m_isBigEndian ) SwapEndian_64(linearOffset);
|
|
|
897 elementsWritten += fwrite(&linearOffset, sizeof(linearOffset), 1, m_indexStream);
|
|
|
898 }
|
|
|
899
|
|
|
900 // return success/failure of write
|
|
|
901 return ( elementsWritten == (size_t)(offsetCount + 1) );
|
|
|
902 }
|
|
|
903
|
|
|
904 // write index data for a single reference to new index file
|
|
|
905 bool BamStandardIndex::WriteReference(const ReferenceIndex& refEntry) {
|
|
|
906 bool refOk = true;
|
|
|
907 refOk &= WriteBins(refEntry.Bins);
|
|
|
908 refOk &= WriteLinearOffsets(refEntry.Offsets);
|
|
|
909 return refOk;
|
|
|
910 }
|
