Fix bug in Tabix handling

Author: jeromekelleher

bio2zarr/vcf_utils.py

@@ -63,7 +63,7 @@ def __post_init__(self):
             assert self.start > 0
         if self.end is not None:
             self.end = int(self.end)
-            assert self.end > 0
+            assert self.end > self.start
 
     def __str__(self):
         s = f"{self.contig}"
@@ -641,15 +641,17 @@ def partition_into_regions(
                 next_contig = self.sequence_names[region_contigs[i + 1]]
                 next_start = region_starts[i + 1]
                 end = next_start - 1  # subtract one since positions are inclusive
+                # print("next_start", next_contig, next_start)
                 if next_contig == contig:  # contig doesn't change
                     regions.append(Region(contig, start, end))
                 else:
                     # contig changes, so need two regions (or possibly more if any
                     # sequences were skipped)
                     regions.append(Region(contig, start))
                     for ri in range(region_contigs[i] + 1, region_contigs[i + 1]):
-                        regions.append(self.sequence_names[ri])
-                    regions.append(Region(next_contig, 1, end))
+                        regions.append(Region(self.sequence_names[ri]))
+                    if end >= 1:
+                        regions.append(Region(next_contig, 1, end))
 
         # Add any sequences at the end that were not skipped
         for ri in range(region_contigs[-1] + 1, len(self.sequence_names)):

tests/data/vcf/multi_contig.vcf.gz

@@ -32,6 +32,7 @@
         ("1kg_2020_chrM.bcf.csi", {"chrM": 23}),
         ("1kg_2020_chr20_annotations.bcf.csi", {"chr20": 21}),
         ("NA12878.prod.chr20snippet.g.vcf.gz.tbi", {"20": 301778}),
+        ("multi_contig.vcf.gz.tbi", {str(j): 933 for j in range(5)}),
     ],
 )
 def test_index_record_count(index_file, expected):
@@ -53,6 +54,7 @@ def test_index_record_count(index_file, expected):
         ("1kg_2020_chrM.bcf.csi", ["chrM:1-"]),
         ("1kg_2020_chr20_annotations.bcf.csi", ["chr20:49153-"]),
         ("NA12878.prod.chr20snippet.g.vcf.gz.tbi", ["20:1-"]),
+        ("multi_contig.vcf.gz.tbi", ["0:1-"] + [str(j) for j in range(1, 5)]),
     ],
 )
 def test_partition_into_one_part(index_file, expected):
@@ -63,6 +65,18 @@ def test_partition_into_one_part(index_file, expected):
     assert [str(r) for r in regions] == expected
 
 
