Generating the VCF input for `stoat vcf`

stoat vcf takes as input a VCF file representing the variants present in a set of samples. These variants must be represented relative to the pangenome graph used by stoat vcf. To produce such a VCF, we recommend using the vg snakemake pipeline, specifically the genotyping pipeline.

Briefly, sequencing reads are mapped with vg giraffe, the coverage is computed with vg pack, and genotypes are predicted with vg call. Minimal commands would look like:

vg giraffe -Z graph.gbz -t {threads} -f sampleXXX.R1.fq.gz -f sampleXXX.R2.fq.gz -o gaf | gzip > sampleXXX.reads.gaf.gz
vg pack -x graph.gbz -g sampleXXX.reads.gaf.gz -Q 5 -t {threads} -o sampleXXX.cov.pack
vg call -t {threads} -k sampleXXX.cov.pack -aA -s sampleXXX -z graph.gbz | bgzip > sampleXXX.vgcall.vcf.gz

The single-sample VCFs are then merged into one VCF (the input for STOAT) using bcftools, see below.

Running the vg snakemake

The easiest way to run this pipeline is to clone the repo and write a new config file for your experiment.

Write a config file

Here is an example config.yaml file with parameters needed by stoat vcf. Other parameters may be necessary depending on your data and computing environment.

graph: 'graph_name'
gbz: "path_to_graph/graph_name.gbz"
dist: "path_to_graph/graph_name.dist"
ref_paths_list: 'reference_path_list.tsv'       # The paths in the graph to be used as references, one per line. Can be found with `vg paths -L`
gt_ref: True                                    # This ensures that the output variants match between samples
sample_tsv: 'sample_list.tsv'                   # A tsv of sample name and read file location(s). see https://github.com/vgteam/vg_snakemake?tab=readme-ov-file#input-fastqcram
gt_min_var_len: 0                               # Genotype all sizes of variants, including SNPs

Run the pipeline

The pipeline can be run as described here, or on an HPC as described here.

Merge the VCFs

The VCFs across several samples are then merged using bcftools merge, merging multi-allelic records by ID (-m id). For example, something like:

bcftools merge -m id sample1.vgcall.vcf.gz sample2.vgcall.vcf.gz sample3.vgcall.vcf.gz -Oz -o merged.vgcall.vcf.gz

If you used the Snakemake pipeline described above, the individual VCF files for each sample should be in path_to_vg_snakemake_dir/results/[sample_name]/[sample_name].[graph_name].gt.minlen0.vcf.gz. To merge them, using the sample list from the config file to get all sample names, you could do:

bcftools merge --file-list <(tail -n +2 [sample_list] | awk  -v OFS='' '{print "path_to_vg_snakemake_dir/results/", $1, "/", $1, ".[graph_name].gt.minlen0.vcf.gz"}') -m id -Oz -o merged_vcf.vcf.gz

Note: before merging, each sample VCF should be indexed with bcftools index -f [sample_vcf]