Four diploid individuals were sequenced at a 20 base pair (bp) locus. Every individual contributes two haplotype sequences that we will call a and b. Thus, 4 diploid individuals give eight sequences in the sample. All sites are biallelic with no missing data.
An extra sequence is shown here from an outgroup. We do not count it in our SFS. It is only here to help with polarizing the unfolded SFS.
The table is laid out similar to a variant call format (VCF) file, with individuals as columns and sites rows.
site outgroup ind1_a ind1_b ind2_a ind2_b ind3_a ind3_b ind4_a ind4_b
---- -------- ------ ------ ------ ------ ------ ------ ------ ------
1 A A A A A A A A A
2 C C C C C C C C C
3 G A G G G G G G G
4 T T T T T T T T T
5 A A A A A A A A A
6 C C C C C T C C C
7 G G G G G G G G G
8 T T T T T T T T T
9 A A G G A A A A A
10 C C C C C C C C C
11 G G G G G G G G G
12 T C T T C T C T T
13 A A A A A A A A A
14 C C C C C C C C C
15 G A A A A A A G G
16 T T T T T T T T T
17 A A A A A A A A A
18 C C C C C C C C C
19 G G T T T T T T T
20 T T T T T T T T T
Here is that same sequence data represented as a simplified VCF file. Notice that we drop any sites that are not segregating in the sample. Our initial outgroup is now the REF column, since we used it as the reference genome for variant calling. The two sequences from an individual are collapsed into genotype.
site REF ALT ind1 ind2 ind3 ind4
---- --- --- ----- ----- ----- -----
3 G A 0/1 0/0 0/0 0/0
6 C T 0/0 0/0 0/1 0/0
9 A G 0/1 0/1 0/0 0/0
12 T C 0/1 0/1 0/1 0/0
15 G A 1/1 1/1 1/1 0/0
19 G T 0/1 1/1 1/1 1/1
Here, 0 means the REF base and 1 means the ALT base — so 0/0 is homozygous reference, 0/1 is a heterozygote, and 1/1 is homozygous alternate. The fourteen sites where every sequence agreed were removed to make our lives easier, but a real VCF can retain information on all sites.
Be careful with REF and ALT. They are labels for the two bases seen at a site, and REF is simply whichever one the reference genome happens to carry. We often try to chose a reference genome that is close to but not a part of the sample, and assume that the reference base is the ancestral allele. That assumption is sometimes wrong, which is why the folded SFS can be useful.
The site frequency spectrum is a histogram, nothing more. For now, we will only think of the unfolded SFS since it is a little easier to count. For each variable site you count how many of the eight sequences carry the ALT allele. If three sequences carry the ALT allele, you add one to the bin for derived count 3. If seven sequences carry the ALT allele, you add one to the bin for derived count 7.
Find the segregating sites
How many segrgating sites (\(S\)) are there in the sample? The VCF where we only kept segregating sites makes this obvious.
\(S =\) ____________ segregating sites, out of \(L = 20\).
Most of the twenty rows are invariant. Those sites are real data, but \(S\) only cares about variation and the SFS only has bins for counts of 1 through \(n-1\).
Which alleles are derived?
At each segregating site, one of the two bases matches the outgroup (REF) and this is encoded as 0. The ALT base is encoded as 1.If we consider the REF ancestral and the ALT derived, we can get the derived count by simply taking the sum of that row.
Count how many of the eight sequences carry the derived base for each segregating site.
Site
Derived count \(i\)
Build the unfolded spectrum
The unfolded spectrum \(\xi_i\) is the number of sites whose derived count is \(i\). With \(n=8\) sequences, \(i\) goes from 1 to 7. Go through your derived counts and tally them.
\(i\) (derived count)
1
2
3
4
5
6
7
\(\xi_i\) (sites)
Sites in bin 1 are singletons — a derived allele seen once, in one sequence. In a large real sample singletons are always the tallest bar, because most mutations are young and most young mutations are still rare.
Fold the spectrum
Maybe your referenece is not a suitable outgroup for polarization or the derived specrum visably has problems (the “U” shape). If you are not confident that the derived SFS is accurate, you can use the minor-allele count instead:
\[\text{minor count} = \min(i,\; n - i)\]
Maybe a site has 7 derived copies has 1 ancestral copy. If we do not know the ancestral state, this could be 7 ancestral and 1 derived. The minor count here is 1. If a site has 3 alleles that match the REF and 5 match the ALT, the minor count is 3. The minor count is always \(\le n/2\).
Site
Derived count \(i\)
Minor count \(\min(i, n-i)\)
The folded spectrum \(\eta_i\) only needs bins up to \(\lfloor n/2 \rfloor = 4\).
\(i\) (minor count)
1
2
3
4
\(\eta_i\) (sites)
Again \(\sum_i \eta_i = S\). Folding does not throw sites away — it puts two bins on top of each other.
Summarize the spectra
\(i\)
1
2
3
4
5
6
7
unfolded \(\xi_i\)
folded \(\eta_i\)
—
—
—
—
Interpret
Your unfolded spectrum has sites out at 6 and 7 copies; your folded spectrum has nothing above
Both still total \(S\), so no information was lost. In which bins did the sites with 7 and 7 derived alleles go?
Sites 3, 6 and 19 all sit in the same folded bin. Two of them are brand-new mutations sitting at one copy. The third is an old allele that has nearly swept through the population. Which is which? What assumption would you have to make to tell them apart?
How can you diagnose if your referece genome is a suitable outgroup for polarization? What would you do if it is not?
The same spectrum, from a VCF
The file sfs-example.vcf holds the same six sites and the same four individuals. The Info field is left blank with a “.”, but if you had other external information about ancestral states, you could store that information here.
library(vcfR)vcf <-read.vcfR("sfs-example.vcf", verbose =FALSE)samples <-colnames(vcf@gt)[-1] # drop the FORMAT columnn <-2*length(samples) # diploids -> chromosomes# Assume the REF allele is ancestral, so the ALT count is the derived countaa <- vcf@fix[, "REF"]# count ALT alleles per site: split each genotype on / or | and count the 1sgt <-extract.gt(vcf, element ="GT")derived <-apply(gt, 1, function(g) sum(unlist(strsplit(g, "[/|]")) =="1"))minor <-pmin(derived, n - derived)data.frame(POS = vcf@fix[, "POS"], REF = vcf@fix[, "REF"], ALT = vcf@fix[, "ALT"],AA = aa, ALT_count = derived, derived_count = derived, minor_count = minor)
POS REF ALT AA ALT_count derived_count minor_count
chr1_3 3 G A G 1 1 1
chr1_6 6 C T C 1 1 1
chr1_9 9 A G A 2 2 2
chr1_12 12 T C T 3 3 3
chr1_15 15 G A G 6 6 2
chr1_19 19 G T G 7 7 1
Critical thing to remember: \(n\) is 8, not 4. The spectrum counts the number of haplotypes, not individuals.
xi <-tabulate(derived, nbins = n -1)eta <-tabulate(minor, nbins = n %/%2)cat("unfolded xi[1..7] =", xi, "\n")
unfolded xi[1..7] = 2 1 1 0 0 1 1
cat("folded eta[1..4] =", eta, "\n")
folded eta[1..4] = 3 2 1 0
par(mfrow =c(1, 2), mar =c(4, 4, 2, 1))barplot(xi, names.arg =1:(n -1), col ="#CC6677",xlab ="derived count", ylab ="sites", main ="Unfolded")barplot(eta, names.arg =1:(n %/%2), col ="#117733",xlab ="minor count", ylab ="sites", main ="Folded")
Compare these to the tables you filled in by hand. Do they match? If so, you can make an SFS.
