This is a test version of Biostars. For the public version, visit https://www.biostars.org.
best practice for diploid variant calling

Hi I'm working on an experiment to benchmark variant calling using both a graph and a reference-based approach. I got the graph part working smoothly, but the reference approach is giving me some trouble... In particular, I need to evaluate calls for both haplotypes of the same individuals and I'm testing three strategies:

  1. merge hap1 and hap2 in a single FASTA with univocal headers, then align
  2. align independently to each haplotype and merge the BAMs (see below for details), then call
  3. TBD – same as point 2. but merge the VCF generated after calling from each BAM independently

I'm not sure which is the more “correct” way to face this problem considering there will always be a compromise to be made. Specifically,

point 1

Everything works just fine down to variant calling but then I'm having problem to set the benchmark with hap.py since the tool requires a unique reference – in my case I picked hap1 as it is the better one – and in order to fix chromosome headers and duplicate calls in the output VCF after variant calling I've done the following:

zcat <output>.vcf.gz | awk -F '\t' '{ OFS="\t"; gsub("_hap[0-9]", "", $1); print }' | bcftools sort | bcftools norm -D -Oz -o <output>-bench.vcf.gz && tabix -p vcf <output>-bench.vcf.gz

point 2

I order to merge calls from BAMs aligned to different references (haplotype of the same individual in this case) I used bam-mergeRef since I'm not sure of a way to do so with samtools merge; the problem with this is that the variant calling step yield an empty VCF...

This is probably the more parsimonious way to do this I guess, so if anyone has input/experience about a possible workaround let me know!

point 3

I haven't thoroughly tested this yet, but my idea is to independently call BAM files aligned to each haplotype, produce the corresponding VCF files and then merge calls; I set up a command to do so but I have never done this before, so also here I'm open to any input/feedback:

bcftools merge --merge all <output>-hap1.vcf.gz <output>-hap2.vcf.gz -Oz -o <output>-merged.vcf.gz

First of all, I'll be mostly interested in knowing which one is the more unbiased way to address this issue; secondly, I'm aware point 1 would require a bit of engineering on the hap.py side to understand what is causing an error as well as that is probably not advised but if someone has done this before any help is appreciated!

Sorry for the long post, I just want to finish saying I'm working with HiFi long reads and aligning them with minimap2, then calling variants with deepvariant.

variant-calling haplotype samtools diploid bcftools
• 6,219 views • link
written by Matteo Ungaro • Matteo Ungaro

It sounds like you are calling variants with respect to a "phased assembly" (e.g. a diploid assembly has both a maternal and paternal haploid genome assembly). are you aligning the sample against itself e.g. the reads used to create the assembly are aligned against the assembly itself?

• 0 views • link
written by cmdcolin cmdcolin

@cmdcolin indeed, you're correct! We do have a parental cell line which was used to build the two phased assemblies with Verkko, then this same cell line has been exposed to different conditions and de novo sequenced.

By aligning the reads used to assemble the two haplotypes as well as those exposed to different experimental conditions, we are looking to better understand mutational patterns taking advantage of the ground truth in our two assemblies.

• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

that is quite a cool experiment. I am not actually sure of the best approach (sorry to disappoint). I think that if you look up "phased assembly variant calling" there are some discussions like in Variant calling analysis on phased assembly

I think that aligning to a graph that is composed of both phased assemblies is, at least theoretically, pretty ideal. the other solutions like you describe above seem to be a bit more complex...you have to consider whether the aligner would understand having extremely similar sequences if you align reads to the concatenation of both the maternal and paternal haplotypes....or you would have to consider the penalties involved for the aligner to align to both the maternal and paternal haplotype separately, which would implicitly align e.g. maternal reads to the paternal haplotype, and vice versa. So, tricky problem!

• 0 views • link
written by cmdcolin cmdcolin

2 answers

Cool experiment! Aligning HiFi reads back to phased Verkko assemblies for mutation benchmarking sounds solid.

For the reference-based approach, I think point 3 (calling variants separately on each haplotype then merging VCFs) is the least biased. It avoids multi-mapping issues from concatenating FASTAs (point 1) and doesn't force BAM merging across differing refs (point 2). Use bcftools merge --merge all, but first normalize contig names (e.g., remove _hap suffixes) to align chromosomes. Watch for overlapping positions—run bcftools norm -m+any afterward to handle multi-allelic sites.

If haplotypes have structural differences, consider lift-over tools like CrossMap for better alignment before merging. This keeps things phased-aware without projecting everything to one hap.

