Does incomplete coverage affect SNPGenie accuracy? And how to determine valid sites with lofreq?

[YangJingqii]

Hi Nelson,
I'm using SNPGenie to analyze pooled NGS data and have two related questions about coverage.
Background:

My sample inclusion criteria: ≥80% of sites with >200x depth
In practice, it's nearly impossible to achieve 100% coverage — some sites are always under-covered or not covered at all

Questions:

Does incomplete or uneven coverage affect the accuracy of SNPGenie's estimates (e.g., π, dN/dS)? If so, what's the recommended way to handle under-covered sites?
I'm using lofreq for SNP calling, which explicitly states that users don't need to filter by depth since it handles this internally. However, this means the output doesn't tell me which sites were effectively analyzed. How should I determine the denominator for per-site diversity calculations? Do I need to independently calculate coverage from BAM files to define valid sites?

Thanks!

[singing-scientist]

Greetings @YangJingqii, thanks for using SNPGenie! These are great questions and I think there is really no standard answer.

Regarding uneven coverage, the (unsatisfying) answer is, it depends. If there is some minimum coverage that you believe gives you enough confidence and power to detect minor alleles, you might simply use that. One rule of thumb (for better or worse) is to use 10x the reciprocal of your desired minimum allele frequency. For example, if you'd like to detect intrahost SNPs (also called iSNVs) down to AF = 5%, then you might consider only sites with 10 * (1 / 0.05) = 200 coverage or more. But is there really a meaningful difference between sites with 199 and 201 coverage? No. Also if there are big regions with disturbingly low (or high) coverage, you'll of course want to investigate that (is there deletion? duplication? mismatches within the PCR primer regions?). In reality this will require exploring your data and seeing what makes the most sense, which is often an iterative process.

I cannot speak on how lofreq works, as I haven't used that caller extensive. I have no doubt that lofreq considers depth in its calculations and filtering, but I'm a bit skeptical that one would NEVER need to filter by depth. However, if I recall it does provide the coverage or depth—this may be a count of the number of reads it decided could be included (i.e., passed its QC) at that site. Unfortunately I cannot speak further on lofreq, and it may be best to contact them with those questions.

Hope this helps!
Chase

[YangJingqii]

Hi Nelson,
Thanks for your previous reply! I'd like to follow up on a question about coverage (proportion of sites covered), not depth.
Background:
My variant calling pipeline only calls SNPs at sites with ≥200x depth. This means sites below 200x will not appear in the VCF at all.
Question about π:
When SNPGenie calculates π, it uses gene length as the denominator. However, if some sites are not called due to insufficient depth, they won't appear in the VCF. SNPGenie will treat these sites as "no variation," while the denominator remains the full gene length. This would underestimate π.
My idea is to independently calculate the number of sites with ≥200x coverage from BAM files as the "effective length," then correct:
π_corrected = π_snpgenie × (gene_length / effective_length)
Does this approach make sense?
Question about πN/πS:
This gets more complicated, since N_sites and S_sites are calculated based on codon structure and cannot be simply corrected by a ratio. For this situation:

Is there a recommended approach to handle this?
Or in practice, if sample coverage is sufficiently high (e.g., ≥80% of sites meeting the threshold), is it acceptable to just use SNPGenie's output directly without additional normalization?

Thanks!

[singing-scientist]

No problem! Ah, now I understand the question better. Yes, if some sites are being excluded (e.g., because of low coverage) you could calculate π across just the sites that pass, and I agree that's a good idea if there is bad or inconsistent coverage. (Of course, you'd still probably want to investigate why, in a high-coverage environment, many sites fail to be covered.)

I think a better (more direct and general) calculation might be π = sum(π) / L, where you sum the π values across the L sites that are being included. For example, simply apply a filter operation in R or Excel and do the necessary summation—depending on how it's being used, SNPGenie might include the depth in the sites_results file (IF it was given by lofreq) which would make such filtering quick and easy after a simple left join. (Note this is the same as the simple mean of π across the subset of sites that "pass".)

For πN or πS the best calculation would be πN = sum(N_diffs) / sum(N_sites) and πS = sum(S_diffs) / sum(S_sites), again limiting to only including the codons that "pass". If codons contain a mix of passing and failing sites, you'll have to decide how to handle those (e.g., exclude codons with even one failing site?).

Let me know if that makes sense.

[YangJingqii]

That's a brilliant solution! I really appreciate it.
I've decided to use the conservative approach: for each codon, I'll calculate the minimum depth across all three nucleotide positions. If the minimum depth < 200x, I'll exclude the entire codon from the calculation.
Thanks again for the clear explanation!

[singing-scientist]

Thanks so much! That all sounds great, and it's heartening that people like you are thinking through these important upstream issues. Please don't hesitate to let me know if you have any further questions!

Chase