1 Purpose of the analysis

This analysis estimates the Weir and Cockerham (1984) fixation indices for two maize populations. The three statistics of interest are:

\(\theta\), interpreted as \(F_{ST}\), which measures genetic differentiation between populations;
\(f\), interpreted as \(F_{IS}\), which measures heterozygote deficiency within populations; and
\(F\), interpreted as \(F_{IT}\), which measures heterozygote deficiency relative to the total population.

The analysis uses genotype data from the Genomes to Fields 2018 dataset. The original VCF contains 312 samples and 1,353,525 variants. Sixty individuals are sampled from the Mo44_PHW65 population and 70 individuals are sampled from the PHN11_PHW65 population. A reproducible random sample of 2,000 variants is taken before quality control.

The final analysis obtained:

Statistic	Estimate
\(\theta\) (\(F_{ST}\))	0.3329271
\(f\) (\(F_{IS}\))	1.0000000
\(F\) (\(F_{IT}\))	1.0000000

2 Understanding the input files

Two files are required.

The VCF file contains the genetic variants and genotypes.
The sample-decoder file connects each VCF sample ID to its genetic line.

The VCF genotype field commonly uses the following coding:

VCF genotype	Biological meaning	Numeric dosage
`0/0`	Two reference alleles	0
`0/1` or `1/0`	One reference and one alternative allele	1
`1/1`	Two alternative alleles	2
`./.`	Missing genotype	`NA`

The dosage is the number of alternative alleles carried by an individual at a SNP.

3 Load the required package

The vcfR package reads and manipulates VCF files. Installation is required only once and is therefore shown but not executed automatically during knitting.

install.packages("vcfR")

library(vcfR)

4 Define the file locations

The following paths assume that the two downloaded files are stored in a folder named data inside the user’s home directory.

data_dir <- path.expand("~/data")

vcf_file <- file.path(
  data_dir,
  "G2F_PHG_minreads1_Mo44_PHW65_MoG_assemblies_14112019_filtered_plusParents.vcf"
)

decoder_file <- file.path(
  data_dir,
  "G2F_PHG_minreads1_Mo44_PHW65_MoG_assemblies_14112019_filtered_plusParents_sampleDecoder.txt"
)

stopifnot(file.exists(vcf_file))
stopifnot(file.exists(decoder_file))

5 Read and prepare the sample decoder

The decoder has no header. Therefore, header = FALSE is necessary. If header = TRUE were used, R would incorrectly treat the first individual as a column name.

decoder <- read.table(
  decoder_file,
  header = FALSE,
  stringsAsFactors = FALSE
)

names(decoder) <- c("VCF_ID", "Genetic_Line")

head(decoder)

##     VCF_ID     Genetic_Line
## 1 GE18N213 PHN11_PHW65_0206
## 2 GE18N080  Mo44_PHW65_0316
## 3 GE18N143 PHN11_PHW65_0042
## 4 GE18N031  Mo44_PHW65_0134
## 5 GE18N070  Mo44_PHW65_0052
## 6 GE18N351   PHW65_MoG_0621

dim(decoder)

## [1] 312   2

The first column is the sample name used in the VCF. The second column is the corresponding genetic line.

6 Create the population variable

The beginning of each genetic-line name identifies its population. In the regular expressions below, ^ means “begins with.”

decoder$Population <- ifelse(
  grepl("^Mo44_PHW65", decoder$Genetic_Line),
  "Mo44_PHW65",
  ifelse(
    grepl("^PHN11_PHW65", decoder$Genetic_Line),
    "PHN11_PHW65",
    ifelse(
      grepl("^PHW65_MoG", decoder$Genetic_Line),
      "PHW65_MoG",
      "Other"
    )
  )
)

head(decoder)

##     VCF_ID     Genetic_Line  Population
## 1 GE18N213 PHN11_PHW65_0206 PHN11_PHW65
## 2 GE18N080  Mo44_PHW65_0316  Mo44_PHW65
## 3 GE18N143 PHN11_PHW65_0042 PHN11_PHW65
## 4 GE18N031  Mo44_PHW65_0134  Mo44_PHW65
## 5 GE18N070  Mo44_PHW65_0052  Mo44_PHW65
## 6 GE18N351   PHW65_MoG_0621   PHW65_MoG

table(decoder$Population)

