Vg call in a specific region

Question

[arunkumarramrsh]

Hi,

I am trying to call variants in a narrow region using vg chunk and vg call. It worked with an older version of vg tools but does not seem to work with the newer one.

In a previous version, this worked for me

vg chunk -x ref_alt.xg -a SRR5647740_paired_wasp_mapped_sorted.gam -g -p 3:3716672-3718990 -c 10 > chunk.vg

for file in chunk_0_3_3716352_3719310.gam; do vg augment -a pileup -t 20 -Z ${file/.gam/_aug.trans} -S ${file/.gam/_aug.support}  chunk.vg $file  > ${file/.gam/_aug.vg} 2>err_vg_aug ; done 

for file in chunk_0_3_3716352_3719310_aug.vg; do vg call -z ${file/_aug.vg/_aug.trans} -s ${file/_aug.vg/_aug.support} $file > ${file/_aug.vg/_calls.vcf}; done

I tried to adapt it to the current release of vg tools and this is what I have tried

vg chunk -x ref_alt.xg -a SRR5647740_paired_wasp_mapped_sorted.gam -g -p 3:3716672-3718990 -c 10 > chunk.vg
vg augment chunk.vg chunk_0_3_3716352_3719310.gam -A aug.gam > aug.vg
vg index aug.vg -x aug.xg
vg pack -x aug.xg -g aug.gam -Q 5 -o aln_aug.pack
vg call aug.xg -k aln_aug.pack > SRR5647740_calls.vcf

Is this correct? Its producing a vcf file but with no variants. The old version for the same input file did. Thanks!

Answer 1

[glenn.hickey]

It's less necessary to break the graph into tiny chunks with vg chunk -p. Here's an example on chunking by chromosome:

https://github.com/vgteam/vg/wiki/Whole-genome-calling-and-genoyping#calling-the-graph-by-first-splitting-into-components

That said, workflows like that should still work. You can use vg depth to verify the coverage in you GAM and pack files. Lack of coverage is the only reason I could think of why call would produce no output. Perhaps the -Q 5 filter in the new workflow is filtering everything out?

[arunkumarramrsh]

Thanks so much! I found out the error had nothing to do with vg tools. It was some odd error encountered when running the program on a HPC. Got it resolved. Thanks!