資料の内容

　この資料では、GxE解析について、オオムギの多環境試験（Multi-environmental Trial: MET）データを例として、実際の解析の流れを説明します。

　利用するデータは、オオムギの遺伝・育種学的研究に長らく用いられてきたSteptoe×Morexの交配から得られた150の倍加半数体（doubled haploid、DH）系統のデータです。利用するデータはRのパッケージagridatに含まれています。

必要なパッケージ

　最初に、解析に必要となるパッケージを読み込みます。

hereは、階層構造をもつフォルダにあるデータファイルの指定に便利です。tdyverseは、様々なユーテリティ機能が備わっているパッケージですが、ここでは、データフレーム（data.frame)構造をもつデータの操作・変形に用います。agridatはデータが含まれているパッケージであり、fields、gplotsは、画像を用いた2次元データの描画に用いられます。

require(here)
require(tidyverse)
require(agridat)
require(fields)
require(gplots)

データの読み込み

agridatから、オオムギDH系統の多環境試験で収集された表現型データを読み込みます。

pheno <- steptoe.morex.pheno
head(pheno)
##       gen   env amylase diapow hddate lodging malt height protein  yield
## 1 Steptoe  MN92    22.7     46  149.5      NA 73.6   84.5    10.5 5.5315
## 2 Steptoe MTi92    30.1     72  178.0      10 76.5     NA    11.2 8.6403
## 3 Steptoe MTd92    26.7     78  165.0      15 74.5   75.5    13.4 5.8990
## 4 Steptoe  ID91    26.2     74  179.0      NA 74.1  111.0    12.1 8.6290
## 5 Steptoe  OR91    19.6     62  191.0      NA 71.5   90.0    11.7 5.3440
## 6 Steptoe  WA91    23.6     54  181.0      NA 73.8  112.0    10.0 6.2700

各変数（phenoの各列）の属性データを確認します。

str(pheno)
## 'data.frame':    2432 obs. of  10 variables:
##  $ gen    : Factor w/ 152 levels "Morex","SM1",..: 152 152 152 152 152 152 152 152 152 152 ...
##  $ env    : Factor w/ 16 levels "ID91","ID92",..: 4 8 6 1 11 15 7 5 9 10 ...
##  $ amylase: num  22.7 30.1 26.7 26.2 19.6 23.6 21 NA NA NA ...
##  $ diapow : int  46 72 78 74 62 54 62 NA NA NA ...
##  $ hddate : num  150 178 165 179 191 ...
##  $ lodging: int  NA 10 15 NA NA NA NA NA 0 50 ...
##  $ malt   : num  73.6 76.5 74.5 74.1 71.5 73.8 70.8 NA NA NA ...
##  $ height : num  84.5 NA 75.5 111 90 112 98 82 77.5 95 ...
##  $ protein: num  10.5 11.2 13.4 12.1 11.7 10 12 NA NA NA ...
##  $ yield  : num  5.53 8.64 5.9 8.63 5.34 ...

　表現型データとして、遺伝子型（系統、gen）、環境（試験地と年の組合せ、env）、および、8形質のデータ（αアミラーゼ amylase、糖化力 diapow、出穂日 hdddate、倒伏率 lodging、モルティング品質 malt、草丈 height、穀粒のタンパク質含量 protein、収量 yield）が含まれています。本資料では収量のデータを例として用います。収量データには、16環境で計測された単位収量（t/ha）が含まれています。

まず、遺伝子型と環境がどのような組合せで計測されたデータか確認してみます。

tab <- table(pheno$gen, pheno$env) # table関数で組合せを数え上げ
tab[1:10, 1:10] # 系統と環境の各組合せが1回ずつ計測されている（最初の10系統、10環境のみ）
##        
##         ID91 ID92 MA92 MN92 MTd91 MTd92 MTi91 MTi92 NY92 ON92
##   Morex    1    1    1    1     1     1     1     1    1    1
##   SM1      1    1    1    1     1     1     1     1    1    1
##   SM10     1    1    1    1     1     1     1     1    1    1
##   SM103    1    1    1    1     1     1     1     1    1    1
##   SM104    1    1    1    1     1     1     1     1    1    1
##   SM105    1    1    1    1     1     1     1     1    1    1
##   SM11     1    1    1    1     1     1     1     1    1    1
##   SM110    1    1    1    1     1     1     1     1    1    1
##   SM112    1    1    1    1     1     1     1     1    1    1
##   SM116    1    1    1    1     1     1     1     1    1    1
dim(tab) # データの次元を調べる（152系統、16環境）
## [1] 152  16
table(tab == 1) # 抜けがないか確認。すべて1である（TRUE）ことがわかる
## 
## TRUE 
## 2432