## 
##  Mo44_PHW65       Other PHN11_PHW65   PHW65_MoG 
##         100           6         142          64

The decoder contains 100 Mo44_PHW65 individuals, 142 PHN11_PHW65 individuals, 64 PHW65_MoG individuals, and 6 other entries. The first two populations are used because they can provide the required sample sizes of 60 and 70.

7 Read the VCF and select the working dataset

Reading the full VCF uses several gigabytes of memory. To avoid repeating that expensive operation, the code first checks whether the reduced 130-individual, 2,000-variant object has already been saved.

filtered_rds <- file.path(
  data_dir,
  "geno_filtered_130individuals_2000SNPs.rds"
)

set.seed(123)

pop1_ids <- decoder$VCF_ID[
  decoder$Population == "Mo44_PHW65"
]

pop2_ids <- decoder$VCF_ID[
  decoder$Population == "PHN11_PHW65"
]

selected_pop1 <- sample(
  pop1_ids,
  size = 60,
  replace = FALSE
)

selected_pop2 <- sample(
  pop2_ids,
  size = 70,
  replace = FALSE
)

selected_ids <- c(selected_pop1, selected_pop2)

if (file.exists(filtered_rds)) {
  geno_filtered <- readRDS(filtered_rds)
} else {
  geno <- read.vcfR(
    vcf_file,
    verbose = TRUE
  )

  n_snps <- nrow(geno@fix)
  n_individuals <- ncol(geno@gt) - 1
  vcf_ids <- colnames(geno@gt)[-1]

  stopifnot(n_snps == 1353525)
  stopifnot(n_individuals == 312)
  stopifnot(sum(decoder$VCF_ID %in% vcf_ids) == 312)

  snp_idx <- sort(
    sample(
      seq_len(nrow(geno@fix)),
      size = 2000,
      replace = FALSE
    )
  )

  geno_filtered <- geno[
    snp_idx,
    c("FORMAT", selected_ids)
  ]

  saveRDS(geno_filtered, filtered_rds)
  rm(geno)
  gc()
}

geno_filtered

## ***** Object of Class vcfR *****
## 130 samples
## 10 CHROMs
## 2,000 variants
## Object size: 2.4 Mb
## 3.284 percent missing data
## *****        *****         *****

length(selected_pop1)

## [1] 60

length(selected_pop2)

## [1] 70

set.seed(123) makes the sampling reproducible. replace = FALSE ensures that an individual or SNP cannot be selected more than once.

8 Retain biallelic SNPs

A biallelic SNP has only a reference allele and one alternative allele. The simple dosage coding 0, 1, and 2 assumes that there are only two alleles.

geno_bi <- geno_filtered[
  is.biallelic(geno_filtered),
]

geno_bi

## ***** Object of Class vcfR *****
## 130 samples
## 10 CHROMs
## 1,891 variants
## Object size: 2.3 Mb
## 3.338 percent missing data
## *****        *****         *****

9 Extract and standardize the genotype field

The genotype field is extracted from the VCF. A phased genotype uses |, whereas an unphased genotype uses /. Replacing | with / gives one consistent representation.

gt <- extract.gt(
  geno_bi,
  element = "GT"
)

gt_clean <- gsub(
  pattern = "\\|",
  replacement = "/",
  x = gt
)

table(gt_clean, useNA = "ifany")

## gt_clean
##    0/0    1/1   <NA> 
## 150715  86909   8206

10 Convert genotypes to the numeric matrix

The numeric matrix is called M_wc. Missing genotypes remain missing rather than being imputed.

M_wc <- matrix(
  NA_real_,
  nrow = nrow(gt_clean),
  ncol = ncol(gt_clean),
  dimnames = dimnames(gt_clean)
)

M_wc[gt_clean == "0/0"] <- 0
M_wc[gt_clean == "0/1" | gt_clean == "1/0"] <- 1
M_wc[gt_clean == "1/1"] <- 2

# Transpose exactly once: individuals become rows and SNPs become columns.
M_wc <- t(M_wc)

dim(M_wc)

## [1]  130 1891

range(M_wc, na.rm = TRUE)

## [1] 0 2

The transpose must be performed exactly once. Before transposition, the SNPs are rows and individuals are columns. After transposition, the 130 individuals are rows and the SNPs are columns.

11 Match the population labels to the matrix rows

population_wc <- decoder$Population[
  match(
    rownames(M_wc),
    decoder$VCF_ID
  )
]

