GWAS Analysis of Flowering Time and Plant Height in a Rice Diversity Panel

Population Structure, Heritability, and Marker-Trait Association

1 1. Introduction

This report runs a GWAS on the Zhao et al. (2011) rice panel: 413 Oryza sativa accessions, ~36,900 SNPs (44K array, MSU6), phenotyped for flowering time at Arkansas and plant height. Given known subpopulation structure (indica, aus, temperate/tropical japonica, aromatic, admixed), the workflow runs PCA, heritability estimation, and a Q+K association test in rrBLUP (Endelman, 2011), plus candidate-gene mapping.

2 2. Methodology

2.1 2.1 Data Import and Sample Matching

Germplasm metadata, the 44K genotype matrix, and the 34-trait phenotype table (ricediversity.org) got pulled in and matched by accession ID – only accessions present in all three were kept.

library(rrBLUP)
library(dplyr)
library(tidyr)
library(ggplot2)

germplasm <- read.csv("D:\\OneDrive - Universitas Islam Indonesia\\College\\R for NGS\\RiceDiversity.44K.germplasm.csv")

genotype <- read.csv("D:\\OneDrive - Universitas Islam Indonesia\\College\\R for NGS\\RiceDiversity.44K.MSU6.Genotypes.csv")

phenotype <- read.table(
  "http://www.ricediversity.org/data/sets/44kgwas/RiceDiversity_44K_Phenotypes_34traits_PLINK.txt",
  header = TRUE, sep = "\t", stringsAsFactors = FALSE, check.names = FALSE
)

idx_ft <- grep("^Flowering.*time.*Arkansas$", trimws(names(phenotype)), ignore.case = TRUE)
idx_ph <- grep("^Plant.*height$", trimws(names(phenotype)), ignore.case = TRUE)

stopifnot(
  "Column for 'Flowering time at Arkansas' not found" = length(idx_ft) == 1,
  "Column for 'Plant height' not found" = length(idx_ph) == 1
)

names(phenotype)[idx_ft] <- "Flowering.time.at.Arkansas"
names(phenotype)[idx_ph] <- "Plant.height"

2.2 2.2 Genotype Encoding and Cleaning

Since these lines are inbred (homozygous), each SNP got recoded biallelically: major allele = -1, minor allele = +1, missing calls (“0”) turned into NA; accessions were intersected across all three tables.

map <- genotype[, c("id", "chr", "position")]
geno_mat <- as.matrix(genotype[, -(1:3)])
rownames(geno_mat) <- genotype$id

convert_geno <- function(snp_row) {
  called <- snp_row[snp_row != "0"]
  alleles <- unique(called)
  if (length(alleles) == 0 || length(alleles) > 2) return(rep(NA, length(snp_row)))
  if (length(alleles) == 1) return(rep(NA, length(snp_row)))
  freq <- sort(table(called), decreasing = TRUE)
  major <- names(freq)[1]; minor <- names(freq)[2]
  ifelse(snp_row == "0", NA, ifelse(snp_row == major, -1, 1))
}

geno_num <- apply(geno_mat, 1, convert_geno)
rownames(geno_num) <- colnames(geno_mat)
colnames(geno_num) <- rownames(geno_mat)

acc_id <- gsub("NSFTV_", "", rownames(geno_num))
rownames(geno_num) <- acc_id

common_id <- Reduce(intersect, list(acc_id, as.character(germplasm$NSFTV.ID), as.character(phenotype$NSFTVID)))

geno_num  <- geno_num[common_id, ]
germ_sub  <- germplasm[match(common_id, germplasm$NSFTV.ID), ]
pheno_sub <- phenotype[match(common_id, phenotype$NSFTVID), ]

stopifnot(all(rownames(geno_num) == germ_sub$NSFTV.ID))
stopifnot(all(rownames(geno_num) == pheno_sub$NSFTVID))

2.3 2.3 Kinship Matrix and Population Structure

Kinship matrix \(K\) came from A.mat() (min.MAF = 0.05, mean-imputed), following VanRaden (2008). Structure was checked via PCA (Jolliffe & Cadima, 2016), projecting accessions onto the components with the most genotypic variance.

2.4 2.4 Narrow-Sense Heritability

Heritability per trait came from the mixed model \(y=\mu+g+e\), solved by REML in mixed.solve() (Endelman, 2011), after dropping missing-phenotype accessions and subsetting \(K\).

2.5 2.5 GWAS Model and Multiple-Testing Correction

