Understanding the HG002 truthset

[KasperH2]

Hi
Sorry for posting here, as this is mostly related to the HG002 truthset, but there doesn't seem to be much activity (and no discussion section) on the genome-data-integration github and it is also somewhat related to Truvari. Hope you can help anyway.

I have been using Truvari with HG002_SVs_Tier1_v0.6.vcf.gz and HG002_SVs_Tier1_v0.6.bed,
However I see some different numbers in the article for the HG002 benchmark
There is 12,745 SVs in the truthset with the filter tag PASS.

”These regions define our Tier 1 benchmark set, which spans 2.51 Gbp and includes 5,262 insertions and 4,095 deletions. These regions exclude 1,837 of the 12,745 SVs because they were within 50 bp of a 20–49-bp indel; they exclude an additional 856 SVs within 50 bp of a candidate SV for which no consensus genotype could be determined; and they exclude an additional 411 calls that were not fully supported by a diploid assembly as the only SV in the region.”

Now, 5262 + 4095 = 9357 variants
But 12745 - (1837 + 856 + 411) = 9641 variants

Not sure if I am missing something and this is intended?
Anyway, using Truvari with --passonly with the HG002_SVs_Tier1_v0.6.bed gives 9641 variants, so I assume this is the correct number.

Also, although only related to the HG002_SVs_Tier1_v0.6.vcf.gz truthset, I have looked at the number of different REPTYPEs contained under the INS and DEL SVTYPEs. Does this table make sense? Can an INS have a SIMPLEDEL as a REPTYPE and a DEL have DUP as REPTYPE for example (these are including all variants):

37412	DEL		36600	INS
1	ARTIFACT		1	ARTIFACT
11867	CONTRAC		90	CONTRAC
12	DUP		27297	DUP
24531	SIMPLEDEL		234	SIMPLEDEL
			7639	SIMPLINS
999	SUBSDEL		78	SUBSDEL
2	SUBSINS		1261	SUBSINS

Just asking, because I don't fully know. But to me it seems odd?

Thanks,
Kasper Holm

[ACEnglish]

Hello,

I'll try my best to answer your questions with the disclaimer that I'm not the official maintainer of the GIAB data.

It appears that the DEL/INS SV counts in the manuscript are incorrect. 9641 is the count of PASS/Tier1 SVs which breaks down to

   4199 DEL
   5442 INS

From my understanding, the REPTYPE isn't a technical definition the same as SVTYPE's DEL/INS as much as it is an annotation. It would appear you have found edge cases where the REPTYPE annotation could be incorrect. However, it may be useful for you to explore how many of these are Tier1 SVs. Generally, I choose not to stray outside of the PASS/Tier1 regions.

$ bcftools view  -i "REPTYPE == 'DUP' & SVTYPE == 'DEL'" Tier1_Only.vcf.gz | grep -v "#" | cut -f1-6
19	10636680	HG2_Ill_svaba_21979	TATATACATATACATATACATATACATATACATATACACATATACATATAC	T	10

So for this example it's a 50bp deletion that looks to be in some kind of simple repeat that may be the cause of the missed annotation.
Furthermore, looking for this SV in the Truvari Raw Data generated as part of our preprint shows that high quality long-read assemblies also discover the SV.

[KasperH2]

Thanks for taking your time to look into this. That clarifies a few things. Yes I assume the algorithms doing the REPTYPEs can be fooled by the type of sequences deleted/inserted. For your information this is what the table for the Tier1 SVs look like:

<style> </style>

5464	del		7281	ins
1738	CONTRAC		19	CONTRAC
1	DUP		2818	DUP
3442	SIMPLEDEL		87	SIMPLEDEL
			3942	SIMPLINS
283	SUBSDEL		34	SUBSDEL
			381	SUBSINS

Used:
cat $VCFFILE | grep -v "^#" | grep "SVTYPE=DEL" | perl -ne 'print "$1\n" if /[;\t]REPTYPE=([^;\t]+)/' | sort | uniq -c
cat $VCFFILE | grep -v "^#" | grep "SVTYPE=INS" | perl -ne 'print "$1\n" if /[;\t]REPTYPE=([^;\t]+)/' | sort | uniq -c

[ACEnglish]

Our slight differences may be due to the region files. I believe that bcftools -R doesn't exactly fit how Truvari defines an SV fitting into an --includebed region.

$ bcftools query -R HG002_SVs_Tier1_v0.6.bed -i "FILTER == 'PASS'" -f "%INFO/SVTYPE\t%INFO/REPTYPE\n" HG002_SVs_Tier1_v0.6.vcf.gz | sort | uniq -c
   1176 DEL	CONTRAC
      1 DEL	DUP
   2870 DEL	SIMPLEDEL
    214 DEL	SUBSDEL
     10 INS	CONTRAC
   1771 INS	DUP
     63 INS	SIMPLEDEL
   3321 INS	SIMPLEINS
     28 INS	SUBSDEL
    251 INS	SUBSINS

[KasperH2]

The difference is because I did not intersect the Tier1s with the bed file, but used all those that are listed as PASS in the HG002_Tier1_SVs_v0.6.vcf (12745 SVs). Intersecting it with bedtools intersect gives 9719 SVs and using the truvari --bed way gives the 9641 SVs that is also used in the paper. The difference, as I understand it, is that bedtools intersect does not require both start and stop intervals to be within one bed region, whereas Truvari does. You have 9705 with bcftools query.

On another note, I tried to do a workaround to test for performance of DUP calling, as I see some cnv callers can find them.
The problem is that the cnv callers will say SVTYPE=DUP, while the truthset says SVTYPE=INS,REPTYPE=DUP, which will be seen as a mismatch by truvari?
As a solution I changed the truthset so that if a variant had REPTYPE=DUP I would set SVTYPE=DUP in it. Got no matches, but I think Truvari can handle this solution? It will compare the SVTYPE=DUP, right?

[ACEnglish]

See developer documentation for how svtype is determined.