主成分分析・クラスタ解析・分類器

はじめに

本講義では、イネの種子形をゲノムワイドマーカー遺伝子型から予測をした論文（Iwata et al. 2015)で使われたデータを用いて、主成分分析、クラスタ解析、分類器を用いた解析について、Rを用いた具体的な解析方法を説明します。

必要なパッケージ

最初に、解析に必要となるパッケージを読み込みます。 hereは階層構造をもつフォルダにあるデータファイルの指定に便利です。 plotlyは、インタラクティブなグラフの作成に、apeはデンドログラムの描画に利用します。 また、kernlabは、各種カーネル法を用いた解析を行うためのパッケージです。 ここでは、サポートベクターマシンを用いた分類器の構築に用います。

require(here)
require(plotly)
require(ape)
require(kernlab)

必要なデータの読み込み

　イネ種子形状表現型データの読み込みを行います。 これは、第1回講義で説明した楕円フーリエ記述子を用いてイネの種子形状を定量化したデータです。 イネ遺伝資源374系統の種子の画像から種子の輪郭を抽出し、それを楕円フーリエ記述子で定量化した後、系統平均を求めたものです。

ef.data <- read.csv(here("data",
                    "FourierDescriptor_datasetC.csv"),
               row.names = 1)
head(ef.data)[,1:8]

##   a1        b1        c1        d1           a2           b2           c2
## 1  1 -2.70e-17  4.80e-17 0.6111694  0.001039338 -0.005764824 -0.009906241
## 3  1 -4.45e-18  9.03e-18 0.5004835  0.002268561 -0.010033800 -0.012319853
## 4  1 -1.95e-17  3.62e-17 0.4762106  0.004475037 -0.010257737 -0.016015754
## 5  1 -2.68e-17  8.04e-17 0.3628859  0.000684247 -0.005984455 -0.019866794
## 6  1  7.11e-18 -5.04e-18 0.6107227  0.003886584  0.002052005  0.001065063
## 7  1 -1.76e-17  3.64e-17 0.4894189 -0.001395340 -0.001280532 -0.013163296
##             d2
## 1  0.010476213
## 3  0.006070962
## 4  0.007118144
## 5 -0.006931079
## 6  0.016307643
## 7  0.003520876

dim(ef.data)

## [1] 374  80

374×80（374系統分の80個のフーリエ記述子）のデータです。なお、最初の3列（\(a_1, b_1, c_1\)）は、輪郭のサイズ、向き、計測開始点の違いに対する標準化により定数（\(a_1 = 1, b_1 = 0, c_1 = 0\)）となっていることに注意しましょう。

　次に、ゲノムワイドマーカー遺伝子型データの読み込みを行います。これは、種子形状を定量化した386系統のゲノムワイド一塩基多型（single nucleotide polymorphisms: SNPs）の遺伝子型データです。

gs.data <- read.csv(here("data",
                    "GenotypeScore_datasetC.csv"),
               row.names = 1)
head(gs.data)[,1:8]

##   id1000001 id1000003 id1000005 id1000007 id1000008 id1000011 id1000013
## 1         2         2         0         2         2         0         0
## 3         0         0         0         0         0         0         0
## 4         0         0         0         0         0         2         0
## 5         0         0         2         2         0         0         2
## 6         0         0         0         0         0         2         0
## 7         2         2         0         2         2         0         0
##   id1000015
## 1         2
## 3         0
## 4         0
## 5         0
## 6         0
## 7         2

dim(gs.data)

## [1]   374 36901

374x36901（374系統分の36,901 SNPsのマーカー遺伝子型）のデータです。なお、数値として記述されていて、参照配列（日本晴 Nipponbare）の全ゲノム配列と同じ多型を保つ場合は0、変異型の対立遺伝子をホモでもつ場合は2となっています。なお、これらイネ遺伝資源は、基本的に純系であるためほぼ全てのSNPがホモ接合であると考えられますが、一部のデータで、0, 2ではない中間的な値をもつ場合があります。これは、欠測補完を行った際に期待値として補完されたか、あるいは、実際に観察されたのが中間値（例えば、ヘテロ接合の場合は1となる）であった場合です。

イネの種子輪郭の再現

では、イネの種子の輪郭を再現してみます。まずは、一次近似を行います。 楕円フーリエ記述子による一次近似は、楕円による近似となります。 ここでは、最初の系統（1行目のデータ）について、一次近似に関わる4変数\(a_1, b_1, c_1, d_1\)から楕円を描いてみる。 その計算式は、 \[ x(\theta) = a_1 \cos \theta + b_1 \sin \theta \\ y(\theta) = c_1 \cos \theta + d_1 \sin \theta \] となります。 ここでは、\(\theta\)を0.01ラジアンずつ変化させて、\(x(\theta)\)と\(y(\theta)\)を計算しています。

