生物科学専門実習IV 集団構造解析、ゲノムワイド関連解析

スクリプトを実行する前に、こちらのファイルをダウンロード、解凍して１箇所のフォルダ（ディレクトリ）にまとめておく。

ライブラリのインストール

一度インストールすれば、２回目以降はこの部分は実行不要。

install.packages("maps") # 地図描画
install.packages("ggplot2") # plot用ライブラリ
install.packages("rrBLUP") # GWAS
install.packages("qqman") # GWAS
install.packages("dplyr") # データフレームの整形など
install.packages("reshape2") #ヒートマップなどに便利

ライブラリの読み込み

以下は、Rを起動のたびに実行します。

library(maps) 
library(rrBLUP) 
library(qqman) 
library(dplyr) 
library(ggplot2)
library(reshape2)

表現型・SNPデータのロード

setwd("****") # ファイルを置いたディレクトリへのパスを***に書く
load("genot_trim.Rdata") # SNP データ 
load("SNP_INFO_250k_new.Rdata") # SNP のゲノム上の位置情報
load("Atwell_pheno199.Rdata") # 表現型データ
accession_list_ <- read.csv("call_method_75_info_Trim199.csv", header=T, fill = T, na.strings = c("", "NA") ) # 各系統の位置情報等
accession_list_ <- accession_list_[order(accession_list_$ecotype_id),] # リストの順番を揃える

表現型・SNPデータの整理

head(pheno199_df) # pheno199_df に表現型が107種類入っている。表現型の情報はファイル「study_1_phenotypes.csv」及び論文（Atwell et al. 2010 Nature）を参照。
pheno__ <- pheno199_df[,2] # 2列目の"Diameter Field" という表現型を選択（ここで入れる数字を変えて表現型を変える）
names(pheno__) <- pheno199_df$accession_id
pheno__ <- pheno__[!is.na(pheno__)] # NAデータを除く
pheno_ <- pheno__[names(pheno__) %in% rownames(genot_trim)] # SNPデータのない系統の表現型は除く

gt <- genot_trim[rownames(genot_trim) %in% names(pheno_),] # SNPデータの数を絞る（表現型データに揃える）

# SNPデータと表現型データの順番を揃える(order関数)
gt <- gt[order(as.numeric(rownames(gt))),]
pheno <- pheno_[order(as.numeric(names(pheno_)))]

# accession listのサイズ も表現型データに揃える
accession_list <- accession_list_[accession_list_$ecotype_id %in% names(pheno),]

gt[,1] # SNPデータを眺めてみよう。0がCol-0型、1が non-Col-0型。

head(accession_list) # 系統リストを眺めてみる。緯度経度、地域の情報

head(SNP_INFO) # SNP情報を眺めてみる。SNP のゲノム上の位置情報

GWAS を行う

# 最終的にSNPデータはgという行列にまとめる
g <- data.frame(SNP_INFO, t(gt))
rownames(g) <- 1:nrow(g)
colnames(g) <- c("marker", "chrom", "pos", rownames(gt))

# kinship matrix の計算（時間がかかるので気長に待つ）
Amat <- A.mat(gt, shrink = T)
colnames(Amat) <- rownames(Amat) <- rownames(gt)

# 最終的に表現型データはpという行列にまとめる
p <- data.frame(rownames(gt), pheno)
colnames(p) <- c("ecotype_id", "y")

# minor allele frequency の計算
p.al1 <- apply(gt + 1, 2, mean, na.rm = T) / 2
maf <- pmin(p.al1, 1 - p.al1)

gwa <- GWAS(p, g, K = Amat, min.MAF = 0.05, plot = F, n.core = 3)
mht <- data.frame(SNP = gwa$marker, CHR = gwa$chrom, BP = gwa$pos, P = 10^(-gwa$y))

# 中身を見てみる
head(mht) 

# gwa plot with Bonferroni and FDR
q.values <- p.adjust(mht$P, method = "BH")
manhattan(mht,suggestiveline=FALSE,genomewideline=-log10(0.05/length(mht[,1])),highlight=mht[q.values<0.05,]$SNP) 

# qq-plot
qq(mht$P)

有意なSNPを出力

# SNP with FDR < 0.05
mht[q.values<0.05,]$SNP

# SNP with -log10(p-value) > 5（少し緩い基準で出してみる）
mht[-log10(mht$P)>5,]

特定の領域のクローズアップを出力

# 5番染色体のみ出力
manhattan(subset(mht, CHR == 5))

# 5番染色体の一部を出力
manhattan(subset(mht, CHR == 5), 
          xlim = c(15070000, 15090000))

集団構造を考慮しないnaive GWAS