rrBLUP::GWAS() fits \(y=X\beta+Zg+e\). Three correction levels were compared – naive (none), Q (3 PCs), Q+K (3 PCs + kinship) – via genomic control \(\lambda=\text{median}(\chi^2)/0.456\) (Devlin & Roeder, 1999; Yu et al., 2006), significance called via Bonferroni \(\alpha=0.05/m\) (Dunn, 1961).

lambda_gc <- function(pvals) {
  chisq <- qchisq(1 - pvals, df = 1)
  median(chisq, na.rm = TRUE) / qchisq(0.5, df = 1)
}

build_gwas_geno <- function(dat, map) {
  data.frame(marker = colnames(dat$geno),
             chrom  = map$chr[match(colnames(dat$geno), map$id)],
             pos    = map$position[match(colnames(dat$geno), map$id)],
             t(dat$geno), check.names = FALSE)
}

build_gwas_pheno <- function(dat, trait_name) {
  df <- data.frame(NSFTVID = names(dat$y), trait = as.numeric(dat$y))
  names(df)[2] <- trait_name
  df
}

run_gwas <- function(dat, map, trait_name, n.PC = 3, use_K = TRUE) {
  geno_df  <- build_gwas_geno(dat, map)
  pheno_df <- build_gwas_pheno(dat, trait_name)
  GWAS(pheno = pheno_df, geno = geno_df, K = if (use_K) dat$K else NULL,
       n.PC = n.PC, min.MAF = 0.05, P3D = TRUE, n.core = 1, plot = FALSE)
}

2.6 2.6 Candidate Gene Mapping

Significant markers got checked against MSU/RAP-DB annotations within a $$100 kb window, roughly this panel’s LD block (Zhao et al., 2011), cross-referenced against known loci for these traits:

Trait	Candidate gene	Chromosome
Flowering time	Hd1	6
Flowering time	Hd3a	6
Flowering time	Ghd7	7
Flowering time	DTH8 / Ghd8	8
Flowering time	OsGI	1
Plant height	SD1	1
Plant height	OsGA20ox family	multiple

3 3. Results

3.1 3.1 Population Structure (PCA)

Amat_out <- A.mat(geno_num, min.MAF = 0.05, impute.method = "mean", return.imputed = TRUE)
K        <- Amat_out$Amat
geno_imp <- Amat_out$imputed

pca <- prcomp(geno_imp, scale. = TRUE)
summary(pca)$importance[, 1:5]

##                             PC1      PC2      PC3      PC4      PC5
## Standard deviation     78.48424 45.97621 31.78910 20.20176 17.78001
## Proportion of Variance  0.35579  0.12209  0.05837  0.02357  0.01826
## Cumulative Proportion   0.35579  0.47788  0.53625  0.55982  0.57808

pca_df <- data.frame(pca$x[, 1:3], Sub.population = germ_sub$Sub.population)

ggplot(pca_df, aes(PC1, PC2, color = Sub.population)) +
  geom_point(size = 2, alpha = 0.8) +
  theme_bw() +
  labs(title = "PCA of the Rice Diversity Panel", x = "PC1", y = "PC2")

The accessions don’t scatter randomly, they clump into pockets lining up with subpopulation. PC1 does most of the splitting (35.6% of variance): temperate/tropical japonica huddle together, indica and aus cluster separately further along the same axis, aromatic and admixed lines land in between (fitting their mixed ancestry). PC2 (12.2% more) mostly pushes aromatic further from japonica. Together the two axes cover almost half the genotypic variance (47.8%) that is a clear sign this isn’t one uniform population, hence PCs and/or kinship belong in the GWAS model.

3.2 3.2 Narrow-Sense Heritability

prep_trait <- function(trait_vec, ids, K, geno) {
  names(trait_vec) <- ids
  idx <- !is.na(trait_vec)
  list(y = trait_vec[idx], K = K[idx, idx], geno = geno[idx, , drop = FALSE])
}

dat_flower <- prep_trait(pheno_sub$Flowering.time.at.Arkansas, rownames(geno_num), K, geno_imp)
dat_height <- prep_trait(pheno_sub$Plant.height,               rownames(geno_num), K, geno_imp)

h2 <- function(dat) {
  fit <- mixed.solve(y = dat$y, K = dat$K)
  c(Vg = fit$Vu, Ve = fit$Ve, h2 = fit$Vu / (fit$Vu + fit$Ve))
}

h2_flower <- h2(dat_flower)
h2_height <- h2(dat_height)

h2_flower

##           Vg           Ve           h2 
## 2.142593e-07 1.594473e+02 1.343762e-09

