バイオメトリックス第7回 ゲノミック選抜（GS）
本チュートリアルでは、Rを用いて GWASとGSのための予測モデル構築を行う方法について説明する。
例として利用するデータ
ここでは、Rを用いてGWASおよびGSモデル構築を行う方法を解説する。 なお、Zhaoら（2011; Nature Communications 2:467）がイネ遺伝資源のGWASに用いた表現型データと、 McCouchら（2016; Nature Communication 7:10532）が用いたHigh Density Rice Array (HDRA)を用いてジェノタイピングされたSNP遺伝子型データ（いずれも、Rice Diversityからダウンロードできる）を、例として用いる。 後者のデータは、70万SNPsのデータのうち、他のマーカーとの冗長性が無く、遺伝子型データの欠測率が低く（< 1%）、かつ、マイナーアリル頻度（minor allele frequency: MAF）が高い（> 5%）、38,769 SNPsのデータを用いる。なお、欠測値については、Beagle 5.1（Blowning et al. 2018）を用いて補完をしてある。
利用するRのパッケージ
SNPデータの読み込み・解析にgastonとpcaMethodsを、GSの予測モデル構築にはBGLRパッケージを用いる。いずれもRのパッケージインストーラで簡単にダウンロード・インストールが可能である。
require(gaston) # for reading a vcf file and performing GWAS
require(glmnet) # for genomic prediction
require(BGLR) # for genomic predictionデータの読み込み
まず、vcf 形式のマーカー遺伝子型データ(McCouch ら 2016) を読み込む。vcf 形式のデータの読み込みには、gastonパッケージに含まれる関数read.vcfを用いる。
# read a vcf file
bm <- read.vcf("HDRA-G6-4-RDP1-RDP2-NIAS.AGCT_MAF005MIS001NR_imputed.vcf.gz") ## ped stats and snps stats have been set.
## 'p' has been set.
## 'mu' and 'sigma' have been set.bm # show the summary of the file## A bed.matrix with 1568 individuals and 38769 markers.
## snps stats are set
## ped stats are setイネ遺伝資源アクセッション1568系統についての、38,769 SNPsの遺伝子型データが含まれている。
bmには、SNP遺伝子型が数値としてコードされたものが含まれている。
# show the first 10 rows and 10 columns of marker score data
as.matrix(bm[1:10, 1:10])## SNP-1.8562. SNP-1.18594. SNP-1.24921. SNP-1.25253.
## IRGC121316@c88dcbba.0 0 0 0 0
## IRGC121578@950ba1f2.0 0 0 0 0
## IRGC121491@0b7d031b.0 0 0 0 2
## IRGC121440@d02dfe4e.0 0 0 0 2
## IRGC121896@89688646.0 0 0 0 0
## IRGC121323@57c71a3b.0 0 0 0 0
## IRGC121251@392016d6.0 0 0 0 0
## IRGC121902@f1be6ab1.0 0 0 0 1
## IRGC124483@9a9187e6.0 0 0 0 0
## IRGC121314@6c2404e4.0 0 0 0 0
## SNP-1.29213. SNP-1.30477. SNP-1.31732. SNP-1.32666.
## IRGC121316@c88dcbba.0 0 2 0 2
## IRGC121578@950ba1f2.0 0 0 0 2
## IRGC121491@0b7d031b.0 0 0 2 2
## IRGC121440@d02dfe4e.0 0 0 2 2
## IRGC121896@89688646.0 0 0 0 2
## IRGC121323@57c71a3b.0 0 0 0 2
## IRGC121251@392016d6.0 0 0 0 2
## IRGC121902@f1be6ab1.0 0 0 1 2
## IRGC124483@9a9187e6.0 0 2 0 2
## IRGC121314@6c2404e4.0 0 0 0 2
## SNP-1.40457. SNP-1.75851.
## IRGC121316@c88dcbba.0 0 0
## IRGC121578@950ba1f2.0 0 0
## IRGC121491@0b7d031b.0 0 0
## IRGC121440@d02dfe4e.0 0 0
## IRGC121896@89688646.0 0 0
## IRGC121323@57c71a3b.0 0 0
## IRGC121251@392016d6.0 0 0
## IRGC121902@f1be6ab1.0 0 0
## IRGC124483@9a9187e6.0 0 0
## IRGC121314@6c2404e4.0 0 0SNP遺伝子型が、“0”と”2”でコードされていることが分かる。また、4列目の上から８番目のように、“1”でコードされているSNP遺伝子型もある。“0”, “1”, “2”は、それぞれ、A1/A1、A1/A2、A2/A2の遺伝子型を表し、“0”と”2”はホモ接合体、“1”はヘテロ接合体である。なお、イネは自殖であるため、ヘテロ接合の頻度は低い。
なお、A1やA2の対立遺伝子が、それぞれ、どの塩基に対応するのかは、SNP情報で確認できる。
# show the first 10 SNPs
head(bm@snps, 10)bm@snpsには、染色体番号（chr）、マーカー名（id）、物理位置（pos）、A1とA2の対立遺伝子（A1, A2)、3つのSNP遺伝子型（A1/A1, A1/A2, A2/A2）の観察数（N0, N1, N2）、非欠測率（callrate）、マイナーアリル頻度（maf）などの情報が含まれている。
系統/品種の属性情報と表現型データの読み込み
次に、おなじ遺伝資源アクセッションについて得られた表現型データを読み込む。
# read phenotypic data
pheno <- read.csv("RiceDiversityPheno.csv", row.names = 1)
str(pheno) # show the data summary## 'data.frame': 388 obs. of 36 variables:
## $ Flowering.time.at.Arkansas : num 75.1 89.5 94.5 87.5 89.1 ...
## $ Flowering.time.at.Faridpur : int 64 66 67 70 73 NA 74 75 55 NA ...
## $ Flowering.time.at.Aberdeen : int 81 83 93 108 101 158 122 81 74 NA ...
## $ FT.ratio.of.Arkansas.Aberdeen : num 0.927 1.078 1.016 0.81 0.882 ...
## $ FT.ratio.of.Faridpur.Aberdeen : num 0.79 0.795 0.72 0.648 0.723 ...
## $ Culm.habit : num 4 7.5 6 3.5 6 3 2.5 6.5 3 NA ...
## $ Leaf.pubescence : int 1 0 1 1 1 1 NA 1 NA NA ...
## $ Flag.leaf.length : num 28.4 39 27.7 30.4 36.9 ...
## $ Flag.leaf.width : num 1.283 1 1.517 0.892 1.75 ...
## $ Awn.presence : int 0 0 0 0 0 1 0 1 1 NA ...
## $ Panicle.number.per.plant : num 3.07 4.05 3.12 3.7 2.86 ...
## $ Plant.height : num 111 144 128 154 148 ...
## $ Panicle.length : num 20.5 26.8 23.5 29 30.9 ...
## $ Primary.panicle.branch.number : num 9.27 9.17 8.67 6.36 11.17 ...
## $ Seed.number.per.panicle : num 4.79 4.44 5.08 4.52 5.54 ...
## $ Florets.per.panicle : num 4.91 4.73 5.15 4.73 5.68 ...
## $ Panicle.fertility : num 0.879 0.75 0.935 0.813 0.871 0.907 0.781 0.887 0.537 NA ...
## $ Seed.length : num 8.06 7.71 8.24 9.71 7.12 ...
## $ Seed.width : num 3.69 2.95 2.93 2.38 3.28 ...
## $ Seed.volume : num 2.59 2.11 2.15 1.94 2.2 ...
## $ Seed.surface.area : num 3.91 3.66 3.71 3.7 3.66 ...
## $ Brown.rice.seed.length : num 5.79 5.6 6.12 7.09 5.14 ...
## $ Brown.rice.seed.width : num 3.11 2.57 2.47 2.04 2.8 ...
## $ Brown.rice.surface.area : num 3.51 3.29 3.33 3.3 3.28 ...
## $ Brown.rice.volume : num 7.74 5.11 5.13 4.09 5.56 ...
## $ Seed.length.width.ratio : num 2.19 2.61 2.81 4.08 2.17 ...
## $ Brown.rice.length.width.ratio : num 1.86 2.18 2.48 3.47 1.83 ...
## $ Seed.color : int 0 0 0 0 0 0 0 0 0 NA ...
## $ Pericarp.color : int 0 0 0 0 0 0 0 0 0 NA ...
## $ Straighthead.suseptability : num 4.83 7.83 NA 8.33 8.17 ...
## $ Blast.resistance : int 8 4 3 5 3 2 2 9 2 1 ...
## $ Amylose.content : num 15.6 23.3 23.1 19.3 23.2 ...
## $ Alkali.spreading.value : num 6.08 5.64 5.53 6.03 5.44 ...
## $ Protein.content : num 8.45 9.2 8 9.6 8.5 ...
## $ Year07Flowering.time.at.Arkansas: num 73 88 105.5 86.5 85.5 ...
## $ Year06Flowering.time.at.Arkansas: num 77.2 91 83.5 88.5 92.7 ...phenoには、36形質（一部は、形質と試験地あるいは年次との組合せ）のデータが含まれている。
ゲノミック選抜
表現型がある系統と無い系統
では、ここから、ゲノミックセレクション(genomic selection: GS)で用いられるゲノミック予測(genomic prediction: GP)のためのモデル構築を行ってみよう。ここでは、まず、glmnetパッケージを用いて、正則化線形回帰を行う。
なお、GSは、GWASでは有意なマーカーが検出されないような形質で特に有効である。ここでは、対象形質を、「個体あたりの稔性のある頴花数」として解析を行う。これを非説明変数（目的変数）\(y\)とする。
# prepare y
y <- pheno$Panicle.number.per.plant * pheno$Florets.per.panicle * pheno$Panicle.fertility
# Fertile florets per plant
names(y) <- rownames(pheno) # naming the elements of y with the rownames
y <- na.omit(y) # remove NA entriesまた、説明変数\(X\)については、SNP遺伝子型のスコアを用いる。bmデータからSNP遺伝子型スコア行列を取り出す。なお、SNP遺伝子型スコアについては、事前に平均0、標準偏差1に標準化しておく。
なお、RStudio Cloudでは、ここで実行が止まってしまうことが分かった。そこで、データ内のSNP数を減らすこととした。10個に1つを使うことで、SNP数を1/10にした。
# extract the marker score matrix from the original data, bm
standardize(bm) <- "mu_sigma"
bm <- bm[, (1:ncol(bm)) %% 10 == 0]
X <- as.matrix(bm)次に、説明変数\(X\)について、表現型データ\(y\)がある系統と、無い系統に分ける。最終的には、前者のデータから後者を予測する。後述するように、こうした予測を行うことによって、表現型データの無い遺伝資源についても、有望な系統のスクリーニングを行うことが可能となる。
X.wp <- X[names(y), ] # X with phenotypes (wp)
X.wop <- X[!(rownames(X) %in% names(y)), ] # X without phenotypes (wop)
c(dim(X.wp), dim(X.wop))## [1] 356 3876 1212 3876表現型データ\(y\)を持つ系統が356系統、持たない系統が1212系統ある。
まずは、glmnetパッケージを用いて正則化線形回帰を行うための準備を行う。glmnetでは、alphaというパラメータがあり、これを0に設定するとリッジ回帰、1に設定するとLASSO、0〜1の中間的な値に設定するとElasticNet(「伸び縮みする網」という意味。網目の細かいリッジと粗いLASSOの間を臨機応変に設定できるので、このような名前がついたと思われる)を用いてモデル構築が行われる。まずは、alphaを0にして、リッジ回帰を行う。
# develop a prediction model with ridge regression
alpha <- 0 # 0 corresponds to ridge regression表現型データのある356系統のSNP遺伝子型データX.wpから、多型の無いマーカーを除く。多型の有無のチェックには、分散の計算を用いている。分散が0でなければ多型があると判断できる。
# select only polymorphic markers
poly <- apply(X.wp, 2, var) > 0
X.wp <- X.wp[, poly]表現型データの無い1212系統のSNP遺伝子型データ（X.wop）についても、同じSNPのセットに揃えておく。
# select the same markers also for lines without phenotypic data
X.wop <- X.wop[, poly]glmnetパッケージのcv.glmnet関数を用いてリッジ回帰のモデルを構築する。cv.glmnet関数は、回帰係数の大きさに対して課すペナルティの大きさを調整するパラメータを、交差検証（cross-validation）によって最適化してくれる。関数名の最初のcvは交差検証を意味している。データの標準化はあらかじめ行ってあるので、ここでは行わない。なお、計算には時間がかかる。
# build a model with ridge regression
model.ridge <- cv.glmnet(X.wp, y, alpha = alpha, standardize = F)次に、得られたモデルをもとに、マーカー遺伝子型データから表現型データを計算する。ただし、ここでは、モデル構築に用いた系統のSNP遺伝子型データを入力にするため、モデルに基づく「予測」（prediction）ではなく、モデルの「あてはめ」（fitting）である。関数名predictに惑わされないよう注意する。
# predict y based on X (but actually this is not prediction)
y.fit <- predict(model.ridge, newx = X.wp, s = "lambda.min") あてはめられた\(y\)と観察値の\(y\)について比較する。両者の相関係数を計算する。
# compare fitted y with observed y
plot(y, y.fit)
abline(0, 1)
cor(y, y.fit)## lambda.min
## [1,] 0.88462なお、入力される説明変数の数が多い場合には、あてはまりの良いモデルを構築するのは容易である。したがって、あてはまりの良さは、モデルの予測能力を反映しない。モデルの予測能力を評価するためには、モデル構築に用いなかったデータに対する「予測」を行い、その精度を評価する必要がある。そこで、本当はモデル構築に用いることができるデータの一部を、モデル構築から除いておき、除いておいたデータを予測して、その予測値を観察された値と比較する。除いておいたデータは、「本当はモデル構築に用いることができるデータ」であり、表現型データ（\(y\)）があるデータである。したがって、予測値を計算して、「本当は観察されていた\(y\)の値」との比較が可能である。
\(n\)-fold交差検証とよばれる方法では、無作為にデータを\(n\)グループ分割し、そのうち1グループをモデル構築から除いておき、残りの\(n-1\)クラスで予測モデルを構築する。そして、除いておいた1グループについて、予測を行い、実測値との比較を行う。では、無作為に1〜\(n\)のIDを系統にふってみよう。ここでは、\(n=3\)とする。無作為にIDを並べ替えるにはsample関数を用いる。
# prepare ids (1 to n.fold) for all lines and sort them
n.fold <- 3
id <- 1:length(y) %% n.fold + 1
id## [1] 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2
## [38] 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## [75] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1
## [112] 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2
## [149] 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## [186] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1
## [223] 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2
## [260] 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## [297] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1
## [334] 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3cv.id <- sample(id)
cv.id## [1] 1 1 2 2 2 3 1 3 1 1 1 1 2 2 1 2 1 3 3 2 2 2 1 2 3 3 3 2 2 2 1 1 2 2 1 2 1
## [38] 3 1 2 3 1 2 2 1 3 2 1 1 1 3 1 3 1 2 1 2 1 2 2 3 3 3 2 1 1 3 1 3 3 1 2 3 1
## [75] 3 2 2 2 3 1 2 2 1 3 2 1 1 1 3 3 2 2 2 2 2 3 1 2 3 2 3 1 1 1 1 2 2 3 1 3 2
## [112] 1 1 3 2 2 2 1 3 1 2 2 1 1 2 3 3 2 1 3 2 3 2 1 3 3 3 3 3 3 2 3 2 3 2 3 2 3
## [149] 3 2 3 1 1 3 2 1 2 2 2 2 1 3 1 3 2 2 3 3 2 3 2 1 2 1 1 1 1 2 2 1 1 3 2 3 3
## [186] 3 3 2 3 3 3 1 1 2 1 1 1 3 1 3 2 2 3 3 3 3 1 3 2 2 1 3 2 3 3 1 1 3 3 2 3 1
## [223] 3 1 3 2 1 1 1 2 1 3 1 1 3 3 3 1 1 1 3 3 2 3 2 2 1 2 3 1 3 1 3 3 1 1 1 1 2
## [260] 2 3 2 2 1 1 1 1 2 2 2 3 1 3 3 3 2 2 1 3 2 1 1 3 2 1 2 3 1 3 2 1 3 1 3 3 1
## [297] 2 2 1 1 1 2 1 1 3 1 2 3 3 1 2 2 1 3 3 2 2 2 2 2 2 3 2 3 2 2 3 2 2 1 3 1 3
## [334] 3 2 1 2 2 3 3 3 1 1 3 2 1 3 3 2 3 1 3 1 3 2 1次に、IDが1番の系統を除いて、モデル構築用データを作成してみる。
# closs validate
y.train <- y[cv.id != 1]
X.wp.train <- X.wp[cv.id != 1, ]データを分割したことにより、多型がなくなるマーカーが出現する。そこでこれらを改めて取り除く。
#' remove monomorphic markers
poly <- apply(X.wp.train, 2, var) > 0
X.wp.train <- X.wp.train[, poly]訓練用データでモデルを構築し、1分割目の予測を行なう。
# build a model
model <- cv.glmnet(X.wp.train, y.train, alpha = alpha, standardize = F)
# predict y of the 1st fold
y.pred <- predict(model, newx = X.wp[cv.id == 1, poly], s = "lambda.min")1分割目について、観察値と予測値を比較する。観察値と予測値の相関係数を精度の指標として計算する。
# compare prediction with observation
plot(y[cv.id == 1], y.pred)
abline(0, 1)
cor(y[cv.id == 1], y.pred)## lambda.min
## [1,] 0.7127171実際には、各分割について、同様の計算を繰り返し、それらをまとめて精度評価する。そこで、以上の操作を、for文を用いて1〜3の分割について繰り返す。
# repeat this analysis for #1 - #5 folds
y.pred <- rep(NA, length(y))
for(i in 1:n.fold) {
print(i)
y.train <- y[id != i]
X.wp.train <- X.wp[id != i, ]
poly <- apply(X.wp.train, 2, var) > 0
X.wp.train <- X.wp.train[, poly]
model <- cv.glmnet(X.wp.train, y.train, alpha = alpha, standardize = F)
y.pred[id == i] <- predict(model, newx = X.wp[id == i, poly], s = "lambda.min")
}## [1] 1
## [1] 2
## [1] 3予測値と観察値を比較する。相関係数も精度の指標として計算する。
# comparison between prediction and observation
plot(y, y.pred)
abline(0, 1)
cor(y, y.pred)## [1] 0.6981772リッジ回帰で推定された各SNPマーカーの効果（回帰係数）を図示してみる。全てのデータを用いてモデルを構築し、それをもとにSNPマーカーの効果を図示する。効果（回帰係数）は、coef関数で取り出せる。
# estimate and draw marker effects
beta.ridge <- coef(model.ridge, s = "lambda.min")[-1] # model.ridge has been obtained above
plot(abs(beta.ridge), main = "Ridge", ylab = "Coefficients", xlab = "Position (bp)")
推定効果をプロットしてみると、特に目立って効果の大きなSNPsがあるわけではなく、多くのSNPがある程度の遺伝効果を持っていることが分かる。これは、GWASの際にも示唆されたように、特に効果の大きな遺伝子が存在しないことによると推察される。こうした形質で、観察値と予測値の相関が0.74程度であることは重要である。なぜなら、例えば、少数の遺伝子に支配される形質は従来のマーカー利用選抜（marker assisted selection: MAS）で選抜でき、GSの利用は必ずしも必要ではない。しかし、GSでは、特に目立った効果を示す遺伝子が存在しない（あるいは、見つからない）形質でも、ゲノムワイドマーカー多型をもとに選抜可能である。収量や品質などの経済的に重要な農業関連形質の多くは、多数の遺伝子に支配される所謂「複雑形質（complex traits）であり、GS利用の意義が大きい。
最後に、表現型データがある356系統から得られた予測モデルを、表現型データの無い1212系統に適用して、後者の表現型値を予測する。
# prediction of phenotypes of genotypes without phenotypic records
y.pred.wop <- predict(model.ridge, newx = X.wop, s = "lambda.min")
conf <- hist(c(y.pred, y.pred.wop), col = "#ff00ff40", border = "#ff00ff")
hist(y.pred.wop, col = "#0000ff40", border = "#0000ff", breaks = conf$breaks, add = T)
上図は、モデル作成に用いられた表現型データをもつ356系統の交差検証による予測値と、表現型データの無い1212系統の予測値のヒストグラムである。ここで、個体あたりの稔性頴花数が多いほうが良いとすると、このような予測から、そうした系統が1212系統にどのくらい含まれているかが分かる。
例えば、この値が16を超える系統が1212系統中70以上あることが分かる。
sum(y.pred.wop>16)## [1] 75最近、多数の遺伝資源についてゲノムデータが解析・公開されるようになってきている。例えば、イネでは遺伝資源3,000系統について全ゲノム配列のリシークエンスが行われ、そのデータが既に公開されている（The 3,000 Rice Genome Project 2014, Gigasicence 3:7; Sun et al. 2017, Nucl. Aci. Res. 45:597-605）。また、今回例として用いているデータも1,500系統以上について70万SNPsの遺伝子型データが解析され、公開されている。こうしたデータを活用することにより、表現型データのスクリーニングが難しいために、死蔵されることが多かった遺伝資源を簡易にスクリーニングできるようになる。したがって、GSを用いることで、有用な遺伝資源を迅速に発見し、育種に効率的に活用する道が開ける（Yu et al. 2016, Nature Plants 2: 16150; Tanaka and Iwata 2017, Theor. Appl. Genet. 131:93-105）。
次に、LASSOを用いた正則化線形回帰を行う。LASSOを用いる場合には、alphaを1にする。なお、alphaを1にする以外は、上で行った計算と全く同じである。
# buid a model with LASSO
alpha = 1 # do the same thing with alpha = 1
y.pred <- rep(NA, length(y))
for(i in 1:n.fold) {
print(i)
y.train <- y[id != i]
X.wp.train <- X.wp[id != i, ]
poly <- apply(X.wp.train, 2, var) > 0
X.wp.train <- X.wp.train[, poly]
model <- cv.glmnet(X.wp.train, y.train, alpha = alpha, standardize = F)
y.pred[id == i] <- predict(model, newx = X.wp[id == i, poly], s = "lambda.min")
}## [1] 1
## [1] 2
## [1] 3plot(y, y.pred)
abline(0,1)
cor(y, y.pred) ## [1] 0.6910018全データを用いて推定されたSNP効果（回帰係数）の絶対値を図示する。
# estimate and draw marker effects
model.lasso <- cv.glmnet(X.wp, y, alpha = alpha, standardize = F)
beta.lasso <- coef(model.lasso, s = "lambda.min")[-1]
plot(abs(beta.lasso), main = "LASSO", ylab = "Coefficients", xlab = "Position (bp)")
リッジ回帰に比べ、効果をもつSNPsの数が少ないことが分かる。また、効果をもつ場合は、リッジ回帰に比べ、効果の大きさが大きいことも分かる。しかし、LASSOを用いた場合でも、複数のSNPsに重み付けされており、この形質に関わる遺伝子の数が少なくないことが示唆される。
リッジ回帰とLASSOで推定されたSNPの重み（回帰係数）について散布図を描いて比較をしてみる。
plot(beta.ridge, beta.lasso)
LASSOでは、リッジ回帰でも効果が大きかったSNPの一部が大きな効果をもっていることがわかる。
最後に、先と同じように、表現型データが無い1212系統の予測値を計算する。
# prediction of phenotypes of lines without phenotypic data
y.pred.wop <- predict(model.lasso, newx = X.wop, s = "lambda.min")
conf <- hist(c(y.pred,y.pred.wop), col = "#ff00ff40", border = "#ff00ff")
hist(y.pred.wop, col = "#0000ff40", border = "#0000ff", breaks = conf$breaks, add = T)
sum(y.pred.wop>16)## [1] 48用いたモデル化手法とモデルが異なるため、リッジ回帰の場合と少し異なる結果が得られる。実際に適用する場合には、複数の手法の中から、それぞれの形質に合った手法を経験的に選んで用いるのが良い。どの形質にも万能な手法は無い（データ科学の世界では、これを「No free lunch定理」とよぶ）。実際のGSにはここでは紹介しなかった手法が多数用いられる。その中でも特にGBLUPとよばれる線形混合モデルとベイズ回帰が重要であるが、時間の都合上解説はしない。
以上でRを用いたGWASとGS予測モデル構築の方法に関する説明は終わりである。ここでは種子の縦横比と個体あたりの稔性のある頴花数の２形質を用いて解析を行ったが、他の形質についても解析を行ってみるといろいろな発見があるだろう。
付録：ゲノム血縁行列に基づくベイズ回帰モデル
BGLRパッケージのBGLR関数を用いて、ゲノム血縁行列に基づくベイズ回帰モデルを構築する。同手法は、線形混合モデルを基にしたGBLUPをベイズの枠組みで実行するものである。
なお、この手法では、系統間の遺伝共分散を規定するカーネルを指定する必要がある。ここでは、ゲノム関係行列を用いる。まずは、表現型データがある系統間の遺伝共分散を規定するゲノム関係行列を計算する。
K.wp <- tcrossprod(X.wp) / ncol(X.wp)つぎに、表現型データがない系統と、表現型データがある系統との間の遺伝共分散を規定する関係行列を計算する。なお、この行列は正方行列ではなく、表現型データのない系統の数を\(N\)、表現型データがある系統の数を\(M\)とすると\(N\times M\)行列となる。
K.wop2wp <- X.wop %*% t(X.wp) / ncol(X.wp)
dim(K.wop2wp)## [1] 1212 356BGLRでは、カーネル行列を含めたモデルの構造をETAという引数で指定する。具体的には、以下のように計算する。なお、計算にはマルコフ連鎖モンテカルロ（MCMC）アルゴリズムが用いられるが、その反復計算の長さを指定するパラメータも指定する。
eta <- list(list(K=K.wp, model="RKHS"))
nIter <- 5000
burnIn <- 1000
model.br <- BGLR(y, ETA = eta, nIter = nIter, burnIn = burnIn, verbose = F)あてはめられた\(y\)と観察値の\(y\)について比較をする。
# compare estimated y (g) with observed y
y.fit <- model.br$yHat
plot(y, y.fit)
abline(0, 1)
cor(y, y.fit)## [1] 0.9176515交差検証を行い、予測精度を評価する。
# repeat this analysis for the first to the fifth folds
y.pred <- rep(NA, length(y))
for(i in 1:n.fold) {
print(i)
y.train <- y
y.train[id == i] <- NA # masking the i-th fold as NA
model.cv <- BGLR(y.train, ETA = eta, nIter = nIter, burnIn = burnIn, verbose = F)
y.pred[id == i] <- model.cv$yHat[id == i]
}## [1] 1
## [1] 2
## [1] 3繰り返して予測をした後、観察値と予測値の散布図の描画と、両者の相関係数の計算を計算する。
# comparison between prediction and observation
plot(y, y.pred)
abline(0, 1)
cor(y, y.pred)## [1] 0.6940705表現型データがある356系統から得られた予測モデルを、表現型データの無い1212系統に適用して、後者の表現型値を予測する。
先ほど計算しておいた関係行列をもとに、予測値を計算できる。
y.wop.pred <- model.br$mu + K.wop2wp %*% solve(K.wp) %*% model.br$ETA[[1]]$uなお、以下のようにするとSNPに対する重みを取り出すこともできる。
beta.br <- t(X.wp / ncol(X.wp)) %*% solve(K.wp) %*% model.br$ETA[[1]]$u
plot(abs(beta.br))
実は、こうして計算されるSNPの重みは、先に説明したリッジ回帰におけるSNPの重み、すなわち、リッジ回帰の回帰係数と強い関係をもつ。
plot(beta.br, beta.ridge)
ほぼ、同じ値が得られていることがわかる。推定方法などは異なるが、互いに強い関係をもつ手法であることが、この結果からも分かる・
なお、こうして得られるSNPの重みから、新たなデータに対する予測値は、以下のようにも計算できる。
y.wop.pred.2 <- model.br$mu + X.wop %*% beta.br両者は、完全に一致する。
plot(y.wop.pred, y.wop.pred.2)
abline(0,1)
最後に、表現型データがある356系統のyの推定値と、表現型データがない1212系統の予測値を比較する。
# prediction of phenotypes of lines without phenotypic data
conf <- hist(c(y.pred,y.pred.wop), col = "#ff00ff40", border = "#ff00ff")
hist(y.pred.wop, col = "#0000ff40", border = "#0000ff", breaks = conf$breaks, add = T)
sum(y.pred.wop>16)## [1] 48ここでも、リッジ回帰の結果とよく似た結果が得られることが分かる。