nservant / HiC-Pro

HiC-Pro: An optimized and flexible pipeline for Hi-C data processing
Other
386 stars 182 forks source link

mapped_2hic_fragments possible memory issue #306

Closed auberginekenobi closed 4 years ago

auberginekenobi commented 4 years ago

Issue: For some of my HiC datasets, mapping reliably finishes but downstream hic processing fails. Version: HiC-Pro 2.11.3-beta

Repro:

time ( HiC-Pro -c config_template_tumor.txt -i ../input -o . -s mapping ) > mapping.log 2>&1
time ( HiC-Pro -c config_template_tumor.txt -i bowtie_results/bwt2 -o . -s proc_hic -s merge_persample -s quality_checks -s build_contact_maps -s ice_norm ) > all_other_steps.log 2>&1
>
--------------------------------------------
Fri Jan 24 14:24:53 PST 2020
Assign alignments to restriction fragments ...
Logs: logs/rawdata/mapped_2hic_fragments.log
make: *** [/datasets/home/11/211/ochapman/bin/HiC-Pro_2.11.3-beta/bin/../scripts//Makefile:152: mapped_2hic_fragments] Error 137

Issue reproduced on HiC-Pro nonparallel mode, using kubernetes pods; 20cpus x 64GB ram and 2cpus x 8G ram

Comments: My limited expertise says this looks like a memory allocation issue with sort, which gets killed according to the logfiles. For 2 samples, on the 20x64 cpu x ram configuration both return a "Killed" message; on 2x8, one now returns a "read failed: <...> cannot allocate memory". Both samples are of comparable size as others which processed successfully. Logs abbreviated below. Thoughts?

Logs:

mapped_2hic_fragments.py -v -S -t 100 -m 100000 -a -f /datasets/home/11/211/ochapman/genomes/HiC-Pro/grch38_1kgmaj.MboI.bed -r bowtie_results/bwt2/rawdata/MB174_grch38_1kgmaj.bwt2pairs.bam -o hic_results/data/rawdata
## overlapMapped2HiCFragments.py
## mappedReadsFile= bowtie_results/bwt2/rawdata/MB174_grch38_1kgmaj.bwt2pairs.bam
## fragmentFile= /datasets/home/11/211/ochapman/genomes/HiC-Pro/grch38_1kgmaj.MboI.bed
## minInsertSize= None
## maxInsertSize= None
## minFragSize= 100
## maxFragSize= 100000
## allOuput= True
## SAM ouput= True
## verbose= True

## Loading Restriction File Intervals ' /datasets/home/11/211/ochapman/genomes/HiC-Pro/grch38_1kgmaj.MboI.bed '...
Warning :  1457166 fragment(s) outside of range and discarded.  5748545  remaining.
load_restriction_fragment function took 59691.402 ms
## Opening SAM/BAM file ' bowtie_results/bwt2/rawdata/MB174_grch38_1kgmaj.bwt2pairs.bam '...
## Classifying Interactions ...
<... Numbers ...>
## Sorting valid interaction file ...
LANG=en; sort -T tmp -k2,2V -k3,3n -k5,5V -k6,6n -o hic_results/data/rawdata/MB174_grch38_1kgmaj.bwt2pairs.validPairs hic_results/data/rawdata/MB174_grch38_1kgmaj.bwt2pairs.validPairs
sort: read failed: hic_results/data/rawdata/MB174_grch38_1kgmaj.bwt2pairs.validPairs: Cannot allocate memory
/HiC-Pro_2.11.3-beta/scripts/mapped_2hic_fragments.py -v -S -t 100 -m 100000 -a -f /datasets/home/11/211/ochapman/genomes/HiC-Pro/grch38_1kgmaj.MboI.bed -r bowtie_results/bwt2/rawdata/MB106_grch38_1kgmaj.bwt2pairs.bam -o hic_results/data/rawdata
## overlapMapped2HiCFragments.py
## mappedReadsFile= bowtie_results/bwt2/rawdata/MB106_grch38_1kgmaj.bwt2pairs.bam
## fragmentFile= /datasets/home/11/211/ochapman/genomes/HiC-Pro/grch38_1kgmaj.MboI.bed
## minInsertSize= None
## maxInsertSize= None
## minFragSize= 100
## maxFragSize= 100000
## allOuput= True
## SAM ouput= True
## verbose= True

## Loading Restriction File Intervals ' /datasets/home/11/211/ochapman/genomes/HiC-Pro/grch38_1kgmaj.MboI.bed '...
Warning :  1457166 fragment(s) outside of range and discarded.  5748545  remaining.
load_restriction_fragment function took 59820.016 ms
## Opening SAM/BAM file ' bowtie_results/bwt2/rawdata/MB106_grch38_1kgmaj.bwt2pairs.bam '...
## Classifying Interactions ...
< ... Numbers ...>
## Sorting valid interaction file ...
LANG=en; sort -T tmp -k2,2V -k3,3n -k5,5V -k6,6n -o hic_results/data/rawdata/MB106_grch38_1kgmaj.bwt2pairs.validPairs hic_results/data/rawdata/MB106_grch38_1kgmaj.bwt2pairs.validPairs
/datasets/home/home-00/11/211/ochapman/bin/HiC-Pro_2.11.3-beta/scripts/hic.inc.sh: line 86:  1169 Killed                  sort -T tmp -k2,2V -k3,3n -k5,5V -k6,6n -o hic_results/data/rawdata/MB106_grch38_1kgmaj.bwt2pairs.validPairs hic_results/data/rawdata/MB106_grch38_1kgmaj.bwt2pairs.validPairs
nservant commented 4 years ago

Thanks for the detailed issue. Indeed it looks like a sort issue ... again. In the mapped_2hic_fragments.sh, could you try to add a -S 80% in the sort command. This should define the amount of memory to use. Note that in this case, the sort command may swap with your TMP_DIR folder.

auberginekenobi commented 4 years ago

Hi Nicolas, Alas, same error:

## Sorting valid interaction file ...
LANG=en; sort -T tmp -S 80% -k2,2V -k3,3n -k5,5V -k6,6n -o hic_results/data/rawdata/MB277_grch38_1kgmaj.bwt2pairs.validPairs hic_results/data/rawdata/MB277_grch38_1kgmaj.bwt2pairs.validPairs
/datasets/home/home-00/11/211/ochapman/bin/HiC-Pro_2.11.3-beta/scripts/hic.inc.sh: line 86:  1173 Killed                  sort -T tmp -S 80% -k2,2V -k3,3n -k5,5V -k6,6n -o hic_results/data/rawdata/MB277_grch38_1kgmaj.bwt2pairs.validPairs hic_results/data/rawdata/MB277_grch38_1kgmaj.bwt2pairs.validPairs

I have been able to get around the issue by feeding HiC-Pro 128G of RAM, but this doesn't seem like a general solution...

nservant commented 4 years ago

In theory, using -T and -S are the best ways to improve unix sort performance. Unfortunately, I do not have any other ideas ... sorry. How many reads do you have on your fastq files ? Note that using HiC-Pro, you can also split reads in chunks ... all chunks will be processed in parallel until the final merge to build the maps. So I would say that you will still have a big sort at the end, but the processing of each chunck should be much faster. N

auberginekenobi commented 4 years ago

Hi Nicolas, I presume that the mapped_2hic_fragments.log file is outputting the number of records to be sorted. If true, 2,160,500,000 reads. I had hoped to avoid using the -p option because I'm running on kubernetes, not an HPC scheduler, and I didn't see documentation for chunks on a non-cluster configuration.

nservant commented 4 years ago

You can still slip the reads into chunks without cluster option to avoid sorting memory issue. The only point is that the read chunks will be run one by one ... which may take a lot of time. Best