Separating indels and snps in vcf

[balsampoplar]

Task at hand: Remove simple indel loci from vcf file.

------ example vcf file excerpt ----------
CHR POS ID REF ALT
1 20293 S_20293 A G
1 22689 S_22689 A -
1 23251 S_23251 - T
--------------------------------------------

$ bcftools filter -e'%TYPE="snp"' in.vcf > indels.vcf
Results in a VCF file containing only header

$ bcftools filter -e'%TYPE="indels"' in.vcf > snp.vcf
Results in a VCF file containing all variants (snps and indels)

It appears that bcftools is not treating the last two sites above as indels, but as snps. I saw the same behvior with --remove-indels filter in vcftools.

Does anyone know what's going on?

[sarvidsson]

What VCF version should that be (and what software called those)? It is not proper VCF4.1 or 4.2 - dashes are not allowed as bases, so no wonder that the tools you tried couldn't identify those as indels.

VCF 4.2 spec, check section 1.4.1 for the REF and ALT fields, as well as section 5.2 for examples with properly formatted indels: http://samtools.github.io/hts-specs/VCFv4.2.pdf

[balsampoplar]

sarvidsson: My VCF files were generated using the Tassel GBS pipeline and the header indicates they conform to 4.0 specification.

[sarvidsson]

Then their conformity is broken: http://www.1000genomes.org/node/101; the allowed contents of the REF and ALT fields are the same for 4.0 as 4.1.

[balsampoplar]

That's interesting. According to the 1000genomes 4.0 spec, the '.' in the ALT column indicates monomorphic site. That's clearly not the same as a '-' which is supposed to indicate an indel, right?

[sarvidsson]

A quick and dirty way to filter these files would be to

grep -v $'\t-\t' in.vcf > no_indels.vcf

in case there are no other tab-separated fields with dashes. Check that with:

grep $'\t-\t' in.vcf > indels.vcf

[sarvidsson]

That's interesting. According to the 1000genomes 4.0 spec, the '.' in the ALT column indicates monomorphic site. That's clearly not the same as a '-' which is supposed to indicate an indel, right?

Exactly. I'd check with the TASSEL people that those records are really indels, perhaps they can fix the VCF output for a future version.

[balsampoplar]

Thanks. One liners are always useful. I have been doing this by first subsetting REF/ALT columns containing dashes in R and then using that list with --positions filter in vcftools to separate out snp and indels in individual vcf files.

Still I would like to figure out what dash versus periods mean in the REF/ALT fields.

Edit: Will post a question to the Tassel list.

[sarvidsson]

Thanks. One liners are always useful. I have been doing this by first subsetting REF/ALT columns containing dashes in R and then using that list with --positions filter in vcftools to separate out snp and indels in individual vcf files.

You don't need to load the whole thing in R for that,

cut -f 3,4 in.vcf | grep '\-$' | cut -f 1 > ref_dash_ID-list
cut -f 3,5 in.vcf | grep '\-$' | cut -f 1 > alt_dash_ID-list

would be much faster

[balsampoplar]

I have been only bringing in the first 5 columns in R, but I like your solution as well. Thanks.