vcfdist v3.0.0

.github/workflows/build.yml → .github/workflows/test.yml

@@ -1,20 +1,24 @@
-name: build
+name: Tests
 
 on:
   push:
-    branches: [ "master" ]
+    branches: [ "master", "dev" ]
   pull_request:
-    branches: [ "master" ]
+    branches: [ "master", "dev" ]
 
 jobs:
-  build:
+  Tests:
 
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-24.04
 
     steps:
-    - uses: actions/checkout@v3
-    - name: install htslib
+    - uses: actions/checkout@v4
+      with:
+        ref: ${{ github.ref }}
+    - name: Install htslib
       run: |
+       sudo apt-get update
+       sudo apt-get install -y libbz2-dev liblzma-dev
        wget https://github.com/samtools/htslib/releases/download/1.17/htslib-1.17.tar.bz2
        tar -xvf htslib-1.17.tar.bz2
        cd htslib-1.17
@@ -23,11 +27,17 @@ jobs:
        make install
        echo LIBRARY_PATH=`pwd`:$LIBRARY_PATH >> $GITHUB_ENV
        echo LD_LIBRARY_PATH=`pwd`:$LD_LIBRARY_PATH >> $GITHUB_ENV
-    - name: make vcfdist
+    - name: Make vcfdist
       run: |
        cp -r htslib-1.17/htslib src
        echo $LIBRARY_PATH
        echo $LD_LIBRARY_PATH
        cd src
        make
        ./vcfdist -v
+    - name: Test vcfdist
+      run: |
+       sudo apt-get install -y libgtest-dev
+       pip install pytest pytest-workflow
+       cd tests
+       pytest -vv

.gitignore

@@ -1,4 +1,5 @@
 # ignore certain file types
+*.a
 .*.swp
 *.o
 *.log
@@ -13,8 +14,15 @@
 *.gz
 *.gprof
 *.png
-src/vcfdist
+*.json
+*.log
+*.aux
 src/test
+__pycache__
+
+# ignore binaries
+src/vcfdist
+tests/unit/build/test_vcfdist
 
 # ignore analysis
 analysis/*/*/

README.md

@@ -1,5 +1,6 @@
 # vcfdist: benchmarking phased variant calls