stopifnot(sum(is.na(population_wc)) == 0)
table(population_wc)

## population_wc
##  Mo44_PHW65 PHN11_PHW65 
##          60          70

The expected table contains 60 Mo44_PHW65 individuals and 70 PHN11_PHW65 individuals.

12 Remove unusable SNPs

First, SNPs with no observed genotype are removed. Next, fixed SNPs are removed. A fixed SNP has allele frequency 0 or 1 and therefore contains no variation.

has_observations <- colSums(!is.na(M_wc)) > 0
M_wc <- M_wc[, has_observations, drop = FALSE]

allele_frequency_wc <- colMeans(
  M_wc,
  na.rm = TRUE
) / 2

keep_polymorphic <- (
  is.finite(allele_frequency_wc) &
    allele_frequency_wc > 0 &
    allele_frequency_wc < 1
)

M_wc <- M_wc[, keep_polymorphic, drop = FALSE]

dim(M_wc)

## [1]  130 1237

table(population_wc)

## population_wc
##  Mo44_PHW65 PHN11_PHW65 
##          60          70

sum(is.na(M_wc))

## [1] 5273

range(M_wc, na.rm = TRUE)

## [1] 0 2

The final matrix contains 130 individuals and 1,237 informative polymorphic SNPs. Missing genotypes are deliberately retained.

13 Why missing values are not imputed here

Mean imputation is useful when constructing a VanRaden genomic relationship matrix because matrix multiplication requires complete numeric data. It is not appropriate for this Weir-Cockerham calculation.

For example, an imputed value of 0.73 is not a biological genotype. It cannot be classified as 0, 1, or 2 and cannot be correctly counted as homozygous or heterozygous. Instead, each SNP is analysed using its observed individuals, and its sample size is calculated after excluding missing values.

14 Weir and Cockerham functions

14.1 Sample-size correction

wc_Nc <- function(n) {
  r <- length(n)
  n_bar <- mean(n)

  (
    r * n_bar -
      sum(n^2) / (r * n_bar)
  ) / (r - 1)
}

14.2 Variance components

wc_components <- function(n, p, h) {
  r <- length(n)
  n_bar <- mean(n)
  n_c <- wc_Nc(n)

  p_bar <- sum(n * p) / sum(n)

  s2 <- sum(
    n * (p - p_bar)^2
  ) / (n_bar * (r - 1))

  h_bar <- sum(n * h) / sum(n)

  a <- (n_bar / n_c) *
    (
      s2 -
        (
          p_bar * (1 - p_bar) -
            ((r - 1) / r) * s2 -
            h_bar / 4
        ) / (n_bar - 1)
    )

  b <- (n_bar / (n_bar - 1)) *
    (
      p_bar * (1 - p_bar) -
        ((r - 1) / r) * s2 -
        ((2 * n_bar - 1) / (4 * n_bar)) * h_bar
    )

  c <- h_bar / 2

  list(
    a = a,
    b = b,
    c = c,
    p_bar = p_bar,
    s2 = s2,
    h_bar = h_bar
  )
}

The three components represent different levels of genetic variation:

a: variation between populations;
b: variation among individuals within populations; and
c: variation between the two alleles within an individual.

14.3 One-SNP function

wc_one_snp <- function(g1, g2, snp_name = NA_character_) {
  n1 <- sum(!is.na(g1))
  n2 <- sum(!is.na(g2))

  if (n1 < 2 || n2 < 2) {
    return(
      data.frame(
        SNP = snp_name,
        n1 = n1,
        n2 = n2,
        p1 = NA_real_,
        p2 = NA_real_,
        h1 = NA_real_,
        h2 = NA_real_,
        a = NA_real_,
        b = NA_real_,
        c = NA_real_,
        theta = NA_real_,
        f = NA_real_,
        F = NA_real_
      )
    )
  }

  p1 <- mean(g1, na.rm = TRUE) / 2
  p2 <- mean(g2, na.rm = TRUE) / 2

  h1 <- mean(g1 == 1, na.rm = TRUE)
  h2 <- mean(g2 == 1, na.rm = TRUE)

  comp <- wc_components(
    n = c(n1, n2),
    p = c(p1, p2),
    h = c(h1, h2)
  )

  a <- comp$a
  b <- comp$b
  c <- comp$c

  denominator_abc <- a + b + c
  denominator_bc <- b + c

  theta <- ifelse(
    denominator_abc == 0,
    NA_real_,
    a / denominator_abc
  )

  f_value <- ifelse(
    denominator_bc == 0,
    NA_real_,
    b / denominator_bc
  )

  F_value <- ifelse(
    denominator_abc == 0,
    NA_real_,
    (a + b) / denominator_abc
  )

  data.frame(
    SNP = snp_name,
    n1 = n1,
    n2 = n2,
    p1 = p1,
    p2 = p2,
    h1 = h1,
    h2 = h2,
    a = a,
    b = b,
    c = c,
    theta = theta,
    f = f_value,
    F = F_value
  )
}

