I am conducting a populations genetics study using SNPs derived from RADseq data (total of 18469 SNPs for 28 individuals and 16 populations). In this context I need to compute pairwise Fst between all my populations. I would like to know which methods would be best to compute the Fst from big SNP datasets and with small sampling sizes for my populations (I have 2 individuals for most of the populations, and some with only 1). If I can also compute other statistics like Qst and Gst, it would be perfect.
I already tried using different R packages (pegas, hierfstat, popgenome), but none of them were very satisfying so I would like to have some opinions on that.