should bcftools norm --multiallelics - respect lexicographical order?

[freeseek]

This is not an issue, but a behavior that maybe could be more optimal.

I have noticed that when splitting multiallelic variants in different VCF files, you don't necessarily end up with variants in the same order. Here is an example:

(echo "##fileformat=VCFv4.1"
echo "##contig=<ID=chr19>"
echo -e "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"
echo -e "chr19\t50359054\t.\tC\tT,A\t.\t.\t.") > in.vcf

Now if I split with bcftools norm --multiallelics -:

$ bcftools norm --multiallelics - --no-version in.vcf
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr19>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
chr19	50359054	.	C	T	.	.	.
chr19	50359054	.	C	A	.	.	.
Lines   total/split/realigned/skipped:	1/1/0/0

But if I further sort the file:

$ bcftools norm --multiallelics - --no-version in.vcf | bcftools sort
Writing to /tmp/bcftools-sort.HmQM3B
Lines   total/split/realigned/skipped:	1/1/0/0
Merging 1 temporary files
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=chr19>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
chr19	50359054	.	C	A	.	.	.
chr19	50359054	.	C	T	.	.	.
Cleaning
Done

The order of the variants has changed. I suppose it would be quite a bit of work to rewrite bcftools norm to sort the variants after splitting them so that you would not require a sort but reporting nevertheless just in case.

I have noticed this as some HTSlib tools (such as IMPUTE5 v1.1.4) do not work if VCFs have different orders in the variants.

[pd3]

Adding this could be quite straightforward by using vcfbuf_t with a new mode SORT_BY_ALLELE mode. This would make all output from norm sorted this way.

[astaric]

+1 for having the output of the norm command sorted.

This is already the case when -a --atom-overlaps . is used, but not when -m --multi-overlaps . is used. There is already code that sorts the output lines, but it only considers ref and pos. Replacing that comparison with the one used by sort command will make the norm command output the lines in the same order as sort command. Patch for the change bellow.

@pd3, not sure I understand you comment. Would you prefer this be a new flag for bcftools norm? If yes, I can try creating a PR that adds such a flag.

diff --git a/vcfnorm.c b/vcfnorm.c
index 05922e37..b0a18cca 100644
--- a/vcfnorm.c
+++ b/vcfnorm.c
@@ -42,6 +42,7 @@ THE SOFTWARE.  */
 #include "abuf.h"
 #include "gff.h"
 #include "regidx.h"
+#include "vcfsort.h"
 
 #define CHECK_REF_EXIT 1
 #define CHECK_REF_WARN 2
@@ -2291,7 +2292,9 @@ static void normalize_line(args_t *args, bcf1_t *line)
         args->lines[i] = bcf_dup(line);
         while ( rbuf_prev(&args->rbuf,&i) )
         {
-            if ( args->lines[i]->rid==args->lines[j]->rid && args->lines[i]->pos > args->lines[j]->pos ) SWAP(bcf1_t*, args->lines[i], args->lines[j]);
+            bcf_unpack(args->lines[i], BCF_UN_STR);
+            bcf_unpack(args->lines[j], BCF_UN_STR);
+            if ( cmp_bcf_pos(&args->lines[i], &args->lines[j]) > 0) SWAP(bcf1_t*, args->lines[i], args->lines[j])
             j = i;
         }
         if ( args->atomize!=SPLIT ) break;
diff --git a/vcfsort.h b/vcfsort.h
new file mode 100644
index 00000000..2dd4dd8c
--- /dev/null
+++ b/vcfsort.h
@@ -0,0 +1,4 @@
+#ifndef __VCFSORT_H__
+#define __VCFSORT_H__
+int cmp_bcf_pos(const void *aptr, const void *bptr);
+#endif

[pd3]

This is now supported via a new option -S, --sort lex for lexicographical order

[astaric]

Thank you!