efm <- as.matrix(ef.data)
ef <- efm[1,] # 最初のデータをとりだす
theta <- seq(0, 2*pi, 0.01)
x <- ef[1] * cos(theta) + ef[2] * sin(theta)
y <- ef[3] * cos(theta) + ef[4] * sin(theta)
plot(x, y, type = "l", xlim = c(-1.2, 1.2), asp = 1)

つぎに、高次のフーリエ記述子も用いて、近似輪郭を再現します。楕円フーリエ記述子の系統平均から得られる近似輪郭は、系統の平均的形状を表します。計算式は、 \[ x(\theta) = \sum_k^N (a_k \cos k\theta + b_k \sin k\theta) \\ y(\theta) = \sum_k^N (c_k \cos i\theta + d_k \sin k\theta) \] です。

x <- 0; y <- 0
for(i in 1:20) {
  x <- x + ef[i * 4 - 3] * cos(i * theta) + ef[i * 4 - 2] * sin(i * theta)
  y <- y + ef[i * 4 - 1] * cos(i * theta) + ef[i * 4    ] * sin(i * theta)
}
plot(x, y, type = "l", xlim = c(-1.2, 1.2), asp = 1)

なお、同計算を簡単にするために、function関数を用いて、自作の関数を作成します。この関数は、楕円フーリエ記述子を入力として、輪郭座標\(x(\theta), y(\theta)\)を出力する関数となります。

ef2coord <- function(ef, theta = seq(0, 2*pi, 0.01)) {
  x <- 0; y <- 0
  for(i in 1:(length(ef)/4)) {
    x <- x + ef[i * 4 - 3] * cos(i * theta) + ef[i * 4 - 2] * sin(i * theta)
    y <- y + ef[i * 4 - 1] * cos(i * theta) + ef[i * 4    ] * sin(i * theta)
  }
  coord <- list(x = x, y = y)
  return(coord)
}

では、自作関数ef2coordを用いて、最初の100系統（1〜100行目）の平均的な種子形状を表す輪郭を描いてみましょう。

par(mfrow = c(10, 10), mar = c(1, 1, 1, 1))
for(i in 1:100) {
  ef <- efm[i,]
  coord <- ef2coord(ef)
  plot(coord$x, coord$y, type = "l", xlim = c(-1.2, 1.2), asp = 1,
       ann = F, axes = F)
}

種子の長幅比に大きな変異があることがわかります。

イネの種子輪郭形状の主成分分析

次に、これらデータをもとに、主成分分析を行ってみましょう。 なお、主成分分析については第2回講義で説明しました。 まず、最初の3列（\(a_1, b_1, c_1\)）のデータは定数なので、主成分分析では除いて解析を行う必要があります。 ここでは、これら3変数を除いた374×77のデータを行列\(X\)とします。

X <- efm[, -(1:3)]
head(X)[,1:8]

##          d1           a2           b2           c2           d2         a3
## 1 0.6111694  0.001039338 -0.005764824 -0.009906241  0.010476213 0.04236411
## 3 0.5004835  0.002268561 -0.010033800 -0.012319853  0.006070962 0.05898568
## 4 0.4762106  0.004475037 -0.010257737 -0.016015754  0.007118144 0.06463415
## 5 0.3628859  0.000684247 -0.005984455 -0.019866794 -0.006931079 0.08003286
## 6 0.6107227  0.003886584  0.002052005  0.001065063  0.016307643 0.04761987
## 7 0.4894189 -0.001395340 -0.001280532 -0.013163296  0.003520876 0.06502755
##             b3           c3
## 1  0.003442847  0.000106244
## 3 -0.003249424 -0.008229188
## 4  0.002837122 -0.010016006
## 5 -0.000073500 -0.006161375
## 6 -0.000371803 -0.003549855
## 7 -0.000604691 -0.004130003

行列\(X\)について、77変数間の分散・共分散行列を計算します。

vx <- cov(X)
head(vx)[,1:6]

##               d1            a2            b2            c2            d2
## d1  8.171798e-03 -9.415175e-07 -1.328381e-05  1.970368e-04  4.234751e-04
## a2 -9.415175e-07  1.942862e-05  1.758903e-06  4.424538e-06 -4.651919e-06
## b2 -1.328381e-05  1.758903e-06  1.371549e-05  1.789137e-05 -8.760429e-07
## c2  1.970368e-04  4.424538e-06  1.789137e-05  6.081415e-05  1.517491e-05
## d2  4.234751e-04 -4.651919e-06 -8.760429e-07  1.517491e-05  6.778285e-05
## a3 -1.282188e-03 -3.901075e-06  3.512568e-06 -2.970499e-05 -6.730872e-05
##               a3
## d1 -1.282188e-03
## a2 -3.901075e-06
## b2  3.512568e-06
## c2 -2.970499e-05
## d2 -6.730872e-05
## a3  2.113467e-04