次に、収量データについて、512行（系統数）×16列（環境数）のかたちに整形する。

まずは、gen、env、yieldを抜き出す。

phe.org<- pheno %>% select(gen, env, yield)
head(phe.org)
##       gen   env  yield
## 1 Steptoe  MN92 5.5315
## 2 Steptoe MTi92 8.6403
## 3 Steptoe MTd92 5.8990
## 4 Steptoe  ID91 8.6290
## 5 Steptoe  OR91 5.3440
## 6 Steptoe  WA91 6.2700

次に、整形を行う。

phe <- phe.org %>% pivot_wider(names_from = "env", values_from = "yield") %>%
  column_to_rownames("gen")
dim(phe) # 152行、16列になっているか確認
## [1] 152  16
head(phe[,1:8]) # 最初の6行8列を表示
##           MN92  MTi92  MTd92  ID91  OR91  WA91 MTi91 MTd91
## Steptoe 5.5315 8.6403 5.8990 8.629 5.344 6.270 4.099 7.066
## Morex   4.1305 4.1760 3.7141 7.951 8.329 5.969 1.738 4.666
## SM1     4.6100 6.5800 5.2260 6.095 6.287 4.389 5.070 2.499
## SM2     5.1590 5.3760 5.3360 7.201 6.008 4.414 5.490 2.810
## SM3     4.7010 6.0720 5.8700 6.014 4.507 5.849 5.768 2.480
## SM4     5.2970 6.3200 4.7980 8.087 6.717 6.492 5.362 3.566

　以降に行う解析に関わるのですが、欠測のないバランスの良いデータかどうかを確認します。

sum(is.na(phe))
## [1] 0

欠測しているデータはなく、バランスの良いデータだと分かります。

　欠測がない場合は、今回の講義で説明したように、列平均・行平均をとることで、遺伝子型や環境の平均を容易に求めることができます。ここでは、環境平均を計算します。

(ye <- apply(phe, 2, mean))
##     MN92    MTi92    MTd92     ID91     OR91     WA91    MTi91    MTd91 
## 4.902026 5.974555 5.013441 7.499632 5.922605 5.506697 5.815270 3.235270 
##     NY92     ON92     WA92     MA92    SKo92    SKg92    SKk92     ID92 
## 5.743118 3.287153 3.556970 6.918011 7.486533 3.732182 5.149829 4.949679

　計算された環境平均を\(x\)軸にして、各遺伝子型の各環境における表現型値をプロットしてみましょう。 1つの遺伝子型のデータは、折れ線でつながったかたちで散布されます。

matplot(ye, t(phe), type = "l", lty = "dotted")

この図をみると、遺伝子型と環境の交互作用パタンをある程度把握することができます。環境間で、遺伝子型の順位が大きく変化しているところがあり、「質的な」交互作用があることが分かります。

環境間の類似性の解析

今回の講義では、GxEがあるときの対応(Bernardo 2002)として、

を説明しました。

このうち、「小さくする」は、ターゲットとしてる環境群をいくつかの小さくて、より均一な部分群に分類する、ということを説明しました。このための解析を行ってみましょう。

　まず、環境間で、全遺伝子型（系統）の表現型の相関係数\(r\)を計算します。また、相関係数\(r\)から、第8回講義で紹介した式を用いて、距離\(d\)に変換します。また、計算された距離行列を、ヒートマップとして表示してみます。

r <- cor(phe, use = "pair")
n <- nrow(r)
d <- 2 * (1 - 1/n) * (1 - r)
image.plot(x = 1:n, y = 1:n, z = d, col = heat.colors(64))

　計算された距離をもとに、クラスタ解析を行ってみます。第3回の講義で説明したように、クラスタ間の距離の定義を変えると、異なる解析結果が得られます。ここでは、クラスタ間の最遠の距離をクラスタ間の距離とする完全連結（complete linkage）法（最遠隣法（furthest neighbor method）ともよばれる）を用います。

