ここでは、Rを用いて各種データ解析法を用いて、多数の遺伝資源から有望な系統を効率的に選び出してみよう。 なお、Zhaoら（2011; Nature Communications 2:467）がイネ遺伝資源のGWASに用いた表現型データと、 McCouchら（2016; Nature Communication 7:10532）が用いたHigh Density Rice Array (HDRA)を用いてジェノタイピングされたSNP遺伝子型データ（いずれも、Rice Diversityからダウンロードできる）を、例として用いる。なお、後者のデータは、70万SNPsのデータのうち、他のマーカーとの冗長性が低く（LD < 0.5）、遺伝子型データの欠測率が低く（< 1%）、かつ、マイナーアリル頻度（minor allele frequency: MAF）が高く（> 5%）、10,144 SNPsのデータを用いる。なお、欠測値については、Beagle 5.1（Blowning et al. 2018）を用いて補完をしてある。
require(tidyverse)
require(plotly)
require(Rtsne)
require(ranger)まず、マーカー遺伝子型データを読み込む。ここでは、RDS形式として保存してあるので、readRDS関数を用いて読み込む。
geno <- readRDS("geno.rds")
geno[1:5,1:5] # show the head of the file## SNP.1.8562. SNP.1.30477. SNP.1.32666. SNP.1.140083.
## NSFTV1@86f75d2b.0 0 0 2 0
## NSFTV3@5ef1be74.0 2 0 2 0
## NSFTV4@81d03b86.0 2 0 2 0
## NSFTV5@5533f406.0 0 0 2 2
## NSFTV6@0d125c0e.0 2 0 2 0
## SNP.1.163649.
## NSFTV1@86f75d2b.0 0
## NSFTV3@5ef1be74.0 0
## NSFTV4@81d03b86.0 0
## NSFTV5@5533f406.0 0
## NSFTV6@0d125c0e.0 0イネ遺伝資源アクセッション388系統についての、10,144 一塩基多型（single nucleatide polymorphisms）SNPsの遺伝子型データが含まれている。行がアクセッション（品種・系統）、列がSNPに対応している。
SNP遺伝子型は、“0”と”2”でコードされている。また、中には”1”でコードされているSNP遺伝子型もある。“0”, “1”, “2”は、それぞれ、A1/A1、A1/A2、A2/A2の遺伝子型を表し、“0”と”2”はホモ接合体、“1”はヘテロ接合体である。なお、イネは自殖であるため、ヘテロ接合の頻度は低い。
次に、おなじ遺伝資源アクセッションについて得られた表現型データを読み込む。
pheno <- read.csv("pheno.csv", row.names = 1, stringsAsFactors = T)
str(pheno) # show the data summary## 'data.frame': 388 obs. of 38 variables:
## $ Accession.Name : Factor w/ 384 levels "1021","27","318",..: 10 12 239 240 19 20 22 26 27 28 ...
## $ Sub.population : Factor w/ 6 levels "ADMIX","AROMATIC",..: 5 4 3 2 3 6 6 5 5 4 ...
## $ 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形質（一部は、形質と試験地あるいは年次との組合せ）のデータが含まれている。
388系統の遺伝的背景を調べるために、マーカー遺伝子型データを用いて、第3回講義で紹介した主成分分析を行う。マーカー遺伝子型データは、上述したように”0”, “1”, “2”でスコア化され、行列のようなかたちでgenoに保存されている。
このデータを取り出し、その主成分分析を行う。ここでは、prcompを用いて解析をする。
res <- prcomp(geno) # use prcomp for pca
pcs <- res$x[,1:4]なお、主成分分析の結果は、サイズが大きいので、主成分得点を抽出した後に削除しておく。
主成分分析の結果を図示する。
acc <- pheno %>% select(Accession.Name, Sub.population)
df <- data.frame(acc, pcs)
plot_ly(data = df, x = ~PC1, y = ~PC2, z = ~PC3,
color = ~Sub.population, text = ~Accession.Name,
type = "scatter3d", mode = "markers")分集団がかたまって散布されていることが分かる。これは、イネの遺伝資源に比較的明瞭な分集団構造があることを示している。なお、分離されたグループ間に主に散布されているADMIXは、分集団間の混合された遺伝背景をもつ系統が含まれるグループである。
次に、t-SNEを用いた解析も行ってみる。
res <- Rtsne(geno, perplexity = 30)df <- data.frame(acc, tSNE1 = res$Y[,1], tSNE2 = res$Y[,2])
plot_ly(data = df, x = ~tSNE1, y = ~tSNE2,
color = ~Sub.population, text = ~ Accession.Name,
type = "scatter", mode = "markers")2次元でもうまく分離できている。この図でも、ADMIXは中間的な位置に散布されていることも分かる。
以下の4形質はいずれも二値形質である。
table(pheno$Leaf.pubescence)##
## 0 1
## 42 235table(pheno$Awn.presence)##
## 0 1
## 272 81table(pheno$Seed.color)##
## 0 1
## 354 7table(pheno$Pericarp.color)##
## 0 1
## 330 31予測モデルを作るためにデータを準備する。
ここでは、葉の毛茸（leaf pubescence）の有無を予測するモデルを作成する。葉の毛茸をyとして抜き出し、ゲノムデータと結合しておく。また、欠測値を除いておく。
y <- as.factor(pheno$Leaf.pubescence)
gp <- data.frame(y, geno) %>% drop_na()random Forestを用いてモデルを作る。
model <- ranger(y ~ ., data = gp)得られた結果を正解と比較する。
t <- table(gp$y, model$predictions)
t##
## 0 1
## 0 25 17
## 1 5 230sum(diag(t)) / sum(t)## [1] 0.9205776正解率は90%以上。しかし、これはモデルの訓練に使ったデータに対する正解率である。たいてい、過大評価となる。
そこで、交差検証を行う。5分割の交差検証を行ってみよう。
まずは、無作為に1〜5のIDを発生させる。
n.fold <- 5
n.sample <- nrow(gp)
cv.id <- sample(seq(n.sample) %% n.fold) + 1そのあと、1分割ずつ検証データ（test）として除いたデータ（train）をもとに訓練をし、検証データの予測を行う。5分割なので、5回繰り返す
y.pred <- rep(NA, n.sample)
for(i in 1:n.fold) {
test <- cv.id == i
train <- !test
model <- ranger(y ~ ., data = gp[train, ])
y.pred[test] <- predict(model, gp[test, ])$predictions
}精度を評価する。
t <- table(gp$y, y.pred)
print(t)## y.pred
## 1 2
## 0 24 18
## 1 3 232accuracy <- sum(diag(t)) / sum(t)
accuracy## [1] 0.924187790%程度の精度である。悪くない。
オリジナルのデータでは、表現型データが無い品種・系統についてもゲノムデータが提供されている。そこで、これら品種・系統についてゲノムから形質の予測を行ってみる。
まずは、ゲノムデータを読み込む。
geno.wp <- readRDS("geno_wp.rds")
dim(geno.wp)## [1] 500 5945500系統のデータがある。
まずは、表現型データがあるすべてのデータ（gp）を用いてモデルを作る。
model <- ranger(y ~ ., data = gp)次に、このモデルを用いて、表現型データが無い品種・系統について、毛茸の有無を予測する。
y.pred.wp <- predict(model, geno.wp)$predictions予測値を確認する。
table(y.pred.wp)## y.pred.wp
## 0 1
## 40 460ほとんどの系統が毛茸をもつことがわかる。実際に観察しなくても、予測ができることがわかる（もちろん、予測の誤りがあることに注意する）。
ここでは、個体あたりの穂数（panicle number per plant）の予測を行うモデルを作成してみる。
穂数をyとして抜き出し、ゲノムデータと結合する。また、欠測データを除いておく。
y <- pheno$Panicle.number.per.plant
gp <- data.frame(y, geno) %>% drop_na()ここでは、random Forestを用いて予測をおこなってみる。
model <- ranger(y ~ ., data = gp)y.pred <- model$predictions
cor(gp$y, y.pred)## [1] 0.82463025分割の交差検証を行う。
n.fold <- 5
n.sample <- nrow(gp)
cv.id <- sample(seq(n.sample) %% n.fold) + 1予測を行ってみる。
y.pred <- rep(NA, n.sample)
for(i in 1:n.fold) {
test <- cv.id == i
train <- !test
model <- ranger(y ~ ., data = gp[train, ])
y.pred[test] <- predict(model, gp[test, ])$predictions
}cor(gp$y, y.pred)## [1] 0.8167404表現型データが無い品種・系統を予測する。まずは、すべてのデータを用いてモデルを作る。
model <- ranger(y ~ ., data = gp)作ったモデルをもとに、表現型データが無い品種・系統を予測する。
y.pred.wp <- predict(model, geno.wp)$predictions予測された値（表現型データが無い品種・系統）と、観察されている品種・系統（訓練データ）の観察値でヒストグラムを描く。
plot_ly(x = y.pred.wp, type = "histogram",
name = "without phenotypes") %>%
add_trace(x = gp$y, type = "histogram",
name = "with phenotypes (obs)")予測された値（表現型データが無い品種・系統）と、観察されている品種・系統（訓練データ）予測値でもヒストグラムを描く。
plot_ly(x = y.pred.wp, type = "histogram",
name = "without phenotypes") %>%
add_trace(x = y.pred, type = "histogram",
name = "with phenotypes (pred)")予測値では、観察値に比べると、値の範囲が狭いことがわかる。
| 形質名 | 説明 | 形質のタイプ |
|---|---|---|
| “Flowering.time.at.Arkansas” | Arkansasにおける開花日（播種後日数） | 連続値 |
| “Flowering.time.at.Faridpur” | Faridpurにおける開花日（播種後日数） | 連続値 |
| “Flowering.time.at.Aberdeen” | Aberdeenにおける開花日（播種後日数） | 連続値 |
| “FT.ratio.of.Arkansas.Aberdeen” | ArkansasとAberdeen間の開花日の比 | 連続値 |
| “FT.ratio.of.Faridpur.Aberdeen” | FaridpurとAberdeen間の開花日の比 | 連続値 |
| “Culm.habit” | 主稈の基部の垂直からの推定平均傾斜角 | 階級値 |
| “Leaf.pubescence” | 葉の毛茸の有無 | 二値 |
| “Flag.leaf.length” | 止め葉の長さ | 連続値 |
| “Flag.leaf.width” | 止め葉の幅 | 連続値 |
| “Awn.presence” | 芒（のぎ）の有無 | 二値 |
| “Panicle.number.per.plant” | 個体あたりの穂数 | 連続値（計数値） |
| “Plant.height” | 草丈 | 連続値 |
| “Panicle.length” | 穂の長さ | 連続値 |
| “Primary.panicle.branch.number” | 一次枝梗の数 | 連続値（計数値） |
| “Seed.number.per.panicle” | 穂あたりの種子数 | 連続値（計数値） |
| “Florets.per.panicle” | 穂あたりの穎花数 | 連続値（計数値） |
| “Panicle.fertility” | 穂の種子の稔性 | 連続値（割合) |
| “Seed.length” | 種子長 | 連続値 |
| “Seed.width” | 種子幅 | 連続値 |
| “Seed.volume” | 種子の体積 | 連続値（計算値） |
| “Seed.surface.area” | 種子の表面積 | 連続値（計算値） |
| “Brown.rice.seed.length” | 玄米長 | 連続値 |
| “Brown.rice.seed.width” | 玄米幅 | 連続値 |
| “Brown.rice.surface.area” | 玄米の表面積 | 連続値（計算値） |
| “Brown.rice.volume” | 玄米の体積 | 連続値（計算値） |
| “Seed.length.width.ratio” | 種子の長幅比 | 連続値（計算値） |
| “Brown.rice.length.width.ratio” | 玄米の長幅比 | 連続値（計算値） |
| “Seed.color” | 種子の色 | 二値 |
| “Pericarp.color” | 果皮色 | 二値 |
| “Straighthead.suseptability” | 青立病感受性 | 階級値 |
| “Blast.resistance” | いもち病抵抗性 | 階級値 |
| “Amylose.content” | アミロース含量 | 連続値 |
| “Alkali.spreading.value” | アルカリ崩壊度 | 階級値 |
| “Protein.content” | タンパク質含量 | 連続値 |
| “Year07Flowering.time.at.Arkansas” | 2007年のArkansasにおける開花日 | 連続値 |
| “Year06Flowering.time.at.Arkansas” | 200６年のArkansasにおける開花日 | 連続値 |
例えば、穂の数と長さを予測してみる。
穂の数のモデル構築と予測
y.1 <- pheno$Panicle.number.per.plant
gp <- data.frame(y.1, geno) %>% drop_na()
model <- ranger(y.1 ~ ., data = gp)
y.pred.wp.1 <- predict(model, geno.wp)$predictions穂の長さのモデル構築と予測
y.2 <- pheno$Panicle.length
gp <- data.frame(y.2, geno) %>% drop_na()
model <- ranger(y.2 ~ ., data = gp)
y.pred.wp.2 <- predict(model, geno.wp)$predictions予測された穂の数と穂の長さの散布図
plot_ly(x = y.pred.wp.1, y = y.pred.wp.2,
type = "scatter", mode = "markers", text = rownames(geno.wp))