第2回講義でも説明したように、分散・共分散行列の固有値分解を行うことで、主成分の固有値と固有ベクトルを計算できます。 まずは、固有値分解を行います。 Rでは、eigenという関数を用いて固有値分解ができます。

eig <- eigen(vx)
str(eig)

## List of 2
##  $ values : num [1:77] 8.51e-03 7.91e-05 5.59e-05 3.85e-05 3.28e-05 ...
##  $ vectors: num [1:77, 1:77] 9.80e-01 1.53e-06 -1.57e-03 2.38e-02 5.11e-02 ...
##  - attr(*, "class")= chr "eigen"

固有値分解の結果得られた固有値を図示します。

barplot(eig$values)

barplot(eig$values[1:10])

第1主成分の固有値が非常に大きく、形状変異のほとんどが同主成分で説明できることが分かります。

次に、寄与率と累積寄与率を計算します。 寄与率とは、全固有値の和分のその主成分の固有値の割合で、各主成分の説明する変動（分散）の割合を表します。 累積寄与率は、寄与率の累積です。 なお、第2回の講義でも説明したように、全固有値の和は、もとの全変数のもつ分散の和に一致します。

cntr <- eig$values / sum(eig$values)
plot(cntr[1:10], ylim = c(0, 1), type = "h")
lines(1:10, cumsum(cntr[1:10]), col = "green")

なお、同じ解析を主成分分析専用の関数prcompを用いて行うこともできます。

pca <- prcomp(X)
summary(pca)

## Importance of components:
##                            PC1      PC2     PC3      PC4      PC5      PC6
## Standard deviation     0.09227 0.008895 0.00748 0.006205 0.005723 0.005081
## Proportion of Variance 0.96307 0.008950 0.00633 0.004350 0.003710 0.002920
## Cumulative Proportion  0.96307 0.972020 0.97835 0.982700 0.986410 0.989330
##                             PC7      PC8      PC9     PC10     PC11     PC12
## Standard deviation     0.003904 0.003532 0.003209 0.003113 0.002715 0.002449
## Proportion of Variance 0.001720 0.001410 0.001160 0.001100 0.000830 0.000680
## Cumulative Proportion  0.991050 0.992460 0.993630 0.994720 0.995560 0.996230
##                            PC13     PC14     PC15     PC16     PC17     PC18
## Standard deviation     0.002083 0.001901 0.001818 0.001527 0.001461 0.001404
## Proportion of Variance 0.000490 0.000410 0.000370 0.000260 0.000240 0.000220
## Cumulative Proportion  0.996730 0.997130 0.997510 0.997770 0.998010 0.998240
##                            PC19     PC20    PC21     PC22     PC23      PC24
## Standard deviation     0.001283 0.001219 0.00111 0.001073 0.001002 0.0009341
## Proportion of Variance 0.000190 0.000170 0.00014 0.000130 0.000110 0.0001000
## Cumulative Proportion  0.998420 0.998590 0.99873 0.998860 0.998970 0.9990700
##                             PC25      PC26      PC27      PC28      PC29
## Standard deviation     0.0008949 0.0007895 0.0007808 0.0007392 0.0007129
## Proportion of Variance 0.0000900 0.0000700 0.0000700 0.0000600 0.0000600
## Cumulative Proportion  0.9991600 0.9992300 0.9993000 0.9993600 0.9994200
##                             PC30      PC31      PC32      PC33      PC34
## Standard deviation     0.0006562 0.0006357 0.0006032 0.0005796 0.0005238
## Proportion of Variance 0.0000500 0.0000500 0.0000400 0.0000400 0.0000300
## Cumulative Proportion  0.9994700 0.9995200 0.9995600 0.9996000 0.9996300
##                             PC35      PC36      PC37      PC38      PC39
## Standard deviation     0.0005132 0.0004908 0.0004799 0.0004749 0.0004569
## Proportion of Variance 0.0000300 0.0000300 0.0000300 0.0000300 0.0000200
## Cumulative Proportion  0.9996600 0.9996800 0.9997100 0.9997400 0.9997600
##                             PC40     PC41      PC42      PC43      PC44
## Standard deviation     0.0004244 0.000393 0.0003748 0.0003554 0.0003496
## Proportion of Variance 0.0000200 0.000020 0.0000200 0.0000100 0.0000100
## Cumulative Proportion  0.9997800 0.999800 0.9998100 0.9998300 0.9998400
##                             PC45      PC46     PC47     PC48      PC49
## Standard deviation     0.0003328 0.0003187 0.000309 0.000299 0.0002951
## Proportion of Variance 0.0000100 0.0000100 0.000010 0.000010 0.0000100
## Cumulative Proportion  0.9998500 0.9998700 0.999880 0.999890 0.9999000
##                             PC50      PC51      PC52     PC53      PC54
## Standard deviation     0.0002824 0.0002768 0.0002594 0.000252 0.0002434
## Proportion of Variance 0.0000100 0.0000100 0.0000100 0.000010 0.0000100
## Cumulative Proportion  0.9999000 0.9999100 0.9999200 0.999930 0.9999300
##                             PC55      PC56      PC57      PC58      PC59
## Standard deviation     0.0002329 0.0002317 0.0002268 0.0002157 0.0002053
## Proportion of Variance 0.0000100 0.0000100 0.0000100 0.0000100 0.0000000
## Cumulative Proportion  0.9999400 0.9999500 0.9999500 0.9999600 0.9999600
##                             PC60      PC61      PC62      PC63     PC64
## Standard deviation     0.0001953 0.0001923 0.0001778 0.0001684 0.000166
## Proportion of Variance 0.0000000 0.0000000 0.0000000 0.0000000 0.000000
## Cumulative Proportion  0.9999700 0.9999700 0.9999800 0.9999800 0.999980
##                             PC65      PC66      PC67      PC68      PC69
## Standard deviation     0.0001507 0.0001482 0.0001426 0.0001299 0.0001236
## Proportion of Variance 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## Cumulative Proportion  0.9999800 0.9999900 0.9999900 0.9999900 0.9999900
##                            PC70      PC71      PC72      PC73      PC74
## Standard deviation     0.000106 0.0001028 9.605e-05 9.366e-05 8.815e-05
## Proportion of Variance 0.000000 0.0000000 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  0.999990 0.9999900 1.000e+00 1.000e+00 1.000e+00
##                             PC75      PC76      PC77
## Standard deviation     8.513e-05 8.328e-05 7.244e-05
## Proportion of Variance 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  1.000e+00 1.000e+00 1.000e+00

