Out of memory when attempting to read BAM header in REDI3 analyze

[void-tiff]

Hello! Thank you so much for your work on this tool, it's been a lot easier to use for my huge dataset than REDI2.

I have a BAM file that's ~30GB that I'm trying to analyze, but it kills the job pretty quickly due to running out of memory trying to read the header:

[E::bam_hdr_read] Out of memory
[ERROR] (<class 'ValueError'>) file does not have a valid header (mode='r') - is it BAM/CRAM format?
[ERROR] Killing job

I have specified 96 threads (running on 48 CPUs), so after looking it up I suspect it might be an issue with the number of threads but I'm not sure. I can run samtools view -H in the command line just fine so I know this isn't an issue with the header itself being corrupted or not being there at all.

If it helps, the full command I'm running is
python -m reditools analyze {input BAM file} -r {trinity fasta} -o {output file} -t 96

It could also be the fact that I'm not using the --window command, but I don't have a good intuition on how big or small that number should be, nor do I know if that would affect the reading of the header.

I'm brand new to HPCs and am learning as I go, so your help and patience is very appreciated!!

[ahanden]

Interesting. We had a user report a similar issue that is being caused by the cache storing the reference fasta file (A patch for this will be included in the next release). You might try adding MD tags to your BAM file (samtools calmd) and then run REDItools without a fasata file.

If you do want to give the window option a try, I usually start with the length of my largest chromosome divided by 10.

[void-tiff]

Thank you so much for your fast response!!
I'll try adding MD tags to the file. Do you know if running without a fasta reference will change the output at all?

[ahanden]

You will get the exact same output, no worries there.

[void-tiff]

Unfortunately, even after adding MD tags and running without a FASTA reference I'm still encountering this bug. samtools calmd and index ran within a couple minutes, but as soon as it got to REDItools it gave me the same error.

[ahanden]

Try reducing the number of threads to something like 10 or 20 and let me know if you still get the error.

[void-tiff]

I have reduced the thread count to 10 and it seems to be working fine!! This is with the md-tagged files and without the fasta.

As a complete side note, do you have a way to estimate about how long a run takes with number of threads/"chromosomes"/file size? If not, that's cool! Just wanted a way to properly estimate how long it will take for when I sbatch my slurm jobs.

[void-tiff]

Is there any kind of progress checking that the program does in case it fails to complete? Besides my verbose run, two of my out files are blank and the output files I specified don't exist. If I run it again, will it pick up where it left off or start over completely?

[ahanden]

I have reduced the thread count to 10 and it seems to be working fine!! This is with the md-tagged files and without the fasta.

Great to hear!

As a complete side note, do you have a way to estimate about how long a run takes with number of threads/"chromosomes"/file size? If not, that's cool! Just wanted a way to properly estimate how long it will take for when I sbatch my slurm jobs.

Runtime depends not only on file size, threads, and chromosomes, but also on the parameter values (and there are a lot of parameters), so it's been difficult coming up with general guidelines. Here's the best I have so far: using options -s 1 -q 25 -bq 25 -l 10 -me 3 -C on a BAM file aligned to hg38 with 49 million mapped reads, I've found the following:

Threads	Runtime (minutes)
1	457
5	63
20	22

(note that this is on the development code a.k.a. v3.7)

Is there any kind of progress checking that the program does in case it fails to complete? Besides my verbose run, two of my out files are blank and the output files I specified don't exist. If I run it again, will it pick up where it left off or start over completely?

There is no resume option, but that's an interesting feature idea that is definitely worth pursuing. The "out files" that are blank, are those the temporary files? It's possible they could be blank if they correspond to chromosomes without any reads (or without edits). The file naming convention doesn't offer any help on what chromosome each temporary file corresponds to though.

I don't think there will be time to add this to release 3.7, but for 3.8, these changes will be a priority.

[void-tiff]

Runtime depends not only on file size, threads, and chromosomes, but also on the parameter values (and there are a lot of parameters), so it's been difficult coming up with general guidelines.

Ah, I see. I have ~500k chromosomes, and in using the defaults for all of those parameters on 20 threads, 48 hours was not enough. Do you have any tips for speeding up the process? I'm assuming this is some kind of user error on my end.

The "out files" that are blank, are those the temporary files? It's possible they could be blank if they correspond to chromosomes without any reads (or without edits).

Not temporary files, I actually couldn't find where those were being put (so I will specify a location in my next run). I meant the slurm out file. The verbose run had a lot of output (obviously), so I knew that that was being run, but without the output files being created incrementally and nothing being written to the slurm out file I'm not sure how far along the run got.

Thanks so much with all your help on this!

[ahanden]

~500,000 chromosomes? I can't say I've talked to anyone trying to run REDItools on something of that magnitude before.

If you determine the path for the temporary files, you can use the verbose output to determine the state of the analysis. Using the verbose output, when you see a message like [INFO] [REGION=chr1] 12345 total reads, it means that the region/chromosome has finished processing. So the temporary file with 'chr1' in the first column is finished.

My next suggestion would be to run REDItools for each chromosome using the -g/--region option. That would let you pick up where you left off in the analysis.