consensus only runs on the first sequence (chrom) of fasta file

[DrYak]

when running:

bcftools consensus --output results/SPFL100652/references/consensus.bcftools.fasta --chain results/SPFL100652/references/consensus.bcftools.chain --fasta-ref references/concat_EPI_ISL_412866_EPI_ISL_1653999.fasta --mark-del - --mask results/SPFL100652/references/coverage_mask_lowcoverage.bed --mask-with n -H A -i "INFO/AD[0]<INFO/AD[*]" results/SPFL100652/references/consensus.bcftools.bcf.gz

The resulting .fasta and .chain files only contain a consensus of the first sequence (the first chrom) - EPI_ISL_412866 (RSV-A) - present in the reference file. The second sequence - EPI_ISL_1653999 (RSV-B) - is missing from the output, even if it is present in the reference fasta, the pileup bcf and the mask bed inputs.
It happens in both bcftools 1.20 (that we currently use) and bcftools 1.22 (the most recent avilable on bioconda).

This ZIP file contains all input and outputs mentionned above: SPFL100652.zip

Am I doing something wrong? Passing the wrong parameters or something?

(This could be overlapping with #2052 maybe ? That issues also mentions some sequence not being present in the output FASTA)

Context:

• This is the last step of the rule in V-pipe that genereates consensus sequences using bcftools.
• These are swab test samples of RSV viruses (sub-group unknown at this stage) and we're aligning the reads against a fasta file with references that includes sequence for both sub-group A and sub-group B RSV.
• The expected behaviour is that we get a consensus with two sequences (potentially one sequence being mostly or entirely n depending on which sub-group was present in the swab sample)
• Currently we only get EPI_ISL_412866 (RSV-A) outputs.

[DrYak]

(Note: merely swaping the order of sequence in the fast reference file causes EPI_ISL_1653999 (RSV-B) only to be output instead and now EPI_ISL_412866 (RSV-A) is missing instead).

[pd3]

Thank you for the issue and the test case, I believe this is now fixed.

[DrYak]

Will test now...

[DrYak]

...both the fix (92c31be) and the work-around (remove blank lines in .fasta for bcftools version 1.22) confirmed to work.

Thank you very much for the very fast handling of our issue!
🙏

[DrYak]

BTW: Do you have an idea when you would be tagging the next version 1.23 of bcftools?

(And thus when we'll get the autobumped version in bioconda?)

[pd3]

It will take a while, likely sometime in autumn / winter