prcompで得られる結果をplot関数でそのまま図示すると、固有値の棒グラフが表示されます。

plot(pca)

次に、主成分得点を計算してみます。

第2回講義でも説明したように、主成分得点は、もとの変数を固有ベクトルの要素で重み付けした一次結合として計算されます。 計算は、もとの変数のある標本の値をベクトルとして表すと、そのベクトルと固有ベクトルの内積として表せます。 ここでは、全ての標本のベクトルを束ねた行列と、全ての主成分の固有ベクトルを束ねた行列の積として、全ての系統の全ての主成分の得点を同時にに計算しています。

X.centered <- scale(X, center = T, scale = F)
z <- X.centered %*% eig$vectors

zが計算された主成分得点であり、各行が各標本、各列が各主成分に対応します。 plot関数を用いて、第1、第2主成分の対散布を描いてみます。

plot(z[, 1:2])

同じ図は、prcompを用いて計算された結果をもとに描くこともできます。

pca <- prcomp(X)
plot(pca$x[, 1:2])

イネのゲノムワイドSNP遺伝子型データの主成分分析

　イネの輪郭形状の主成分分析と同様の解析は、ゲノムワイドSNPの遺伝子型データについても適用できます。 まずは、同データを行列として取り出してみましょう。

G <- as.matrix(gs.data)
dim(G)

## [1]   374 36901

rm(gs.data) # メモリを圧迫するので削除しておく

374x36901の行列として取り出すことができます。

注意したいことは、このデータの場合は、上で行ったように、分散・共分散行列を計算してから主成分分析を行うと非常に時間とコンピュータ上のメモリを使います。場合によっては、PC上では計算できないかもしれません。

これは、変数の数が36,901あるためです（上の楕円フーリエ記述子の例では77でした）。 したがって、分散・共分散行列の次元は、\(36901\times36901\)となります。この計算には時間とメモリを要します。

そこで、別の方法を用いて主成分分析を行います。 第2回の講義で説明したように、距離行列か内積行列が分かれば、そのスペクトル分解から主成分分析と同様の解析ができます。 ここでは、内積行列（各系統のSNP遺伝子型をベクトルとして扱い、全系統間で内積を計算した行列）を計算して、 それをもとに主成分分析を行います。

まずは、SNP遺伝子型のスコア（0-2）について、scale関数を用いて、SNPごとに平均0になるように基準化します。 次に、tcrossprod関数を用いて内積行列を計算しています。

G.centered <- scale(G, center = T, scale = F)
GG <- tcrossprod(G.centered)
dim(GG)

## [1] 374 374

rm(G.centered)

系統数\(\times\)系統数の行列となっていることが分かります。

次に、この内積行列の固有値分解（スペクトル分解）を行います。

eig.G <- eigen(GG)
str(eig.G)

## List of 2
##  $ values : num [1:374] 4090900 845143 670540 256766 151349 ...
##  $ vectors: num [1:374, 1:374] 0.05213 -0.07137 -0.05728 0.00636 -0.05457 ...
##  - attr(*, "class")= chr "eigen"

