-r parameter is not taken when "AC=" is presented in vcf INFO

[oushujun]

Hello @jporton,

I used -r to specify using the REF as the ancestral state, and found SoFoS was not taking this parameter for some occasions. This seems to be related to the INFO column of the vcf file where there is AC (allele count) information. When removing it or replacing the AC information, it takes the -r parameter again.

This test.vcf.zip can reproduce the issue:

sofos -n 10 -u -r -a 1.0 -b 1.0 test.vcf 
#SoFoS v2.0
#date=2021-12-18T16:33:38-0600
#epoch=1639866818307
#path=test.vcf
#alpha=1
#beta=1
#size=10
#folded=0
#refalt=1
#
## [0.0006s elapsed] 60 sites processed
Number,Prior,Observed,Posterior
0,5.4545454545454959,0,0.69853641197197858
1,5.4545454545454959,1,1.3048755954406486
2,5.4545454545454959,1,1.8488491162524154
3,5.4545454545454959,1,2.4702060322842092
4,5.4545454545454959,0,3.2635682519285871
5,5.4545454545454959,1,4.2688676722197485
6,5.4545454545454959,0,5.4985717806044025
7,5.4545454545454959,0,6.9696393737318694
8,5.4545454545454959,0,8.7375634335076597
9,5.4545454545454959,0,10.95262220018968
10,5.4545454545454959,56,13.986700131869155

sed 's/AC=/AA=/' test.vcf | sofos -n 10 -u -r -a 1.0 -b 1.0 - 
#SoFoS v2.0
#date=2021-12-18T16:34:02-0600
#epoch=1639866842418
#path=-
#alpha=1
#beta=1
#size=10
#folded=0
#refalt=1
#
[W::vcf_parse] INFO 'AA' is not defined in the header, assuming Type=String
## [0.0006s elapsed] 60 sites processed
Number,Prior,Observed,Posterior
0,5.4545454545455394,0,16.557384357674806
1,5.4545454545455394,55,18.025665060544522
2,5.4545454545455394,2,12.323574756722911
3,5.4545454545455394,1,6.569310900606677
4,5.4545454545455394,0,2.9275361416968226
5,5.4545454545455394,0,1.1301463637682869
6,5.4545454545455394,0,0.41433400150816324
7,5.4545454545455394,0,0.22850857449687428
8,5.4545454545455394,0,0.30356657899338751
9,5.4545454545455394,1,0.54636205470466404
10,5.4545454545455394,1,0.97361120928371347

Any ideas? Thanks!

Shujun

[jporton]

It appears that your AN and AF tags are off. Your AN should be 20 not 2, and your AF is 1.0 saying the frequency of the alternative is 1.
It appears that in your file AN, AC, and AF are being calculated in your vcf as if allele of 0 is a missing allele.

[oushujun]

@jporton Thank you for helping to check the vcf file. Did you mean AN, AC, AF values will be used and not recalculated even they are misspecified and disregard the -r parameter?