#時間がかかるので気長に待つ
naive_p <- numeric(0)
for(i in 1:length(g[,1])){
  naive_p[i] <- anova(lm(p$y ~ unlist(g[i,4:length(g[1,])])))$Pr[1]
}
naive_mht <- data.frame(SNP = gwa$marker, CHR = gwa$chrom, BP = gwa$pos, P = naive_p)
naive_mht <- na.omit(naive_mht)
manhattan(naive_mht,suggestiveline=FALSE,genomewideline=-log10(0.05/length(mht[,1])),highlight=mht[q.values<0.05,]$SNP) # マンハッタンプロット

# qq-plot (普通のGWASとどう違うだろうか)
qq(naive_mht$P)

集団構造：主成分分析（PCA）

pca_result <- prcomp(gt, scale. = TRUE) #少し待つ
# 結果の要約を表示
# "Proportion of Variance" は、その軸が全体の情報の何%を説明しているかを示す
summary(pca_result)$importance[, 1:5]
pca_scores <- as.data.frame(pca_result$x)

pca_data <- cbind(pca_scores, Population = accession_list$region)
ggplot(pca_data, aes(x = PC1, y = PC2, color = Population)) +
  geom_point(size = 3, alpha = 0.8) +          
  theme_minimal() +                           
  labs(title = "PCA",
       x = "PC1",
       y = "PC2") +
  theme(legend.position = "right")

集団構造：ペアワイズなFst

calculate_pairwise_fst_matrix <- function(geno_mat, pop_vec) {
  populations <- unique(pop_vec)
  n_pops <- length(populations)
  fst_matrix <- matrix(0, nrow = n_pops, ncol = n_pops)
  rownames(fst_matrix) <- populations
  colnames(fst_matrix) <- populations
  
  pairs <- combn(populations, 2)
  for (i in 1:ncol(pairs)) {
    pop1 <- pairs[1, i]
    pop2 <- pairs[2, i]
    
    # 2集団のデータを抽出
    sub_geno1 <- geno_mat[pop_vec == pop1, , drop = FALSE]
    sub_geno2 <- geno_mat[pop_vec == pop2, , drop = FALSE]
    
    # 各集団のアレル頻度 (p1, p2)
    p1 <- colMeans(sub_geno1, na.rm = TRUE)
    p2 <- colMeans(sub_geno2, na.rm = TRUE)
    
    # Hs (集団内ヘテロ接合度の平均)
    # Hs = (2*p1*q1 + 2*p2*q2) / 2
    Hs_per_snp <- ( (2 * p1 * (1 - p1)) + (2 * p2 * (1 - p2)) ) / 2
    
    # Ht (全体のアレル頻度に基づくヘテロ接合度)
    # p_bar = (p1 + p2) / 2
    p_bar <- (p1 + p2) / 2
    Ht_per_snp <- 2 * p_bar * (1 - p_bar)
    
    # --- Hudsonの推定値によるゲノム平均Fst ---
    # 分子(Ht - Hs)の総和 / 分母(Ht)の総和
    numerator_sum <- sum(Ht_per_snp - Hs_per_snp, na.rm = TRUE)
    denominator_sum <- sum(Ht_per_snp, na.rm = TRUE)
    
    if (denominator_sum == 0) {
      fst_val <- 0
    } else {
      fst_val <- numerator_sum / denominator_sum
    }
    # 行列に値を格納 (対称行列にする)
    fst_matrix[pop1, pop2] <- fst_val
    fst_matrix[pop2, pop1] <- fst_val
  }
  return(fst_matrix)
}

fst_mat <- calculate_pairwise_fst_matrix(gt, accession_list$region) # Fst計算関数を実行

# 行列をggplotで扱える形式(Long format)に変換
fst_melted <- melt(fst_mat)
colnames(fst_melted) <- c("Pop1", "Pop2", "Fst")

# ヒートマップ描画
ggplot(fst_melted, aes(x = Pop1, y = Pop2, fill = Fst)) +
  geom_tile(color = "white") + 
  
  geom_text(aes(label = sprintf("%.3f", Fst)), color = "black", size = 4) + # 数値をタイルの中に表示
  scale_fill_gradient(low = "white", high = "firebrick", name = "Fst") +   # 色の設定
  theme_minimal() + # デザイン調整
  labs(title = "Pairwise Genome-wide Fst Heatmap", x = NULL, y = NULL) +
  coord_fixed() + 
  
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 12, face = "bold"),
    axis.text.y = element_text(size = 12, face = "bold")
  )

集団構造と表現型の関係

pheno_region <- data.frame(pheno, Population = accession_list$region)

ggplot(pheno_region, aes(x = Population, y = pheno)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.5) +
    theme_light() +
  scale_fill_brewer(palette = "Pastel1") +
  theme(axis.text.x = element_text(size = 10,    
                                   angle = 45,  
                                   hjust = 1))