先ほどと同様に、固有値から寄与率を計算することができます。

barplot(eig.G$values[1:10] / sum(eig.G$values))

寄与率と累積寄与率を図で表してみます。

cntr <- eig.G$values / sum(eig.G$values)
plot(cntr[1:10], ylim = c(0, 1), type = "h")
lines(1:10, cumsum(cntr[1:10]), col = "green")

第2回講義で説明したように、固有ベクトルと固有値の平方根の積を計算することで、主成分得点の計算ができます。

z.G <- eig.G$vectors[,1:100] %*% diag(sqrt(eig.G$values[1:100]))

計算された結果について、図示をしてみましょう。

plot(z.G[, 1:2], asp = 1)

明瞭なクラスタがあるように見えます。

\(k\)平均法を用いた非階層クラスタリング

そこで、非階層クラスタリング法の一つである\(k\)平均法を用いてクラスタ解析を行ってみます。 \(k\)平均法については、第3回講義で説明をしました。 kmeansという関数を用いると実行することができます。なお、k平均法は初期値依存性が高いため、100個の異なる初期値から繰り返して計算するようにしています。また、分類を行うクラスタ数は\(k=5\)としてあります。

res <- kmeans(z.G, centers = 5, nstart = 100)
str(res)

## List of 9
##  $ cluster     : int [1:374] 2 5 1 3 1 4 4 2 2 5 ...
##  $ centers     : num [1:5, 1:100] -121.5746 100.3248 0.0356 73.1081 -134.674 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "1" "2" "3" "4" ...
##   .. ..$ : NULL
##  $ totss       : num 8514372
##  $ withinss    : num [1:5] 508918 375307 207925 770146 1043227
##  $ tot.withinss: num 2905523
##  $ betweenss   : num 5608849
##  $ size        : int [1:5] 55 104 19 109 87
##  $ iter        : int 3
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

kmeans関数から得られた結果resのうちの、clusterというオブジェクト内に、 各系統について、分類されたクラスタ（5個のクラスタ）のIDが格納されています。 \(k\)平均法による分類結果を、主成分分析の結果を色付けして確認してみましょう。

plot(z.G[, 1:2], col = res$cluster, asp = 1)

同じ色の点（同じクラスタに分類された系統）がクラスタを形成していることが分かります。

第2、3主成分についても確認してみましょう。

plot(z.G[, 2:3], col = res$cluster, asp = 1)

plotlyパッケージの関数plot_lyを用いて3次元のグラフを描いてみます。

df <- data.frame(z.G[, 1:3], res$cluster)
colnames(df) <- c("PC1", "PC2", "PC3", "cluster")
plot_ly(data = df, x = ~PC1, y = ~PC2, z = ~PC3,
        color = ~cluster, type = "scatter3d", mode = "markers")

## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

2次元よりも、より直感的に主成分得点とクラスタの関係を把握できることがわかります。

なお、楕円フーリエ記述子の主成分分析の結果についても、 ゲノムワイドSNPの遺伝子型をもとにしたクラスタ解析の結果で色付けして表示してみましょう。

plot(z, col = res$cluster)

特に第1主成分において、クラスタ毎にある一定の傾向があることが分かります。

箱ひげ図を用いて、クラスタごとの主成分得点の分布を確認してみましょう。

boxplot(z[,1] ~ res$cluster)

boxplot(z[,2] ~ res$cluster)

とくに、第1主成分で、クラスタ間に主成分得点に違いがあることが分かります。

イネの種子形状の主成分が説明する変動の視覚化

では、第1主成分とは、何を評価する主成分なのでしょうか？ 主成分の意味を把握するためには、通常、biplotとよばれるグラフを描きます。 Rではprcompの結果をもとに、biplot関数を用いて描くことができます。

biplot(pca)

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length =
## arrow.len): zero-length arrow is of indeterminate angle and so skipped

変数\(d_1\)が第1主成分に大きな重みをもっていることが分かります。 \(d_1\)は、一次近似楕円の短軸の長さを表す変数で、第1主成分が種子の長幅比を評価していることが推察されます。

なお、主成分分析では、ある特定の主成分得点をもつ場合の、もとの変量の値を計算することができます。 この計算法について、第2回の講義で説明しました。

実際にR上で計算するには次のようにします。計算に必要なのは、固有値と固有ベクトルの情報です。 固有値は、対応する主成分の得点の分散を表します。 そこで、第1主成分の得点が、平均（主成分得点の平均は0）\(\pm 2\)倍の標準偏差の値をとるときの、 楕円フーリエ記述子を逆に求めてみましょう。

z.u <- 0 + 2 * sqrt(eig$values[1]) # 標準偏差なので平方根。平均は0。
z.l <- 0 - 2 * sqrt(eig$values[1])
x.u <- z.u * eig$vectors[,1]
x.l <- z.l * eig$vectors[,1]

