Segmentation fault during f5c call-methylation

[csaw2606]

Dear developers,
We are new to ONT sequencing and are running into a segmentation fault error during modification calling and would appreciate any advice you could give us!

f5c call-methylation -b sample.filtered.sorted.bam -g reference_genome.fa -r sample.fastq --slow5 sample.blow5 -t 32 -B 5M --cuda --progress
[meth_main::INFO]�[1;34m Default methylation tsv output format is changed from f5c v0.7 onwards to match latest nanopolish output. Set --meth-out-version=1 to fall back to the old format.�[0m
[init_core::INFO]�[1;34m builtin DNA R10 nucleotide model loaded�[0m
[init_core::INFO]�[1;34m builtin DNA R10 cpg model loaded�[0m
[meth_main::83.235*1.00] 512 Entries (1.3M bases) loaded
[pthread_processor::83.912*1.20] 512 Entries (1.3M bases) processed
..
..
[pthread_processor::46548.537*11.78] 512 Entries (0.9M bases) processed
[meth_main::46548.611*11.78] 512 Entries (0.9M bases) loaded
[sig_handler::ERROR]�[1;31m I regret to inform that a segmentation fault occurred. But at least it is better than a wrong answer.�[0m

The f5c version is v1.6.
The blow5 file was created using slow5tools merge
The output of slow5tools quickcheck is:

[main] cmd: slow5tools quickcheck sample.blow5
[main] real time = 0.005 sec | CPU time = 0.003 sec | peak RAM = 0.003 GB

Thank you for your help!

[hasindu2008]

1. Where did you get these two flags --cuda --progress. f5c does not have thouse options: https://hasindu2008.github.io/f5c/docs/commands. But this is unlikely to be a cause for the segfault

2. Is this a very data large dataset and/or publicaly available - asking because then I could reproduce the problem on my side.

3. Could you please run the tail command on the methylation called output to see around which genomic coordinates caused this issue?

[csaw2606]

Hello,
Thank you for your quick response!

1. We saw it in some other publication but deleted it now from the script
2. We are running the command with a file that has ~470GB and a second that has ~625GB
3. tail meth.tsv

KI270337.1      -       1020    1020    31c8b093-adcb-424f-afd9-ee7f8b3ffcae  -0.95    -200.87 -199.92 1       1       TTATGTACCGTACATAA
KI270337.1      -       558     558     3f519c30-72ff-4ff5-913f-a0fdaba1fc38  0.40     -251.03 -251.43 1       1       TTATATACCGTACATAA
KI270337.1      +       558     558     71d42f38-7ed9-400f-ba66-4959e87249ab  0.26     -199.34 -199.59 1       1       TTATATACCGTACATAA
KI270337.1      +       732     732     71d42f38-7ed9-400f-ba66-4959e87249ab  0.51     -163.33 -163.84 1       1       ACTGTACACGAAATATC
KI270337.1      +       558     558     97d4ce41-c9a5-47d5-9ee6-951f306efc6c  -0.96    -161.16 -160.20 1       1       TTATATACCGTACATAA
KI270337.1      +       732     732     97d4ce41-c9a5-47d5-9ee6-951f306efc6c  -0.21    -136.13 -135.92 1       1       ACTGTACACGAAATATC
KI270337.1      +       1020    1020    97d4ce41-c9a5-47d5-9ee6-951f306efc6c  -2.39    -146.53 -144.14 1       1       TTATGTACCGTACATAA
KI270337.1      +       558     558     25d7d772-5e48-4001-a714-195ea839bab1  -0.39    -193.88 -193.49 1       1       TTATATACCGTACATAA
KI270337.1      +       732     732     25d7d772-5e48-4001-a714-195ea839bab1  0.71     -154.28 -154.99 1       1       ACTGTACACGAAATATC
KI270337.1      +       1020    1020    25d7d772-5e48-4001-a714-195ea839bab1  0.33     -161.35 -161.68 1       1       TTATGTACCGTACATAA

[hasindu2008]

Could you please run as:

f5c call-methylation -b sample.filtered.sorted.bam -g reference_genome.fa -r sample.fastq --slow5 sample.blow5 -t 32 -B 5M -w KI270337.1

giving the -w to make f5c only process KI270337.1, then see if the same segmentation fault occurs?

I am trying to narrow down from which reference contig it is coming from, which we can use to further narrow down to find which set of reads.

[csaw2606]

The output of the command is:
f5c call-methylation -b sample.filtered.sorted.bam -g reference_genome.fa -r sample.fastq --slow5 sample.blow5 -t 32 -B 5M -w KI270337.1

[meth_main::INFO] Default methylation tsv output format is changed from f5c v0.7 onwards to match latest nanopolish output. Set --meth-out-version=1 to fall back to the old format.
[init_core] Iterating over region: KI270337.1

[init_core::INFO] builtin DNA R10 nucleotide model loaded
[init_core::INFO] builtin DNA R10 cpg model loaded
[meth_main::60.454*0.99] 512 Entries (0.9M bases) loaded
[meth_main::61.031*1.04] 109 Entries (0.2M bases) loaded
[pthread_processor::61.500*1.09] 512 Entries (0.9M bases) processed
[pthread_processor::61.702*1.12] 109 Entries (0.2M bases) processed
[meth_main] skipped unmapped: 0, skipped secondary: 0, skipped low_mapq: 0
[meth_main] total entries: 621, qc fail: 9, could not calibrate: 10, no alignment: 8, bad reads: 7
[meth_main] total bases: 1.1 Mbases
[meth_main] Data loading time: 1.153 sec
[meth_main]     - bam load time: 0.015 sec
[meth_main]     - fasta load time: 0.857 sec
[meth_main]     - slow5 load time: 0.281 sec
[meth_main] Data processing time: 1.248 sec
[meth_main] Data output time: 0.001 sec
[main] Version: 1.6
[main] CMD: f5c call-methylation -b sample.filtered.sorted.bam -g reference_genome.fa -r sample.fastq --slow5 sample.blow5 -t 32 -B 5M -w KI270337.1
[main] Real time: 63.329 sec; CPU time: 70.495 sec; Peak RAM: 8.912 GB

[hasindu2008]

Oh it has run. Could you try a few reference contigs just after this one in the BAM file and see if it gives the segfault?

[csaw2606]

For f5c call-methylation -b sample.filtered.sorted.bam -g reference_genome.fa -r sample.fastq --slow5 sample.blow5 -t 32 -B 5M -w KI270394.1 it gives an error:

[meth_main::INFO] Default methylation tsv output format is changed from f5c v0.7 onwards to match latest nanopolish output. Set --meth-out-version=1 to fall back to the old format.
[init_core] Iterating over region: KI270394.1

[init_core::INFO] builtin DNA R10 nucleotide model loaded
[init_core::INFO] builtin DNA R10 cpg model loaded
chromosome      strand  start   end     read_name       log_lik_ratio   log_lik_methylated     log_lik_unmethylated    num_calling_strands     num_motifs    sequence
[meth_main::59.684*1.00] 0 Entries (0.0M bases) loaded
[pthread_processor::59.685*1.00] 0 Entries (0.0M bases) processed
[meth_main] skipped unmapped: 0, skipped secondary: 0, skipped low_mapq: 0
[meth_main] total entries: 0, qc fail: 0, could not calibrate: 0, no alignment: 0, bad reads: 0
[meth_main] total bases: 0.0 Mbases
[meth_main] Data loading time: 0.001 sec
[meth_main]     - bam load time: 0.000 sec
[meth_main]     - fasta load time: 0.000 sec
[meth_main]     - slow5 load time: 0.001 sec
[meth_main] Data processing time: 0.001 sec
[meth_main] Data output time: 0.000 sec
[meth_main::ERROR] all 0 out of 0 reads failed. Check the input files.

[hasindu2008]

That has failed because there were no mapped reads to that reference contig, rather than a segmentation fault.

Let us get 4000 mappings before and after this last reference contig successfully processed and then run f5c individually on each transcript as follows:

samtools view  sample.filtered.sorted.bam | grep "KI270337.1" -A4000 -B4000 | cut -f 3 | sort -u > regions.txt
while read reg; do f5c call-methylation -b sample.filtered.sorted.bam -g reference_genome.fa -r sample.fastq --slow5 sample.blow5 -t 32 -B 5M -w $reg >$reg.txt 2> $reg.log; done < regions.txt

Then we can see which log file would have "segmentation fault".

If there is no segmentation fault, then we could iterate through all the reference contigs one by one like, which will take a bit of time.

cut -f 1 reference_genome.fa.fai | while read reg; do f5c call-methylation -b sample.filtered.sorted.bam -g reference_genome.fa -r sample.fastq --slow5 sample.blow5 -t 32 -B 5M -w $reg >$reg.txt 2> $reg.log; done 

Then we can see which log file would have "segmentation fault".

[hasindu2008]

Another question I wanted to ask is, if you rerun the original f5c command for the whole data set a couple of times, does this segmentation fault reproducibly?

[csaw2606]

Regarding the segmentation fault reproducibilty: we tried several times (with slightly varying parameters) but the error always occured when the output methylation file had exactly 25.1GB. So, it seems to always stop at the same point.

We are still not through looking at every contig. For some that we checked already, we get different warnings and errors:

KI270396.1

[meth_main::WARNING] 1 out of 5 reads missing in FAST5/SLOW5.

KI270411.1

 [meth_main::WARNING] 1 out of 6 reads missing in FAST5/SLOW5.

KI270528.1

 [meth_main::WARNING] 2 out of 2 reads failed. Double check --pore and --rna options.
 [meth_main::WARNING] 2 out of 2 reads missing in FAST5/SLOW5.
 [meth_main::ERROR] all 2 out of 2 reads failed. Check the input files.

KI270322.1

 [meth_main::WARNING] 1 out of 8 reads missing in FAST5/SLOW5.

KI270424.1

[meth_main::WARNING] 1 out of 1 reads failed. Double check --pore and --rna options.
[meth_main::ERROR] all 1 out of 1 reads failed. Check the input files.

KI270735.1

[init_core] Iterating over region: KI270735.1 Error: Invalid json-rpc request

[hasindu2008]

That invalid json-rpc request is likely to be coming from something else - f5c doesn't use json or rpc or print such error.

That 25GB observation, does it happen on other samples too or just this sample?

Is this dataset sharable or publically so I could run on my side?

[csaw2606]

Hello,
we are currently working on making the data accessible. Since this is not so easy for us with such huge data files we are trying to find a solution with our computational department. In the meantime, we try to narrow the files down from which the error comes from.

Yes, the error always occurs when the methylation file has 25.1GB

[hasindu2008]

Would you be able to download this large blow5 file and the corresponding fastq from https://gentechgp.github.io/gtgseq/docs/data.html#na24385-hg002-promethion-data-20x
then see if the error still happens at 25GB?

This is a dataset that I have tested on so many computers with f5c. If the error still comes, then we can conclude it is not relevant to the dataset, instead specific to your system.

[hasindu2008]

Also, are you using the precompiled binaries that we provide on GitHub or did you compile yourself or installed through conda?

[csaw2606]

Hello,
unfortunately we are experiencing issues with the server that we are using to share data. Would it also be a possible option to send you the data on a USB (from central Europe, depending a bit on where you are based)?

Thank you for your effort!

[hasindu2008]

Before that, I think it it is worthwhile checking if this problem is relevant to the dataset or not. Could you try what i suggested in the messages before - downloading a large well tested blow5 file and seeing if the problem can be reproduced.

Also what is the binary download source?