h2_height

##           Vg           Ve           h2 
## 4.803884e-07 4.448929e+02 1.079784e-09

Here the numbers don’t match expectations. Both traits come back with heritability basically at zero (\(h^2\approx1\times10^{-9}\)), almost all variance in the residual (\(V_e\) = 159.5, 444.9), next to nothing genetic (\(V_g\approx2\)-\(5\times10^{-7}\)). Taken at face value that means these traits are almost purely environmental, flatly contradicting prior reports on this dataset (\(h^2\) usually 0.6-0.9).

3.3 3.3 GWAS Model Comparison and Selection

gwas_naive <- run_gwas(dat_flower, map, "Flowering.time.at.Arkansas", n.PC = 0, use_K = FALSE)

## [1] "GWAS for trait: Flowering.time.at.Arkansas"
## [1] "Variance components estimated. Testing markers."

gwas_Q     <- run_gwas(dat_flower, map, "Flowering.time.at.Arkansas", n.PC = 3, use_K = FALSE)

## [1] "GWAS for trait: Flowering.time.at.Arkansas"
## [1] "Variance components estimated. Testing markers."

gwas_QK    <- run_gwas(dat_flower, map, "Flowering.time.at.Arkansas", n.PC = 3, use_K = TRUE)

## [1] "GWAS for trait: Flowering.time.at.Arkansas"
## [1] "Variance components estimated. Testing markers."

data.frame(
  model  = c("naive", "Q (PCs only)", "Q+K"),
  lambda = c(lambda_gc(10^-gwas_naive$Flowering.time.at.Arkansas),
             lambda_gc(10^-gwas_Q$Flowering.time.at.Arkansas),
             lambda_gc(10^-gwas_QK$Flowering.time.at.Arkansas))
)

##          model    lambda
## 1        naive 0.7848352
## 2 Q (PCs only) 0.8031715
## 3          Q+K 0.8031715

For flowering time, \(\lambda\) climbs from 0.785 (naive) to 0.803 with three PCs, then stalls then Q and Q+K both land on 0.803, kinship adds nothing measurable. None crack 1, so there’s no structure-inflation; the test runs slightly conservative if anything. Basically the three PCs already do most of the correction work for this trait. Q+K stayed as the final model anyway, since kinship still guards against finer cryptic relatedness PCs alone can miss.

3.4 3.4 Final GWAS and Manhattan Plots

gwas_flower <- run_gwas(dat_flower, map, "Flowering.time.at.Arkansas", n.PC = 3, use_K = TRUE)

## [1] "GWAS for trait: Flowering.time.at.Arkansas"
## [1] "Variance components estimated. Testing markers."

gwas_height <- run_gwas(dat_height, map, "Plant.height",               n.PC = 3, use_K = TRUE)

## [1] "GWAS for trait: Plant.height"
## [1] "Variance components estimated. Testing markers."

manhattan_plot <- function(gwas_result, trait, title) {
  df <- gwas_result
  names(df)[names(df) == trait] <- "score"
  data_cum <- df %>% group_by(chrom) %>% summarise(max_pos = max(pos), .groups = "drop") %>%
    mutate(pos_add = lag(cumsum(max_pos), default = 0)) %>% select(chrom, pos_add)
  df <- df %>% left_join(data_cum, by = "chrom") %>% mutate(pos_cum = pos + pos_add)
  axis_df <- df %>% group_by(chrom) %>% summarise(center = mean(pos_cum), .groups = "drop")
  bonferroni <- -log10(0.05 / nrow(df))
  ggplot(df, aes(x = pos_cum, y = score, color = factor(chrom %% 2))) +
    geom_point(alpha = 0.7, size = 1) +
    scale_color_manual(values = c("steelblue", "grey40")) +
    scale_x_continuous(labels = axis_df$chrom, breaks = axis_df$center) +
    geom_hline(yintercept = bonferroni, linetype = "dashed", color = "red") +
    labs(title = title, x = "Chromosome", y = expression(-log[10](italic(p)))) +
    theme_bw() + theme(legend.position = "none")
}

manhattan_plot(gwas_flower, "Flowering.time.at.Arkansas", "GWAS: Flowering Time at Arkansas")

manhattan_plot(gwas_height, "Plant.height", "GWAS: Plant Height")

The two Manhattan plots don’t agree at all. Flowering time has one point on chromosome 6 that breaks away and crosses the red Bonferroni line, everything else well below. Plant height shows nothing near that threshold that is just noise left to right. Only flowering time turns up a real hit here; plant height doesn’t, at least not at this sample size. The Flowering.time.at.Arkansas / Plant.height column is already the \(-\log_{10}(p)\) score per marker.