clust.method <- "complete"
tre <- hclust(as.dist(d), method = clust.method)
plot(tre)

これにより、環境間の類似性を視覚的に確認できます。

安定性解析（Joing Regression Analysis, a.k.a. Finlay-Wilkinsonモデル）

今回の講義で説明したGxEがあるときの対応(Bernardo 2002)のうち「利用する」は、 統計モデルを用いてGxEのパターンを明らかにして、GxEパターンをうまく利用することで、 特定の環境においてパフォーマンスが最良となる遺伝子型を選び出す方法です。

そうした方法について、実際に解析を行いながら見ていきましょう。

　まず、この後行う安定性解析のために、ここでは、整形する前のデータ（phe.org）を使います。

head(phe.org)
##       gen   env  yield
## 1 Steptoe  MN92 5.5315
## 2 Steptoe MTi92 8.6403
## 3 Steptoe MTd92 5.8990
## 4 Steptoe  ID91 8.6290
## 5 Steptoe  OR91 5.3440
## 6 Steptoe  WA91 6.2700

　まず、各環境における全遺伝子型の平均を求めます。これを環境の指標とするのが、Joint regression analysis (JRA)です。この環境指標（各遺伝子型における全遺伝子型平均）を、\(x\)という変数名で、もとの縦長フォーマットのデータに新たな列として加えます。

ye <- apply(phe, 2, mean)
ye
##     MN92    MTi92    MTd92     ID91     OR91     WA91    MTi91    MTd91 
## 4.902026 5.974555 5.013441 7.499632 5.922605 5.506697 5.815270 3.235270 
##     NY92     ON92     WA92     MA92    SKo92    SKg92    SKk92     ID92 
## 5.743118 3.287153 3.556970 6.918011 7.486533 3.732182 5.149829 4.949679

yeの各要素には環境のIDがついているので、このIDを用いて並べ替えます。なお、phe.org$envのとおりに並べ替えたいのですが、これらを文字列とすることで、yeを意図するに並べ替えることができます。

phe2 <- data.frame(phe.org, x = ye[as.character(phe.org$env)])
head(phe2)
##       gen   env  yield        x
## 1 Steptoe  MN92 5.5315 4.902026
## 2 Steptoe MTi92 8.6403 5.974555
## 3 Steptoe MTd92 5.8990 5.013441
## 4 Steptoe  ID91 8.6290 7.499632
## 5 Steptoe  OR91 5.3440 5.922605
## 6 Steptoe  WA91 6.2700 5.506697

　次に、準備したデータに次のモデルを当てはめて、モデルパラメータの推定を行います。

\[ y_{ij} = \mu + g_i + b_i x_j + e_{ij} \]

model <- lm(yield ~ gen * x, data = phe2)

　得られた回帰パラメータをもとに、求めたいパラメータ（各品種の環境指数に対する回帰係数）を計算します。また、切片も品種ごとに異なるので、この情報も抽出します。最後に抽出されたパラメータをもとに、各品種の各環境における表現型を環境指数に回帰した際の回帰直線を描きます。

コードは少し複雑なので、結果に注目して下さい。

b <- c(coef(model)[153], coef(model)[154:304] + coef(model)[153])
names(b) <- levels(phe2$gen)
hist(b)

a <- c(coef(model)[1], coef(model)[2:152] + coef(model)[1])
matplot(ye, t(phe), type = "p", pch = 1, col = "gray")
col <- rep("blue", length(b))
names(col) <- names(b)
col[b < quantile(b, 3/4)] <- "green"
col[b < quantile(b, 1/4)] <- "red"
for(i in 1:length(b)) {
  abline(a[i], b[i], col = col[i])
}

　上の図では、回帰直線の傾きに応じて異なる色で回帰直線を描いています。 第8回の講義で説明したように、これら回帰係数の推定値によって（0に近い、1に近い、1より大きい） 各品種の環境応答の特徴づけができます。

\(b_i\)が - 0のとき 遺伝子型iのパフォーマンスは、いずれの環境でも変わらない - 1のとき 遺伝子型iのパフォーマンスは、他の遺伝子型の平均的パフォーマンスと変わらない - 1より大きいとき 良環境では平均より良く、悪環境では平均より悪い - 1より小さいとき 良環境では平均より悪く、悪環境では平均より良い

