1 Introduction

Population structure determines how genetic diversity is distributed among populations, among individuals within populations, and within individuals. This study quantified these hierarchical components of genetic variation in two maize populations using the Weir and Cockerham (1984) method-of-moments estimators. Genotypic data were obtained from the Genomes to Fields 2018 dataset. A reproducible sample comprising 60 individuals from the Mo44_PHW65 population and 70 individuals from the PHN11_PHW65 population was evaluated. From 2,000 randomly selected variants, quality control retained 1,237 biallelic, polymorphic SNPs. The multilocus estimates were \(F_{ST}=0.3329\), \(F_{IS}=1.0000\), and \(F_{IT}=1.0000\). The \(F_{ST}\) estimate indicates substantial differentiation between the two maize populations. The values of \(F_{IS}\) and \(F_{IT}\) reflect the absence of observed heterozygotes in the retained marker set, consistent with the highly inbred nature of the evaluated parental lines. Overall, genetic variation was distributed strongly between populations and among individuals within populations, whereas the estimated within-individual component was zero.

Genetic diversity in a structured population can be partitioned into three hierarchical components: variation among populations, variation among individuals within populations, and variation between alleles within individuals. The fixation indices estimated by the Weir and Cockerham (1984) method provide a formal framework for quantifying these components.

The variance components are:

  • \(a\): variation attributable to differences among populations;
  • \(b\): variation attributable to differences among individuals within populations; and
  • \(c\): variation attributable to differences between alleles within individuals.

The corresponding fixation indices are

\[F_{ST}=\theta=\frac{a}{a+b+c},\]

\[F_{IS}=f=\frac{b}{b+c},\]

and

\[F_{IT}=F=\frac{a+b}{a+b+c}.\]

Thus, \(F_{ST}\) describes genetic differentiation among populations, \(F_{IS}\) describes heterozygote deficiency within populations, and \(F_{IT}\) describes heterozygote deficiency relative to the total population.

2 Objective

The objective was to quantify genetic variation:

  1. among the Mo44_PHW65 and PHN11_PHW65 maize populations;
  2. among individuals within these populations; and
  3. within individuals.

This was accomplished by estimating the Weir-Cockerham variance components \(a\), \(b\), and \(c\) and their associated multilocus fixation indices.

3 Dataset and software

The Genomes to Fields 2018 VCF contained 312 samples, 1,353,525 variants across 10 chromosomes, and approximately 4.26% missing genotype calls. A sample decoder linked the VCF identifiers to their genetic lines. Genotypes 0/0, 0/1, 1/0, and 1/1 were coded as alternative-allele dosages 0, 1, 1, and 2, respectively. Missing genotypes remained NA and were excluded locus by locus.

library(vcfR)

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), file.exists(decoder_file))

4 Population identification

decoder <- read.table(decoder_file, header=FALSE, stringsAsFactors=FALSE)
names(decoder) <- c("VCF_ID", "Genetic_Line")

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")))

population_inventory <- as.data.frame(table(decoder$Population))
names(population_inventory) <- c("Population", "Number_of_individuals")
population_inventory
##    Population Number_of_individuals
## 1  Mo44_PHW65                   100
## 2       Other                     6
## 3 PHN11_PHW65                   142
## 4   PHW65_MoG                    64

The decoder contained 100 Mo44_PHW65 individuals, 142 PHN11_PHW65 individuals, 64 PHW65_MoG individuals, and 6 other entries. The first two populations provided the required samples of 60 and 70.

5 Reproducible sampling

