beta2: segfault in utgcns and bug in Snakemake

[ptrebert]

Hi,

I successfully processed two samples with the beta2 release (installed via bioconda), but for a third one, I am running into the following issue:

-- Opening output FASTA file '../7-consensus/packages/part054.fasta'.
--
-- Computing consensus for b=0 to e=4294967295 with errorRate 0.2000 (max 0.4000) and minimum overlap 2500
--

generateTemplateStitch()-- COPY READ read #0 62 (len=2510432 to 0-2510432)

Failed with 'Segmentation fault'; backtrace (libbacktrace):
./packages/part054.sh: line 13:  1325 Segmentation fault      /gpfs/project/projects/medbioinf/projects/hgsvc/sig_chrY/run_folder/.snakemake/conda/450628c96885b8740733903fdcdb9b04/lib/verkko/bin/utgcns -V -V -V -threads 8 -import ../7-consensus/packages/part054.cnspack -A ../7-consensus/packages/part054.fasta -C 2 -norealign -maxcoverage 50 -e 0.20 -l 2500 -edlib

I am not sure if the above simply triggered the below error in Snakemake (known, and supposedly fixed, see snakemake/snakemake#823), or if that is indicating another problem in the pipeline:

mkdir -p packages

cat > ./packages/part054.sh <<EOF
#!/bin/sh
set -e

/gpfs/project/projects/medbioinf/projects/hgsvc/sig_chrY/run_folder/.snakemake/conda/450628c96885b8740733903fdcdb9b04/lib/verkko/bin/utgcns \\
    -V -V -V \\
    -threads 8 \\
    -import ../7-consensus/packages/part054.cnspack \\
    -A ../7-consensus/packages/part054.fasta \\
    -C 2 -norealign \\
    -maxcoverage 50 \\
    -e 0.20 \\
    -l 2500 \\
    -edlib
EOF

chmod +x ./packages/part054.sh

./packages/part054.sh > ../7-consensus/packages/part054.err 2>&1
        
        (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)
    cluster_jobid: 2049544.hpc-batch

Error executing rule generateConsensus on cluster (jobid: 69, external: 2049544.hpc-batch, jobscript: /gpfs/project/projects/medbioinf/projects/hgsvc/sig_chrY/run_folder/output/hybrid/verkko/AFR-MSL-SL72-M_HG03579.HIFIRW.ONTUL.na.wg/.snakemake/tmp.i2kert9m/snakejob.gener
Job failed, going on with independent jobs.
Exiting because a job execution failed. Look above for error message
BUG: Out of jobs ready to be started, but not all files built yet. Please check https://github.com/snakemake/snakemake/issues/823 for more information.
Remaining jobs:
 - generateConsensus: 7-consensus/packages/part054.fasta
 - combineConsensus: 7-consensus/unitig-popped.fasta
 - verkko: assembly.homopolymer-compressed.gfa, assembly.homopolymer-compressed.layout, assembly.hifi-coverage.csv, assembly.ont-coverage.csv, assembly.fasta
The code used to generate one or several output files has changed:
    To inspect which output files have changes, run 'snakemake --list-code-changes'.
    To trigger a re-run, use 'snakemake -R $(snakemake --list-code-changes)'.
The input used to generate one or several output files has changed:
    To inspect which output files have changes, run 'snakemake --list-input-changes'.
    To trigger a re-run, use 'snakemake -R $(snakemake --list-input-changes)'.
The params used to generate one or several output files has changed:
    To inspect which output files have changes, run 'snakemake --list-params-changes'.
    To trigger a re-run, use 'snakemake -R $(snakemake --list-params-changes)'.

I can share the input data via Globus if needed.

Best,
Peter

[skoren]

Did it try the job twice, if so I think snakemake worked correctly despite the weird error reporting. It stopped since the output wasn't successful. If you can share the offending input, 7-consensus/packages/part054.cnspack, either to globus or following the instruction on the FAQ here: https://canu.readthedocs.io/en/latest/faq.html#how-can-i-send-data-to-you

[ptrebert]

yes, the run was restarted once - I read that as a "should no longer happen" error, but ok, the segfault seems more critical for now.
I shared the data via Globus:

T2T:/Verkko_test-data/github-issue53/

Globus reported a few connection problems, but the file seems complete (size-wise, at least); just ping me if it has nevertheless been corrupted during the transfer. Thanks!

[skoren]

Cute, the issue here is canu's consensus has a max read limit of 2,097,151bp (2^21 bits) and your read is too long. We can either skip too long reads on load or allow canu to store longer reads.

[ptrebert]

"640 kB ought to be enough for anybody."

Quite unexpected - I assume that limit can only be changed when compiling Canu?

[skoren]

There's actually no reason to enforce/use that limit in consensus. It affects overlap and sequence data structures but ONT reads are never loaded into those. If you built verkko from source, can you try the following patch in the canu code:

--- a/src/utgcns/unitigConsensus.C
+++ b/src/utgcns/unitigConsensus.C
@@ -298,7 +298,7 @@ unitigConsensus::generateTemplateStitch(void) {
   char        *fragment = seq->getBases();
   uint32       readLen  = seq->length();

-  uint32       tigmax = AS_MAX_READLEN;  //  Must be at least AS_MAX_READLEN, else resizeArray() could fail
+  uint32       tigmax = std::max(readLen, AS_MAX_READLEN);  //  Must be at least AS_MAX_READLEN, else resizeArray() could fail
   uint32       tiglen = 0;
   char        *tigseq = NULL;

It is running on my end so if you can't change the code, I can also share the final fasta assuming it runs to completion.

[ptrebert]

It is running on my end so if you can't change the code, I can also share the final fasta assuming it runs to completion.

As a right now solution, that would be much appreciated. If that problem occurs with just one sample, then I would not have to go through the process of deploying a source build of Verkko on our cluster.

[skoren]

I put the finished fasta on globus under T2T:/Verkko_test_data_results/github-issue53/

[ptrebert]

Thanks a lot, it all worked as expected with one exception: the hifi-corrected.fasta is not being dumped (not a problem for me right now).
I assume the above patch in unitigConsensus will then be included in the next (major) Verkko release as bug fix?

[skoren]

Yes it will. In the meantime, you can check/trim any read over 2mb in your ONT datasets to avoid hitting it. In the case above, the read was the start of the contig and we trim back to at least 2x coverage so most of it was erased anyway from the final consensus.

[ptrebert]

Great, and thanks for the advice.