3.5 3.5 Significant SNPs and Candidate Gene

get_sig_snps <- function(gwas_result, trait) {
  bonferroni <- -log10(0.05 / nrow(gwas_result))
  gwas_result %>% rename(score = all_of(trait)) %>% filter(score > bonferroni) %>% arrange(desc(score))
}

sig_flower <- get_sig_snps(gwas_flower, "Flowering.time.at.Arkansas")
sig_height <- get_sig_snps(gwas_height, "Plant.height")
sig_flower

##              marker chrom     pos    score
## id6006123 id6006123     6 9652192 5.946863

sig_height

## [1] marker chrom  pos    score 
## <0 rows> (or 0-length row.names)

boxplot_top_snp <- function(dat, gwas_sig, marker_col = "marker") {
  if (nrow(gwas_sig) == 0) { message("No SNP passed Bonferroni."); return(invisible(NULL)) }
  top_marker <- gwas_sig[[marker_col]][1]
  geno_top   <- dat$geno[, top_marker]
  box_df <- data.frame(trait = dat$y,
                        genotype = factor(ifelse(geno_top < 0, "Major allele", "Minor allele"),
                                           levels = c("Major allele", "Minor allele")))
  ggplot(box_df, aes(x = genotype, y = trait, fill = genotype)) +
    geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.15, alpha = 0.5) +
    theme_bw() + labs(title = paste("Allele effect at", top_marker), x = "Genotype class", y = "Phenotype")
}

boxplot_top_snp(dat_flower, sig_flower)

annotate_window <- function(sig_snps, flank_bp = 100000) {
  sig_snps %>% mutate(
    window_start = pmax(0, pos - flank_bp), window_end = pos + flank_bp,
    rapdb_link = paste0("https://rapdb.dna.affrc.go.jp/viewer/gbrowse2/irgsp1/?name=Chr",
                         chrom, ":", window_start, "..", window_end))
}

annotate_window(sig_flower)

##              marker chrom     pos    score window_start window_end
## id6006123 id6006123     6 9652192 5.946863      9552192    9752192
##                                                                                 rapdb_link
## id6006123 https://rapdb.dna.affrc.go.jp/viewer/gbrowse2/irgsp1/?name=Chr6:9552192..9752192

That one flowering-time SNP is id6006123 (chr6, position 9,652,192 bp, score = 5.95). The boxplot backs it up major and minor allele accessions form two visibly different phenotype distributions, not one blob, as expected of a real marker. Nothing hit that bar for plant height. The $$100 kb window (chr6: 9,552,192-9,752,192) doesn’t directly overlap a known flowering-time gene, but Hd1 sits only ~215 kb away (chr6: ~9,335,337-9,337,606; Yano et al., 2000)

4 4. Conclusion

Population structure here is substantial (PC1+PC2 \(\approx\) 48% of genotypic variance), and Q+K comes out mildly conservative (\(\lambda<1\)) – one credible flowering-time locus near Hd1 on chromosome 6, nothing for plant height. Main caveat: heritability. Both traits show implausible near-zero \(h^2\), likely a REML boundary/scaling issue rather than real absence of additive variance, worth re-fitting before trusting GWAS power. The SNP hit looks solid.

5 5. Reference

Devlin, B., & Roeder, K. (1999). Genomic control for association studies. Biometrics, 55(4), 997-1004.

Dunn, O. J. (1961). Multiple comparisons among means. Journal of the American Statistical Association, 56(293), 52-64.

Endelman, J. B. (2011). Ridge regression and other kernels for genomic selection with R package rrBLUP. The Plant Genome, 4(3), 250-255.

Jolliffe, I. T., & Cadima, J. (2016). Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A, 374(2065), 20150202.

VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. Journal of Dairy Science, 91(11), 4414-4423.

Yano, M., Katayose, Y., Ashikari, M., et al. (2000). Hd1, a major photoperiod sensitivity quantitative trait locus in rice, is closely related to the Arabidopsis flowering time gene CONSTANS. The Plant Cell, 12(12), 2473-2483.

Yu, J., Pressoir, G., Briggs, W. H., et al. (2006). A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nature Genetics, 38(2), 203-208.

Zhao, K., Tung, C.-W., Eizenga, G. C., et al. (2011). Genome-wide association mapping reveals a rich genetic architecture of complex traits in Oryza sativa. Nature Communications, 2, 467.