bcftools call Doesn't Seem to Understand Negated Sample Argument

[fbschmtt]

Hi bcftools Team,
thank you very much for this very useful tool!
I am currently trying to run the sequence of mpileup | call using a sample level grouping. I want to use all present samples in the original BAM file, for mpileup I used -s "^notinhere", is there a more elegant solution? I tried running with -s without any specified sample list, however this seemed to result in messed up argument parsing.
Also, once I pipe the output into bcftools call, the -s "^notinhere" strategy doesn't seem to work, resulting in the following error:
"Warning: No such sample in the VCF: ^notinhere"
If I manually specify all the samples in the original bam file it seems to work. Following is the full command (quoting / unquoting the sample argument doesn't seem to have an effect):
bcftools mpileup -Ou -s "^notinhere" -f .../genome.fa STAR_out/Aligned.sortedByCoord.out.bam | bcftools call -s "^notinhere" -mv -Oz -o calls_samples.vcf.gz

Version Used: 1.7, though I also tested with 1.21 and got the same errors.

Best

[pd3]

I don't understand what are you trying to do. The negation is used to exclude all but the listed samples. From your example it seems the sample "notinhere" is not present in the BAM file, therefore it cannot be excluded and the program does the correct thing when complaining.

If the samples option is not given, the program does per sample calling, if @RG is present in the BAM.

Probably best if you could provide a small test case for us to understand the problem.

[fbschmtt]

Thank you for the reply.
Is there an easy way to generate test data? As I work on human data sharing my original data isn't possible.

My understanding was that there can be @rg as well as @sm and these can be used to group reads differently, i.e. multiple sequencing runs from one sample / patient. I got this impression from this issue #1018

I think the parsing isn't working because the error message still contains the negation, i.e. it says "^notinhere" not in the file and not "notinhere".
For mpileup the ^notinhere seems to work properly, the idea is if there are sample a,b,c in the file and I say I want ^d it should pick all samples, as I didn't find an option to set the grouping to sample instead of read group level. Or is the default behavior to group by sample anyways?

[pd3]

There is the script https://github.com/pd3/mpileup-tests/blob/main/misc/create-bam-test which you could use to create a small bam slice and anonymize the data, possibly with some manual editing of the RG and SM tags

[fbschmtt]

Sorry for the delay in response. I now managed to create some synthetic data to show the behavior I mean.
While testing I noticed that my approach was unnecessary and pretty much corresponds to the default behavior without any sample arguments. So for me this doesn't have any priority.

However, the sample argument appears to behave strangely for bcftools call, unless I don't understand what it is supposed to do. As I understand it, it is basically supposed to filter samples? So filtering samples at the mpileup or call stage should lead to the same results? Using the following command, the sample option for mpielup works without issue, so this would be an easy workaround.
bcftools mpileup -Ou -f nod2.fa neat_simulated_snps_nod2.bam | bcftools call -mv -Ov | less

When specifying different sample arguments to call, the output appears to not work correctly. Here are some arguments which produce unexpected outputs:

• When using -s "sample1" on call, the vcf file still contains a sample2 column, although this appears to not contain correct info (It is different from not specifying -s at all).
• When using -s "sample2" on call, the vcf file only contains a sample2 column, no sample1.
• When using -s "^sample1" on call, the command fails, as it doesn't appear to parse the negation correctly. Putting the argument to mpileup works as expected and only returns sample2.

Thanks for working on this.

test_data_call_sample_argument.zip

[pd3]

Thank you for the bug report and the test case, this is now fixed. Please try it out and let me know if you encounter any problems

[fbschmtt]

Thank you for the quick fix. Everthing is working as expected now.