表現型の地図上の分布

pheno_map <- data.frame(pheno, Lon = accession_list$longitude, Lat = accession_list$latitude)
world_map <- map_data("world")
ggplot() +
  geom_polygon(data = world_map,
               aes(x = long, y = lat, group = group),
               fill = "gray90",
               color = "gray60",
               size = 0.2) +  
  geom_point(data = pheno_map, 
             aes(x = Lon, y = Lat, color = pheno), 
             size = 3,          
             alpha = 0.8) +     
  scale_color_viridis_c(option = "C", name = "Phenotype\nValue") +

  #coord_quickmap() +
  coord_quickmap(xlim = c(-10, 120), ylim = c(10, 80)) + # 地図上の表示箇所を絞りたい場合はここで

  theme_bw() + 
  labs(title = "Geographical Distribution of Phenotypes",
       x = "Longitude", 
       y = "Latitude") +
  
  theme(panel.grid = element_blank())

ゲノムワイドなFstのプロット

calculate_fst_each_snp <- function(geno_mat, pop_vec) {
 
  populations <- unique(pop_vec)
  n_total <- nrow(geno_mat)
  
  p_total <- colMeans(geno_mat, na.rm = TRUE)
  Ht <- 2 * p_total * (1 - p_total)
  
  Hs_sum <- numeric(ncol(geno_mat)) # 0で初期化
  
  for (p in populations) {
    # その集団のデータを抽出
    sub_geno <- geno_mat[pop_vec == p, , drop = FALSE]
    n_sub <- nrow(sub_geno)
    # その集団内の頻度とヘテロ接合度
    p_sub <- colMeans(sub_geno, na.rm = TRUE)
    H_sub <- 2 * p_sub * (1 - p_sub)
    # 重み付けして加算 (サンプルサイズ比率を掛ける)
    Hs_sum <- Hs_sum + H_sub * (n_sub / n_total)
  }
  Hs <- Hs_sum
  Fst <- (Ht - Hs) / Ht
  
  # Htが0（全個体が同じ値）の場合、NaNになるので0に置換
  Fst[is.nan(Fst)] <- 0
  # 計算上の誤差でマイナスになる場合は0にする
  Fst[Fst < 0] <- 0
  return(Fst)
}

fst_values <- calculate_fst_each_snp(gt, accession_list$region)

df_plot <- cbind(SNP_INFO, Fst = fst_values)
df_plot <- df_plot %>%
  group_by(Chr) %>%
  summarise(chr_len = max(Pos)) %>%
  mutate(tot = cumsum(chr_len) - chr_len) %>%
  select(-chr_len) %>%
  left_join(df_plot, ., by = c("Chr" = "Chr")) %>%
  arrange(Chr, Pos) %>%
  mutate(BPcum = Pos + tot) 

axisdf <- df_plot %>%  # X軸のラベルを表示する位置 (各染色体の中央) を計算
  group_by(Chr) %>% 
  summarize(center = (max(BPcum) + min(BPcum)) / 2)


ggplot(df_plot, aes(x = BPcum, y = Fst)) +
  geom_point(aes(color = as.factor(Chr)), alpha = 0.8, size = 1.5) +
  scale_color_manual(values = rep(c("grey60", "black"), 5)) +
  scale_x_continuous(label = axisdf$Chr, breaks = axisdf$center) +
  scale_y_continuous(limits = c(0.1, 1.05), expand = c(0, 0)) + # 0.1以下はデータ数が多いので省略
  theme_bw() +
  theme(
    legend.position = "none",           
    panel.border = element_blank(),     
    panel.grid.major.x = element_blank(), 
    panel.grid.minor.x = element_blank(),
    axis.line.y = element_line(color = "black"), 
    axis.line.x = element_line(color = "black")  
  ) +
   labs(title = "Genome-wide Distribution of Global Fst",
       x = "Chromosome",
       y = "Global Fst") +
  # Fstが高い上位1%のラインを引く
  geom_hline(yintercept = quantile(df_plot$Fst, 0.99, na.rm=TRUE), 
             linetype = "dashed", color = "red")
## 警告が出るが気にしない

生物科学専門実習IV 集団構造解析、ゲノムワイド関連解析

2025-12-17

ライブラリのインストール

ライブラリの読み込み

表現型・SNPデータのロード

表現型・SNPデータの整理

GWAS を行う

有意なSNPを出力

特定の領域のクローズアップを出力

集団構造を考慮しないnaive GWAS

集団構造：主成分分析（PCA）

集団構造：ペアワイズなFst

集団構造と表現型の関係

表現型の地図上の分布

ゲノムワイドなFstのプロット