We've moved!
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!

How to get GenotypeGVCFs to not emit sites with no data

Hello,

I am working with RADseq data and the vast majority of my reference genome thus has no data for any individual. But when I run GenotypeGVCFs, I am getting a huge VCF file because a missing data genotype is listed for every position of the reference (see below). My question is how can I get GenotypeVCFs to emit only confident homozygous ref and confident variant sites and not emit sites with no data in any individual?

Thanks,
Ben

chr1 1 . N . . MLEAC=.;MLEAF=. GT:AD:DP:RGQ ./.:0,0:0:0 ./.:0,0
:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0
:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0
:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0
chr1 2 . N . . MLEAC=.;MLEAF=. GT:AD:DP:RGQ ./.:0,0:0:0 ./.:0,0
:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0
:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0
:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0 ./.:0,0:0:0

Answers

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭
    edited October 2016

    @evansb
    Hi Ben,

    What kind of data are you working with? Can you tell us how you pre-processed the data?

    Thanks,
    Sheila

    Edit: I don't think there is a way to do this with GenotypeGVCFs, but you can use SelectVariants. I asked my questions because it is odd to have so many no-call sites in your final VCF.

  • evansbevansb Member

    Hi Sheila.

    With RADseq data, we expect no data at all for most of the genome and then high coverage near rare cutting restriction sites (SbfI sites, in my case). So the distribution of no-call sites is what I expect, I just wish there were a flag in GenotypeGVCF that would allow us to leave out the sites with no data. Anyhow, thank you for answering my question - my next step is to filter as you point out.

    Thanks!
    Ben

  • SheilaSheila Broad InstituteMember, Broadie ✭✭✭✭✭

    @evansb
    Hi Ben,

    Ah, I see. Did you run the whole Best Practices pipeline on the entire genome? You can use -L to only run on intervals you are interested in. Have a look at this article as well.
    I hope this helps!

    -Sheila

Sign In or Register to comment.