なお、これだけでは、元の変数の値になりません。 なぜなら、主成分分析を行う前に、楕円フーリエ記述子について平均0になるように基準化してあるためです。 その際に、減算に使われた楕円フーリエ記述子の平均を、上で得られたx.uやx.lに加えてやる必要があります。 なお、scale関数を用いて基準化した場合には、基準化された結果の中に、平均の情報が保存されています。

それを取り出してみます。

str(X.centered)

##  num [1:374, 1:77] 0.10563 -0.00505 -0.02933 -0.14265 0.10519 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:374] "1" "3" "4" "5" ...
##   ..$ : chr [1:77] "d1" "a2" "b2" "c2" ...
##  - attr(*, "scaled:center")= Named num [1:77] 0.50554 0.00249 -0.00493 -0.00956 0.00389 ...
##   ..- attr(*, "names")= chr [1:77] "d1" "a2" "b2" "c2" ...

x.mean <- attr(X.centered, "scaled:center")

次に、この値をx.u, x.lに加えます。しかし、これだけでも輪郭を描くのに十分な情報が得られていません。 それは、定数であった\(a_1, b_1, c_1\)の情報が欠けているためです。 そこで、これらは定数としてベクトルの先頭に付け加えます。

ef.u <- c(1, 0, 0, x.mean + x.u)
ef.m <- c(1, 0, 0, x.mean      )
ef.l <- c(1, 0, 0, x.mean + x.l)

先ほど作成しておいた関数を用いて、特定の主成分得点をもつ輪郭座標を計算してみます。 これは、上で作成したef.u、ef.m、ef.lをもとに座標変換することで実行できます。

coord.u <- ef2coord(ef.u)
coord.m <- ef2coord(ef.m)
coord.l <- ef2coord(ef.l)

それでは、図として描いてみましょう。

plot(coord.u, type = "l", xlim = c(-1.2, 1.2), asp = 1, col = "green")
lines(coord.m, col = "gray")
lines(coord.l, col = "orange")

緑色の輪郭が、第1主成分得点が\(0 + 2 * \sigma\)（\(\sigma\)は主成分得点の標準偏差）のときに対応する輪郭、 緑色の輪郭が、第1主成分得点が平均（0）のときに対応する輪郭、 オレンジ色の輪郭が、第1主成分得点が\(0 - 2 * \sigma\)のときに対応する輪郭である。

同じようにして第2主成分についても、図を描いてみましょう。

z.u <- 0 + 2 * sqrt(eig$values[2]) # 標準偏差なので平方根。平均は0。
z.l <- 0 - 2 * sqrt(eig$values[2])
x.u <- z.u * eig$vectors[,2]
x.l <- z.l * eig$vectors[,2]
ef.u <- c(1, 0, 0, x.mean + x.u)
ef.l <- c(1, 0, 0, x.mean + x.l)
coord.u <- ef2coord(ef.u)
coord.l <- ef2coord(ef.l)
plot(coord.u, type = "l", xlim = c(-1.2, 1.2), asp = 1, col = "green")
lines(coord.m, col = "gray")
lines(coord.l, col = "orange")

さきほどと同じく、緑色の輪郭が、第2主成分得点が\(0 + 2 * \sigma\)（\(\sigma\)は主成分得点の標準偏差）のときに対応する輪郭、緑色の輪郭が、第2主成分得点が平均（0）のときに対応する輪郭、オレンジ色の輪郭が、第2主成分得点が\(0 - 2 * \sigma\)のときに対応する輪郭である。

イネのゲノムワイドSNP遺伝子型データに見られたクラスタを解釈する

なお、先ほどの解析で、ゲノムワイドSNPの多型の変異は、あるクラスタを形成しており、それは、楕円フーリエ記述子の主成分として評価された種子形状の変異との関連もあることが示されました。

ここからは、このクラスタが何を意味するものであるかを推察していきたいと思います。

まずは、各系統がもつパスポート情報を読み込みます。

passport <- read.csv(here("data", "RiceDiversityLine.csv"), stringsAsFactors = T)
head(passport)

##   GSOR.ID        IRGC.ID NSFTV.ID Accession.Name Country.of.origin  Latitude
## 1  301001 To be assigned        1       Agostano             Italy 41.871940
## 2  301003         117636        3  Ai-Chiao-Hong             China 27.902527
## 3  301004         117601        4       NSF-TV 4             India 22.903081
## 4  301005         117641        5       NSF-TV 5             India 30.472664
## 5  301006         117603        6       ARC 7229             India 22.903081
## 6  301007 To be assigned        7          Arias         Indonesia -0.789275
##   Longitude Sub.population     PC1     PC2     PC3     PC4
## 1  12.56738            TEJ -0.0486  0.0030  0.0752 -0.0076
## 2 116.87256            IND  0.0672 -0.0733  0.0094 -0.0005
## 3  87.12158            AUS  0.0544  0.0681 -0.0062 -0.0369
## 4  75.34424       AROMATIC -0.0073  0.0224 -0.0121  0.2602
## 5  87.12158            AUS  0.0509  0.0655 -0.0058 -0.0378
## 6 113.92133            TRJ -0.0293 -0.0027 -0.0677 -0.0085

