Greetings, @YangJingqii ! Thanks for the questions.

1. snpgenie_within_group calculates the π values for each codon, whereas the bootstrap script performs the statistical test of πN=πS much more efficiently. If I did it all over again I'd probably use Python and combine it all into one, but unfortunately this two-step process is what I have. Regarding the error, please feel free to provide a minimum reproducible example so I can address it.
2. snpgenie.pl doesn't itself bootstrap because it's a computationally intensive second step (see 1) but also because it's often not necessary with intrahost samples. The reason is that the sample rather than the codon is the unit of observation. Thus, because each sample has both a πN and πS, you can do standard statistics like a paired Wilcoxon test of πN = πS. I'm not sure what the second part of the question is asking, but using trees can be both more powerful when done well, but also assumes the tree is correct and therefore can also be more subject to error. Usually with intrahost samples there's so little variation as to make trees worthless.
3. πN/πS and mean dN/dS are the same basic calculation, except sometimes a correction like Jukes-Cantor is applied for the former. Such correction is not necessary for the very low levels of divergence seen in intrahost samples. If you're instead asking about dN and dS vs. the reference sequence, this is no longer all pairwise comparisons, but instead all sequence reads vs. the reference sequence. It's a measure of divergence (from the REF) rather than diversity.
4. I don't understand the question, what is the basis — what is meant by 'basis'?

Thanks!
Chase

snpgenie_within_group calculates the π values for each codon,I want to know if we can calculate piN/piS value for each codon. However, I'm uncertain whether the piN/piS value should be 0 or undefined when S_sites = 0, N_diffs = 0, and N_sites = 3 for a codon.

@YangJingqii thank you, it is possible to calculate πN/πS for individual codons, but there can be substantial uncertainty, and in the extreme many codons will have πS=0 or πS undefined. In the instance you show, πS is undefined, and πN/πS is undefined as well. Depending on whether you're just offering descriptive statistics or trying to test a specific hypothesis, other approaches with πN/πS (i.e. identifying a priori meaningful groups of codons which can be bootstrapped) or even dN/dS may be necessary.

Thank you for your response very much!
We know the formula for calculating diversity is
.
If I only want to calculate the diversity of variants with a frequency lower than 50%, for example, at a certain site where the reference is C (count 470), one alt is A (count 15), and another alt is T (count 800), in the case of variants with a frequency below 50%, T would not be included. If I directly delete this row, will snpgenie mistakenly count these 800 reads as reference C, resulting in an incorrect calculation? If I want to calculate the diversity of only those variants with a frequency below 50%, can snpgenie be used? What form should I use for this(In the second table, the count of T is subtracted from the coverage. In the third table, the count of T and C is subtracted from the coverage. )? Or is it that snpgenie can only calculate the diversity of all variants at this site and cannot specify just a subset (since reference reads are always included, and the reference reads are often greater than 50%)? If there is only one low-frequency variant at a site, is the diversity for that site considered to be 1? Or does it have to be calculated together with the reference?However, in the paper "Within-host genetic diversity of SARS-CoV-2 lineages in unvaccinated and vaccinated individuals", snpgenie was used to calculate the diversity of variants below 50%, which confuses me. Could you please clarify? Thank you very much!

Hello @YangJingqii !

Regarding the possibility of ignoring certain variant frequencies for calculating π, this is no longer π. In other words, π is a property of a site, not a property of particular variants. (The one exception is excluding certain variants or frequencies which you suspect to be sequencing errors, e.g., AF<1%.) What you are suggesting is to calculate the mean number of pairwise differences among only a subset of variants, which is not customary and would probably be difficult to interpret. It's always possible to invent a new metric, but I'm not sure of your goal or what the meaning of the resultant metric would be.

Regarding the paper you cite, I am not a coauthor so I do not know what the authors mean by only calculating variation among variants lower than 50%. Looking at the paper, I don't see a statement like that in the main text. Could you point out which part you are referring to?

Chase