This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!
homo-ref/snp representation in vcf, Unified Genotyper, snp filtering, homo-ref GQ?
This is my first post/question on this forum. I should say I am relatively new to bioinformatics (self-teaching for the last four months), so my understanding is somewhat on a beginner's level. I tried to see if anyone posted the same/similar question, but I apologize if this is a duplicate.
First, this is what I am trying to accomplish:
- I am analyzing several lines of whole-genome sequence data from two closely-related plant species to a single reference genome. My goal is to create a vcf file with reference and non-reference calls for all these lines. I have successfully done this, but have questions about the output/results and options for cleaning up/filtering the data. Note that I am equally interested in variant/snp sites and reference/non-variant sites because I want to use these data for genome-level population genetic analyses.
Here is my work flow:
I am aligning whole genome sequence data from 19 lines to a single reference line. These 19 lines represent two closely related species, one of which is quite variable among populations. The data is all from paired-end illumina sequencing, mostly 2x75 and 2x100 bp reads, and a few lines with 2x35 bp reads. Average coverage ranges from 5x to 30x. I am using bwa with a q30 cutoff for mapping quality, i.e. I throw out reads that map lower than q30.
I align each line separately using bwa. I dedupe each bam separately and create a dedupped, sorted bam for each line. Following dedupping, I merge all the bams together with RG information preserved. Then I run local realignment on this merged bam. Following realignment, I run the Unified Genotyper to identify variants. (There is not a known set of snps for these species for quality recalibration.) Because I want to identify regions/sites of the genome that are not variable as well as snps, I have been specifying that the Genotyper output all callable sites, i.e. snps AND reference-identical sites. I have been using a quality threshold of 40 to 50.
Finally, I have two observations/questions about the output I am getting:
1) Sites that are statistically reference have quals in the 40s and 50s, i.e. relatively low. This makes sense given that the qual score is related to the probability of a site being statistically non-reference. What seem to be high quality snps have qual scores in the 100s and even over 1000 in somey cases. Many of this sequence data is high coverage, so that is not surprising. However, some of the snps have qual scores just under 100 and or just over. These lower-qual snps tend to be in regions where the surrounding sites, which appear to be mostly identical to the reference based on visually skimming in IGV, have qual scores below 40. These sections tend to have messier data, i.e. several singleton snps (note that this could be due to sequencing errors, or this could represent biological reality if we have sequenced, for example, duplicated regions that are aligning to the same position). The problem is that this creates "skippy" sections in the vcf file, i.e. sections where there are floater snps without any flanking non-variant calls. This is a problem, because I want to use this data for population genetic inference (e.g. this could be upward biasing our estimates of pi). My understanding is that the qual value is a combination of overall callability and the probability that a site is non-reference. So it makes sense that this would happen. What I want to be able to do is generate a vcf file that has quality reference calls and quality snp calls, and that does not contain these "floater snps." However, just increasing my qual threshold to 100 will not produce what I want, because all my reference sites will get filtered out. It would be great if I could somehow filter out these floater snps all together.
So, is there a way to adjust the calling parameters or filter the vcf to achieve what I want?
I do not see an obvious parameter in the Unified Genotyper, but perhaps I am missing something.
I have also read through some of the newer features (e.g. filtering, select variants), and do not see an obvious way here either. It seems I would need to filter based on # of continuous neighboring sites (or something to that effect), which I think would require custom code since this is not a tagged item in the vcf.
2) For the sites with homo ref calls, there is not a GQ value. Is this just a simple function of depth? Or is there a way to specifically request a GQ for homo-ref sites?