GT are all ./. in SVDSS VCF

[yangao07]

Hi, I ran SVDSS-v2.1.0 on HiFi BAMs with the following commands, the output VCFs all look good except the GT fields are all ./., is this intended?
Thanks!

   clu_weight=2
   min_sv_len=30
   SVDSS smooth --reference $ref --threads $threads --bam $bam > $smooth_bam
   samtools index -@ $threads $smooth_bam
   SVDSS search --threads $threads --index $ref_idx --bam $smooth_bam > $out_sfs
   SVDSS call --reference $ref --bam  $smooth_bam --sfs $out_sfs \
               --min-sv-length $min_sv_len --min-cluster-weight ${clu_weight} --threads $threads > $out_vcf

[ldenti]

Hi, yes that's intended. Currently we just call variations and we rely on external genotyper for getting the GT information.

In v2.1.0 we suggest to use kanpig. Adding it to your script should be straightforward. Here the command we use:

SVDSS/run_svdss

Line 197 in 0787a23

/usr/bin/time -vo ${wd}/genotyping.time ${kanpig} gt --input $wd/variations.vcf.gz --reads $BAM --reference $FA | bcftools sort -Oz >$wd/variations-gt.vcf.gz

Let me know if you need more info or you have any other question.

[yangao07]

I see. Thanks!
I should use the smooth.bam instead of raw bam for kanpig, right?

[ldenti]

Never thought about it. I always used the original BAM file since I was assuming that removing SNPs and small indels from the alignments can impact the genotyping (kanpig should expect them in the alignments).

However, from a small test I did on some example data I have, I got 3 SVs (out of 733) genotyped differently using the two BAMs. By inspecting the three loci in IGV, the genotypes computed from the smoothed BAM seems the correct one.

Since kanpig should be fast, if you have time to test both BAMs and provide me the genotyped VCFs, I can check this more. (in this case, I could also need the reads - if they are publicly available, I could also run everything by myself though)

[yangao07]

I am running it on both bam files.
Not finished on all datasets yet, but based on truvari, the two vcf files are very similar.
Will update the final result here later.

[yangao07]

Hi, just want to update that: across different coverage 15~40X, and on both hifi and ONT data, the precision/recall of raw and smooth.bam are all very similar.

One side quesion: I noticed clear drop of precision on ONT data when coverage increase, e.g., 15 to 40X.
Should I use different parameters for different coverage? Like --min-cluster-weight?

[yangao07]

Here are the truvari result for 15X:

cat ../../15X/svdss/svdss.Truvari_v5.3/refine.variant_summary.json
{
    "TP-base": 24502,
    "TP-comp": 23314,
    "FP": 2215,
    "FN": 3686,
    "precision": 0.9132359277684202,
    "recall": 0.8692351355186604,
    "f1": 0.890692445767286,
    "base cnt": 28188,
    "comp cnt": 25529
}

and 40X:

cat svdss.Truvari_v5.3/refine.variant_summary.json
{
    "TP-base": 24784,
    "TP-comp": 23834,
    "FP": 5592,
    "FN": 3404,
    "precision": 0.8099639774349215,
    "recall": 0.8792393926493544,
    "f1": 0.8431811683541761,
    "base cnt": 28188,
    "comp cnt": 29426
}

Looks like there are much more SVs called for 40X.

[ldenti]

Keep in mind that SVDSS is not optimized for ONT. We still need time to test it and optimize it (there are some heuristics that we believe do not fully work with higher error rates - what chemistry are you using?).

So you can use it but you might get lower accuracy. Just for reference, can you also share the accuracy on the HiFi data you mentioned?

[yangao07]

I am using ONT R10 datasets: HG002_PAW70337, 2025-01, probably the latest datasets.

For HiFi 15X:

cat 15X/svdss/svdss.Truvari_v5.3/refine.variant_summary.json
{
    "TP-base": 23291,
    "TP-comp": 21436,
    "FP": 928,
    "FN": 4897,
    "precision": 0.9585047397603291,
    "recall": 0.8262735915992621,
    "f1": 0.887490776833189,
    "base cnt": 28188,
    "comp cnt": 22364
}

and HiFi 40X:

cat 40X/svdss/svdss.Truvari_v5.3/refine.variant_summary.json
{
    "TP-base": 24817,
    "TP-comp": 23085,
    "FP": 1337,
    "FN": 3371,
    "precision": 0.9452542789288346,
    "recall": 0.8804101035901802,
    "f1": 0.9116806195041478,
    "base cnt": 28188,
    "comp cnt": 24422
}

I am using HG002 T2T v1.1 benchmark set.