Clarification on VCF Output

[HLHsieh]

Hi there,

I noticed differences in the VCF outputs between v1.5 and v1.5.3 when processing the same locus on the same sample, and I’m looking for some clarity on how to interpret the REF and ALT fields in v1.5.3. In v1.5, the REF and ALT fields seem to provide the full repeat sequences. However, in v1.5.3 the fields have changed and I am wondering whether I can get the similar sequences output in v1.5.3. I would prefer to continue using v1.5.3, since v1.5 occasionally throws an IndexError: list index out of range.

Additionally, I observed substantial differences in the reported results between v1.5 and v1.5.3. Could these discrepancies be due to the updated GMM clustering post-processing introduced in v1.5.2?

In v1.5.3,

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	.
chr1	1435798	.	T		.	PASS	RUL_REF=10	GT:DP:AD	0:15:13

In v1.5,

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	.
chr1	1435798	.	GGCGCGGAGCGGCGCGGAGC	GGCGCTCGGCCTGGCCCTGAGCTTGCGGCTGGCGCTGGCGCGGAGCGGCGCGAGCGCGGTGAGTGCGGCGGGCGGCCGGGCCGGG	.	PASS	END=1435818;MOTIF=GGCGCGGAGC;COPIES=2.0	GT:DP:AL:ALR:AC:ACR:AD:ALT_MOTIF	1/1:15:17.9/17.9:2-48/2-48:3.0/3.0:0.3-8.0/0.3-8.0:14/14:GGCGCT(10);GGCCG(1);GCGGCGG(1)

Any suggestions would be appreciated.

Best,
Hsin

[readmanchiu]

Hi @HLHsieh
Thanks for using Straglr.
The VCF format updated is trying to following the VCFv4.5 (to the best of my interpretation)
Having said that, the ALT sequence(s) should be output. But your example seems to be a homozygous loci with reference-sized alleles, therefore the ALT sequences are not output in the VCF
I can probably add a custom field to output the reference repeat sequence if desired.
From 1.5 to 1.5.3, the major change is the clustering. If you look at the copy_number and size columns in the tsv output, they should probably be the same, but the genotype and allele columns may be different.
Let me know if there is something you think is not correct since the version change. I will have a look.

[davidlougheed]

Hello,

I have a similar question with Straglr. I am trying to benchmark various STR callers with Truvari using the GIAB TR benchmark, but Straglr outputs only symbolic alleles (<CNV:TR>). With symbolic alleles, there are some metrics we do not get (e.g., TruScore) from Truvari.

I also get an error when trying to parse the output VCF from this version of Straglr; since AD has 2 values but the AD FORMAT is set to R ("The field has one value for each possible allele, including the reference" - from spec).

[readmanchiu]

I agree that the symbolic alleles specified in vcf4.5 are quite lacking in info, but all the relevant information should be reported in the INFO fields. We don't report any score in Straglr, you won't be able to find in any of the outputs.
The AD has should report a comma-separated text, as it does in the vcf output of the test data. Please post the output line where AD is set to R for me to have a look.

[davidlougheed]

the line in question I'm referring to is

##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Total read depth for each allele">

I believe R is not the right number here, as I think AD in Straglr is per-genotype not per-allele.

[readmanchiu]

sorry @davidlougheed, not sure why you think AD is per genotype, as it reports the read depth of each allele as stated in the description. Can you give an example?

[davidlougheed]

Here's an example:

chr1    29745   .       A       <CNV:TR>,<CNV:TR>       .       PASS    RUS_REF=CCCACG;SVLEN=25,25;RN=1,1;RUS=CCCCG,CCCCG;RUC=2.6,6.4;CIRUC=0.0,0.0,0.0,0.0     GT:DP:AD        1/2:18:6,12

note that in the spec it says "The field has one value for each possible allele, including the reference" - which would mean AD would have three entries here with Number=R. the tricky part here is that because SVs are kind of hacked onto the VCF spec, and the VCF spec is not very precise, I believe this means the read support for each possible value (as in, reads that support the reference vs reads that support each alt.)

[readmanchiu]

Thanks @davidlougheed
The intention of the "AD" fields here are exactly meant for read support of each detected allele. Originally I thought using AD will save creating extra custom fields. But it seems like AD has pretty strict definition so I will probably create a custom field to avoid the contention.

[davidlougheed]

yeah, I've ran into this before as well!