fst: Haplotype-based Fst

In ytutsunomiya/GHap: Genome-Wide Haplotyping

Haplotype-based Fst

Description

Multi-alleleic Fst computed using block summary statistics generated from ghap.blockstats.

Usage

 ghap.fst(blockstats.pop1,
          blockstats.pop2,
          blockstats.tot)

Arguments

blockstats.pop1	A data.frame containing block statistics computed on population 1.
blockstats.pop2	A data.frame containing block statistics computed on population 2.
blockstats.tot	A data.frame containing block statistics computed on population 1 + population 2.

Details

This function calculates Fst (Nei, 1973) based on the formula for multi-allelic markers:

Fst = (Ht - Hs) / Ht

where Ht is the total gene diversity (i.e., expected heterozygosity in the population) and Hs is the subpopulation gene diversity (i.e., the average expected heterozygosity in the subpopulations).

Value

The function returns a data.frame with the following columns:

BLOCK	Block alias.
CHR	Chromosome name.
BP1	Block start position.
BP2	Block end position.
EXP.H.pop1	Expected heterozygosity in population 1.
EXP.H.pop2	Expected heterozygosity in population 2.
EXP.H.tot	Expected heterozygosity in the total population.
FST	Fst value.

Author(s)

Yuri Tani Utsunomiya <[e-mail]>
Marco Milanesi <[e-mail]>

References

M. Nei. Analysis of Gene Diversity in Subdivided Populations. PNAS. 1973. 70, 3321-3323.

Examples

# #### DO NOT RUN IF NOT NECESSARY ###
# 
# # Copy phase data in the current working directory
# exfiles <- ghap.makefile(dataset = "example",
#                          format = "phase",
#                          verbose = TRUE)
# file.copy(from = exfiles, to = "./")
# 
# # Load data
# phase <- ghap.loadphase("example")
# 
# # Generate blocks of 5 markers
# blocks <- ghap.blockgen(phase, windowsize = 5,
#                         slide = 5, unit = "marker")
# 
# # Haplotyping
# ghap.haplotyping(phase = phase, blocks = blocks, outfile = "example",
#                  binary = T, ncores = 1)
# 
# # Load haplotype genotypes using prefix
# haplo <- ghap.loadhaplo("example")
# 
# ### RUN ###
# 
# # Compute block statistics for population 1
# ids <- which(haplo$pop == "Pure1")
# haplo <- ghap.subset(haplo, ids = ids,
#                      variants = haplo$allele.in,
#                      index = TRUE)
# hapstats1 <- ghap.hapstats(haplo)
# blockstats1 <- ghap.blockstats(hapstats1)
# 
# # Compute block statistics for population 2
# ids <- which(haplo$pop == "Pure2")
# haplo <- ghap.subset(haplo, ids = ids,
#                      variants = haplo$allele.in,
#                      index = TRUE)
# hapstats2 <- ghap.hapstats(haplo)
# blockstats2 <- ghap.blockstats(hapstats2)
# 
# # Compute block statistics for combined populations
# ids <- which(haplo$pop %in% c("Pure1","Pure2"))
# haplo <- ghap.subset(haplo, ids = ids,
#                      variants = haplo$allele.in,
#                      index = TRUE)
# hapstats12 <- ghap.hapstats(haplo)
# blockstats12 <- ghap.blockstats(hapstats12)
# 
# # Compute FST
# fst <- ghap.fst(blockstats1, blockstats2, blockstats12)
# ghap.manhattan(data = fst, chr = "CHR", bp = "BP1",
#                y = "FST", type = "h")

ytutsunomiya/GHap documentation built on Oct. 15, 2024, 5:27 p.m.