For each SNP, n1 and n2 are the numbers of observed individuals, p1 and p2 are the allele frequencies, and h1 and h2 are the observed heterozygote frequencies.

15 Test the calculation on one SNP

g1_test <- M_wc[
  population_wc == "Mo44_PHW65",
  1
]

g2_test <- M_wc[
  population_wc == "PHN11_PHW65",
  1
]

test_result <- wc_one_snp(
  g1 = g1_test,
  g2 = g2_test,
  snp_name = colnames(M_wc)[1]
)

test_result

##          SNP n1 n2        p1        p2 h1 h2           a         b c
## 1 S1_1202006 60 70 0.4166667 0.5285714  0  0 0.002389124 0.2502046 0
##         theta f F
## 1 0.009458364 1 1

For the observed analysis, the first SNP was S1_1202006. Its allele frequencies were approximately 0.4167 and 0.5286, and its \(\theta\) estimate was approximately 0.00946.

16 Analyse every retained SNP

results_list <- vector(
  mode = "list",
  length = ncol(M_wc)
)

for (j in seq_len(ncol(M_wc))) {
  g1 <- M_wc[
    population_wc == "Mo44_PHW65",
    j
  ]

  g2 <- M_wc[
    population_wc == "PHN11_PHW65",
    j
  ]

  results_list[[j]] <- wc_one_snp(
    g1 = g1,
    g2 = g2,
    snp_name = colnames(M_wc)[j]
  )
}

snp_results <- do.call(
  rbind,
  results_list
)

head(snp_results)

##          SNP n1 n2        p1        p2 h1 h2             a         b c
## 1 S1_1202006 60 70 0.4166667 0.5285714  0  0  0.0023891236 0.2502046 0
## 2 S1_2871916 60 70 0.5666667 0.0000000  0  0  0.1587741815 0.1151042 0
## 3 S1_2888789 60 70 0.5666667 0.5285714  0  0 -0.0031647268 0.2513765 0
## 4 S1_3275604 60 70 0.4333333 0.4714286  0  0 -0.0031647268 0.2513765 0
## 5 S1_3443770 60 70 1.0000000 0.5000000  0  0  0.1228841146 0.1367188 0
## 6 S1_3776359 60 70 0.5833333 0.4857143  0  0  0.0008873432 0.2505394 0
##          theta f F
## 1  0.009458364 1 1
## 2  0.579725205 1 1
## 3 -0.012750108 1 1
## 4 -0.012750108 1 1
## 5  0.473354232 1 1
## 6  0.003529231 1 1

dim(snp_results)

## [1] 1237   13

summary(snp_results$theta)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.      NAs 
## -0.08896  0.04204  0.34417  0.29589  0.47858  1.00000        2

The analysis produced 1,237 SNP rows and 13 columns. Two SNPs had undefined \(\theta\) values because their relevant denominator was zero. A small number of negative per-SNP estimates can occur through sampling variation and do not automatically indicate a programming error.

17 Calculate the final multilocus estimates

The correct multilocus estimator is not the arithmetic mean of the per-SNP ratios. The variance components must be added across SNPs before the final ratios are calculated.

A_total <- sum(snp_results$a, na.rm = TRUE)
B_total <- sum(snp_results$b, na.rm = TRUE)
C_total <- sum(snp_results$c, na.rm = TRUE)

theta_total <- A_total / (
  A_total + B_total + C_total
)

f_total <- B_total / (
  B_total + C_total
)

F_total <- (
  A_total + B_total
) / (
  A_total + B_total + C_total
)

final_results <- data.frame(
  Statistic = c(
    "theta (FST)",
    "f (FIS)",
    "F (FIT)"
  ),
  Estimate = c(
    theta_total,
    f_total,
    F_total
  )
)