+def test_tabix_multi_chrom_bug():
+    index_file = "multi_contig.vcf.gz.tbi"
+    vcf_path = data_path / (".".join(list(index_file.split("."))[:-1]))
+    indexed_vcf = vcf_utils.IndexedVcf(vcf_path, data_path / index_file)
+    regions = indexed_vcf.partition_into_regions(num_parts=10)
+    # An earlier version of the code returned this, i.e. with a duplicate
+    # for 4 with end coord of 0
+    # ["0:1-", "1", "2", "3", "4:1-0", "4:1-"]
+    expected = ["0:1-", "1", "2", "3", "4:1-"]
+    assert [str(r) for r in regions] == expected
+
+
 class TestCsiIndex:
     @pytest.mark.parametrize(
         "filename",

tests/data/vcf/multi_contig.vcf.gz.tbi

@@ -1,13 +1,21 @@
 all: 1kg_2020_chr20.bcf.csi \
 	1kg_2020_chr20.vcf.gz.tbi \
+	1kg_2020_chr20.vcf.gz.csi \
 	1kg_2020_chr20_annotations.vcf.gz.tbi \
+	1kg_2020_chr20_annotations.vcf.gz.csi \
 	1kg_2020_chr20_annotations.bcf.csi\
 	1kg_2020_others.bcf.csi \
 	1kg_2020_others.vcf.gz.tbi \
+	1kg_2020_others.vcf.gz.csi \
+	1kg_2020_freebayes.vcf.gz.tbi \
+	1kg_2020_freebayes.vcf.gz.csi \
+	1kg_2020_freebayes.bcf.csi \
 	1kg_p3_all_chr1.bcf.csi \
 	1kg_p3_all_chr1.vcf.gz.tbi\
+	1kg_p3_all_chr1.vcf.gz.csi\
 	1kg_p1_all_chr6.bcf.csi\
-	1kg_p1_all_chr6.vcf.gz.tbi
+	1kg_p1_all_chr6.vcf.gz.tbi \
+	1kg_p1_all_chr6.vcf.gz.csi
 
 .PRECIOUS: %.bcf
 .PRECIOUS: %.vcf.gz
@@ -17,6 +25,7 @@ all: 1kg_2020_chr20.bcf.csi \
 1KG_2020_GENOTYPES_URL="${1KG_2020_BASE}chr20.recalibrated_variants.vcf.gz"
 1KG_2020_ANNOTATIONS_URL="${1KG_2020_BASE}chr20.recalibrated_variants.annotated.vcf.gz"
 1KG_2020_OTHERS_URL="${1KG_2020_BASE}others.recalibrated_variants.vcf.gz"
+1KG_2020_FREEBAYES="http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20200515_EBI_Freebayescalls/ALL.chrX.freebayes.20200518.snps_indels.NYhc.GRCh38.vcf.gz"
 
 # 1000 genomes phase 3
 1KG_P3_ALL_URL=http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/release/20190312_biallelic_SNV_and_INDEL/ALL.chr1.shapeit2_integrated_snvindels_v2a_27022019.GRCh38.phased.vcf.gz
@@ -29,7 +38,7 @@ all: 1kg_2020_chr20.bcf.csi \
 %.vcf.gz.tbi: %.vcf.gz
 	tabix $<
 
-%.bcf.csi: %.bcf
+%.csi: %
 	bcftools index $<
 
 %.bcf: %.vcf.gz
@@ -41,6 +50,9 @@ all: 1kg_2020_chr20.bcf.csi \
 1kg_2020_others.vcf.gz:
 	bcftools view ${1KG_2020_OTHERS_URL} | head -n 10000 | bcftools view -O z -o $@
 
+1kg_2020_freebayes.vcf.gz:
+	bcftools view ${1KG_2020_FREEBAYES} | head -n 10000 | bcftools view -O z -o $@
+
 1kg_2020_chr20_annotations.vcf.gz:
 	bcftools view ${1KG_2020_ANNOTATIONS_URL} | head -n 100000 | bcftools view -O z -o $@
 

tests/test_vcf_utils.py

@@ -12,6 +12,8 @@
 @click.argument("vcfs", nargs=-1)
 @click.option("-p", "--worker-processes", type=int, default=1)
 @click.option("-f", "--force", is_flag=True, default=False)
+# TODO add options for verbose and to force the use of a given
+# index file
 def cli(vcfs, worker_processes, force):
     data_path = pathlib.Path("validation-data")
     if len(vcfs) == 0:

validation-data/Makefile

@@ -0,0 +1,30 @@
+# Development script to automate generating multi-contig
+# VCFs of various lenghts for testing the indexing code.
+# TODO should move to a "scripts" directory or something
+import sys
+import click
+
+
+def write_header(num_contigs):
+    click.echo("##fileformat=VCFv4.2")
+    click.echo(f"##source={' '.join(sys.argv)}")
+    click.echo('##FILTER=<ID=PASS,Description="All filters passed">')
+    for contig in range(num_contigs):
+        click.echo(f"##contig=<ID={contig}>")
+    header = "\t".join(["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO"])
+    click.echo(header)
+
+
+@click.command
+@click.argument("contigs", type=int)
+@click.argument("rows-per-contig", type=int)
+def cli(contigs, rows_per_contig):
+    write_header(contigs)
+    for j in range(contigs):
+        for k in range(rows_per_contig):
+            pos = str(k + 1)
+            click.echo("\t".join([str(j), pos, "A", "T", ".", ".", ".", "."]))
+
+
+if __name__ == "__main__":
+    cli()