The random seed ensured that the same individuals and variants could be recovered. Individuals and markers were sampled without replacement. The reduced VCF was saved so the complete 1.7 GB file would not have to be loaded repeatedly.

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, 60, replace=FALSE)
selected_pop2 <- sample(pop2_ids, 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)
  vcf_ids <- colnames(geno@gt)[-1]
  stopifnot(nrow(geno@fix) == 1353525)
  stopifnot(ncol(geno@gt) - 1 == 312)
  stopifnot(sum(decoder$VCF_ID %in% vcf_ids) == 312)

  snp_idx <- sort(sample(seq_len(nrow(geno@fix)), 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
## *****        *****         *****

6 Marker filtering and genotype coding

Only biallelic variants were retained. Phased separators were standardized, and genotypes were converted into a numeric individual-by-SNP dosage matrix.

geno_bi <- geno_filtered[is.biallelic(geno_filtered), ]
gt <- extract.gt(geno_bi, element="GT")
gt_clean <- gsub("\\|", "/", gt)

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
M_wc <- t(M_wc)

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

Markers without an observed genotype were discarded. Monomorphic loci with allele frequency 0 or 1 were also removed. Missing genotypes were not imputed because the Weir-Cockerham estimator uses the observed sample size, allele frequency, and heterozygote frequency at each locus.

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]

quality_control_summary <- data.frame(
  Quantity=c("Individuals", "Initial sampled variants",
             "Biallelic variants", "Polymorphic SNPs analysed",
             "Missing genotype observations retained"),
  Value=c(nrow(M_wc), 2000, nrow(geno_bi@fix), ncol(M_wc), sum(is.na(M_wc)))
)
quality_control_summary
##                                 Quantity Value
## 1                            Individuals   130
## 2               Initial sampled variants  2000
## 3                     Biallelic variants  1891
## 4              Polymorphic SNPs analysed  1237
## 5 Missing genotype observations retained  5273

After quality control, 130 individuals and 1,237 informative polymorphic SNPs remained.

7 Estimation of variance components

For each locus, \(n_i\) denotes observed sample size, \(p_i\) denotes alternative-allele frequency, and \(h_i\) denotes observed heterozygote frequency in population \(i\). The sample-size correction is

\[n_c=\frac{r\bar n-\sum_i n_i^2/(r\bar n)}{r-1}.\]

The weighted allele frequency, among-population allele-frequency variance, and weighted heterozygosity are

\[\bar p=\frac{\sum_i n_i p_i}{\sum_i n_i},\qquad s^2=\frac{\sum_i n_i(p_i-\bar p)^2}{\bar n(r-1)},\qquad \bar h=\frac{\sum_i n_i h_i}{\sum_i n_i}.\]

wc_Nc <- function(n) {
  r <- length(n)
  n_bar <- mean(n)
  (r*n_bar - sum(n^2)/(r*n_bar))/(r-1)
}

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)
}

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, p2=NA,
      h1=NA, h2=NA, a=NA, b=NA, c=NA, theta=NA, f=NA, F=NA))
  }

  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(c(n1,n2), c(p1,p2), c(h1,h2))
  a <- comp$a; b <- comp$b; c <- comp$c
  den_abc <- a+b+c
  den_bc <- b+c

  data.frame(
    SNP=snp_name, n1=n1, n2=n2, p1=p1, p2=p2, h1=h1, h2=h2,
    a=a, b=b, c=c,
    theta=ifelse(den_abc == 0, NA, a/den_abc),
    f=ifelse(den_bc == 0, NA, b/den_bc),
    F=ifelse(den_abc == 0, NA, (a+b)/den_abc)
  )
}

8 Genome-wide and multilocus estimation

results_list <- vector("list", 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, g2, 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

Multilocus indices were calculated as ratios of summed variance components, not as arithmetic means of locus-specific ratios:

\[\hat\theta=\frac{\sum_j a_j}{\sum_j(a_j+b_j+c_j)},\quad \hat f=\frac{\sum_j b_j}{\sum_j(b_j+c_j)},\quad \hat F=\frac{\sum_j(a_j+b_j)}{\sum_j(a_j+b_j+c_j)}.\]

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)

component_results <- data.frame(
  Source=c("Among populations",
           "Among individuals within populations",
           "Within individuals"),
  Component=c("A = sum(a)", "B = sum(b)", "C = sum(c)"),
  Estimate=c(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)
)

component_results
##                                 Source  Component  Estimate
## 1                    Among populations A = sum(a)  87.41466
## 2 Among individuals within populations B = sum(b) 175.14927
## 3                   Within individuals C = sum(c)   0.00000
final_results
##     Statistic  Estimate
## 1 theta (FST) 0.3329271
## 2     f (FIS) 1.0000000
## 3     F (FIT) 1.0000000

9 Results

9.1 Locus-specific differentiation

The 1,237 locus-specific \(F_{ST}\) estimates ranged from approximately -0.089 to 1.000, with a median of 0.344. Two loci had undefined estimates because their total estimated variance was zero. Negative locus-specific method-of-moments estimates can arise through sampling variation and were retained when accumulating the variance components.

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
data.frame(
  Quantity=c("Undefined FST", "Undefined FIS", "Undefined FIT"),
  Number=c(sum(is.na(snp_results$theta)),
           sum(is.na(snp_results$f)),
           sum(is.na(snp_results$F)))
)
##        Quantity Number
## 1 Undefined FST      2
## 2 Undefined FIS      3
## 3 Undefined FIT      2

