No read at coverage 1 and below

[jdalino]

Describe the bug
No reads generated in paired-end mode at coverage=1. The output BAM is empty (0 reads). Same for coverage below 1 (0.5, etc.)

To Reproduce
Run NEAT read-simulator in paired-end mode with coverage 1 or less:

reference: GRCh38.fa
read_len: 151
ploidy: 2
coverage: 1
paired_ended: true
fragment_mean: 300
fragment_st_dev: 30
produce_bam: true
threads: 4

Expected behavior
A BAM with 1x

Desktop (please complete the following information):

• OS: Ubuntu 22.04 (AWS EC2)
• NEAT version: latest from GitHub

Possible cause
An issue in the number of reads calculation:
number_reads = number_reads_per_layer * (options.coverage // 2)
Should probably be something like:
number_reads = int(number_reads_per_layer * (options.coverage / 2))

[joshfactorial]

We didn't really account for fractional coverage

[joshfactorial]

I'll take a look

[keshav-gandhi]

Hello @jdalino! I recently performed similar tests on the p-arm of chromosome 22. I will try to replicate your bug.

[joshfactorial]

I was looking into an unrelated (I thought) issue and ran into this bug. @jdalino you are correct, we never enforce coverage as being an integer in the code (and Python doesn't enforce types anyway). So even a number like "5.0" for coverage ends up throwing an indexing error because we are trying to index up to "number_reads." I see your point about the fractional coverage. If we were to allow a 0.5 coverage, then we would want to do a straight division rather than an integer division, but maybe // is the wrong choice since it also precludes a coverage = 1. I think maybe the fix is to use ceil along with regular division (instead of int as in the suggestion) to account for low coverage and small contigs. I will test that out.

[keshav-gandhi]

Replicated the same bug with coverage = 1 and 0 < coverage < 1.