このパスポートデータのうち、Sub.populationというのが、イネの系統の分類群（分集団）の情報です。 なお、パスポートデータの並びは、ゲノムワイドSNP遺伝子型データや、楕円フーリエ記述子データと異なるものになっています。 そこで、まずは、Sub.populationのデータを取り出し、それをゲノムワイドSNP遺伝子型データと対応付けます。

subpop <- passport$Sub.population
names(subpop) <- passport$NSFTV.ID
head(subpop)

##        1        3        4        5        6        7 
##      TEJ      IND      AUS AROMATIC      AUS      TRJ 
## Levels: ADMIX AROMATIC AUS IND TEJ TRJ

ゲノムワイドSNP遺伝子型データの行名を示してみます。

head(rownames(G))

## [1] "1" "3" "4" "5" "6" "7"

これをもとに、subpopの情報から必要な系統分を抜き出し、並び替えれば良さそうです。

subpop.gs <- subpop[rownames(ef.data)]
head(subpop.gs)

##        1        3        4        5        6        7 
##      TEJ      IND      AUS AROMATIC      AUS      TRJ 
## Levels: ADMIX AROMATIC AUS IND TEJ TRJ

できました。

では、抽出された分集団のデータsubpop.gsと、先ほど、非階層クラスタリングを\(k\)平均法で行った結果を table関数を用いて分割表を作成して比較してみましょう。

table(subpop.gs, res$cluster)

##           
## subpop.gs   1  2  3  4  5
##   ADMIX     3 17  5 24  7
##   AROMATIC  0  0 14  0  0
##   AUS      52  0  0  0  0
##   IND       0  0  0  0 80
##   TEJ       0 87  0  0  0
##   TRJ       0  0  0 85  0

先ほど\(k\)平均法を用いて得られたクラスタが、各系統の所属する分集団によく一致していることが分かります。 なお、ADMIXというグループは、これら分集団間の混合（例えば、異なる分集団からの系統を両親や最近の祖先にもつ） から生成された系統を集めたグループですが、所属する系統が5つのクラスタにわたって所属していることがわかります。

イネのゲノムワイドSNP遺伝子型データの階層的クラスタリング

まず、ゲノムワイドSNPの遺伝子型データに基づき、階層的クラスタリングを行ってみます。 これは、第3回講義の内容です。系統間の距離をSNPの遺伝子型スコアをもとにしたユークリッド距離として計算します。

なお、次のコードを用いると計算に時間を要します。

#G.dist <- dist(G)

そこで、第2回の講義で説明した式を用いて、内積行列から距離行列を計算します。

dij <- sqrt(matrix(diag(GG), nrow(GG), ncol(GG)) + matrix(diag(GG), nrow(GG), ncol(GG), byrow = T) - 2 * GG)
G.dist <- as.dist(dij)
rm(dij)

次に、この距離情報をもとにして、完全連結（complete linkage、最遠隣（furthest neighbor））法を用いて、階層的クラスタリングを行います。

hc <- hclust(G.dist)
plot(hc, hang = -1, cex = 0.5)

図を見やすくするため、また、デンドログラム（樹形図）の末端ノードに分集団に応じた色付けをするために、 apeパッケージのas.phylo関数と、それに応じたplot関数を用いてデンドログラムを描いてみます。

plot(as.phylo(hc), cex = 0.5, tip.color = as.numeric(subpop.gs))

なお、群内分散と群間分散の比を最大にするようにクラスタを形成していくWard法を用いて階層的クラスタ解析を行い、 その結果を図示してみます。

hc.ward <- hclust(G.dist, method = "ward.D2")
plot(as.phylo(hc.ward), cex = 0.5, tip.color = as.numeric(subpop.gs))

分集団が非常に明瞭に分かれていることが分かります。

イネのゲノムワイドSNP遺伝子型データに基づく分類器の学習

最後に、機械学習法の一つであるサポートベクターマシンを用いて、教師付き学習による分類を行ってみます。 その目的は、ADMIX以外の分集団で学習を行い、ADMIXの系統を学習した分集団に割り振ることです。 最終的には、各分集団への所属確率を計算し、それをもとに、ADMIXの各系統が、 どのような遺伝的背景をもつのか推察できるようにしたいと思います。

まず、ゲノムワイドSNP遺伝子型データを入力として、分集団に分類するモデルを構築してみましょう。 モデル構築には、kernlabのksvm関数を用います。 ここで、kernelとしてrbrdotが指定されていますが、これは動径基底関数（radial basis function: RBF）とよばれる非線形のカーネルです。

