Hello,
I am using the vg Giraffe haplotype-sampling workflow to genotype SVs in a minigraph cactus pangenome graph from shortread sequencing data, with the goal of estimating allele frequencies for the fixed catalogue of SVs (from the pangenome graph) across many samples. I understand the Giraffe mapping/haplotype-sampling portion of the workflow, but I am less certain about the recommended use of vg pack and vg call afterward.
Previously, I mapped reads against the full GBZ and used vg pack, then vg call full.gbz -k sample.pack -a -A -z ...
My understanding of these options is:
-a: genotype every snarl, including reference calls
-A: genotype all snarls, including nested child snarls, I needed this to run vcfbub later
-z: restrict genotypes to traversals represented in the GBZ
This was convenient because I could get a genotype for essentially every snarl and then merge calls across samples to calculate allele frequencies. However, I found that this approach produced incorrect genotypes for some SVs for which I already knew the genotype independently.
I then switched to Giraffe haplotype sampling. I generated a sample specific GBZ, mapped the reads to that graph, ran vg pack against the sampled graph, and then ran vg call -z against the sampled graph. For the previously problematic variants, this produced the correct genotypes, including for some relatively complex SVs that other short read genotyping approaches had difficulty with.
The issue is that the sampled GBZ is a subgraph of the full graph, so it does not necessarily contain every snarl/allele in my original pangenome SV catalogue. Consequently, the VCF generated from the sampled graph does not contain a genotype for every catalogue variant.
I also tried using the gam produced by mapping to the sampled graph to run vg pack against the full graph, followed by vg call -a against the full graph. Since the sampled graph is a subgraph of the original graph, I thought this might allow me to recover genotypes for the complete catalogue. However, I still do not obtain calls for every snarl in the full graph. My main questions are:
What is the recommended workflow for using haplotype sampling for mapping while still obtaining genotypes for a fixed catalogue of snarls/SVs from the original full graph?
Is it appropriate to map to the sampled GBZ and then run vg pack/vg call against the original full graph, given that the sampled graph is a subgraph and the alignments should therefore be valid on the original graph? If so, should vg call -a against the full graph theoretically emit a genotype for every catalogue snarl?
More generally, is there a recommended way to distinguish:
-a variant that was excluded from the sampled graph because haplotype sampling considered it unlikely to be present -a variant that is represented but has evidence for the reference allele
- and a variant for which there is not enough read evidence to make a genotype call?
For population allele-frequency estimation, I want to avoid incorrectly treating “not emitted by vg call” as equivalent to 0/0 if the absence of the record instead reflects the sampled graph or insufficient genotype evidence.
Thank you!
0 answers
No answers yet.
Log in to answer this question.