Folded SFS?

[JoshuaDanielRubin]

Good morning,

(I don't know if this repo is actively maintained so no worries if I dont hear back)

Thanks for this tool, it is very useful. I was just wondering about something with the SFS. I first polarized my VCF, and then did vcfdo sfs to create an sfs based on the newly populated AA fields. However I find that my SFS is folded, i.e. for a population vcf of 10 diploid individuals, the sfs is length 10 (10 bins). I dont understand why, as I would have thought it should be 20 if it were unfolded?

I need an unfolded sfs for my further processing, so hopefully you can help me clarify what is going on.

Best,
Joshua

[andrewparkermorgan]

The short answer to your question is that the SFS are unfolded but assume a haploid genome. (In the original use case, I was working with pseudo-diploid genotype calls from biological samples that were actually derived from mixtures of haploid genomes.) Changing that behavior would not be too difficult, if it is of interest, but would take me some time.

Also note that the SFS produced by this tool is not just counting alleles -- it performs draws from a binomial distribution for which the expected proportion is determined by either the genotype, the allelic balance in the read counts or the genotype posterior (FORMAT/GP). If you want simple allele counts you should use something different.

[JoshuaDanielRubin]

Ah, I see. Thank you for the clarification.

Apropos of nothing, I don't suppose you might have a suggestion about the easiest way to get unfolded SFS for human samples (simple allele counts)? I have found it surprisingly difficult to find a tool which can do this.

[andrewparkermorgan]

It seems like there would be a "standard" tool but I think there are some complications just below the surface, in particular how to handle missing data. The script I wrote for myself (before writing vcfdo) is here: https://gist.github.com/andrewparkermorgan/e705fa4ccaa7c6f77efae2de93899a1a

It should work as advertised for unfolded SFS in diploids. It includes the option to produce a folded SFS but that is not well tested, so I would not trust it. Missing data is handled as follows. If there are n individuals (assume all diploid), the size of the unfolded SFS should be 2n + 1. Suppose m < n individuals have non-missing genotypes, and that the observed number of derived alleles is k. Let q = int( [ k/(2m) ] * 2 n ) and give credit to the q-th bin in the unfolded SFS. That is, assume that genotypes are missing at random so the allele frequency among the non-missing individuals is the same as in the whole sample. It's probably a reasonable approach if the amount of missing data is small.

More complicated but maybe more principled is a tool like easySFS which handles missing data by "projecting" the SFS to a smaller sample size. See that repo for details.

[JoshuaDanielRubin]

Understood! thanks so much for the information and for creating this tool, it's super useful.