なお、このあとのSVMによる解析は、入力にGをそのまま用いると、メモリ不足からかRStudio Cloudでの計算が異常終了してしまうようでした。 そこで、代わりに上位100主成分の主成分得点を入力として用いることにします。

model <- ksvm(x = z.G, y = subpop.gs, kernel = "rbfdot", scale = F)

次に、作成したモデルをもとに、predict関数を用いて予測をしてみます。 ただ、ここでは本当は予測ではありません。 モデル構築に用いたデータを入力にしているのでモデルの当てはめ（fitting）です。

subpop.pred <- predict(model, z.G)
table(subpop.pred, subpop.gs)

##            subpop.gs
## subpop.pred ADMIX AROMATIC AUS IND TEJ TRJ
##    ADMIX       53        0   0   0   0   1
##    AROMATIC     0       14   0   0   0   0
##    AUS          0        0  52   0   0   0
##    IND          0        0   0  80   0   0
##    TEJ          3        0   0   0  87   0
##    TRJ          0        0   0   0   0  84

その結果、ADMIXという分類群が別にできてしまっていることが分かります。

そこで、ADMIXに属する系統を訓練用データ（機械学習でモデルを作るためのデータ）から除いてモデルを作り、 そのモデルにADMIXを当てはめて、その遺伝的背景を予測してみます。

まず、ADMIXを訓練データから除きます。 入力であるゲノムワイドSNP遺伝子型データだけからではなく、 出力である分集団データについても、ADMIXを除いたデータを作成します。

selector <- subpop.gs != "ADMIX"
z.G.woadm <- z.G[selector, ]
subpop.gs.woadm <- as.factor(as.character(subpop.gs[selector])) # 一度文字列にしてから、新たにfactorにしている

次に、ADMIXだけのデータを準備します。

z.G.adm <- z.G[!selector, ]

このデータは、モデルを訓練した後で用います。

では、分類モデルを訓練データから作成してみましょう。

model <- ksvm(x = z.G.woadm, y = subpop.gs.woadm, kernel = "rbfdot", scale = F)

このデータに、訓練データを当てはめてみます。 これも、モデルを訓練したデータをそのままモデルに当てはめているので、 予測ではありません。

subpop.pred <- predict(model, z.G.woadm)
table(subpop.pred, subpop.gs.woadm)

##            subpop.gs.woadm
## subpop.pred AROMATIC AUS IND TEJ TRJ
##    AROMATIC       14   0   0   0   0
##    AUS             0  52   0   0   0
##    IND             0   0  80   0   0
##    TEJ             0   0   0  87   0
##    TRJ             0   0   0   0  85

最後に、ADMIXのデータを訓練データから得られたモデルに当てはめて、予測をしてみます。

adm.pred <- predict(model, z.G.adm)
head(adm.pred)

## [1] TRJ AUS TRJ TEJ TEJ IND
## Levels: AROMATIC AUS IND TEJ TRJ

table(adm.pred)

## adm.pred
## AROMATIC      AUS      IND      TEJ      TRJ 
##        0        1       15       17       23

その結果、予測はされたのですが、ADMIXの各系統が1つ分集団に分類されてしまい、 その遺伝的背景（実際には、混合された系統）を推察することができません。

ksvm関数は、実は、分類確率を計算するモデルを構築することができます。 そこで、改めて、分類確率の計算を可能にするモデルを構築してみます。 prob.modelというオプションをT(TRUE)にするだけです。

model <- ksvm(x = z.G.woadm, y = subpop.gs.woadm, kernel = "rbfdot", prob.model = T)

構築されたモデルに、ADMIXの系統のデータを当てはめてみましょう。

adm.pred <- predict(model, z.G.adm, type = "probabilities")
head(adm.pred)

##        AROMATIC        AUS        IND        TEJ        TRJ
## [1,] 0.02388834 0.07680274 0.04874128 0.17441726 0.67615038
## [2,] 0.04647561 0.30852647 0.52221657 0.06185483 0.06092652
## [3,] 0.02743097 0.19823776 0.60272421 0.06064153 0.11096553
## [4,] 0.07794615 0.14472813 0.16371629 0.34389459 0.26971484
## [5,] 0.05676639 0.21826995 0.19920226 0.15699817 0.36876323
## [6,] 0.05435664 0.29826162 0.50596836 0.07245852 0.06895485

今度は、確率として計算されていることが分かります。 先ほどの分類では、各系統が、この確率が最も大きくなる分集団に分類されていました。

最後に、ADMIXの各系統の遺伝的背景を表すと考えられる、分集団への所属確率を図として描いてみましょう。

barplot(t(adm.pred), col = rainbow(5), horiz = T)

このような図を描くことで、各系統のもつ遺伝的背景を、視覚的に把握できます。