Generating the VCF input for `stoat vcf` - Pa-Tou/stoat GitHub Wiki

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 and include the path through the graph of each allele. To produce such a VCF, we recommend using the vg snakemake pipeline, specifically the genotyping pipeline, or our branch of PanGenie.

Genotyping with vg call

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]

Genotyping with PanGenie

For full instructions, see the PanGenie documentation.

Pre-processing the graph

The pangenome graph must be pre-processed before running PanGenie using this pipeline. The VCF file can be built using vg deconstruct --all-snarls -P and vcfbub -l 0 -r 100000. The gfa must contain W-lines (see wiki).

Running PanGenie

If the graph contains haploid samples, these will need to be converted to diploid. This can be done with bcftools, and sed to convert 1/1 to 1|1 as PanGenie expects.

bcftools +fixploidy -O v {input.vcf} -- -f 2 | sed 's/\([0-9*]\)\/\([0-9]*\)/\\1|\\2/g' > {output.vcf}

Then run PanGenie's indexing step:

PanGenie-index -v {input.vcf} -r {input.fasta} -o {pangenie_prefix}

And genotype:

PanGenie -f {pangenie_prefix} -i {reads} -s {sample}

As with vg call, the sample VCFs can then be merged.