feat: avoid freebayes runnning forever on high coverage low-complexity regions (and parallelize more smoothly)

[dlaehnemann]

We regularly run into situations where freebayes does not seem to parallelize particularly well, with some genomic regions taking forever to process, while all other genomic regions have long finished.

I think this usually comes down to the situation that is described in this issue on the freebayes repository. Namely, we are using the snakemake wrapper for freebayes, which in turn uses the fasta_generate_regions.py script to generate equally sized stretches of reference genome sequence to process in parallel. There are stretches of low complexity genome, in our current case an only-G repeat in this region on chromosome 5 of the human genome (I think it was this), which can accumulate a lot of coverage, even in targeted sequencing data. This will lead to freebayes on such regions running for a very long time on a single CPU (as in days). And if the snakemake workflow reserved a lot more cores for a freebayes job, say 32, all the other CPUs reserved (31 in this case) will sit idle once all other genomic regions have been done. And if you kill the job to avoid wasting the reserved CPU time, you loose all the work that was done by this freebayes job completely.

Thus, I see three possible things we could address here:

1. Set a limit on the number of alleles to consider at any site.
   This is recommended for performance tuning and can be done via the --use-best-n-alleles command line argument. The default here is all, which I guess can easily balloon in low-complexity regions. So I am assuming that even a rather permissive value would be fine here, to just avoid crazy out there numbers. But where exactly this lies, I am unsure -- could be that 20 is good, or 12, or 8, or.... I'll try to get a test case out of our above example and test run times with different values.
2. Improve the generation of regions to process in parallel, most importantly making them contain an equal amount of sequencing data (as opposed to reference genome stretches of equal length).
   Both the issue already mentioned above and another issue on long freebayes run times due to ultra high coverage mention a clear solution for this: Using the coverage_to_regions.py script to get more even coverage across genome regions that are processed. We could for example add this functionality to the snakemake wrapper for freebayes, adding a params: number_of_regions= configuration analogous to the existing params: chunksize= for fasta_generate_regions.py. Alternatively, we could create a separate snakemake wrapper to generate such regions beforehand and/or using different tools than mentioned in those issues.
3. Avoid the loss of all work done if we nevertheless hit a hard to process region.
   This is really annoying, especially for whole genome sequencing data, where freebayes can run for days, even with heavy parallelization. For me, it would thus make sense to split the execution of freebayes into distinct jobs even for a single sample group, so that snakemake can take care of the parallelization (instead of the parallel-based freebayes-parallel script). Thus, each job would run on a single thread and could run independently from all other jobs, and we would have to have a gather-rule to aggregate the resulting VCF/BCF output. For a compute cluster environment with a batch submission system, it would then make sense to group jobs per sample, as they share the same input files, or to group them with the gather/aggregate rule. This is especially important for the case where input and output data are transferred into and out of a compute node via a storage plugin.

[dlaehnemann]

Two more possibilities I see, are to:

4. Use some kind of reference masking of low complexity regions during mapping of reads, so no reads get mapped there.
5. Use some kind of low complexity region BED file to exclude regions from freebayes analysis.

[johanneskoester]

Plan: use option 1 (after testing), and expose it via the config file under params/freebayes.

[dlaehnemann]

For some background on the discussion that we had that led to this decision:

Option 1. with the --use-best-n-alleels flag is the easiest to implement and we will try to use our current case to test for a good value to use, to make runtime reasonable.

Options 4. and 5., the masking of low complexity regions, will not usually work in non-human (or at least non-model) species, as such information will not be available. Also, these might not work well (or even be needed any more) with pangenomes.

Option 2., smoothing out coverage across the determined regions for parallelization, is probably useful in general, but we don't expect it to solve the issue at hand: very pointed (a stretch of a few hundred bp) very high coverage (~30000X) of a low complexity region (almost only Gs). This region could not really sanely be broken down any further.

Option 3. is more involved, as freebayes-parallel takes care to remove duplicate VCF entries at region borders and replicating this behaviour in a full-blown snakemake handling of the situation will probably be more complicated.

[dlaehnemann]

That said, I think that Option 3. might actually still be on the table. Counter to what we thought in our meeting, I don't think that freebayes-parallel actually takes care to not loose variants that might span region borders. It only removes duplicates that might arise there, using vcfuniq from the vcflib, and takes care of proper sorting. But from what I can see from the fasta_generate_regions.py script code, it does not generate overlapping regions:
https://github.com/freebayes/freebayes/blob/a9525ebdf14d6ad8d5e6f432889655539432b450/scripts/fasta_generate_regions.py#L12-L37

Thus, I think that we might actually loose MNVs, insertions, deletions and complex events that exactly span the border from one region to the next. This will be very rare, but I think that doesn't mean it won't happen; and it does mean it will be hard to catch this issue, as we would probably only spot this if we were expecting a particular variant from this category at such a site as a positive control.

Thus, instead of loosing important functionality that we thought was there in the script, I think we might actually eliminate (some very rare) false negatives if we properly implemented this ourselves, with overlapping regions and a sane merging of them. But, as we discussed, this will be quite some work to get sorted out.

[johanneskoester]

Oh, interesting! Indeed, then we could as well do the parallelization ourselves at some point.