-![build](https://github.com/timd1/vcfdist/actions/workflows/build.yml/badge.svg)
+![Tests](https://github.com/timd1/vcfdist/actions/workflows/test.yml/badge.svg)
+![Downloads](https://anaconda.org/bioconda/vcfdist/badges/downloads.svg)
 [![DOI](https://zenodo.org/badge/472945373.svg)](https://zenodo.org/badge/latestdoi/472945373)
 
 ## Overview
@@ -14,8 +15,7 @@
 ## Introduction
 vcfdist is a distance-based **germline variant calling evaluation tool** that:
 - simultaneously evaluates **SNPS, INDELs, and SVs** up to 10Kb
-- **requires local phasing** information for truth and query variants
-- can **standardize** query and truth VCF **variant representations**
+- requires **local phasing information** for truth and query variants
 - can evaluate **flip** and **switch phasing errors**
 
 vcfdist provides more accurate SNP, INDEL, and SV precision-recall curves than previous work, particularly when complex variants are involved.
@@ -33,8 +33,8 @@ Please cite the following works if you use vcfdist:
 
 <pre>
 @article{dunn2023vcfdist,
-  author={Dunn, Tim and Narayanasamy, Satish},
   title={vcfdist: Accurately benchmarking phased small variant calls in human genomes},
+  author={Dunn, Tim and Narayanasamy, Satish},
   journal={Nature Communications},
   year={2023},
   volume={14},
@@ -49,18 +49,21 @@ Please cite the following works if you use vcfdist:
 
 <details>
 <summary>
-<a href="https://doi.org/10.1101/2024.01.23.575922" target="_blank"><b>[bioRxiv]</b> Jointly benchmarking small and structural variant calls with vcfdist</a>
+<a href="https://doi.org/10.1186/s13059-024-03394-5" target="_blank"><b>[Genome Biology]</b> Jointly benchmarking small and structural variant calls with vcfdist</a>
 </summary>
 
 <pre>
 @article{dunn2024vcfdist,
-  author={Dunn, Tim and Zook, Justin M and Holt, James M and Narayanasamy, Satish},
   title={Jointly benchmarking small and structural variant calls with vcfdist},
-  journal={bioRxiv},
+  author={Dunn, Tim and Zook, Justin M and Holt, James M and Narayanasamy, Satish},
+  journal={Genome Biology},
+  volume={25},
+  number={1},
+  pages={253},
   year={2024},
-  publisher={Cold Spring Harbor Laboratory},
-  doi={10.1101/2024.01.23.575922},
-  URL={https://doi.org/10.1101/2024.01.23.575922}
+  publisher={Springer},
+  doi={10.1186/s13059-024-03394-5},
+  URL={https://doi.org/10.1186/s13059-024-03394-5}
 }
 </pre>
 </details>
@@ -85,6 +88,7 @@ sudo docker run -it timd1/vcfdist:latest vcfdist --help
 ### Option 3: GitHub source
 vcfdist is developed for Linux and its only dependencies are GCC v9.1+ and HTSlib v1.17. If you don't have HTSlib already, please set it up as follows:
 ```bash
+> sudo apt install wget bzip2 g++ zlib1g-dev libbz2-dev liblzma-dev
 > wget https://github.com/samtools/htslib/releases/download/1.17/htslib-1.17.tar.bz2
 > tar -xvf htslib-1.17.tar.bz2
 > cd htslib-1.17

analysis-v2/phab_cm/1_phab_eval_vcfs.sh

@@ -2,59 +2,60 @@
 
 source globals.sh
 
-# # truvari phab VCF normalization v4.2.1
-# source ~/software/Truvari-4.2.1/venv3.10/bin/activate
-# mkdir -p $out/phab
-# for i in "${!query_names[@]}"; do
-#     echo "phab: mafft for '${query_names[i]}'"
-#     truvari phab \
-#         -r $data/giab-tr-v4.20/GIABTR.HG002.benchmark.chr20.bed \
-#         -b $data/${truth_name}-${truth_version}/split/${truth_name}.most.vcf.gz \
-#         -c $data/${query_names[i]}-${query_versions[i]}/split/${query_names[i]}.most.vcf.gz \
-#         --bSamples HG002 \
-#         --cSamples HG002 \
-#         -f $data/refs/$ref_name \
-#         -t 64 \
-#         -o $out/phab/${query_names[i]}.vcf.gz \
-#         2> $out/phab/${query_names[i]}.log
-# done
-# deactivate
+# truvari phab VCF normalization v4.2.1
+source ~/software/Truvari-4.2.1/venv3.10/bin/activate
+mkdir -p $out/phab
+for i in "${!query_names[@]}"; do
+    echo "phab: mafft for '${query_names[i]}'"
+    truvari phab \
+        -r $data/giab-tr-v4.20/GIABTR.HG002.benchmark.chr20.bed \
+        -b $data/${truth_name}-${truth_version}/split/${truth_name}.most.vcf.gz \
+        -c $data/${query_names[i]}-${query_versions[i]}/split/${query_names[i]}.most.vcf.gz \
+        --bSamples HG002 \
+        --cSamples HG002 \
+        -f $data/refs/$ref_name \
+        -t 64 \
+        -o $out/phab/${query_names[i]}.vcf.gz \
+        2> $out/phab/${query_names[i]}.log
+done
+deactivate
 
-# # split phab normalized VCF into truth/query
-# in_samples=(
-#     "HG002"
-#     "p:HG002"
-# )
-# out_samples=(
-#     "TRUTH"
-#     "QUERY"
-# )
-# for s in "${out_samples[@]}"; do 
-#     echo $s > "${s}.txt" 
-# done
-# for i in "${!query_names[@]}"; do
-#     for s in "${!in_samples[@]}"; do
-#         bcftools view -c1 \
-#                 -s ${in_samples[s]} \
-#                 $out/phab/${query_names[i]}.vcf.gz \
-#             | bcftools reheader -s ${out_samples[s]}.txt \
-#             | sed "s/\//|/g" \
-#             | bcftools norm -a -m-any \
-#             | bcftools view -Oz \
-#             &> $out/phab/${query_names[i]}.${out_samples[s]}.vcf.gz
-#         tabix -f -p vcf $out/phab/${query_names[i]}.${out_samples[s]}.vcf.gz
-#     done
-# done
+# split phab normalized VCF into truth/query
+in_samples=(
+    "HG002"
+    "p:HG002"
+)
+out_samples=(
+    "TRUTH"
+    "QUERY"
+)
+for s in "${out_samples[@]}"; do 
+    echo $s > "${s}.txt" 
+done
+for i in "${!query_names[@]}"; do
+    for s in "${!in_samples[@]}"; do
+        bcftools view -c1 \
+                -s ${in_samples[s]} \
+                $out/phab/${query_names[i]}.vcf.gz \
+            | bcftools reheader -s ${out_samples[s]}.txt \
+            | sed "s/\//|/g" \
+            | bcftools norm -a -m-any \
+            | bcftools view -Oz \
+            &> $out/phab/${query_names[i]}.${out_samples[s]}.vcf.gz
+        tabix -f -p vcf $out/phab/${query_names[i]}.${out_samples[s]}.vcf.gz
+    done
+done
 
 # vcfdist evaluation
 mkdir -p $out/phab-vcfdist
 for i in "${!query_names[@]}"; do
     echo "phab-vcfdist: evaluating '${query_names[i]}'"
-    $timer -v /home/timdunn/vcfdist/src/vcfdist \
+    $timer -v /data-share/timdunn/vcfdist/src/vcfdist \
         $out/phab/${query_names[i]}.QUERY.vcf.gz \
         $out/phab/${query_names[i]}.TRUTH.vcf.gz \
         $data/refs/$ref_name \
         --bed $data/giab-tr-v4.20/GIABTR.HG002.benchmark.chr20.bed \
+        -t 64 -r 64 \
         -l 1000 \
         -p $out/phab-vcfdist/${query_names[i]}. \
         2> $out/phab-vcfdist/${query_names[i]}.log
@@ -75,6 +76,7 @@ for i in "${!query_names[@]}"; do
         -m ga4gh \
         --ref-overlap \
         --all-records \
+        --vcf-score-field=QUAL \
         -o $out/phab-vcfeval/${query_names[i]} \
         2> $out/phab-vcfeval/${query_names[i]}.log
     gunzip $out/phab-vcfeval/${query_names[i]}/output.vcf.gz

analysis-v2/phab_cm/2_gen_cm.py

@@ -6,6 +6,7 @@
 phab = True
 ds = "pav"
 colors = {"hprc": "Greens", "pav": "Reds", "giab-tr": "Blues"}
+names = {"t2t-q100": "Q100-dipcall", "hprc": "hifiasm-dipcall", "pav": "Q100-PAV", "giab-tr": "hifiasm-GIAB-TR"}
 if phab: phab_prefix = "phab-"
 else: phab_prefix = ""
 truvari_prefix = f"./{phab_prefix}truvari/{ds}/result_"
@@ -166,7 +167,7 @@ def match(record1, record2):
 for size_idx in range(SZ_DIMS):
     fig, ax = plt.subplots(figsize=(2,4))
     ax.matshow(np.log(counts[size_idx] + 0.1), cmap=colors[ds])
-    plt.title(f"{ds.upper()} {sizes[size_idx]} Confusion Matrix", fontsize=7)
+    plt.title(f"{names[ds]}\n{sizes[size_idx]} Confusion Matrix", fontsize=7)
     plt.ylabel("vcfdist", fontsize=7)
     ax.set_yticks(list(range(VD_DIMS-1)))
     ax.set_yticklabels(vd_cats[1:], fontsize=5)

analysis-v2/phab_cm/3_sum_cm.py

@@ -25,7 +25,7 @@
 skips = np.zeros((2,2,2), dtype=int)
 
 def record_string(record):
-    return f"{record.CHROM}:{record.POS} {record.REF} {record.ALT[0]} {record.genotype('QUERY')['GT']}"
+    return f"{record.CHROM}:{record.POS}\t{record.REF}\t{record.ALT[0]}\t{record.genotype('QUERY')['GT']}"
 
 def get_size(record):
     ref = record.REF
@@ -169,9 +169,6 @@ def hap2(var):
 
         loc, ref, alt, gt = var.split()
 
-        if vd_type and not ve_type and not tv_type:
-            print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk}")
-
         if vd_type and ve_type and tv_type: # all tools compare
             i = vd_type-1
             j = (ve_type-1)*2 + tv_type-1
@@ -181,29 +178,29 @@ def hap2(var):
                     am += 1
                     # print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
                 else:
-                    # print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
+                    # print(f"truvari TP\t{ds}\t{var}\tVD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
                     pass
             if i == 0 and j == 2:
                 # vcfeval considers two INSs at same location
-                # print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
                 pass
+                # print(f"vcfeval TP\t{ds}\t{var}\tVD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
             if i == 0 and j == 3:
                 pass
-                # print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
+                # print(f"vcfdist FP\t{ds}\t{var}\tVD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
             if i == 1 and j == 0:
+                # print(f"vcfdist TP\t{ds}\t{var}\tVD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
                 pass
-                # print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
             if i == 1 and j == 1:
                 if vd_bk == "lm":
                     pp += 1
                 elif vd_bk == "gm" and ve_bk == "am" and tv_bk == "am":
                     am_diff += 1
                 else:
-                    # print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
+                    # print(f"vcfeval FP\t{ds}\t{var}\tVD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
                     pass
             if i == 1 and j == 2:
+                print(f"truvari FP\t{ds}\t{var}\tVD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
                 pass
-                # print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
 
         elif size == SZ_SNP and vd_type and ve_type and not tv_type: # Truvari skips SNPs
             counts[size][vd_type-1][(ve_type-1)*2] += 1
@@ -215,9 +212,6 @@ def hap2(var):
             #     dels += 1
             pass
 
-        elif ve_type and tv_type and not vd_type:
-            print(f"{ds} {var} VD={cats[vd_type]}:{vd_bk} VE={cats[ve_type]}:{ve_bk} TV={cats[tv_type]}:{tv_bk}")
-
 
 sum_counts = counts[SZ_INDEL][:][:] + counts[SZ_SV][:][:]
 print("[truvari-only TP] more lenient allele match:", am)

analysis-v2/phab_cm/globals.sh

@@ -1,12 +1,12 @@
 #!/bin/bash
 
-data="/home/timdunn/vcfdist/data"
-out="/home/timdunn/vcfdist/analysis-v2/phab_cm"
-parallel="/home/timdunn/parallel-20211222/src/parallel"
-vcfdist="/home/timdunn/vcfdist/src/vcfdist"
-rtg="/home/timdunn/software/rtg-tools-3.12.1/rtg"
-tabix="/home/timdunn/software/htslib-1.16/tabix"
-bgzip="/home/timdunn/software/htslib-1.16/bgzip"
+data="/data-share/timdunn/vcfdist/data"
+out="/data-share/timdunn/vcfdist/analysis-v2/phab_cm"
+parallel="/data-share/timdunn/parallel-20211222/src/parallel"
+vcfdist="/data-share/timdunn/vcfdist/src/vcfdist"
+rtg="/data-share/timdunn/software/rtg-tools-3.12.1/rtg"
+tabix="/data-share/timdunn/software/htslib-1.16/tabix"
+bgzip="/data-share/timdunn/software/htslib-1.16/bgzip"
 timer="/usr/bin/time"
 
 # define reference FASTA