A_total

## [1] 87.41466

B_total

## [1] 175.1493

C_total

## [1] 0

final_results

##     Statistic  Estimate
## 1 theta (FST) 0.3329271
## 2     f (FIS) 1.0000000
## 3     F (FIT) 1.0000000

For the completed analysis:

\[ \sum a = 87.41466,\qquad \sum b = 175.1493,\qquad \sum c = 0. \]

Therefore,

\[ F_{ST} = \frac{87.41466}{87.41466+175.1493+0} = 0.3329271. \]

18 Graphs

par(mfrow = c(1, 2))

plot(
  snp_results$theta,
  type = "l",
  col = "darkblue",
  xlab = "SNP number",
  ylab = expression(theta~"("*F[ST]*")"),
  main = expression("Per-SNP "*theta~"("*F[ST]*")")
)

abline(
  h = theta_total,
  col = "red",
  lwd = 2,
  lty = 2
)

hist(
  snp_results$theta,
  breaks = 30,
  col = "lightblue",
  border = "white",
  xlab = expression(theta~"("*F[ST]*")"),
  main = expression("Distribution of per-SNP "*theta)
)

abline(
  v = theta_total,
  col = "red",
  lwd = 2,
  lty = 2
)

par(mfrow = c(1, 1))

The blue line and histogram describe the variation in \(\theta\) across SNPs. The red dashed line shows the final multilocus estimate.

19 Verify the absence of heterozygotes

sum(M_wc == 1, na.rm = TRUE)

## [1] 0

sum(snp_results$h1 > 0, na.rm = TRUE)

## [1] 0

sum(snp_results$h2 > 0, na.rm = TRUE)

## [1] 0

In the completed analysis, the heterozygote counts were zero. Consequently, \(C_{total}=0\). This is biologically plausible because the dataset represents highly inbred maize parental lines.

20 Interpretation of the final results

The final multilocus estimate was \(F_{ST}=0.3329271\). This indicates strong genetic differentiation between the sampled Mo44_PHW65 and PHN11_PHW65 populations. Approximately 33.3% of the estimated genetic variation was associated with differences between the two population groups.

The estimates \(F_{IS}=1\) and \(F_{IT}=1\) resulted from the absence of observed heterozygotes in the retained markers. Because \(c=\bar h/2\), zero observed heterozygosity gives \(c=0\). Therefore,

\[ F_{IS}=\frac{B}{B+C}=1 \]

and

\[ F_{IT}=\frac{A+B}{A+B+C}=1. \]

These values should be interpreted in the context of the highly inbred maize lines and the genotype coding verified from the VCF.

21 Save the results

results_dir <- file.path(data_dir, "results")
dir.create(results_dir, showWarnings = FALSE)

write.csv(
  snp_results,
  file.path(results_dir, "WC_per_SNP_results.csv"),
  row.names = FALSE
)

write.csv(
  final_results,
  file.path(results_dir, "WC_final_statistics.csv"),
  row.names = FALSE
)

saveRDS(
  list(
    M_wc = M_wc,
    population_wc = population_wc,
    snp_results = snp_results,
    final_results = final_results,
    selected_pop1 = selected_pop1,
    selected_pop2 = selected_pop2
  ),
  file.path(results_dir, "WC_complete_analysis.rds")
)

list.files(results_dir, full.names = TRUE)

## [1] "/home/arginine/data/results/WC_complete_analysis.rds"
## [2] "/home/arginine/data/results/WC_final_statistics.csv" 
## [3] "/home/arginine/data/results/WC_per_SNP_results.csv"

22 Final conclusion

The workflow successfully converted the original VCF data into an individual-by-SNP dosage matrix, retained biallelic and polymorphic markers, preserved missing observations for unbiased per-SNP sample-size calculations, and estimated the Weir-Cockerham variance components. After quality control, 1,237 SNPs from 130 individuals were analysed. The final multilocus estimate of \(F_{ST}=0.3329\) showed substantial genetic differentiation between the two selected maize populations. The estimates \(F_{IS}=1\) and \(F_{IT}=1\) reflected the absence of observed heterozygotes among the retained markers.

Weir and Cockerham Analysis of Two Maize Populations

Rhoda Ayanfeoluwa Ojo

21 September 2026