Limitation on the number of alternative alleles ?

[Wenfei-Xian]

Dear VG team

I want to use VG to genotype the copy number of the short tandem repeats (ATCATCATC) based on short reads. In the alt field, each allele represents one specific copy number, such as: ATC,ATCATC,ATCATCATC

Here is the input vcf to construct the graph:

fileformat=VCFv4.1

fileDate=20141110

source=mutatrix population genome simulator

seed=1415643582

reference=NA

phasing=true

commandline=mutatrix --dry-run -s 0.05 -i 0.01 -p 2 x.fa

INFO=

INFO=

INFO=

INFO=

INFO=

INFO=

FORMAT=

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT zz

random_ATC_4 5001 random_ATC_4_5001 ATCATCATCATC ATC,ATCATC,ATCATCATC,ATCATCATCATCATC,ATCATCATCATCATCATC,ATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC,ATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATCATC . . . . .

Here is the result. ATC(1) which mean the simulated reads carry one copy of ATC. GT = VG (call) output file. ATC(2) which mean the simulated reads carry two copies of ATC. GT = VG (call) output file. ATC(3) which mean the simulated reads carry three copies of ATC. GT = VG (call) output file.

It performs well on the low number of copy number. However, two things confuse me.

1, reference allele. Such as ATC(4). in the input vcf file, ref field is ATCATCATCATC 4 copies. when I map short reads which carry 4 copies of ATC to the graph, ideal GT should be 0/0, but GT in VG call is 4/7.

2, When the short reads which carry more than 17 copies, the GT always 16/16. Does it mean that the maximum number of alternative alleles is 16 ?

vg construct -S -a -r reference.fa -v copynumber.vcf.gz > copynumber.vg vg autoindex --workflow giraffe -r reference.fa -v copynumber.vcf.gz -p reference.fa vg snarls -a copynumber.vg > copynumber.snarls vg giraffe -Z reference.fa.giraffe.gbz -m reference.fa.min -d reference.fa.dist -f short_reads.ATC10copy.1.fq -f short_reads.ATC10copy.2.fq > short_reads.ATC10copy.gam vg pack -x copynumber.vg -g short_reads.ATC10copy.gam -o short_reads.ATC10copy.gam.pack vg call copynumber.vg -k short_reads.ATC10copy.gam.pack -r copynumber.snarls -v copynumber.vcf.gz -s zz > short_reads.ATC10copy.genotype.vcf

Many thanks

Best, Wenfei

[glennhickey]

What is the coverage and length and error model of your reads? I'd guess:

1. (wrong genotypes): this could occur if the coverage is fairly low and a couple of reads somehow match the higher repeat by accident
2. (max genotype 16/16): this could be due to your reads being too short for Giraffe to map across the higher-order alleles. I guess it's also possible that it's a funciton of Giraffe's seeding paramter values.

[Wenfei-Xian]

Thanks for your reply !

reads were simulated by wgsim. Coverage about 100X. reads length = 100bp

When the number of alterative path less than 16 , everything go well ... (no wrong genotype)

Which paramter should I add to avoid the limitation (max genotype 16/16) ?

[Wenfei-Xian]

The problem solved ! When I use manual index instread of autoindex.

vg construct -r random.$1.$2.fa -t 20 -v random.$1.$2.fa.bed.vcf.gz -a > random.$1.$2.fa.vg vg index -x random.$1.$2.fa.xg -g random.$1.$2.fa.vg.gcsa -k 16 -t 100 random.$1.$2.fa.vg

vg gbwt -p -g random.$1.$2.fa.gg -o random.$1.$2.fa.gbwt -v random.$1.$2.fa.bed.vcf.gz -x random.$1.$2.fa.xg vg gbwt -p -g random.$1.$2.fa.sampled.gg -o random.$1.$2.fa.sampled.gbwt -x random.$1.$2.fa.xg -n 64 -l random.$1.$2.fa.gbwt vg minimizer -t 20 -p -i random.$1.$2.fa.min -g random.$1.$2.fa.sampled.gbwt -G random.$1.$2.fa.sampled.gg vg snarls random.$1.$2.fa.xg >random.$1.$2.fa.snarls vg index -j random.$1.$2.fa.dist -s random.$1.$2.fa.snarls random.$1.$2.fa.vg