となります。こうした特徴づけの結果に基づき、対象とする環境（環境指数の大小）に適した遺伝子型を選ぶことができます。

AMMIモデル

　遺伝子型と環境のもつ交互作用を主成分分析と同様の方法でモデル化するのが、AMMI（(additive main effects and multiplicative interaction)モデルです。では、AMMIモデルを適用して、GxEのパターンの解析を行ってみましょう。

　まずは、遺伝子型平均と環境平均を求めます。行平均が前者、列平均が後者です。

Y <- as.matrix(phe)
yg <- apply(Y, 1, mean)
ye <- apply(Y, 2, mean)

　次に、遺伝子型と環境交互作用（およびノイズ、誤差）が含まれた行列\(\boldsymbol Z\)を求めます。この方法は、今回の講義で説明したとおりです。

yge <- sweep(sweep(Y, 2, ye), 1, yg) + mean(Y)

　この行列には、上述したようにGxEだけでなく、計測における誤差（ノイズ）も含まれています。この行列から、できるだけGxEの真のパターン（シグナル）を取り出したいと考えるとき、そのための一つの方法は、特異値分解による行列の主な変動の抽出です。GxEの行列\(\boldsymbol Z\)をもとに、遺伝子型あるいは環境に関する主成分分析を行うった場合と同じ結果が得られます。

　なお、AMMIのモデルは、以下のように表されます。 \[ y_{ij} = \mu + g_i + t_j + \sum_k^n \lambda_k u_{ki} v_{kj} + e_{ij} \] ここで、\(g_i\)は遺伝子型\(i\)の効果、\(t_j\)は環境\(j\)の効果、そして、\(\sum_k^n \lambda_k u_{ki} v_{kj}\)が、特異値分解で求められる部分で、\(\lambda_k\)が\(\boldsymbol Z\)の\(k\)番目の固有値、\(\boldsymbol u_k = [u_{k1}, u_{k2}, ..., u_{kI}]^T\)は\(\boldsymbol Z\)の\(k\)番目の左側固有ベクトル、\(\boldsymbol v_k = [v_{k1}, v_{k2}, ..., v_{kJ}]^T\)は、\(\boldsymbol Z\)の\(k\)番目の右側固有ベクトルです。

　まずは特異値分解\(\boldsymbol {Z = U\Lambda V}^T\)を行います。ここで、\(\boldsymbol U\)は、左側固有ベクトル\(\boldsymbol u_i\)を列として束ねた行列、\(\boldsymbol V\)は、左側固有ベクトル\(\boldsymbol v_j\)を列として束ねた行列である。また、\(\boldsymbol \Lambda = diag(\lambda_1, \lambda_2, \dots, \lambda_n)\)です。

　では、svd関数を用いて固有値分解を行ってみましょう。

sge <- svd(yge)

　次に、右側固有ベクトル、左側固有ベクトルを用いて、遺伝子型と環境のスコアを計算します。

zg <- sge$u %*% diag(sqrt(sge$d))
ze <- sge$v %*% diag(sqrt(sge$d))
rownames(zg) <- names(yg)
rownames(ze) <- names(ye)
colnames(zg) <- colnames(ze) <- paste0("IPC", 1:length(sge$d))

計算されたスコアの散布図を描いてみます。

plot(zg[,1:2], xlim = range(c(zg[,1], ze[,1])),
     ylim = range(c(zg[,2], ze[,2])), col = col[names(yg)], 
     xlab = "IPC1", ylab = "IPC2")
text(ze[,1:2], col = "black", labels = rownames(ze))
abline(h = 0)
abline(v = 0)

　環境を番号で散布しています。 原点からみて、ある環境と同じ向きにある遺伝子型（○で表されている）は、その環境において収量が増加するようなGxEをもちます。

　次に、遺伝子型の平均収量と交互作用の第1主成分（IPCA1）の散布図を描いてみます。

plot(yg, zg[,1], ylim = range(c(zg[,1], ze[,1])), col = col[names(yg)], 
     xlab = "Genotypic and environmental mean", ylab = "IPCA1")
text(ye, ze[,1], col = "black", labels = rownames(ze))
abline(h = 0)
abline(v = mean(yg))