For hap.py benchmarking, build a diploid truth set by combining phased variants from your assemblies, or use one hap as ref but flag potential biases in results.

Graph methods like PGGB (from the pangenome toolkit) are theoretically best for diploid handling—glad that's working smoothly. Check PacBio's HiFi variant calling guide for DeepVariant tweaks: align with pbmm2, phase with WhatsHap post-calling.

No perfect compromise, but point 3 minimizes reference bias. If empty VCFs persist in point 2, debug the bam-mergeRef step—might need a combined ref FASTA for DeepVariant.

Kevin

• 0 views • link
written by Kevin Blighe Kevin Blighe

Kevin thanks I was experimenting with the following since I already pursued most of the approaches you mentioned and, in practice, they are either leading to corrupted outputs or generating some form of bias.

Anyway, I attempted something similar to this

Use bcftools merge --merge all, but first normalize contig names (e.g., remove _hap suffixes) to align chromosomes. Watch for overlapping positions—run bcftools norm -m+any afterward to handle multi-allelic sites.

but instead of -m+any I used -m-any. For some reason, I'm now incurring in the following error prior to the step of bcftools merge:

The REF prefixes differ: CCCAA vs A (5,1) Failed to merge alleles at chr1:3870 in clone0_prt-output-split-hap2.vcf.gz

since the hap1 has been aligned and called against the haplotype 1 of this genome, which is causing the problem above when attempting to merge while using different references.

Is there a way to fix it? I was looking into bcftools norm -f [FASTA1] [VCF1] > [VCF1_norm] and bcftools norm -f [FASTA2] [VCF2] > [VCF2_norm] to left-align indels before merging, but the error persists...

Another option would be to call variants (rerunning DV) for haplotype 2 using haplotype 1 as reference – if DV doesn't complain about potential discrepancies between the BAM and the reference – to get, ultimately, the same set of coordinates before merging.

Let me know, thanks again for the constructive advices!

• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

I found this small tool which aligns to both assemblies, and chooses the best: https://github.com/jheinz27/diploidinator (used by https://github.com/jheinz27/breakinator)

it does not solve all your issues but it is a step towards aligning to a diploid assembly

• 0 views • link
written by cmdcolin cmdcolin

@cmdcolin amazing suggestion! I run diploidinator then attempted to convert PAF to BAM with rustybam; however, shortly after I hit a wall... unfortunately, but obviously, I couldn't merge BAM from with two different coordinate systems since they have duplicated headers.

It appears diploidinator has a way to handle this, and for variant calling I'm in a need of a final BAM, but when attempting conversion to SAM (then piped into BAM with samtools) using rustybam that's when I have issues with duplicated headers...

• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

indeed, you need to treat the two haplotypes as completely independent and there is no sensible to "merge" the results as they will have very different coordinate systems. it is not just the "duplicate headers" the two haplotypes have completely different coordinate systems, different chromosome lengths, etc. the only way to sensibly compare the results is to have an alignment between the two assemblies (e.g. directly align then with minimap hap1.fa hap2.fa > mapping.paf or alternatively map both back to the grch38 reference minimap2 hg38.fa hap1.fa >hap1_vs_hg38.paf and minimap2 hg38.fa hap2 > hap2_vs_hg38.paf) and then you can liftover results back onto the reference genome after hopefully getting more 'accurate' mapping on the actual haplotype genomes. alternative, can go full graph based alignment, by putting hg38, hap1, and hap2 into a graph, and then align and variant call direct on the graph

• 0 views • link
written by cmdcolin cmdcolin

@cmdcolin much appreciated once more! I already pursued the graph-based approach successfully; in fact, the whole point of this analysis is to benchmark a personalized pangenome graph versus T2T references! I guess I will need to do a liftover somehow at some point.

What are your suggestions:

  1. should I lift the PAF from hap1 aligned on hap2 (or vice-versa, is there a difference in this case)
  2. should I lift the VCF to a common reference and then merge

I prefer not to follow with an “external” reference but keep it to the haplotype space of this sample (its T2T genome) to have as little deviation as possible when assessing its specific genetic variants.

P. S. in case I go for point 1. after I have the PAF alignment what can I do in practice? I used CrossMap to liftover in the past — which seems to be working on BAM as well. The end result is expected to be a single BAM file I can use for variant calling. The problem with this tool is that it needs a chain.gz to lift coordinates...

• 0 views • link
written by Matteo Ungaro • Matteo Ungaro

Log in to answer this question.

More posts like this