RでGWASと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を、GWASにはgastonパッケージを、GSの予測モデル構築にはBGLRパッケージを用いる。いずれもRのパッケージインストーラで簡単にダウンロード・インストールが可能である。
require(gaston) # for reading a vcf file and performing GWAS
#require(pcaMethods) # for pca
require(glmnet) # for genomic prediction
require(BGLR) # for genomic predictionGWASの実行
データの読み込み
まず、vcf 形式のマーカー遺伝子型データ(McCouch ら 2016) を読み込む。vcf 形式のデータの読み込みには、gastonパッケージに含まれる関数read.vcfを用いる。
# read a vcf file
bm <- read.vcf("data/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("data/RiceDiversityPheno4GWASGS.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形質（一部は、形質と試験地あるいは年次との組合せ）のデータが含まれている。
ここでは、形質例として種子長を用いる。種子長のデータを抜き出した後、各要素（各系統に対応）に名前を付ける。その後、na.omit関数を用いて欠測値を除く。
y <- pheno[, "Seed.length"] # select seed length as the target trait
names(y) <- rownames(pheno) # naming the elements of y
y <- na.omit(y) # remove missing elements解析に必要なマーカー遺伝子型データの抽出
実は、全ての系統について表現型データが得られているわけではない。そこで、bmから、表現型データをもつ系統のSNPデータを抜き出す。
# select samples with ids contained in rownames(pheno)
bm.wp <- bm[bm@ped$id %in% names(y),]
bm.wp # show the summary of the selected data## A bed.matrix with 361 individuals and 38769 markers.
## snps stats are set
## There is one monomorphic SNP
## ped stats are set356品種/系統のマーカー遺伝子型が抜き出される。
最後に、表現型データ\(y\)を、 SNPデータbm.wpと同じ順番で並び替える。
y <- y[bm.wp@ped$id] # reorder the elements of y as same as the bm@ped$id分集団構造解析
次に、356系統の遺伝的背景を調べるために、マーカー遺伝子型データに基づく主成分分析を行う。マーカー遺伝子型データは、上述したように“0”, “1”, “2”でスコア化され、行列のようなかたちでbm内に保存されている。
このデータを取り出し、その主成分分析を行う。ここでは、pcaMethodsパッケージのpca関数を利用する予定であったが、RStudio Cloudでは計算に長い時間を要するため、prcompを用いることとした。なお、pca関数では、欠測値がある場合でも、主成分分析を行うことができる。
# pca <- pca(as.matrix(bm.wp), "ppca", nPcs = 4, completeObs = F) # obtain the first 4 PCs
#pcs <- scores(pca) # obtain the PC scores
#rm(pca) # because the size of the object "pca" is big, remove it here
res <- prcomp(as.matrix(bm.wp)) # use prcomp for pca
pcs <- res$x[,1:4]
rm(res) # remove the PCA result because the object is large in size.なお、主成分分析の結果は、サイズが大きいので、主成分得点を抽出した後に削除しておく。
主成分分析の結果を図示する。
op <- par(mfrow = c(1,2))
plot(pcs[,1:2]) # a scatter plot graph between pc1 and pc2
plot(pcs[,3:4]) # a scatter plot graph between pc3 and pc4
par(op)遺伝的背景と形質の表現型の間に関連があると、GWASで偽陽性を生じる（実際、GWASでは、遺伝的背景の違いによる偽陽性が大きな問題なる）。ここでは、それを視覚化して確認する。
op <- par(mfrow = c(1,2))
symbols(pcs[,1:2], circle = y - min(y), inches = 0.1)
symbols(pcs[,3:4], circle = y - min(y), inches = 0.1)
par(op)遺伝的背景との関係は少なからずありそうである。
GWASの準備
ここからは、GWASを行う。まずは、そのための準備を行う。
まず、マイナーアリル頻度（MAF）が5%以上のSNPsだけを選ぶ。
bm.wp <- select.snps(bm.wp, maf > 0.05) # filter with MAF
bm.wp # show the data summary## A bed.matrix with 361 individuals and 36112 markers.
## snps stats are set
## ped stats are setこれにより、36,112のSNPマーカーがGWAS解析用に残る。
次に、遺伝的背景を説明するために、ゲノム関係行列を計算しておく。この行列を含めて（具体的には、この行列に規定される共分散構造をもつ多変量正規分布に従う効果を含める、詳細は省く）解析することで、偽陽性を抑え込むことができる。ゲノム関係行列を計算する前に、SNPスコア行列であるbm.wpを標準化（standardize）しておくしておく必要がある。ここでは、平均0、標準偏差1となるように標準化する（なお、自分で標準化を行わない場合は、平均\(2p\)、標準偏差はHardy-Weinberg平衡に従うときの期待値\(2p(1-p)\)で標準化される。ここで、\(p\)はA2対立遺伝子の頻度である）。
standardize(bm.wp) <- "mu_sigma"
grm <- GRM(bm.wp)GWASの実行
準備しておいたデータを用いて、GWASを実行する。
gwa <- association.test(bm.wp, Y = y, method = "lmm", response = "quantitative", test = "lrt", eigenK = eigen(grm), p = 4)gasonのassociation.test関数は計算が速く、短い時間で計算が終わる。
ここで、簡単に関数の引数について説明しておく。
bm.wpは解析に用いるSNPデータ、Y = yは、解析に用いる表現型データを指定している。method = “lmm”は、線形混合モデルをあてはめることを意味する。通常、lmmを指定して解析する。response = “quantitative”は、対象としている形質が量的形質であることを意味する。質的な二値形質の場合は、response = “binary”とする。また、test = “lrt”は、効果の検定法を表し、尤度比検定（Likelihood Ratio Test: LRT)を指定している。この他に、“wald”, “score”を選ぶことができる。これら2つの手法は、LRTよりも計算速度は速いが、計算に時間が問題にならないときには“lrt”を指定するとよい。eigenK = eigen(grm)、および、p = 4は、いずれも集団構造による偽陽性を抑えるための効果についてのオプションであり、前者は、ゲノム関係行列（の固有値分解の結果）を、後者は、遺伝的背景を説明するための主成分を上位何主成分までモデルに含めるかを指定している。ここでは、上位4主成分を遺伝的背景を表す変数としてモデルに組み込んでいる。
では、先頭の10 SNPsについての結果をみてみる。
head(gwa, 10)一番重要なのは、最後の列の\(p\)であり、これは、各SNPの効果の有意性検定の\(p\)値をあらわす。この値が小さいほど、効果が有意であることを意味する。
\(p\)値が小さい順に並べて、先頭10 SNPsについて結果をみてみる。
gwa[order(gwa$p)[1:10], ]第3染色体に高度に有意なSNPがあることが分かる。
GWASの結果を図示する。
次に、マンハッタンプロット（Manhattan plot）と呼ばれる図を描く。縦軸は、\(-log_{10}(p)\)を表す。\(p\)値が小さいと、\(-log_{10}(p)\)は大きくなる。
manhattan(gwa, pch = 20)
なお、非常に多くのSNPについて検定を行っているので、多重検定により生じる偽陽性をコントロールする必要がある。
多重検定に対して\(p\)値を補正して、偽陽性をコントロールする方法の一つとして、BenjaminiとHochberg (1995)が提案した偽発見率（false discovery rate: FDR）がある。これは、検出された陽性の結果のうち、偽陽性が含まれる割合の期待値を計算して、それに基づき偽陽性をコントロールする方法である。ここでは、偽陽性率（FDR）が5%以下になるように、\(p\)値の閾値を調整する。
fdr <- p.adjust(gwa$p, method = "BH")
gwa[fdr < 0.05, ]第3染色体と第5染色体に、有意なSNPが有ることが分かる。
この結果を、図示して確認する。有意なSNPを緑色で表す。
col <- rep("black", nrow(gwa))
col[gwa$chr %% 2 == 1] <- "gray50"
col[fdr < 0.05] <- "green"
manhattan(gwa, pch = 20, col = col)
分集団構造（遺伝的背景）による影響
GWASに用いる集団に分集団構造がある場合には、分集団構造が存在することにより偽陽性が生じやすくなる。
ここでは、主成分得点をGWASのモデルに組み込まない場合、あるいは、ゲノム関係行列をGWASのモデルに組み込まない場合について、陽性となるSNPがどのように増加するかを確認する。なお、ゲノム上の多くのSNPは、対象とする形質の表現型に関わっていない場合がほとんどである。そのような場合、形質の表現型に関わらないSNPの\(p\)値は一様分布に従うと期待される。そこで、\(p\)値が一様分布に従うという仮定のもとでQ-Q（Quantile-Quantile）プロットを描き、一様分布から逸脱するSNPを確認してみる。
まずは、上述したモデルでQ-Qプロットを描く。
qqplot.pvalues(gwa$p, pch = 20)
一部、逸脱するSNPがあるが、これらは表現型を支配する遺伝子の中、あるいは、近傍にあり、一様分布に従わないSNPと考えられる。
一方、主成分を一つも組み込まないモデルでGWASを行う。これは、p = 0$*とすればよい。
gwa.wopc <- association.test(bm.wp, Y = y, method = "lmm", response = "quantitative", test = "lrt", eigenK = eigen(grm), p = 0)Q-Qプロットで、\(p\)値を図示する。
qqplot.pvalues(gwa.wopc$p, pch = 20)
結果は、主成分を含むモデルとほとんど変わらない。つまり、この形質では、主成分の効果をモデルに入れても入れなくても結果は大きく変わらないと考えられる。
次に、ゲノム関係行列で説明される効果を含まないモデルでGWASを実施する。method = “lm”、test = “wald”とする。
gwa.wogrm <- association.test(bm.wp, Y = y, method = "lm", response = "quantitative", test = "wald")Q-Qプロットで、\(p\)値を図示する。
qqplot.pvalues(gwa.wogrm$p, pch = 20)
ほとんど全てのSNPで、\(p\)が一様分布から逸脱していることがわかる。このように、遺伝的背景を無視して、単純にSNPの遺伝子型と、形質の表現型の関連を解析すると、偽陽性が非常に多くなることが分かる。
ゲノミック選抜
表現型がある系統と無い系統
では、ここから、ゲノミックセレクション(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に標準化しておく。
# extract the marker score matrix from the original data, bm
standardize(bm) <- "mu_sigma"
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(nrow(X.wp), nrow(X.wop))## [1] 356 1212表現型データ\(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)## 1
## [1,] 0.9439858なお、入力される説明変数の数が多い場合には、あてはまりの良いモデルを構築するのは容易である。したがって、あてはまりの良さは、モデルの予測能力を反映しない。モデルの予測能力を評価するためには、モデル構築に用いなかったデータに対する「予測」を行い、その精度を評価する必要がある。そこで、本当はモデル構築に用いることができるデータの一部を、モデル構築から除いておき、除いておいたデータを予測して、その予測値を観察された値と比較する。除いておいたデータは、「本当はモデル構築に用いることができるデータ」であり、表現型データ（\(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] 2 2 3 1 2 3 2 3 2 1 3 2 3 1 3 3 1 1 2 1 1 2 1 3 3 3 1 1 2 3 3 2 3 1 3 1 1
## [38] 1 2 3 3 2 2 2 2 1 2 1 3 2 3 2 1 2 1 3 1 3 1 3 3 1 2 3 2 2 1 1 2 3 1 1 3 3
## [75] 2 3 2 1 3 1 3 3 1 3 3 1 2 2 1 2 2 1 3 1 3 1 3 1 2 3 3 2 2 1 2 2 2 3 2 2 3
## [112] 2 3 3 3 1 2 3 2 3 3 3 1 2 2 1 1 1 2 1 2 1 1 1 1 2 3 3 3 2 1 1 3 3 1 1 3 1
## [149] 3 3 3 3 1 1 3 1 3 3 1 2 2 2 3 1 3 2 2 2 3 1 2 2 1 1 3 2 1 3 1 1 2 2 3 3 3
## [186] 2 2 3 2 2 1 1 2 2 1 2 2 1 1 3 1 1 2 2 3 3 3 2 3 3 2 3 1 2 2 1 1 2 3 2 1 1
## [223] 2 2 2 1 2 3 1 2 1 2 2 3 3 2 3 1 2 3 1 3 1 1 3 2 1 2 2 1 2 2 2 1 3 1 3 3 2
## [260] 2 2 1 1 2 3 2 1 2 1 1 2 3 3 1 3 2 1 3 2 3 3 1 1 2 1 1 3 1 3 3 3 1 2 1 1 2
## [297] 1 1 2 3 1 3 1 3 2 2 2 3 1 2 2 3 1 2 1 3 1 3 3 2 2 3 2 3 3 2 3 1 3 1 3 1 3
## [334] 1 2 3 2 2 3 2 1 1 2 3 3 3 1 1 2 1 3 3 1 1 2 2次に、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)## 1
## [1,] 0.7440969実際には、各分割について、同様の計算を繰り返し、それらをまとめて精度評価する。そこで、以上の操作を、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.7055297リッジ回帰で推定された各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] 78最近、多数の遺伝資源についてゲノムデータが解析・公開されるようになってきている。例えば、イネでは遺伝資源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.6767193全データを用いて推定された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] 53用いたモデル化手法とモデルが異なるため、リッジ回帰の場合と少し異なる結果が得られる。実際に適用する場合には、複数の手法の中から、それぞれの形質に合った手法を経験的に選んで用いるのが良い。どの形質にも万能な手法は無い（データ科学の世界では、これを「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, saveAt = "temp/", 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.9379417交差検証を行い、予測精度を評価する。
# 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, saveAt = "temp/", 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.7066438表現型データがある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] 53ここでも、リッジ回帰の結果とよく似た結果が得られることが分かる。