Add infrastructure for validating split files

Author: jeromekelleher

bio2zarr/vcf.py

@@ -199,7 +199,6 @@ def make_field_def(name, vcf_type, vcf_number):
 
 
 def scan_vcf(path, target_num_partitions):
-    logger.debug(f"Scanning {path}")
     with vcf_utils.IndexedVcf(path) as indexed_vcf:
         vcf = indexed_vcf.vcf
         filters = [
@@ -238,6 +237,8 @@ def scan_vcf(path, target_num_partitions):
             pass
 
         regions = indexed_vcf.partition_into_regions(num_parts=target_num_partitions)
+        logger.info(
+            f"Split {path} into {len(regions)} regions (target={target_num_partitions})")
         for region in regions:
             metadata.partitions.append(
                 VcfPartition(

validation-data/Makefile

@@ -15,9 +15,14 @@ all: 1kg_2020_chr20.bcf.csi \
 	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.csi
+	1kg_p1_all_chr6.vcf.gz.csi\
+	1kg_2020_chr20.vcf.gz.2.split\
+	1kg_2020_others.vcf.gz.2.split \
+	1kg_2020_others.vcf.gz.5.split \
+	1kg_2020_freebayes.vcf.gz.5.split \
 
 .PRECIOUS: %.bcf
+.PRECIOUS: %.split
 .PRECIOUS: %.vcf.gz
 
 # Modern 1000 genomes data
@@ -38,6 +43,12 @@ all: 1kg_2020_chr20.bcf.csi \
 %.vcf.gz.tbi: %.vcf.gz
 	tabix $<
 
+%.2.split: %
+	./split.sh $< 2
+
+%.5.split: %
+	./split.sh $< 5
+
 %.csi: %
 	bcftools index $<
 

validation-data/split.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+VCF=$1
+NUM_PARTS=$2
+
+regions=$( PYTHONPATH="../" python3 -m bio2zarr vcf-partition $VCF -n $NUM_PARTS )
+dest=$VCF.$NUM_PARTS.split
+echo $dest
+rm -fR $dest
+mkdir $dest
+# TODO figure out how to do this with xargs
+for r in $regions; do
+	f=$dest/$r.vcf.gz
+	bcftools view $VCF -r $r -Oz > $f
+	bcftools index $f
+done

validation.py

@@ -12,6 +12,7 @@
 # the original unsplit vs split files in there following some
 # naming conventions.
 
+
 @click.command
 @click.argument("vcfs", nargs=-1)
 @click.option("-p", "--worker-processes", type=int, default=1)
@@ -21,17 +22,27 @@
 def cli(vcfs, worker_processes, force):
     data_path = pathlib.Path("validation-data")
     if len(vcfs) == 0:
-        vcfs = list(data_path.glob("*.vcf.gz")) + list(data_path.glob("*.bcf"))
+        vcfs = (
+            list(data_path.glob("*.vcf.gz"))
+            + list(data_path.glob("*.bcf"))
+            + list(data_path.glob("*.split"))
+        )
     else:
         vcfs = [pathlib.Path(f) for f in vcfs]
     tmp_path = pathlib.Path("validation-tmp")
     tmp_path.mkdir(exist_ok=True)
     for f in vcfs:
-        print(f)
+        print("Validate", f)
+        if f.is_dir():
+            files = list(f.glob("*.vcf.gz")) + list(f.glob("*.bcf"))
+            source_file = f.with_suffix("").with_suffix("")
+        else:
+            files = [f]
+            source_file = f
         exploded = tmp_path / (f.name + ".exploded")
         if force or not exploded.exists():
             vcf.explode(
-                [f],
+                files,
                 exploded,
                 worker_processes=worker_processes,
                 show_progress=True,
@@ -51,7 +62,7 @@ def cli(vcfs, worker_processes, force):
                 show_progress=True,
             )
 
-        vcf.validate(f, zarr, show_progress=True)
+        vcf.validate(source_file, zarr, show_progress=True)
 
 
 if __name__ == "__main__":