How can I bypass 'out of memory' problem using a HISAT2 for human genome indexing?

Question

I am trying to perform a HISAT2 indexing using a GCF_000001405.39_GRCh38.p13_genomic.gtf and GCF_000001405.39_GRCh38.p13_genomic.fna files downloaded from https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.39 using the procedure described in Pertea et al. (https://www.ncbi.nlm.nih.gov/pubmed/27560171), e.g. first I create extract_splice_sites.py and extract_exons.py and then my code looks like this

#!/bin/bash
#SBATCH --time=02:30:00
#SBATCH --account=def-myusername
#SBATCH --mem=160000M
#SBATCH [email protected]
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
module load hisat2/2.1.0

hisat2-build --ss imam.ss --exon imam.exon GCF_000001405.39_GRCh38.p13_genomic.fna imam002

I am trying to receive a set of imam002*.ht2 files but every time the procedure ends with a Ran out of memory message. From here https://github.com/griffithlab/rnaseq_tutorial/wiki/Indexing I learned that I must allocate 160Gb of RAM if I want to index the entire genome. Even when I do that and submit my script via sbatch on a HPC I still got this

slurmstepd: error: Detected 1 oom-kill event(s) in step 30827017.batch cgroup. 

Just in case here's the window that I have when I don't sbatch my script, I still have an error

Inactive Modules:
  1) openmpi/2.1.1

The following have been reloaded with a version change:
  1) gcccore/.5.4.0 => gcccore/.7.3.0     2) icc/.2016.4.258 => icc/.2018.3.222     3) ifort/.2016.4.258 => ifort/.2018.3.222     4) imkl/11.3.4.258 => imkl/2018.3.222     5) intel/2016.4 => intel/2018.3

Settings:
  Output files: "imam002.*.ht2"
  Line rate: 7 (line is 128 bytes)
  Lines per side: 1 (side is 128 bytes)
  Offset rate: 4 (one in 16)
  FTable chars: 10
  Strings: unpacked
  Local offset rate: 3 (one in 8)
  Local fTable chars: 6
  Local sequence length: 57344
  Local sequence overlap between two consecutive indexes: 1024
  Endianness: little
  Actual local endianness: little
  Sanity checking: disabled
  Assertions: disabled
  Random seed: 0
  Sizeofs: void*:8, int:4, long:8, size_t:8
Input files DNA, FASTA:
  GCF_000001405.39_GRCh38.p13_genomic.fna
Reading reference sizes
  Time reading reference sizes: 00:00:23
Calculating joined length
Writing header
Reserving space for joined string
Joining reference sequences
  Time to join reference sequences: 00:00:13
  Time to read SNPs and splice sites: 00:00:15
  Ran out of memory; automatically trying more memory-economical parameters.
./hisat2_index.sh: line 24:   943 Killed                  hisat2-build --ss imam.ss --exon imam.exon GCF_000001405.39_GRCh38.p13_genomic.fna imam002

When I use the same script for another genome (not human) - I don't have this problem.

Is it really a problem of RAM and, if it is, how is usually human genome gets indexed by HISAT2? Is it even possible?

Answer 1

You haven't allocated the requisite 160GB of RAM required. Or rather, you allocated 160GB of memory, which is what you were told to do, but they most likely meant 160GiB.

Computers use binary, not decimal, so a gigabyte is 1024 megabytes. This is the same annoying mix-up people run into when their hard drives don't have the space the marketing department says it does. By requesting 160000M in your job script, you actually requested a few gigabytes less than what your source says you need.

$$ 160000 \text{ megabytes} \times \frac{1 \text{ gigabyte}}{1024 \text{ megabytes}} = 156.25 \text{ gigabytes} $$

$$ 160 \text{ gigabytes} \times \frac{1024 \text{ megabytes}}{1 \text{ gigabyte}} = 163840 \text{ megabytes} $$

If "160GB" is what you need, then 163840M is the correct value for the job script.

I would honestly suggest cutting the Gordian knot and not dealing with this "gigabytes vs. gibibytes" ordeal by simply using the G suffix.¹

#SBATCH --mem=160G