Merging overlap variant sites between STACKS .vcf and bcftools call .vcf

[Taryn Cornell]

Hello,

I have used STACKS v2.68 ref_map.pl to find SNPs in a dataset of ~179 samples. My RADseq data is single-end nextRAD. Now, I am hoping to combine this dataset with 23 whole genome short read samples. I used BCFtools mpileup and call to identify SNPs in the 23 whole genomes.  

I want to find overlapping biallelic variant sites. I attempted to do so using the following (simplified) pipeline: bcftools index --> bcftools isec --> bcftools query --> bcftools merge

This yielded an overlap of 20 variants sites between the two datasets. I began with 280,162 (RADseq) and 3,134,304 (WGS). The RADseq data was called using default parameters, in addition to --min-mapq 20 --min-gt-depth 5 and --R 0.75. Both the datasets were filtered for site missingness (<0.1). 

I wanted to know if the .vcf output from STACKS is compatible with this pipeline I described? I used have used --ordered-export and not used it (sorting using BCFtools instead). 

If it is not, are there other recommendations for merging datasets? I have tried using BCFtools to call variants from my RADseq data, but I have much higher missingness per individual and site.

Any input into this situation and these programs is much appreciated. I am a Masters student learning with these data, thank you! 

Taryn