講義の内容
この講義では、QTL解析とGxE解析について、オオムギの多環境試験（Multi-environmental Trial: MET）データを例として、 実際の解析の流れを説明します。
利用するデータは、オオムギの遺伝・育種学的研究に長らく用いられてきたSteptoe×Morexの交配から得られた 150の倍加半数体（doubled haploid、DH）系統のデータです。
利用するデータは下記よりダウンロードできるが、Rのパッケージagridatにも含まれています。 本講義では、予めダンロードして、csv形式で保存しておいたデータを用います。
遺伝子型データ： http://www.genenetwork.org/genotypes/SXM.geno（2020.12.19現在、接続できない）
表現型データ： http://wheat.pw.usda.gov/ggpages/SxM/phenotypes.html
遺伝子型データには7つの染色体全体に分布する218のマーカーのデータが含まれます。 また、表現型データとしては8形質のデータ（αアミラーゼ、糖化力、出穂日、倒伏率、草丈、穀粒のタンパク質含量、収量）が含まれています。 本講義では収量のデータを例として用います。収量データには、16環境で計測された単位収量（t/ha）が含まれています。
必要なパッケージ
最初に、解析に必要となるパッケージを読み込みます。
hereは、階層構造をもつフォルダにあるデータファイルの指定に便利です。 qtlは、QTL解析のためのパッケージです。 tdyverseは、様々なユーテリティ機能が備わっているパッケージですが、ここでは、データフレーム（data.frame)構造をもつデータの操作・変形に用います。 fields, gplotsは、画像を用いた2次元データの描画に用いられます。
require(here)
require(qtl)
require(tidyverse)
require(fields)
require(gplots)疑似乱数の種のセット
以下の解析では、無作為並べ替え検定が何度も行われますが、この講義を受講しているときには、 皆が同じ結果が得られるほうが説明がしやすいので、最初に疑似乱数の種を特定の値に設定しておきます。
set.seed(12345)QTL解析
まずは、QTL解析用に準備しておいたcsv形式のデータファイルを読み込みます。 geno_bc.csvは、マーカー遺伝子型データ、yield.csvは、収量の表現型データです。
pheno.file <- "pheno/yield.csv" # ここを変えれば、いろいろな形質のQTL解析ができる。
cross <- read.cross(format = "csvs", genfile = "geno_bc.csv", phefile = pheno.file)## --Read the following data:
## 150 individuals
## 495 markers
## 17 phenotypes## Warning in summary.cross(cross): Some markers at the same position on chr 2; use
## jittermap().## --Cross type: bcデータを読み込むと何やら警告がでます。これは、第2染色体の同じ位置に2つのマーカーが座乗しているためです。
染色体上の位置が全く同じマーカー遺伝子座がある場合には、事前にどちらかを除いてもいいですが、 jittermap関数を用いて、同じ位置に座乗するマーカー間に小さな隙間をつくることもできます。ここでは、後者の方法を用います。
cross <- jittermap(cross)次に、plot関数を使って様々なグラフ（連鎖地図、欠測値、表現型データのヒストグラム）を表示してみる。
まずは、連鎖地図を表示します。
plotMap(cross)
つぎに、欠測しているマーカー遺伝子型データを表示します。
plotMissing(cross)
次に、表現型データのヒストグラムを描きます。
plotPheno(cross)
変なヒストグラムが描かれますが、これは、yieldデータの最初にあるIDのヒストグラムです。 2列目のデータ以降は、環境での収量データとなります。ここでは、環境１（2列目）の収量のデータを表示してみます。
env.id <- "Env01"
plotPheno(cross, pheno.col = env.id)
なお、授業でも説明したように、区間地図化法（interval mapping）では、区間の両側マーカー遺伝子型をもとに、 区間内のマーカーが座乗していない位置のマーカー遺伝子型の確率をもとめてQTLの位置と効果の推定に用いました。 qtlパッケージの関数を使ってQTL解析を行うためには、あらかじめこの確率を計算しておく必要があります。 ここでは、1 cMごとに遺伝子型の確率を計算しておきます。
interval <- 1 # 1 cMに設定。この値を変えると計算する密度を変えられる
cross <- calc.genoprob(cross, step = interval) # 遺伝子型の確率を計算
head(cross$geno[[1]]$prob[,,1])[, 1:10] # 第1染色体の最初の6サンプル・10座のAAである確率の表示## Hor5 MWG938 loc1 loc2 MWG835A
## [1,] 9.999993e-01 1.000000e+00 9.999670e-01 9.999869e-01 9.999999e-01
## [2,] 7.001991e-07 9.943122e-09 3.301754e-05 1.305333e-05 5.737951e-08
## [3,] 9.999993e-01 1.000000e+00 9.999670e-01 9.999869e-01 9.999999e-01
## [4,] 7.001991e-07 9.943122e-09 3.301754e-05 1.305333e-05 5.737951e-08
## [5,] 7.001991e-07 9.943122e-09 3.301754e-05 1.305334e-05 5.738018e-08
## [6,] 9.999993e-01 1.000000e+00 9.999670e-01 9.999869e-01 9.999999e-01
## loc3 loc4 loc5 loc6 MWG036A
## [1,] 0.9997121014 0.9995821704 0.9996520788 9.999219e-01 9.999998e-01
## [2,] 0.0002878986 0.0004178296 0.0003479212 7.814519e-05 1.965739e-07
## [3,] 0.9997121014 0.9995821704 0.9996520788 9.999219e-01 9.999998e-01
## [4,] 0.0002878986 0.0004178296 0.0003479212 7.814519e-05 1.965739e-07
## [5,] 0.0002879029 0.0004178387 0.0003479349 7.816369e-05 2.160309e-07
## [6,] 0.9997121014 0.9995821704 0.9996520788 9.999219e-01 9.999998e-01次に、1cMごとに計算された確率を図として表示してみます。ここでは、7つある染色体のうちの1番目を表示します。
chr.id <- 1 # ここを変えると1-7番染色体の情報を見ることができる。
image.plot(t(cross$geno[[chr.id]]$prob[,,1]), col = cm.colors(64), horizontal = F)
まずは、単純インターバルマッピング（simple interval mapping、SIM）を行ってみましょう。 ここでは、16環境で計測された収量データのうち、16番目（17列目、列名Env16）を解析してみます。
env.id <- "Env16" # この値を変えると他の環境についてQTL解析ができる
out.em <- scanone(cross, pheno.col = env.id, method = "em")
plot(out.em)
次に、コンポジットインターバルマッピング（composite interval mapping、CIM）を 行って解析を行ってみます。CIMでは、マーカー共変量の数（ここでは4に設定）と、解析の窓の大きさ （その窓の中には共変量マーカーを配置しない。ここでは20（cM）に設定した）を指定して解析を行う必要があります。
n.covar <- 4 # 共変量の数
window.size <- 20 # マーカーを配置しない窓のサイズ
outcim.em <- cim(cross, pheno.col = env.id, method = "em",
n.marcovar = n.covar, window = window.size)
plot(outcim.em)
QTLの存在が示唆されるLOD値のピークがあるが、どの値をしきい値としてQTLの有無を判断すればよいのでしょうか？ ここでは、無作為並べ替え検定（random permutation test）を行ってしきい値を決定する方法を用います。 表現型データを無作為に並び替えると、本来は検出されていたQTLも検出されなくなると考えられます。 それでも、確率的なゆらぎによりLOD値のピークがでてくることがあります。 ただし、無作為に並び替えているので、これらLOD値のピークは、偽のQTLとなります。 そこで、無作為に並び替えることで、LOD値の経験的な帰無分布（帰無仮説（QTLが無い、QTLの効果が0）が正しいときに従うLOD値の分布） をもとめ、この分布を用いて検定を行います。 ここでは時間の節約のために100回の繰り返しで計算していますが、実際には1000回くらいは繰り返してください。
opermcim.em <- cim(cross, pheno.col = env.id, method = "em",
n.marcovar = n.covar, window = window.size,
n.perm = 100)経験的な帰無分布が得られるので、その上側5%点のLOD値をしきい値とします。
summary(opermcim.em, alpha = 0.05)## LOD thresholds (100 permutations)
## [,1]
## 5% 3.76では、QTL解析の結果を描いて、そこにしきい値の水平線を引いてみましょう。
plot(outcim.em)
abline(h = summary(opermcim.em, alpha = 0.05), col = "orange")
オレンジ色の線がしきい値を表し。LOD値のピークがこれより上に出ている場合は、QTLが検出されたと考えます。
なお、以下のようにすると、この閾値のもとで有意なマーカーを列挙できます。
summary(outcim.em, perms = opermcim.em, alpha = 0.05)## chr pos lod
## Ppd-H1 2 35.5 14.8
## Dfr 3 54.4 21.7第2、第3染色体の35.5 cM, 54.4 cMの位置にQTLが検出されます。
このうち、第2染色体に検出されるPpd-H1は、日長反応に関わる遺伝子です。 https://science.sciencemag.org/content/310/5750/1031
第３染色体の遺伝子については、米国の研究者（METの一員だった）に問い合わせてみましたが、 おそらくまだクローニングはされておらず、遺伝子は分かっていないのではないかということでした。
なお、オオムギとコムギの収量関連のQTLや遺伝子については、こちらに総説があります。 https://link.springer.com/article/10.1007/s00122-017-2880-x
次に、検出されたQTLの効果の推定をします。
qtl <- summary(outcim.em, perms = opermcim.em, alpha = 0.05)
qm <- makeqtl(cross, chr = qtl$chr, pos = qtl$pos, what = "prob")
res <- fitqtl(cross, qtl = qm, get.ests = T, method = "hk", pheno.col = env.id)
summary(res)##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 150
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1 + Q2
##
## df SS MS LOD %var Pvalue(Chi2) Pvalue(F)
## Model 2 39.55545 19.777727 26.26003 55.34539 0 0
## Error 147 31.91473 0.217107
## Total 149 71.47019
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value Pvalue(Chi2) Pvalue(F)
## 2@35.5 1 13.93 11.80 19.50 64.18 0 3.21e-13 ***
## 3@54.4 1 22.90 17.62 32.04 105.46 0 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 3.22468 0.03825 84.299
## 2@35.5 -0.61158 0.07634 -8.011
## 3@54.4 -0.78679 0.07661 -10.269第2, 3染色体に検出されたQTLが説明する変動は、それぞれ、全変動の19.5%, 32.0%です。 また、3つのQTL全体で、全表現型変異の55.3%が説明されることが分かります。 さらに、これらQTLのBB遺伝子型のAA遺伝子型に対する効果は、それぞれ、-0.6, -0.8であることがわかります。
次に、マーカー遺伝子型と表現型の関係を図示してみます。 コードはすこしややこしいので、結果のほうに注目してもらえばよいかと思います。
mrk <- find.marker(cross, qtl$chr, qtl$pos)
mg <- NULL
for(i in 1:nrow(qtl)) {
mg <- cbind(mg, cross$geno[[qtl$chr[i]]]$data[,mrk[i]])
}
colnames(mg) <- rownames(mrk)
rownames(mg) <- cross$pheno$id
y <- cross$pheno[, env.id]
for(i in 1:ncol(mg)) {
boxplot(y ~ mg[, i], main = colnames(mg)[i])
}
cmg <- mg[,1]
for(i in 2:ncol(mg)) {
cmg <- paste(cmg, mg[,i], sep = "_")
}
boxplot(y ~ cmg, main = paste(colnames(mg), collapse = "_"))
2つのQTLで、BB遺伝子型（Morexの遺伝子型）になると、収量が大きく低下することが分かります。
では、16環境のそれぞれでCIMを用いたQTL解析を行ってみましょう。
for(i in 2:ncol(cross$pheno)) {
outcim.em <- cim(cross, pheno.col = i, method = "em",
n.marcovar = n.covar, window = window.size)
plot(outcim.em, main = colnames(cross$pheno)[i])
}
QTL解析の結果は、環境間で大きく違います。これは、GxE（遺伝子型と環境の交互作用）によるものです。 では、ここから、第8回講義で行ったいくつかの方法を用いて、GxEに注目して同データの解析を進めていきましょう。
遺伝子型と環境の交互作用パタンの把握
まずは、表現型データを読み込みます。これは、先ほどQTL解析のために読み込んだ表現型ファイルと同じものです。
phe <- read.csv(pheno.file, row.names = 1)
head(phe)## Env01 Env02 Env03 Env04 Env05 Env06 Env07 Env08 Env09 Env10 Env11 Env12
## 1 4.610 6.829 2.56 3.878 7.232 9.2048 3.8938 5.2779 4.798 6.580 5.226 6.095
## 2 5.159 5.785 2.07 4.554 6.326 7.3174 3.2735 5.5314 4.794 5.376 5.336 7.201
## 3 4.701 6.102 3.18 3.016 6.533 7.4983 2.7902 5.0852 4.946 6.072 5.870 6.014
## 4 5.297 6.022 2.93 2.419 7.469 6.8908 4.1540 5.4993 5.131 6.320 4.798 8.087
## 5 4.883 6.716 3.07 2.952 5.997 6.6689 2.5215 4.8926 5.299 6.915 5.929 7.980
## 6 4.681 5.692 2.65 2.968 6.525 8.6451 4.0002 5.5517 3.679 6.608 5.827 5.927
## Env13 Env14 Env15 Env16
## 1 6.287 4.389 5.070 2.499
## 2 6.008 4.414 5.490 2.810
## 3 4.507 5.849 5.768 2.480
## 4 6.717 6.492 5.362 3.566
## 5 7.874 4.932 6.389 3.639
## 6 6.279 4.671 5.095 2.723以降に行う解析にも関わるのですが、欠測のないバランスの良いデータかどうかを確認します。
sum(is.na(phe))## [1] 0欠測しているデータはなく、バランスの良いデータだと分かります。
欠速がない場合は、これも第8回の講義で説明したように、列平均・行平均をとることで、 遺伝子型や環境の平均を容易に求めることができます。ここでは、環境平均を計算します。
(ye <- apply(phe, 2, mean))## Env01 Env02 Env03 Env04 Env05 Env06 Env07 Env08
## 4.902973 5.742707 3.285333 3.555353 6.917340 7.469584 3.734297 5.151366
## Env09 Env10 Env11 Env12 Env13 Env14 Env15 Env16
## 4.934673 5.968773 5.016200 7.489093 5.910420 5.498527 5.853893 3.200193計算された環境平均を\(x\)軸にして、各遺伝子型の各環境における表現型値をプロットしてみましょう。 1つの遺伝子型のデータは、折れ線でつながったかたちで散布されます。
matplot(ye, t(phe), type = "l", lty = "dotted")
この図をみると、遺伝子型と環境の交互作用パタンをある程度把握することができます。 環境間で、遺伝子型の順位が大きく変化しているところがあり、「質的な」交互作用があることが分かります。
環境間の類似性の解析
第8回講義では、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)
距離行列のヒートマップや、クラスタ解析の結果の樹形図よりEnv07, Env08は他とは少し離れた環境であることが分かります。
次に、樹形図のならび順（左から右）に、QTL解析を行い、その結果を示します。
env <- tre$labels[tre$order]
for(i in 1:length(env)) {
outcim.em <- cim(cross, pheno.col = env[i], method = "em",
n.marcovar = n.covar, window = window.size)
plot(outcim.em, main = env[i])
}
Env04, 07, 08, 13では、第3染色体上のQTLが検出されていないことが分かります。 16環境におけるQTL解析の結果を、ヒートマップを用いて、まとめて図示してみます。
env <- tre$labels[tre$order]
lod <- NULL
for(i in 1:length(env)) {
outcim.em <- cim(cross, pheno.col = env[i], method = "em",
n.marcovar = n.covar, window = window.size)
lod <- cbind(lod, outcim.em$lod)
}
colnames(lod) <- env
rownames(lod) <- paste0(outcim.em$chr, ".", round(outcim.em$pos))
heatmap.2(t(lod), trace = "row", dendrogram = "none", Rowv = F,
Colv = F, hline = 3, tracecol = "black", linecol = "cyan",
key = F, keysize = 0.2)
第8回講義（p30）で紹介したようにQTLと環境にも交互作用があり、効果の大きなQTLでも検出されたり、されなかったりすることが分かります。第2染色体のQTLは、先に述べたように日長反応の遺伝子と考えられますが、これは、開花のタイミングが変化することにより、収量に影響を与えていると考えられます。
安定性解析（Joing Regression Analysis, a.k.a. Finlay-Wilkinsonモデル）
第8回講義では、GxEがあるときの対応(Bernardo 2002)のうち「利用する」は、 統計モデルを用いてGxEのパターンを明らかにして、GxEパターンをうまく利用することで、 特定の環境においてパフォーマンスが最良となる遺伝子型を選び出す方法です。
そうした方法について、実際に解析を行いながら見ていきましょう。
まず、この後行う安定性解析のために、データを整形します。 複数環境が複数列にわたっていた横長のフォーマットから、 全ての環境を一つの列にまとめた縦長のフォーマットにします。 しかし、どの環境のデータか分かるように、環境のID（横長フォーマットの列名）を表す列も準備しています。
phe.long <- phe %>% tibble::rowid_to_column("geno") %>%
mutate(geno = as.factor(geno)) %>%
pivot_longer(-geno, names_to = "env", values_to = "yield")
head(phe.long)## # A tibble: 6 x 3
## geno env yield
## <fct> <chr> <dbl>
## 1 1 Env01 4.61
## 2 1 Env02 6.83
## 3 1 Env03 2.56
## 4 1 Env04 3.88
## 5 1 Env05 7.23
## 6 1 Env06 9.20次に、各環境における全遺伝子型の平均を求めます。これを環境の指標とするのが、Joint regression analysis (JRA)です。 この環境指標（各遺伝子型における全遺伝子型平均）を、\(x\)という変数名で、縦長フォーマットにしたデータに新たな列として加えます。
ye <- apply(phe, 2, mean)
phe2 <- data.frame(phe.long, x = (ye[phe.long$env]))
head(phe2)## geno env yield x
## 1 1 Env01 4.6100 4.902973
## 2 1 Env02 6.8290 5.742707
## 3 1 Env03 2.5600 3.285333
## 4 1 Env04 3.8780 3.555353
## 5 1 Env05 7.2320 6.917340
## 6 1 Env06 9.2048 7.469584次に、準備したデータに次のモデルを当てはめて、モデルパラメータの推定を行います。
\[ y_{ij} = \mu + g_i + b_i x_j + e_{ij} \]
model <- lm(yield ~ geno * x, data = phe2)得られた回帰パラメータをもとに、求めたいパラメータ（各品種の環境指数に対する回帰係数）を計算します。 また、切片も品種ごとに異なるので、この情報も抽出します。 最後に抽出されたパラメータをもとに、各品種の各環境における表現型を環境指数に回帰した際の回帰直線を描きます。 コードは少し複雑なので、結果に注目して下さい。
b <- c(coef(model)[151], coef(model)[152:300] + coef(model)[151])
hist(b)
a <- c(coef(model)[1], coef(model)[2:150] + coef(model)[1])
matplot(ye, t(phe), type = "p", pch = 1, col = "gray")
col <- rep("blue", length(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より小さいとき 良環境では平均より悪く、悪環境では平均より良い
となります。こうした特徴づけの結果に基づき、対象とする環境（環境指数の大小）に適した遺伝子型を選ぶことができます。
Finlay-Wilkinsonモデルに基づくQTL解析
では、次に、品種の環境応答の特徴づけにも用いられる\(b_i\)について、QTL解析を行ってみましょう。 上の計算で得られた\(b\)を、表現型データcross$phenoに新たな変数として付け加えます。
cross$pheno$b <- b\(b_i\)について、CIMを用いたQTL解析を行ってみます。
outcim.em <- cim(cross, pheno.col = "b", method = "em",
n.marcovar = n.covar, window = window.size)
opermcim.em <- cim(cross, pheno.col = "b", method = "em",
n.marcovar = n.covar, window = window.size, n.perm = 100)
plot(outcim.em, main = "b_i of Joint Regression Analysis")
abline(h = summary(opermcim.em, alpha = 0.05))
(qtl <- summary(outcim.em, perms = opermcim.em, alpha = 0.05))## chr pos lod
## c2.loc23 2 23.0 4.20
## snp_0728 3 57.1 7.35
## snp_0803 6 70.2 5.48第2、3、6染色体に\(b_i\)に関するQTLが検出されます（無作為並び替え検定の回数が少ないので、 検出の結果が人によって、解析する回によって、異なっている可能性が高いです）。
次に、検出されたQTLの効果の推定をします。
qtl <- summary(outcim.em, perms = opermcim.em, alpha = 0.05)
qm <- makeqtl(cross, chr = qtl$chr, pos = qtl$pos, what = "prob")
res <- fitqtl(cross, qtl = qm, get.ests = T, method = "hk", pheno.col = "b")
summary(res)##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 150
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1 + Q2 + Q3
##
## df SS MS LOD %var Pvalue(Chi2) Pvalue(F)
## Model 3 1.258258 0.41941947 12.26777 31.38331 3.293588e-12 6.314504e-12
## Error 146 2.751065 0.01884291
## Total 149 4.009324
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value Pvalue(Chi2) Pvalue(F)
## 2@23.0 1 0.3297 3.687 8.224 17.50 0 4.94e-05 ***
## 3@57.1 1 0.5840 6.270 14.567 30.99 0 1.21e-07 ***
## 6@70.2 1 0.4543 4.978 11.331 24.11 0 2.41e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 1.00758 0.01127 89.413
## 2@23.0 0.09734 0.02327 4.183
## 3@57.1 -0.12642 0.02271 -5.567
## 6@70.2 0.11115 0.02264 4.9102, 3, 6番染色体に座乗するQTLは、BB型のAA型に対する効果は、0.10, -013, 0.11です。また、全てのQTLをあわせると\(b_i\)のもつ変動の31.4%が説明されることが分かります。
QTLに最寄りのマーカーの遺伝子型を調べ、それらの遺伝子型と\(b_i\)の関係を箱ひげ図で図示してみます。
mrk <- find.marker(cross, qtl$chr, qtl$pos)
mg <- NULL
for(i in 1:nrow(qtl)) {
mg <- cbind(mg, cross$geno[[qtl$chr[i]]]$data[,mrk[i]])
}
colnames(mg) <- rownames(mrk)
rownames(mg) <- cross$pheno$id
mg <- na.omit(mg)
y <- cross$pheno$b[-attr(mg, "na.action")]
for(i in 1:ncol(mg)) {
boxplot(y ~ mg[, i], main = colnames(mg)[i])
}
cmg <- mg[,1]
for(i in 2:ncol(mg)) {
cmg <- paste(cmg, mg[,i], sep = "_")
}
boxplot(y ~ cmg, main = paste(colnames(mg), collapse = "_"))
abline(h = 1)
例えば、\(b_i1\)を大きくするには、2, 3, 6番染色体のQTLの遺伝子型が、BB-AA-BBとするとよいことが分かります。QTLに近接するDNAマーカーで選抜することで、望ましいGxEのパターンをもった遺伝子型を選抜できると期待されます。
なお、回帰直線の切片\(a_i\)についても、QTL解析を行うことができます。
cross$pheno$a <- a
outcim.em <- cim(cross, pheno.col = "a", method = "em",
n.marcovar = n.covar, window = window.size)
opermcim.em <- cim(cross, pheno.col = "a", method = "em",
n.marcovar = n.covar, window = window.size, n.perm = 100)
summary(opermcim.em, alpha = 0.05)## LOD thresholds (100 permutations)
## [,1]
## 5% 3.54plot(outcim.em, main = "a_i of Joint Regression Analysis")
abline(h = summary(opermcim.em, alpha = 0.05))
summary(outcim.em, perms = opermcim.em, alpha = 0.05)## chr pos lod
## DAK595A 2 26.4 5.31第2染色体にQTLが検出されます。
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\)を求めます。この方法は、第8回講義（p12）で説明したとおりです。
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,
xlab = "IPC1", ylab = "IPC2")
text(ze[,1:2], col = "black", labels = substr(rownames(ze), 4, 5))
abline(h = 0)
abline(v = 0)
環境を番号で散布しています。 原点からみて、ある環境と同じ向きにある遺伝子型（○で表されている）は、その環境において収量が増加するようなGxEをもちます。
次に、遺伝子型の平均収量と交互作用の第1主成分（IPCA1）の散布図を描いてみます。
plot(yg, zg[,1], ylim = range(c(zg[,1], ze[,1])), col = col,
xlab = "Genotypic and environmental mean", ylab = "IPCA1")
text(ye, ze[,1], col = "black", labels = substr(rownames(ze), 4, 5))
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.3492697(ad <- -1 / max(ze[,1]))## [1] -0.5747335この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番になれる遺伝子型です。
AMMIモデルのQTL解析
最後に、遺伝子型の全環境における平均、および、交互作用の第1, 2主成分について、QTL解析を行ってみます。
まず、全環境の平均についてQTL解析を行います。
cross$pheno$yg <- yg
outcim.em <- cim(cross, pheno.col = "yg", method = "em",
n.marcovar = n.covar, window = window.size)
opermcim.em <- cim(cross, pheno.col = "yg", method = "em",
n.marcovar = n.covar, window = window.size, n.perm = 100)
plot(outcim.em, main = "yg")
abline(h = summary(opermcim.em, alpha = 0.05))
summary(outcim.em, perms = opermcim.em, alpha = 0.05)## chr pos lod
## ABC325 3 55.7 17.6次に、検出されたQTLの効果の推定をします。
qtl <- summary(outcim.em, perms = opermcim.em, alpha = 0.05)
qm <- makeqtl(cross, chr = qtl$chr, pos = qtl$pos, what = "prob")
res <- fitqtl(cross, qtl = qm, get.ests = T, method = "hk", pheno.col = "b")
summary(res)##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 150
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1
##
## df SS MS LOD %var Pvalue(Chi2) Pvalue(F)
## Model 1 0.5041846 0.50418457 4.37744 12.5753 7.126957e-06 8.492139e-06
## Error 148 3.5051391 0.02368337
## Total 149 4.0093236
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 1.00543 0.01262 79.667
## 3@55.7 -0.11646 0.02524 -4.6143番染色体に座乗するQTLは、BB型のAA型に対する効果は、-011です。次に、QTLに最寄りのマーカーの遺伝子型を調べ、それらの遺伝子型と\(b_i\)の関係を箱ひげ図で図示してみます。
mrk <- find.marker(cross, qtl$chr, qtl$pos)
mg <- NULL
for(i in 1:nrow(qtl)) {
mg <- cbind(mg, cross$geno[[qtl$chr[i]]]$data[,mrk[i]])
}
colnames(mg) <- rownames(mrk)
rownames(mg) <- cross$pheno$id
y <- cross$pheno$yg
for(i in 1:ncol(mg)) {
boxplot(y ~ mg[, i], main = colnames(mg)[i])
}
if(ncol(mg) > 1) {
cmg <- mg[,1]
for(i in 2:ncol(mg)) {
cmg <- paste(cmg, mg[,i], sep = "_")
}
boxplot(y ~ cmg, main = paste(colnames(mg), collapse = "_"))
}次に、交互作用の第1主成分についてもQTL解析します。
cross$pheno$ipca1 <- zg[,1]
outcim.em <- cim(cross, pheno.col = "ipca1", method = "em",
n.marcovar = n.covar, window = window.size)
opermcim.em <- cim(cross, pheno.col = "ipca1", method = "em",
n.marcovar = n.covar, window = window.size, n.perm = 100)
plot(outcim.em, main = "ICPA1")
abline(h = summary(opermcim.em, alpha = 0.05))
summary(outcim.em, perms = opermcim.em, alpha = 0.05)## chr pos lod
## ABG398 3 53.7 14.11
## c7.loc110 7 110.0 4.37次に、検出されたQTLの効果の推定をします。
qtl <- summary(outcim.em, perms = opermcim.em, alpha = 0.05)
qm <- makeqtl(cross, chr = qtl$chr, pos = qtl$pos, what = "prob")
res <- fitqtl(cross, qtl = qm, get.ests = T, method = "hk", pheno.col = "b")
summary(res)##
## fitqtl summary
##
## Method: Haley-Knott regression
## Model: normal phenotype
## Number of observations : 150
##
## Full model result
## ----------------------------------
## Model formula: y ~ Q1 + Q2
##
## df SS MS LOD %var Pvalue(Chi2) Pvalue(F)
## Model 2 0.4096808 0.20484041 3.510878 10.2182 0.0003084055 0.0003625265
## Error 147 3.5996428 0.02448737
## Total 149 4.0093236
##
##
## Drop one QTL at a time ANOVA table:
## ----------------------------------
## df Type III SS LOD %var F value Pvalue(Chi2) Pvalue(F)
## 3@53.7 1 0.40527 3.47504 10.1082 16.5503 0.000 7.7e-05 ***
## 7@110.0 1 0.00803 0.07258 0.2003 0.3279 0.563 0.568
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Estimated effects:
## -----------------
## est SE t
## Intercept 1.00580 0.01286 78.227
## 3@53.7 -0.10463 0.02572 -4.068
## 7@110.0 0.01500 0.02619 0.573第3染色体のQTLは、BB型のAA型に対する効果は、-0.10であることが分かります。
mrk <- find.marker(cross, qtl$chr, qtl$pos)
mg <- NULL
for(i in 1:nrow(qtl)) {
mg <- cbind(mg, cross$geno[[qtl$chr[i]]]$data[,mrk[i]])
}
colnames(mg) <- rownames(mrk)
rownames(mg) <- cross$pheno$id
mg <- na.omit(mg)
y <- cross$pheno$ipca1[-attr(mg, "na.action")]
for(i in 1:ncol(mg)) {
boxplot(y ~ mg[, i], main = colnames(mg)[i])
}
if(ncol(mg) > 1) {
cmg <- mg[,1]
for(i in 2:ncol(mg)) {
cmg <- paste(cmg, mg[,i], sep = "_")
}
boxplot(y ~ cmg, main = paste(colnames(mg), collapse = "_"))
}
abline(h = 0)
この結果から、IPCA1を大きくするには、第3, 7染色体のQTLの遺伝子型をAA-BBとし、逆にIPCA2をおおきくするには、BB-AAとすれば良いことが分かります。
AMMIモデルを用いることで、各遺伝子型の様々な環境におけるポテンシャルを推察できるだけでなく、 GxEパターンをうまく利用できるようにQTLの遺伝子型をデザインできます。
最後に、全ての結果をまとめてプロットしてみましょう。
lod <- NULL
for(i in 2:ncol(cross$pheno)) {
outcim.em <- cim(cross, pheno.col = i, method = "em",
n.marcovar = n.covar, window = window.size)
lod <- cbind(lod, outcim.em$lod)
}
colnames(lod) <- colnames(cross$pheno)[-1]
rownames(lod) <- paste0(outcim.em$chr, ".", round(outcim.em$pos))
heatmap.2(t(lod), trace = "row", dendrogram = "none", Rowv = F,
Colv = F, hline = 3, tracecol = "black", linecol = "cyan",
key = F, keysize = 0.2)