62 if(
ifread(indexFile, magic, 4) != 4)
70 if (magic[0] !=
'B' || magic[1] !=
'A' || magic[2] !=
'I' || magic[3] != 1)
79 if(
ifread(indexFile, &n_ref, 4) != 4)
89 for(
int refIndex = 0; refIndex < n_ref; refIndex++)
95 if(
ifread(indexFile, &(ref->n_bin), 4) != 4)
112 ref->bins.resize(ref->n_bin + 1);
115 for(
int binIndex = 0; binIndex < ref->n_bin; binIndex++)
120 if(
ifread(indexFile, &(binNumber), 4) != 4)
130 Bin* binPtr = &(ref->bins[binIndex]);
131 binPtr->bin = binNumber;
134 if(
ifread(indexFile, &(binPtr->n_chunk), 4) != 4)
144 uint32_t sizeOfChunkList = binPtr->n_chunk *
sizeof(
Chunk);
145 binPtr->chunks = (
Chunk*)malloc(sizeOfChunkList);
146 if(
ifread(indexFile, binPtr->chunks, sizeOfChunkList) != sizeOfChunkList)
155 if(binPtr->bin != MAX_NUM_BINS)
157 for(
int i = 0; i < binPtr->n_chunk; i++)
159 if(binPtr->chunks[i].chunk_beg < ref->minChunkOffset)
161 ref->minChunkOffset = binPtr->chunks[i].chunk_beg;
163 if(binPtr->chunks[i].chunk_end > ref->maxChunkOffset)
165 ref->maxChunkOffset = binPtr->chunks[i].chunk_end;
167 if(binPtr->chunks[i].chunk_end > maxOverallOffset)
169 maxOverallOffset = binPtr->chunks[i].chunk_end;
177 ref->n_mapped = binPtr->chunks[binPtr->n_chunk - 1].chunk_beg;
178 ref->n_unmapped = binPtr->chunks[binPtr->n_chunk - 1].chunk_end;
183 if(
ifread(indexFile, &(ref->n_intv), 4) != 4)
193 uint32_t linearIndexSize = ref->n_intv *
sizeof(uint64_t);
194 ref->ioffsets = (uint64_t*)malloc(linearIndexSize);
195 if(
ifread(indexFile, ref->ioffsets, linearIndexSize) != linearIndexSize)
204 int32_t numUnmapped = 0;
205 if(
ifread(indexFile, &numUnmapped,
sizeof(int32_t)) ==
sizeof(int32_t))
207 myUnMappedNumReads = numUnmapped;
225 if((start >= end) && (end != -1))
227 std::cerr <<
"Warning, requesting region where start <= end, so "
228 <<
"no values will be returned.\n";
239 refChunk.chunk_beg = getMaxOffset();
242 refChunk.chunk_end = Chunk::MAX_CHUNK_VALUE;
243 return(chunkList.insert(refChunk));
246 if((refID < 0) || (refID >= n_ref))
249 std::cerr <<
"Warning, requesting refID is out of range, so "
250 <<
"no values will be returned.\n";
262 if(ref->maxChunkOffset == 0)
268 refChunk.chunk_beg = ref->minChunkOffset;
269 refChunk.chunk_end = ref->maxChunkOffset;
270 return(chunkList.insert(refChunk));
280 end = MAX_POSITION + 1;
285 uint64_t minOffset = 0;
286 getMinOffsetFromLinearIndex(refID, start, minOffset);
288 bool binInRangeMap[MAX_NUM_BINS+1];
290 getBinsForRegion(start, end, binInRangeMap);
293 for(
int i = 0; i < ref->n_bin; ++i)
295 const Bin* bin = &(ref->bins[i]);
296 if(binInRangeMap[bin->bin] ==
false)
303 for(
int chunkIndex = 0; chunkIndex < bin->n_chunk; chunkIndex++)
309 if(bin->chunks[chunkIndex].chunk_end < minOffset)
314 if(!chunkList.insert(bin->chunks[chunkIndex]))
317 std::cerr <<
"Warning, Failed to add a chunk, so "
318 <<
"no values will be returned.\n";
326 return(chunkList.mergeOverlapping());
400 std::cout <<
"BAM Index: " << std::endl;
401 std::cout <<
"# Reference Sequences: " << n_ref << std::endl;
403 unsigned int startRef = 0;
404 unsigned int endRef = myRefs.size() - 1;
405 std::vector<Reference> refsToProcess;
414 for(
unsigned int i = startRef; i <= endRef; ++i)
416 std::cout << std::dec
417 <<
"\tReference ID: " << std::setw(4) << i
418 <<
"; #Bins: "<< std::setw(6) << myRefs[i].n_bin
419 <<
"; #Linear Index Entries: "
420 << std::setw(6) << myRefs[i].n_intv
421 <<
"; Min Chunk Offset: "
422 << std::setw(18) << std::hex << std::showbase << myRefs[i].minChunkOffset
423 <<
"; Max Chunk Offset: "
424 << std::setw(18) << myRefs[i].maxChunkOffset
427 if(myRefs[i].n_mapped != Reference::UNKNOWN_MAP_INFO)
429 std::cout <<
"; " << myRefs[i].n_mapped <<
" Mapped Reads";
431 if(myRefs[i].n_mapped != Reference::UNKNOWN_MAP_INFO)
433 std::cout <<
"; " << myRefs[i].n_unmapped <<
" Unmapped Reads";
435 std::cout << std::endl;
440 std::vector<Bin>::iterator binIter;
441 for(binIter = myRefs[i].bins.begin();
442 binIter != myRefs[i].bins.end();
445 Bin* binPtr = &(*binIter);
446 if(binPtr->bin == Bin::NOT_USED_BIN)
452 std::cout <<
"\t\t\tBin Name: " << binPtr->bin << std::endl;
453 std::cout <<
"\t\t\t# Chunks: " << binPtr->n_chunk << std::endl;
454 std::cout << std::hex << std::showbase;
456 for(
int chunkIndex = 0; chunkIndex < binPtr->n_chunk;
462 if((binPtr->bin != MAX_NUM_BINS) ||
463 (chunkIndex != (binPtr->n_chunk - 1)))
465 std::cout <<
"\t\t\t\tchunk_beg: "
466 << binPtr->chunks[chunkIndex].chunk_beg
468 std::cout <<
"\t\t\t\tchunk_end: "
469 << binPtr->chunks[chunkIndex].chunk_end
474 std::cout << std::dec;
477 for(
int linearIndex = 0; linearIndex < myRefs[i].n_intv;
480 if(myRefs[i].ioffsets[linearIndex] != 0)
482 std::cout <<
"\t\t\tLinearIndex["
483 << std::dec << linearIndex <<
"] Offset: "
484 << std::hex << myRefs[i].ioffsets[linearIndex]
bool getChunksForRegion(int32_t refID, int32_t start, int32_t end, SortedChunkList &chunkList)
Get the list of chunks associated with this region.