　散布図における○の色は、さきほどのFinlay-Wilkinsonモデルにおいて\(b_i\)が1より大きい遺伝子型（青）、\(b_i\)が1に近い遺伝子型（緑）、\(b_i\)が1より小さい遺伝子型（赤）で塗り分けられています。 \(b_i\)の値は、遺伝子型平均とも関連があることがわかります。また、交互作用の第1主成分（IPCA1）は、\(b_i\)によって評価されているGxEのパターンとは異なるパターンを捉えているということも分かります。

　この図は、以下に説明するようで、より戦略的な利用が可能です。同図で表されているのは、\(g_i, t_j\)、および、\(\lambda_1 u_{1i} v_{1j}\)なので、AMMIの1次のモデル（AMMI-1とよばれる）です。すなわち、遺伝子型\(i\)の環境\(j\)における表現型の期待値\(\hat y_{ij}\)は、次のように表されます。 \[ \hat y_{ij} = \mu + g_i + t_j + \lambda_1 u_{1i} v_{1j} \]

　このとき、遺伝子型\(i\)と\(i'\)の環境\(j\)における優劣を考えます。まず、遺伝子型\(i\)と\(i'\)の期待との差は、 \[ \hat y_{ij} - \hat y_{i'j} = (g_i - g_{i'}) + \lambda_1 v_{1j} (u_{1i} - u_{1i'}) \] となります。

　この関係をもとにすると、ある遺伝子型が別の遺伝子型に対して勝るかどうか、ということを図で表すことできます。まず、以下の2つの数値（直線の傾き）を計算します。

(au <- -1 / min(ze[,1]))
## [1] 0.3319154
(ad <- -1 / max(ze[,1]))
## [1] -0.6287173

　この2つの傾きをもつ直線を、対象とする遺伝子型の座標点を通るように引くと、その遺伝子型が勝つことができる遺伝子型を図上でしめすことができます。

plot(yg, zg[,1], ylim = range(c(zg[,1], ze[,1])))
points(ye, ze[,1], pch = 17)
abline(h = 0)
abline(v = mean(yg))
ygm <- yg[which.max(yg)]
zgm <- zg[which.max(yg)]
lines(x = c(ygm, 0), y = c(zgm, zgm - ygm * au), lty = 3)
lines(x = c(ygm, 0), y = c(zgm, zgm - ygm * ad), lty = 3)
text(ygm, zgm, labels = "A")

　この図で、右端の点Aは、収量の全環境平均が最も大きかった遺伝子型（以下、遺伝子型Aとよぶ）に対応します。この点Aから、上述した直線を引き、注目している点の左側で、この2つの直線にはさまれている領域に散布されている遺伝子型は、モデル \[ \hat y_{ij} = \mu + g_i + t_j + \lambda_1 u_{1i} v_{1j} \]　 のもとでは、遺伝子型Aの収量を超えることができません。 　一方、その外側にある遺伝子型は、環境が変わると、GxEの効果も利用して、遺伝子型Aに勝つことができます。

　こうした、GxEのもとでの遺伝子型の優劣をより詳細に図示するには、上述した勝ち負けの関係を全ての遺伝子型間で確認し、環境が変われば1番になれる遺伝子型を探し出します。

　次のスクリプトは、勝ち負けを確認して、いずれかの遺伝子型に負けてしまう場合には、shadowというフラグをTRUEにしています。

bu <- zg[,1] - au * yg
bd <- zg[,1] - ad * yg
shadow <- rep(F, length(bu))
for(i in 1:length(yg)) {
  for(j in 1:length(bu)) {
     if(bu[i] > bu[j] & bd[i] < bd[j]) {
       shadow[i] <- T
     }
  }
}

　では、結果を図示してみます。

col <- rep("gray", length(yg))
col[!shadow] <- "blue"
plot(yg, zg[,1], ylim = range(c(zg[,1], ze[,1])), col = col)
points(ye, ze[,1], pch = 17)
abline(h = 0)
for(i in 1:length(yg)) {
  if(!shadow[i]) {
    ygm <- yg[i]
    zgm <- zg[i,1]
    lines(x = c(ygm, 0), y = c(zgm, zgm - ygm * au), lty = 3)
    lines(x = c(ygm, 0), y = c(zgm, zgm - ygm * ad), lty = 3)
  }
}

　青い点であらわされている遺伝子型は、どこかの環境でGxEを利用して1番になれる遺伝子型です。