SNP calling from RNA-seq for population genetics analysis

Hi,

My final goal is to call SNP, from RNA-seq data, for population genetics analysis but I have some problems.

For the population genetics analysis I need a vcf file with all the samples together (that in each locus has the calling info for each sample). My first attempt was to follow the GATK best practices for RNA-seq (Per-Sample Variant Calling and Filtering ) and then do the merge of vcf files but it doesn't work ( the merged vcf file doesn't have all the info per locus that the population genetics program needs).

I also tried with the EMIT_ALL_CONFIDENT_SITES but the clustered SNPs filter doesn't work.

I thought of two possible solutions:

1) Do the SNP calling with all the samples together. Although in the GATK's guide says clearly: "At the moment, we do not recommend applying the GVCF-based workflow to RNAseq data because although there is no obvious obstacle to doing so, we have not validated that configuration. Therefore, we cannot guarantee the quality of results that this would produce."

2) Do the SNP calling per sample but with and without the EMIT_ALL_CONFIDENT_SITES and then merge the results the two vcf files for each sample. In this way I would have the filtered SNPs and the reference call in one vcf file.

What do you recommend me? A third option?

thanks in advance!,

Facundo

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @fagire
    Hi Facundo,

    What do you mean by "the merged vcf file doesn't have all the info per locus that the population genetics program needs"?

    Thanks,
    Sheila

  • fagirefagire Member

    Hi Sheila,

    With that sentence I mean:

    The problem is that unless your vcf's have a line for every single locus, samples that match consensus at a position will not have an entry for that position, so if you have a SNP in one sample at a given position, other samples will have no vcf entry for that position, and you won't know if your other samples really match consensus there, or if they have low coverage there, and a variant couldn't be safely called. (https://www.biostars.org/p/72696/)

    And the population genetic program treats these locus as unknown.

    Thanks,

    Facundo

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @fagire
    Hi Facundo,

    Thanks for clarifying. Yes, the GVCF method is supposed to address your issue. However, you are correct we have not validated it on RNA-seq data yet, although, we have had users experiment with it. You are welcome to try it out on your own data and let us know how it goes. (There should be no technical difficulties). The other option is to run HaplotypeCaller on all samples together.

    Your two options you stated above are correct. Please note though because we have not validated them yet, we cannot guarantee the best results.

    -Sheila

Sign In or Register to comment.