9.2 Partition of genetic variation

Source of variation Component Estimate
Among populations \(A=\sum a_j\) 87.41466
Among individuals within populations \(B=\sum b_j\) 175.1493
Within individuals \(C=\sum c_j\) 0.0000

The nonzero \(A\) component demonstrates appreciable allele-frequency variation between populations. The larger \(B\) component shows that considerable genetic variation also occurred among individuals within populations. The within-individual component was zero because no heterozygous genotype was observed.

9.3 Multilocus fixation indices

Fixation index Estimate Interpretation
\(F_{ST}\) 0.3329271 Differentiation among populations
\(F_{IS}\) 1.0000000 Heterozygote deficiency within populations
\(F_{IT}\) 1.0000000 Heterozygote deficiency relative to the total population

The final among-population differentiation was

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

Thus, approximately 33.3% of the estimated total genetic variation was attributable to differentiation between the sampled populations. The remaining estimated variation was principally associated with differences among individuals within populations.

9.4 Graphical representation

par(mfrow=c(1,2))
plot(snp_results$theta, type="l", col="darkblue",
     xlab="SNP number", ylab=expression(theta~"("*F[ST]*")"),
     main=expression("Locus-specific "*theta))
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 locus-specific "*theta))
abline(v=theta_total, col="red", lwd=2, lty=2)

par(mfrow=c(1,1))

The locus-specific estimates varied across the genome, indicating that differentiation was not uniform across markers. The red dashed line represents the multilocus estimate.

9.5 Within-individual variation

heterozygosity_check <- data.frame(
  Quantity=c("Observed dosage-1 genotypes",
             "SNPs with heterozygotes in Mo44_PHW65",
             "SNPs with heterozygotes in PHN11_PHW65"),
  Number=c(sum(M_wc == 1, na.rm=TRUE),
           sum(snp_results$h1 > 0, na.rm=TRUE),
           sum(snp_results$h2 > 0, na.rm=TRUE))
)
heterozygosity_check
##                                 Quantity Number
## 1            Observed dosage-1 genotypes      0
## 2  SNPs with heterozygotes in Mo44_PHW65      0
## 3 SNPs with heterozygotes in PHN11_PHW65      0

The absence of dosage-1 genotypes confirmed that observed heterozygosity was zero. Consequently, \(c=\bar h/2=0\) and both \(F_{IS}\) and \(F_{IT}\) equalled 1. These values indicate complete observed homozygosity in the marker set and are biologically consistent with highly inbred parental material.

10 Interpretation

The multilocus \(F_{ST}\) estimate of 0.3329 revealed pronounced genetic differentiation between Mo44_PHW65 and PHN11_PHW65. Their allele-frequency profiles were therefore substantially different, and population membership accounted for an important component of total genetic variation.

The component associated with individuals within populations was larger than the among-population component. Thus, strong population differentiation coexisted with appreciable diversity among individuals belonging to the same population. A high \(F_{ST}\) does not imply genetic uniformity within populations; variation was evident at both hierarchical levels.

The within-individual component was zero because no heterozygotes were observed. Accordingly, \(F_{IS}\) and \(F_{IT}\) reached their upper boundary. These estimates describe the analysed inbred material and marker subset and should not be generalized without reference to the biological nature of the samples.

Missing genotypes were handled using observed locus-specific sample sizes rather than mean-dosage imputation. This preserved the biological distinction among homozygous and heterozygous genotypes. The final multilocus estimates were calculated from summed variance components, maintaining the weighting inherent in the Weir-Cockerham method.

11 Conclusion

The Weir-Cockerham analysis partitioned genetic variation at three hierarchical levels. The among-population component was substantial, yielding \(F_{ST}=0.3329\) and demonstrating strong differentiation between Mo44_PHW65 and PHN11_PHW65. Genetic variation was also evident among individuals within populations, as indicated by the magnitude of the \(B\) component. The within-individual component was zero because the retained markers contained no observed heterozygotes; therefore, \(F_{IS}=F_{IT}=1\). The analysed sample was consequently characterized by marked population structure, meaningful variation among individuals within populations, and complete observed homozygosity within individuals.

12 References

Lima, D. C., Castro Aviles, A., Alpers, R. T., et al. (2023). 2018-2019 field seasons of the Maize Genomes to Fields G x E project. BMC Genomic Data, 24, 29.

Weir, B. S., and Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution, 38(6), 1358-1370.