This is a test version of Biostars. For the public version, visit https://www.biostars.org.
empty vcf from vg call

Hi

I am getting an empty vcf file when I run below commands

vg construct -t 10 -r /data/public-data/hg38/hg38_v2_Homo_sapiens_assembly38.fa -v /data/public-data/1000genomes/ALL.wgs.shapeit2_integrated_snvindels_v2a.GRCh38.27022019.sites_chr_1:Y.vcf >1000g.vg

vg prune -k 12 -e 2 -t 16 -p 1000g.vg >1000g.pruned.vg

#### store the graph in the xg/gcsa index pair
vg index -t 10 -p -b /data/users/philge/tmp/ -x 1000g.xg -g 1000g.gcsa -k 16 1000g.pruned.vg

####align a read to the indexed version of the graph
vg map -t 10 -b /data/users/philge/genomegraph/vg/dbsnp/sorted_A019.bam -x 1000g.xg -g 1000g.gcsa > ./alignment/sorted_A019_1000g.gam

####Compute the read support from the gam (with ignoring qualitiy < 5)
vg pack -t 10 -x 1000g.xg -g ./alignment/sorted_A019_1000g.gam -Q 5 -o ./alignment/sorted_A019_1000g.pack

####Generate a VCF from the support
vg call -t 10 -r ./1000g.snarls ./1000g.xg -k ./alignment/sorted_A019_1000g.pack > ./alignment/sorted_A019_1000g_no_aug.vcf

All my input files for vg call contain data. But output vcf is empty. Any suggestions on what could be the reason?

Thanks
Philge

vg
• 2,112 views • link
updated by kingcohn • kingcohn • written by philgev2 • philgev2

any suggestions, please? I am stuck with the analysis.

• 0 views • link
written by philgev2 • philgev2

2 answers

The commands look okay (though vg snarls seems to be missing). Assuming your snarls file is fine, the other reason for an empty VCF could be not enough mappings found in the pack file. You can use vg depth -k -b 10000 to verify the coverage in your pack file (and with -g to verify the coverage in your GAM).

• 0 views • link
written by glenn.hickey glenn.hickey

I have a similar issue, did you resolve this?

• 0 views • link
written by kingcohn • kingcohn

Log in to answer this question.

More posts like this