calculating sigma for Depth of Coverage cutoffs

sagipolanisagipolani Posts: 41Member

Hi,

I would like to filter variants based on max DP. I understand that in order to define the max DP cutoof there is a need to calculate the Sigma, as stated in GATK's best practices:

"The maximum DP (depth) filter only applies to whole genome data, where the probability of a site having exactly N reads given an average coverage of M is a well-behaved function. First principles suggest this should be a binomial sampling but in practice it is more a Gaussian distribution. Regardless, the DP threshold should be set a 5 or 6 sigma from the mean coverage across all samples, so that the DP > X threshold eliminates sites with excessive coverage caused by alignment artifacts."

Hence, how can one calculate the sigma, and what is it exactly?

Thanks!

Sagi

Best Answer

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 7,123Administrator, GATK Developer admin
    Answer ✓

    That sounds right to me. As always, make sure to examine your results carefully to make sure that this cutoff is appropriate for your data.

    Geraldine Van der Auwera, PhD

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 7,123Administrator, GATK Developer admin
  • sagipolanisagipolani Posts: 41Member

    Thanks!

    What tool would you recommend for this?

  • sagipolanisagipolani Posts: 41Member

    OK, so I found a nice tool for this: QualiMao (http://qualimap.bioinfo.cipf.es/).

    Now, just to be sure, if I have a mean genome coverage of ~10X and a standard deviation (sigma) of 110, so does that mean that the max depth cutoff for variant calling should be: DP(max) = 10 + 5*110 = 560?

    I'd appreciate if you can confirm this.

    Thanks!

    Sagi

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 7,123Administrator, GATK Developer admin
    Answer ✓

    That sounds right to me. As always, make sure to examine your results carefully to make sure that this cutoff is appropriate for your data.

    Geraldine Van der Auwera, PhD

  • sagipolanisagipolani Posts: 41Member

    Thanks Geraldine!

    Helpful as always :)

    Best,

    Sagi

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 7,123Administrator, GATK Developer admin

    I try :)

    Thanks for linking to that tool, looks like it could be useful to others. From what we can tell it's essentially a visualization engine for Picard tools, so you could get all those metrics using Picard from command line. But sometimes it is nice to have everything all plotted out for you!

    Geraldine Van der Auwera, PhD

  • sagipolanisagipolani Posts: 41Member

    Yes it's pretty neat. Even Though I've seen that there are some inconsistencies between the results generated via qualimap and those generated via samtools depth. I've addressed this issue via qualimap's developers and I am awaiting a reply. Once I get an answer I'll post it back here just in case so that everyone may benefit from it.

    Best,

    Sagi

Sign In or Register to comment.