개념설명은 https://www.slideshare.net/secondmath/fst-selection-index 을 참고하기 바란다.

Fst

Load package & Data Load

library(hierfstat);library(data.table);library(knitr);library(DT)
fst.file="https://raw.githubusercontent.com/jinseob2kim/jinseob2kim.github.io/master/fstexample.txt"
a=fread(fst.file)

Read example file. 7 pop & 289 SNPs (PER3 gene)

datatable(a)                            

Basic stat: original Fst

gg=basic.stats(a)
perloc1=gg$perloc                       # per locus statistics
datatable(perloc1) %>% formatRound(1:ncol(perloc1),3)
fstloc1=perloc1$Fst                     # per locus Fst
all1=gg$overall                         # overall locus statistics
kable(t(all1),caption = "overall locus statistics",digits = 3)
overall locus statistics
Ho Hs Ht Dst Htp Dstp Fst Fstp Fis Dest
0.282 0.307 0.332 0.024 0.336 0.029 0.074 0.085 0.081 0.041
fst1=all1[7]                            # overall locus Fst
fst1
##   Fst 
## 0.074

Weir & Cockerham’s theta

gg2=wc(a)
fstloc2=gg2$per.loc$FST                # per locus fst
kable(t(head(fstloc2)),digits=3,caption = "First 6 Fst")
First 6 Fst
1 2 3 4 5 6
0.052 0.052 0.064 0.045 0.061 0.011
fst2=gg2$FST                           # overall locus: mean
fst2
## [1] 0.07931026

Compare

library(ggplot2)
f=data.table(fstloc1,fstloc2)
ggplot(f,aes(fstloc1,fstloc2))+geom_point()+xlab("Wright's Fst")+ylab("Cockerham's theta")+annotate(geom="text", x=0.02, y=0.4,label= paste("r =",round(cor(fstloc1,fstloc2),3)),size=5)

Selection Index: EHH, iHS, xpEHH

Read example file, chr=12: 1424 SNPs & 280 haplotype

library(rehh)
make.example.files()
datatable(fread("bta12_cgu.hap")[1:100],rownames = F)