Repository navigation
Extreme memory consumption when reading certain SRA records? #31
Description
I am reading sequence data from SRA records by first downloading the SRA record with the prefetch command and then iterating through the file using the C++ interface (version 2.10.8), i.e.:
ngs::ReadCollection run("DRR001375");
const size_t num_read = run.getReadCount(ngs::Read::all);
ngs::ReadIterator run_iter = ngs::ReadIterator( run.getReadRange ( 1, num_read, ngs::Read::all ) );
size_t read_count = 0;
while( run_iter.nextRead() ){
++read_count;
while( run_iter.nextFragment() ){
const string seq = run_iter.getFragmentBases().toString();
// Process the read sequence ...
process_read(seq);
}
}
In general, this approach seems to work well. However, I have noticed that for some SRA records (like ERR191522), there is (a) significant memory consumption and (b) a dramatic slow-down when iterating through the file. The following plot shows the speed (in reads per second) and memory consumption (from /proc/meminfo, reported as a fraction of total system memory):
Other SRA records seem to be fine. For DRR001375, the following graph shows fairly constant speed and memory usage:
Is there a way to read SRA records, like ERR191522, without the large memory consumption? If not, is there a way to identify SRA records (in advance) that will exhibit this behavior (as the available RAM on on cluster instances can easily be exhausted while processing a single SRA record).
Activity
durbrow commentedon Jul 27, 2020
The big difference between these two runs is ERR191522 contains alignments and DRR001375 contains only unaligned reads.
In the unaligned case, the reads are stored together in one record. The access pattern in your code snippet will be in the same order as the data is laid out in the file.
In the aligned case, the reads are stored in the order they align on the reference. The two mate pairs are not stored together and might be far apart in the file. Your code snippet will reconstruct the whole read which requires random access I/O to the file. Once upon a time, this wasn't a big problem, but these days big disks are structured as bucket stores with random I/O being very expensive. To try to ameliorate this, we have tuned our underlying database technology to cache aggressively. If your usage does NOT require both mate pairs, it could be sped up.
kwrodarmer commentedon Jul 27, 2020
We are creating an internal ticket to look into this behavior. Thanks for the report and the graphs!
jgans commentedon Jul 27, 2020
For the particular problem we're trying to solve, we do not need to access mate pairs together (or in any particular order). Is there a code example that illustrates how to access read data without the need for aggressive caching or random access I/O?
kwrodarmer commentedon Jul 27, 2020
When a run that is aligned has been detected, the reads will be available first through their primary alignments. If these are retrieved as individual reads in absence of quality scores, they will generate the best performance. The remaining reads will be those that did not align, and they can be fetched separately.
You start with a
ReadCollection.You can then ask how many primary alignments it has with
ReadCollection::getAlignmentCount(Alignment::primaryAlignment).If there are alignments, obtain them via
ReadCollection::getAlignments(Alignment::primaryAlignment)and iterate across them.The aligned read is obtained via
Alignment::getAlignedFragmentBases(). You may need to compensate for orientation.Once finished, look for unaligned reads via
ReadCollection::getReadCount(Read::unaligned), and similarly obtain an iterator on them viaReadCollection::getReads(Read::unaligned).The ReadIterator is an extension of FragmentIterator, and should be treated as such. Otherwise, you will tend to assemble paired reads which will cause issues.
Walk the reads using
FragmentIterator::nextFragment()and obtain bases withFragment::getFragmentBases().It's also useful to avoid converting a
StringRefinto astd::stringfor everything - do so only when necessary. Try to useStringRefto access data and size rather than always NUL-terminating.Okay, I think that's it - off the top of my head. If I left anything out, @durbrow please correct me!
jgans commentedon Jul 27, 2020
When counting the number of sequences in the primary alignment and unaligned reads, I get a slightly different sequence count from the value obtained by iterating over
getReadCount(ngs::Read::all)"raw" reads. Is there another source of read data in an SRA record?For example, in SRR10742149:
getAlignments(ngs::Alignment::primaryAlignment)ReadIterator( run.getReads(ngs::Read::unaligned) )Total number of sequences = 83490707
However, there are 83491358 sequences in 41745679 reads, using
ngs::ReadIterator( run.getReadRange ( 1, getReadCount(ngs::Read::all), ngs::Read::all ) ).kwrodarmer commentedon Jul 27, 2020
Yes, I probably should have mentioned that.
As you see,
ngs::Read::allgives aligned, partially aligned, and unaligned reads. Looking at the primary alignments gets you thengs::Read::alignedcategory, andngs::Read::unalignedgets you the cases where neither mate is aligned.For the cases where one mate is aligned but the other not, these are partially aligned and we don't have a nice, clean category for them. If you ask for all but iterate as fragments, this may be what you need.
jgans commentedon Jul 28, 2020
Thank you for suggesting this strategy! It does indeed provide a significant reduction in memory consumed and a significant speedup as well. As shown below, the memory consumption required to read through ERR191522 is an order-of-magnitude less than the previous, naive strategy, and the run-time is a factor of 14 faster (16 min vs 225 min):
To facilitate comparison with the previous read pairs/sec scale, this graph shows the number of aligned sequences/2 per second. The memory consumption is now the "change in memory consumption" (starting from the beginning of program execution). For my particular application, this performance improvement is large enough to justify ignoring the partially aligned read pairs.
kwrodarmer commentedon Jul 28, 2020
Great to hear!
One of the things that is problematic is assembling the mates. If you can take the mates in isolation, you'll have a better day.
jgans commentedon Aug 13, 2020
While the above strategy of first loading primary alignments and then loading unaligned reads works with many (most?) SRA records (like the above mentioned ERR191522), there appear to be some SRA records for which this strategy generates an error.
For example, all of the 2007231 primary alignments in ERR634825 can be read, but the following C++ exception is thrown when attempting iterate through the 846553 unaligned reads:
VCursorAddColumn failed: '(INSDC:dna:text)CMP_READ' rc = RC(rcVDB,rcCursor,rcUpdating,rcColumn,rcNotFound). An incomplete list of other examples of SRA records that yield this same error include: ERR634636, ERR634688, ERR634695, and ERR634696. All of these records appear to be single read, AB Solid4 experiments from the same study.durbrow commentedon Aug 14, 2020
Ah yes, aligned colorspace! They are broken, and unfortunately they can't (*) be fixed. Is it crucial for you to have these reads? Although I don't have a count, I do know there aren't very many these in the SRA.
* the people who decide such things decided it wasn't worth the expense
jgans commentedon Aug 14, 2020
These reads are not crucial. However, is it possible to interrogate an
ngs::ReadCollectionobject to detect aligned colorspace data prior to loading reads? I could then usengs::ReadIterator(ngs::Read::all)to iterate over all of the reads in these SRA records (instead of reading primary alignments and unaligned reads).As it stands, I can identify records that can successfully load 100% of the
ngs::Alignment::primaryAlignmentdata but fail on the firstngs::Read::unalignedread. The main downside of this approach would be the wast of time reading all of the primary alignments.