バイオスタティスティクス基礎論 第4回講義テキスト
Hiroyoshi Iwata
2022-05-06
階層的クラスタ解析
多数の対象について、それらのもつ多次元の特徴をもとに「似たもの」どうしをグループ（クラスタ cluster）に分類すると便利なことがあります。例えば、DNA多型のデータに基づき遺伝資源に含まれる品種や系統をグループ分けできれば、遺伝資源のもつ形質の変異について、グループの情報をもとに整理し、体系化することができます。 前回の講義でもお話ししたように、多数のサンプルがもつ多数の特徴の変異について、データを眺めるだけで把握するのは困難です。主成分分析では、多数の特徴を低次元の変数で表現することでデータのもつ変異の要約を試みました。クラスタ解析では、多数のデータを少数のグループにまとめることで、データのもつ変異の要約を試みます。今回の講義では、まず、多数のデータを階層的にグループに分類する階層的クラスタ解析について概説します。
今回の講義では、前回までと同様にイネのデータ（Zhao et al. 2011）を用いて説明を進めていきます。今回の講義では、品種・系統データ（RiceDiversityLine.csv）、表現型データ（RiceDiversityPheno.csv）、マーカー遺伝子型データ（RiceDiversityGeno.csv）の3つのデータを用います。いずれも、Rice Diversityのwebページ http://www.ricediversity.org/data/index.cfm からダウンロードして整理したデータです。前回も説明したようにマーカー遺伝子型データは、ソフトウエアfastPHASE（Scheet and Stephens 2006）を用いて欠測値の補間を行ってあります。 まずは、前回と同様に、3種類のデータを読み込んで、それらを結合してみましょう。
# this data set was analyzed in Zhao 2011 (Nature Communications 2:467)
line <- read.csv("RiceDiversityLine.csv", stringsAsFactors = T)
pheno <- read.csv("RiceDiversityPheno.csv", stringsAsFactors = T)
geno <- read.csv("RiceDiversityGeno.csv", stringsAsFactors = T)
line.pheno <- merge(line, pheno, by.x = "NSFTV.ID", by.y = "NSFTVID")
alldata <- merge(line.pheno, geno, by.x = "NSFTV.ID", by.y = "NSFTVID")
rownames(alldata) <- alldata$NSFTV.ID最初に、DNAマーカー（1,311 SNPs）に見られた変異に基づいて、374品種・系統をクラスタに分類してみましょう。まずは、そのためのデータを準備します。
#### analysis of marker data
data.mk <- alldata[, 50:ncol(alldata)]
subpop <- alldata$Sub.population
dim(data.mk)## [1] 374 1311クラスタ解析には様々な方法がありますが、ここではまず1つの方法でクラスタ解析を行ってみます。 まずは、DNAマーカーのデータをもとに、品種・系統間の距離を計算します。
# calculate Euclid distance
d <- dist(data.mk)
head(d)## [1] 54.47141 53.08033 44.70547 52.82571 45.40700 44.36904as.matrix(d)[1:6,1:6]## 1 3 4 5 6 7
## 1 0.00000 54.47141 53.08033 44.70547 52.82571 45.40700
## 3 54.47141 0.00000 37.53194 46.79940 37.68502 49.82169
## 4 53.08033 37.53194 0.00000 44.38481 17.58133 46.49073
## 5 44.70547 46.79940 44.38481 0.00000 43.85254 42.87989
## 6 52.82571 37.68502 17.58133 43.85254 0.00000 46.69070
## 7 45.40700 49.82169 46.49073 42.87989 46.69070 0.00000なお、関数distの返す値は行列（matrix）形式でなく、距離行列特有の形式になっていることに注意して下さい。したがって、例えば、最初の6品種間の総当たりの距離を6×6行列で表示したい場合には、上記のように関数as.matrixで距離行列特有の形式からmatrix形式に変換する必要があります。 では、クラスタ解析を行ってみましょう。
tre <- hclust(d)
tre##
## Call:
## hclust(d = d)
##
## Cluster method : complete
## Distance : euclidean
## Number of objects: 374Callには、回帰分析などと同様に、実行したコマンドがそのまま表示されます。また、Cluster methodにはクラスタ解析の方法（クラスタ間の距離の定義）、Distanceには距離の計算法が表示されます。また、Number of objectsは分類を行った対象（ここでは、品種・系統）の数です。 クラスタ解析の結果を樹形図（dendrogram）で表示してみましょう。
# draw dendrogram
plot(tre)
Fig.1 Dendrogram of 374 varieties / lines obtained based on marker genotype data
Figure 1は関数hclustで得られた結果をそのまま樹形図にしたものです。パッケージapeを用いると、様々な表現様式で樹形図を描くことができます。そのためには、まず関数hclustで得られた結果をパッケージapeで定義されているphyloとよばれるクラスに変換する必要があります。
# convert to a phylo object defined in the ape package
phy <- ape::as.phylo(tre)では、phyloクラスに変換された結果をプロットしてみましょう。
# plot as a phylo object
plot(phy)
Fig.2 Dendrogram converted to the phylo class of the package ape
Figure 2は、品種・系統数が多いこともあり、非常に見にくい図になっています。各品種・系統の遺伝的背景（所属する分集団）と樹形図での位置の関係を枝の色で確認できるようにして、少し見やすい図に描き換えてみましょう。
head(phy$edge)## [,1] [,2]
## [1,] 375 376
## [2,] 376 380
## [3,] 380 236
## [4,] 380 392
## [5,] 392 209
## [6,] 392 334head(subpop[phy$edge[,2]], 10)## [1] <NA> <NA> ADMIX <NA> ADMIX ADMIX <NA> <NA> <NA> <NA>
## Levels: ADMIX AROMATIC AUS IND TEJ TRJcol <- as.numeric(subpop[phy$edge[,2]])
edge.col <- ifelse(is.na(col), "gray", col)
plot(phy, edge.color = edge.col, show.tip.label = F)
Fig.3 Dendrogram colored for each population group of varieties and lines
Figure 3を見ると、同じ分集団に含まれる品種・系統が同じクラスタに含まれる傾向が確認でき、品種・系統のもつ遺伝的背景の違いがクラスタ解析の結果によく反映していることが分かります。 パッケージapeのphyloクラスは、様々な表現の仕方で樹形図を描くことができます。異なるタイプの樹形図を試してみましょう。
# different types of dendrogram
#pdf("fig4.pdf", width = 10, height = 10)
op <- par(mfrow = c(2, 2), mar = rep(0, 4))
plot(phy, edge.color = edge.col, type = "phylogram", show.tip.label = F)
plot(phy, edge.color = edge.col, type = "cladogram", show.tip.label = F)
plot(phy, edge.color = edge.col, type = "fan", show.tip.label = F)
plot(phy, edge.color = edge.col, type = "unrooted", show.tip.label = F)
Fig.4 Various styles of dendrograms drawn using package ape
par(op)
#dev.off()Figure 4に出力されている図は、同じクラスタ解析の結果を異なる様式の樹形図で描いたものです。様式が異なると受ける印象も分かりやすさも異なります。品種・系統の遺伝的関係を大域的に把握したい場合には、4番目の“unrooted”タイプの樹形図が最も目的に適っているのではないかと思われます。 クラスタ解析の結果についてパッケージapeを利用して図にする手順は、途中にphyloクラスへの変換などを必要とするため、少し面倒です。そこで、一連の作業を自作の関数として定義して、クラスタ解析の結果の図示を簡略化してみましょう。
# create an own function
myplot <- function(tre, subpop, type = "unrooted", ...) {
phy <- as.phylo(tre)
col <- as.numeric(subpop[phy$edge[,2]])
edge.col <- ifelse(is.na(col), "gray", col)
plot(phy, edge.color = edge.col, type = type, show.tip.label = F, ...)
}では、自作の関数myplotを使って樹形図を描いてみましょう。
# use the function
d <- dist(data.mk)
tre <- hclust(d)
myplot(tre, subpop)
距離の定義
クラスタ解析では、サンプル間やクラスタ間の距離を計算し、計算された距離に基づいてクラスタリングを行います。したがって、距離の定義が異なると異なる結果が得られることになります。ここでは、サンプル間やクラスタ間の距離の定義について解説します。
まずは、サンプル間の距離についてです。サンプル間の距離を計算するのに、様々な定義があります。まずは、異なる定義の距離に基づいて樹形図を描いてみましょう。
# try different methods for calculating distance
#pdf("fig5.pdf", width = 10, height = 10)
op <- par(mfrow = c(2, 2), mar = rep(0, 4))
d <- dist(data.mk, method = "euclidean") # default method
myplot(hclust(d), subpop)
d <- dist(data.mk, method = "manhattan")
myplot(hclust(d), subpop)
d <- dist(data.mk, method = "minkowski", p = 1.5)
myplot(hclust(d), subpop)
d <- as.dist(1 - cor(t(data.mk)))
myplot(hclust(d), subpop)
Fig.5 Dendrograms calculated based on different definitions of distances between samples
par(op)
#dev.off()今回のデータでは距離の定義が異なっても樹形図のトポロジー（topology）は大きく変わりませんが、データによっては距離の定義が大きく影響する場合があります。
上で用いたサンプル間の距離について、その定義を以下に示します。なお、各サンプルが\(q\)個の特徴で記述されており、\(i\)番目のサンプルのデータベクトルを\(\mathbf{x}_i=(x_{i1},...,x_{iq} )^T\)、\(j\)番目のサンプルのデータベクトルを\(\mathbf{x}_j=(x_{j1},...,x_{jq} )^T\)と表すこととします。このとき、サンプル\(i, j\)間の距離\(d(\mathbf{x}_i,\mathbf{x}_j)\)は、以下のように定義されます。
- ユークリッド（Euclidian）距離 \[ d(\mathbf{x}_i,\mathbf{x}_j) = \sqrt{\sum_{k=1}^q (x_{ik} - x_{jk})^2} \]
- マンハッタン（Manhattan）距離 \[ d(\mathbf{x}_i,\mathbf{x}_j) = \sum_{k=1}^q \left|x_{ik} - x_{jk}\right| \]
- ミンコフスキー（Minkowski）距離 \[ d(\mathbf{x}_i,\mathbf{x}_j) = \left(\sum_{k=1}^q \left|x_{ik} - x_{jk}\right|^p\right)^{1/p} \]
- 相関係数に基づく距離 \[ d(\mathbf{x}_i,\mathbf{x}_j) = r_{ij} = 1 - \frac{\sum_{k=1}^q (x_{ik} - \bar x_i)(x_{jk} - \bar x_j)} {\sqrt{\sum_{k=1}^q (x_{ik} - \bar x_i)^2}\sqrt{\sum_{k=1}^q (x_{jk} - \bar x_j)^2}} \] ここで \(\bar x_i = \frac{1}{n} \sum_{k=1}^q x_{ik}\)
マンハッタン距離は、ニューヨーク市のManhattanのような正方形に区分された市街地を移動する場合の距離というのがその名の由来です。そのような市街地では、例えば、地点\((0,0)\)から地点\((2,3)\)に移動する場合に、建物があるために斜めに移動（ユークリッド距離\(\sqrt{13}\)することができず、道に沿って移動（マンハッタン距離5）する必要があるためです。ミンコフスキー距離は、ユークリッド距離とマンハッタン距離の一般化されたかたちです。\(p=1\)のときはマンハッタン距離、\(p=2\)のときはユークリッド距離に一致します。
相関係数に基づく距離では、「変数間ではなくサンプル間で」相関係数\(r\)を計算して、それを\(1\)から減じたものを距離とします。相関が\(r=1\)の場合は距離\(1-r=0\)、相関が\(r=0\)のときは距離\(1-r=1\)、相関が\(r=−1\)のときには距離が\(1=r=2\)となります。すなわち、相関係数に基づく距離では最大値が2となります。なお、遺伝子間で発現パターンの類似性からクラスタ解析を行う場合には、\(1\)から相関係数\(r\)を減ずる代わりに相関係数の絶対値\(|r|\)を減ずる場合があります。この場合、相関が\(r=−1\)または\(r=1\)のとき（\(|r|=1\)）には距離\(1-|r|=0\)、相関係数がすなわち\(r=|r|=0\)のときに最も距離が遠くなり\(1-|r|=1\)となります。
関数distでは、次のような距離も計算できます。今回のデータには不向きであったので利用しませんでしたが、解析するデータの性質によっては、以下に紹介する距離が適切な場合もあります。
- チェビシェフ（Chebyshev）距離 （関数distでmethod=“maximum”を指定） \[d(\mathbf{x}_i,\mathbf{x}_j)=\max_k\left(|x_{ik}-x_{jk}|\right)\]
- キャンベラ（Canberra）距離 （関数distでmethod=“canberra”を指定） \[d(\mathbf{x}_i,\mathbf{x}_j)=\sum_{k=1}^q\frac{|x_{ik}-x_{jk}|}{|x_{ik}|+|x_{jk}|}\]
- ハミング距離（Hamming）距離 （関数distでmethod=“binary”を指定） \[d(\mathbf{x}_i,\mathbf{x}_j)=\sum_{k=1}^q(1-\delta_{x_{ik},x_{jk}})\] ここで、 \[ \delta_{a, b}= \begin{cases} 1&\text{if } a=b\\ 0 &\text{else} \end{cases} \]
チェビシェフ距離は\(q\)個の特徴のうち最も異なっている1つの特徴の違いだけに基づく距離です。この距離は、ミンコフスキー距離の\(p\rightarrow\infty\)の極限となっています。ハミング距離は情報科学でよく用いられる距離で、同じ長さの数列について、同じ位置の値を比較したときに一致しない位置の数を数えあげたものです。ハミング距離を用いるようなデータでは、\(x_{ik}\)は、連続値ではなく、離散値\((0,1)\)である場合がほとんどです。
ここまではサンプル間の距離の定義について説明してきました。階層的クラスタ解析では、距離の近いサンプルどうしを１つのクラスタにまとめながら、さらに、サンプルとクラスタ、または、クラスタどうしを、上位の階層のクラスタとしてまとめあげていきます。したがって、サンプル間の距離だけでなく、サンプルとクラスタ、または、クラスタ間の距離を定義しておく必要があります。
ここでは、まず、クラスタ間距離の様々な定義に基づいて樹形図を描いてみます。関数hclustでは、クラスタ間の距離の計算方法（定義）をオプションmethodで指定することができます。
#pdf("fig6.pdf", width = 10, height = 10)
d <- dist(data.mk)
op <- par(mfrow = c(2, 3), mar = rep(0, 4))
tre <- hclust(d, method = "complete") # default method
myplot(tre, subpop)
tre <- hclust(d, method = "single")
myplot(tre, subpop)
tre <- hclust(d, method = "average")
myplot(tre, subpop)
tre <- hclust(d, method = "median")
myplot(tre, subpop)
tre <- hclust(d, method = "centroid")
myplot(tre, subpop)
tre <- hclust(d, method = "ward.D2")
myplot(tre, subpop)
Fig.6 Dendrogram based on various inter-cluster distance definitions
par(op)
#dev.off()Figure 6を見ると、クラスタ間距離の定義の違いは、サンプル間距離の定義の違いと異なり、樹形図のトポロジーが大きく変化することが分かります。中には、枝長が負の値になりおかしな樹形図になっている場合もあります（左下、下中央）。また、クラスタ間の違いが非常に強調される場合もあります（右下）。この中からどの手法を選択するかは難しい問題ですが、多くの場合、既知の情報と大きく矛盾が無いものが選ばれます。例えば、ここでは、品種・系統が所属している分集団と矛盾が小さいものを選ぶとよいでしょう。
関数hclustで指定できるクラスタ間の距離の定義を示します。サンプル間の距離\(d(\mathbf{x}_i,\mathbf{x}_j)\)に基づき、クラスタAとBの距離を\(d_{AB}\)は以下のように計算されます。
- 最長距離法（完全連結法） （関数hclustでmethod=“complete”を指定） \[ d_{AB}=\underset{i\in A,j\in B}{\max} d(\mathbf{x}_i,\mathbf{x}_j) \]
- 最短距離法（単連結法） （関数hclustでmethod=“single”を指定） \[ d_{AB}=\underset{i\in A,j\in B}{\min} d(\mathbf{x}_i,\mathbf{x}_j) \]
- 平均距離法 （関数hclustでmethod=“average”を指定） \[ d_{AB}=\frac{1}{n_A n_B} \sum_{i\in A}\sum_{j\in B}d(\mathbf{x}_i,\mathbf{x}_j) \] ここで、\(n_A, n_B\)は、クラスタA, Bに含まれるサンプル数を表す。
以下の3つの定義では、クラスタA, Bが融合して新しくクラスタCができるときに、新しいクラスタCとA,B以外のクラスタO間の距離\(d_{CO}\)を次のように定義する。なお、クラスタAとBの距離を\(d_{AB}\)、クラスタAとOの距離を\(d_{AO}\)、クラスタBとOの距離を\(d_{BO}\)、クラスタA、B、Oに含まれるサンプルの数を\(n_A\), \(n_B\), \(n_O\)と表す。
- 重心法 （関数hclustでmethod=“centroid”を指定） \[ d_{CO}^2=\frac{n_A}{n_A+n_B} d_{AO}^2+\frac{n_B}{n_A+n_B} d_{BO}^2-\frac{n_A n_B}{(n_A+n_B)^2} d_{AB}^2 \]
- メディアン法 （関数hclustでmethod=”median”を指定） \[ d_{CO}=\frac{1}{2} d_{AO}+\frac{1}{2} d_{BO}-\frac{1}{4} d_{AB} \]
- ウォード（Ward）法 （関数hclustでmethod=”ward.D2”を指定） \[ d_{CO}^2=\frac{n_A+n_O}{n_A+n_B+n_O} d_{AO}^2 +\frac{n_B+n_O}{n_A+n_B+n_O} d_{BO}^2-\frac{n_O}{n_A+n_B+n_O} d_{AB}^2 \]
Figure 6において分集団との対応が明瞭と思われる2つの手法について、もう少し詳細に比較してみましょう。
# focus on two clustering methods
op <- par(mfrow = c(1, 2), mar = rep(0, 4))
d <- dist(data.mk)
tre <- hclust(d, method = "complete")
myplot(tre, subpop, type = "phylogram")
tre <- hclust(d, method = "ward.D2")
myplot(tre, subpop, type = "phylogram")
Fig.7 Dendrograms of the longest distance method (left) and Ward method (right)
par(op)多次元データの両側からのクラスタ解析
ここまでは、DNAマーカーデータをもとに品種や系統をクラスタに分類しました。品種や系統のクラスタ解析は、DNAマーカーデータだけでなく、形質データをもとにしても行うことができます。また、全く同じデータについて、品種・系統ではなく、形質を分類する対象とみなして、品種・系統間で似たような変異のパターンをもつ形質どうしを同じクラスタに分類することもできます。ここでは、このようなアプローチについて説明を行います。
まずは、形質データを準備しましょう。全データ（alldata）から、このような解析に適さない形質を除き、形質データを抜き出しましょう。
# preparation of data
required.traits <- c("Flowering.time.at.Arkansas",
"Flowering.time.at.Faridpur", "Flowering.time.at.Aberdeen",
"Culm.habit", "Flag.leaf.length", "Flag.leaf.width",
"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","Brown.rice.seed.length", "Brown.rice.seed.width",
"Straighthead.suseptability","Blast.resistance",
"Amylose.content", "Alkali.spreading.value", "Protein.content")
data.tr <- alldata[, required.traits]
missing <- apply(is.na(data.tr), 1, sum) > 0
data.tr <- data.tr[!missing, ]
subpop.tr <- alldata$Sub.population[!missing]形質データは、形質によって変動の大きさ（分散）が異なります。このデータをそのまま用いると、分散の大きな形質は距離の計算に大きな影響を与え、分散の小さな形質は距離の計算への寄与が小さくなります。そこで、全ての形質について、分散1に基準化しておきます。
# scaling
data.tr <- scale(data.tr)では、形質データについて、品種・系統を分類の対象としたクラスタ解析と、形質を分類の対象としたクラスタ解析を行ってみましょう。
# perform clusetering for both varieties and traits
d <- dist(data.tr)
tre.var <- hclust(d, method = "ward.D2")
d <- dist(t(data.tr))
tre.tra <- hclust(d, "ward.D2")
op <- par(mfrow = c(1, 2))
myplot(tre.var, subpop.tr, type = "phylogram")
plot(tre.tra, cex = 0.5)
Fig.8 Cluster analysis based on trait data. Dendrograms showing the relationship between varieties and lines (left) and the relationship between traits (right)
par(op)Figure 8の右側の樹形図から、植物体のサイズに関わる形質（Plant.height, Panicle.length, Flag.leaf.length）のは互いに関係が強いことが分かります。また、そのクラスタの近くに、3つの環境で計測された開花のタイミング（Flowering.time.at.*****）が位置していることも分かります。また、止め葉の幅（Flag.leaf.width）は、他のサイズ関連形質と異なり、穂の特徴を表す形質との関連が強いことも分かります。このように、多次元データはどちら側からもクラスタ解析を行うことができます。このことを覚えておくと、同じデータを少し違った視点から眺めることができるでしょう。
なお、上述した解析は、関数heatmapを用いてより視覚的に結果を表示できます。
# perform clustring from both sides
#pdf("fig9.pdf")
heatmap(data.tr, margins = c(12,2))
Fig.9 Display of cluster analysis results and heatmap of trait data
#dev.off()ヒートマップを描く関数として、gplotsパッケージに含まれるheatmap.2という関数もあります（他にもたくさんの関数があると思います）。こちらを用いてヒートマップを描いてみます。オプションの与え方が少し異なりますが、見栄えの良い図が描けます。
#pdf("fig9-2.pdf", height = 12)
gplots::heatmap.2(data.tr, margins = c(12,2), col=redgreen(256), trace = "none", lhei = c(2,10), cexRow = 0.3)
Fig.9-2 Display heat map of trait data using heatmap.2 function
#dev.off()以下のようにすると、別に行ったクラスタ解析の結果を反映させることができます。先ほど、同データについて行ったクラスタ解析の結果をヒートマップ表示に反映させてみましょう。
# perform clustering with appropriate methods
#pdf("fig10.pdf")
heatmap(data.tr, Rowv = as.dendrogram(tre.var),
Colv = as.dendrogram(tre.tra),
RowSideColors = as.character(as.numeric(subpop.tr)),
labRow = subpop.tr,
margins = c(12, 2))
Fig.10 The result of changing the cluster analysis method in Figure 9 to Ward’s method
#dev.off()こちらについても、先程と同じくheatmap.2関数を用いて描くこともできます。
# this part is again optional
#pdf("fig10-2.pdf", height = 12)
heatmap.2(data.tr, Rowv = as.dendrogram(tre.var),
Colv = as.dendrogram(tre.tra),
RowSideColors = as.character(as.numeric(subpop.tr)),
labRow = subpop.tr,
margins = c(12,2), col=redgreen(256), trace = "none", lhei = c(2,10), cexRow = 0.3)
Fig.10-2 The result of changing the cluster analysis method of Figure 9-2 to Ward’s method
#dev.off()関数heatmapでは、必ずしも縦横同じデータでクラスタ解析をする必要はありません。例えば、行側について、DNAマーカーデータを用いたクラスタ解析の結果をあてはめることもできます。
# perform clustering based on marker genotypes for determining row order
data.mk2 <- data.mk[!missing, ]
d <- dist(data.mk2)
tre.mrk <- hclust(d, method = "ward.D2")
#pdf("fig11.pdf")
heatmap(data.tr, Rowv = as.dendrogram(tre.mrk),
Colv = as.dendrogram(tre.tra),
RowSideColors = as.character(as.numeric(subpop.tr)),
labRow = subpop.tr,
margins = c(12, 2))
Fig.11 Relationship between the results of cluster analysis based on genetic marker data and traits
#dev.off()こちらについても、先程と同じくheatmap.2関数を用いて描くこともできます。
# this part is optional
#pdf("fig11-2.pdf", height = 12)
heatmap.2(data.tr, Rowv = as.dendrogram(tre.mrk),
Colv = as.dendrogram(tre.tra),
RowSideColors = as.character(as.numeric(subpop.tr)),
labRow = subpop.tr,
margins = c(12,2), col=redgreen(256), trace = "none", lhei = c(2,10), cexRow = 0.3)
Fig.11-2 Relationship between the results of cluster analysis and traits based on genetic marker data expressed using heatmap.2 function
#dev.off()階層的クラスタ解析に基づく分類
サンプル間の類似性についてあるところで線引きし、サンプルを離散的にグループに分類したい場合があります。ここでは、階層的クラスタ解析の結果に基づいてサンプルを決められた数のクラスタに分類する方法について説明します。
DNAマーカーデータに基づく階層的クラスタ解析の結果に基づき、品種・系統を5つのグループに分類してみましょう。5というのは、品種・系統の所属する分集団の数に合わせた数字です。階層的クラスタ解析の結果から、離散的なグループを求めるには関数cutreeを用います。
# classify samples with the cutree function
d <- dist(data.mk)
tre <- hclust(d, method = "ward.D2")
cluster.id <- cutree(tre, k = 5)
head(cluster.id, 10)## 1 3 4 5 6 7 8 9 10 11
## 1 2 3 4 3 5 5 1 1 2クラスタ解析に基づき5つのグループに分類した結果を図示してみましょう。
# visualize the result
op <- par(mfrow = c(1,2), mar = rep(0, 4))
myplot(tre, cluster.id, type = "phylogram")
myplot(tre, subpop, type = "phylogram", direction = "leftwards")
Fig.12 Relationship between classification by cluster analysis (left) and divided groups (right)
par(op)クラスタ解析に基づく分類と分集団の関係を、クロス集計表を作成して確認してみましょう。
# obtain the cross table of classification
table(cluster.id, subpop)## subpop
## cluster.id ADMIX AROMATIC AUS IND TEJ TRJ
## 1 5 0 0 0 84 0
## 2 14 0 0 80 0 0
## 3 1 0 52 0 0 0
## 4 2 14 0 0 0 0
## 5 34 0 0 0 3 85両者は、インディカ（IND）の3品種・系統を除いて、非常によく一致していることが分かります。これは、分集団構造そのものがDNAマーカーデータに基づいて推定されたためだと考えられます。なお、複数の分集団の混合（ADMIX）と推定されている品種については様々なグループに分類されていることも分かります。
階層的クラスタ解析による分類の結果を主成分軸上で確認してみましょう。
# visualize the result in the PCA space
pca <- prcomp(data.mk)
op <- par(mfrow = c(1,2))
plot(pca$x[,1:2], pch = cluster.id, col = as.numeric(subpop))
plot(pca$x[,3:4], pch = cluster.id, col = as.numeric(subpop))
Fig.13 Classification by cluster analysis and relationship between subgroups
par(op)非階層的クラスタリング
ある決められたグループ数に分類する場合には、階層的に分類を行う必要は必ずしもありません。ここでは、非階層的クラスタ解析手法の一つであるk-平均（k-means）法を紹介します。
先ほどと同じデータについて、関数kmeansを用いて5つのグループへの分類を行ってみましょう。
# kmeans clustering
kms <- kmeans(data.mk, centers = 5)
#kms
head(kms$cluster, 20)## 1 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 20 22 24
## 4 5 3 1 3 1 1 2 4 5 1 3 4 1 5 3 3 1 1 1k-平均法では、以下のようなアルゴリズムで決められたグループ数への分類を行います。
- k個のクラスタ中心として、k個のサンプルを無作為に選びだす
- すべてのデータ点とk個のクラスタ中心間の距離を求め、各データ点を中心（重心を中心とする）が最も近いクラスタに分類する
- 形成されたクラスタの中心（重心）を更新する
- クラスタの中心（重心）が変化しなくなるまで、2-3を繰り返す
k-平均法では、最初に無作為に選ばれるサンプルによって結果が変化することがあります。実際に、同じデータで解析を繰り返して、結果のバラツキを確認してみましょう。
# repeat kmeans clustering
for(i in 1:5) {
kms <- kmeans(data.mk, centers = 5)
print(table(kms$cluster, subpop))
}## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 10 0 0 79 0 0
## 2 3 14 0 0 0 0
## 3 23 0 0 0 1 85
## 4 17 0 0 0 86 0
## 5 3 0 52 1 0 0
## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 10 0 0 79 0 0
## 2 3 0 52 1 0 0
## 3 23 0 0 0 1 85
## 4 3 14 0 0 0 0
## 5 17 0 0 0 86 0
## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 4 0 0 25 0 0
## 2 10 0 0 24 0 0
## 3 1 0 0 31 0 0
## 4 40 14 0 0 87 85
## 5 1 0 52 0 0 0
## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 14 0 0 10 0 0
## 2 1 0 0 70 0 0
## 3 16 0 0 0 86 0
## 4 24 14 0 0 1 85
## 5 1 0 52 0 0 0
## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 10 0 0 0 1 31
## 2 13 0 0 0 0 54
## 3 3 14 0 0 0 0
## 4 17 0 0 0 86 0
## 5 13 0 52 80 0 0実行するごとに、異なる結果が得られます。これは、上述したアルゴリズムが初期値依存性が高いことによります。そこで、異なる初期値から何度か繰り返し計算を行ってみます。ここでは、50の異なる初期値をもとに計算をします。
# start from multiple sets of initial points
for(i in 1:5) {
kms <- kmeans(data.mk, centers = 5, nstart = 50)
print(table(kms$cluster, subpop))
}## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 23 0 0 0 1 85
## 2 3 14 0 0 0 0
## 3 17 0 0 0 86 0
## 4 1 0 52 0 0 0
## 5 12 0 0 80 0 0
## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 12 0 0 80 0 0
## 2 23 0 0 0 1 85
## 3 3 14 0 0 0 0
## 4 17 0 0 0 86 0
## 5 1 0 52 0 0 0
## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 1 0 52 0 0 0
## 2 17 0 0 0 86 0
## 3 3 14 0 0 0 0
## 4 23 0 0 0 1 85
## 5 12 0 0 80 0 0
## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 12 0 0 80 0 0
## 2 23 0 0 0 1 85
## 3 17 0 0 0 86 0
## 4 3 14 0 0 0 0
## 5 1 0 52 0 0 0
## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 1 0 52 0 0 0
## 2 23 0 0 0 1 85
## 3 17 0 0 0 86 0
## 4 3 14 0 0 0 0
## 5 12 0 0 80 0 0先ほどと異なり、結果が安定していることが分かると思います。なお、各サンプルが分類されるグループの「番号」は異なる解析の間で異なってきますが、これは5つのグループに任意につけられている番号なので特に問題はありません。
では、k-平均法で分類されたグループと、階層的クラスタ解析で分類されたグループ、および、品種・系統が所属する分集団の関係について、クロス集計表を作成して確認してみましょう。
# compare results
print("k-means vs. subpopulations")## [1] "k-means vs. subpopulations"table(kms$cluster, subpop)## subpop
## ADMIX AROMATIC AUS IND TEJ TRJ
## 1 1 0 52 0 0 0
## 2 23 0 0 0 1 85
## 3 17 0 0 0 86 0
## 4 3 14 0 0 0 0
## 5 12 0 0 80 0 0print("hierarchical clustering vs. subpopulations")## [1] "hierarchical clustering vs. subpopulations"table(cluster.id, subpop)## subpop
## cluster.id ADMIX AROMATIC AUS IND TEJ TRJ
## 1 5 0 0 0 84 0
## 2 14 0 0 80 0 0
## 3 1 0 52 0 0 0
## 4 2 14 0 0 0 0
## 5 34 0 0 0 3 85print("k-means vs. hierarchical clustering")## [1] "k-means vs. hierarchical clustering"table(kms$cluster, cluster.id)## cluster.id
## 1 2 3 4 5
## 1 0 0 53 0 0
## 2 0 1 0 0 108
## 3 88 1 0 0 14
## 4 1 0 0 16 0
## 5 0 92 0 0 0クロス集計表を作成してみると、ADMIXとIND以外の分集団に所属する品種・系統については、k-平均法と階層的クラスタ解析で同じように分類されていることが分かります。3つ目のクロス集計表を見ると、両手法の分類結果はほぼ一致していますが、一部違いが見られます。これは、分集団がADMIX（混合）となっている品種・系統の分類が両手法で異なることが主な原因です。
では、k-平均法と階層的クラスタ解析による分類の結果を、主成分軸上にプロットすることにより確認してみましょう。
# match id between the results of kmeans and hclust
convert.table <- apply(table(kms$cluster, cluster.id), 1, which.max)
convert.table## 1 2 3 4 5
## 3 5 1 4 2cluster.id.kms <- convert.table[kms$cluster]
# draw graph based on the converted id
#pdf("fig14.pdf", width = 8, height = 8)
op <- par(mfrow = c(2,2))
plot(pca$x[,1:2], pch = cluster.id, col = as.numeric(subpop), main = "hclust")
plot(pca$x[,3:4], pch = cluster.id, col = as.numeric(subpop), main = "hclust")
plot(pca$x[,1:2], pch = cluster.id.kms, col = as.numeric(subpop), main = "kmeans")
plot(pca$x[,3:4], pch = cluster.id.kms, col = as.numeric(subpop), main = "kmeans")
Fig.14 Classification by hierarchical cluster analysis (top) and k-means (bottom) Relationship between the scores and the principal component scores
par(op)
#dev.off()適切なグループ数の決定
ここまで分類するグループ数を分集団の数に合わせて5としてきました。では、この5というグループ数が本当に適切な数かどうかを確認するには、どのようにするとよいのでしょうか。適切なグループ数を決めるための1つの方法として、様々な数のグループに分類して、そのときの群（グループ）内平方和（Within groups sum of squares）の減少の程度を基準に決めるという方法があります。
分散分析の際に説明したのと同じように、全平方和は群間平方和と群内平方和に分割されます。したがって、グループ数が1のときは、全平方和が群内平方和となります。その後、2, 3, 4…とグループ数が増えていくと、群間平方和が大きくなり、群内平方和は小さくなって行きます。最終的にグループ数がサンプル数に一致すると、群内平方和は0となります。したがって、群内平方和を最小化するというルールではグループ数が常にサンプル数となってしまい意味がありません。そこで、主成分分析において主成分の数を決めたときのルールと同じように、群内平方和の減少が、急な変化からなだらかな変化に変わる点を見つけ、それを採用するグループ数とします。
では、実際にグループの数を1〜10に変化させて群内平方和を計算し、その減少の様子を図示してみましょう。
# visualize within-group and between-groups sum of squares
n <- nrow(data.mk)
wss <- rep(NA, 10)
wss[1] <- (n - 1) * sum(apply(data.mk, 2, var))
for(i in 2:10) {
print(i)
res <- kmeans(data.mk, centers = i, nstart = 50)
wss[i] <- sum(res$withinss)
}## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10plot(1:10, wss, type = "b", xlab = "Number of groups",
ylab = "Within groups sum of squares")
Fig.15 Change in the within group variation when the number of clusters is 1 to 10 by k-means method
Figure 15を見ると、クラスタ数が5を過ぎたあたりから群内分散の減少が直線的になるのが分かります。この図からも、5という数が適切なグループ数であると考えられます。
分類が曖昧なサンプルの検出
ここまでは、各サンプルを必ず1つのグループに分類してきました。すると、あるグループに分類されたサンプルの中にも、明確にそのグループに分類されたものと、「かろうじて」そのグループに分類されたものが存在することになります。分類の確かさを明示するためにも、後者のように分類が曖昧なサンプルを何らかの基準で検出できると便利です。ここでは、shadow value（Everitt and Hothorn 2011, An introduction to applied multivariate analysis with R. Springer）という統計量をもとに分類の曖昧さを評価する方法を紹介します。
kmsでは、k-平均法で求められた各グループの重心の位置が、計算されています。ここでは、各サンプルからこれらグループの重心までの距離を計算し、最も近いグループまでの距離と、次に近いグループまでの距離を計算し、その違いをもとに曖昧さを評価してみます。
まずは、パッケージfieldsに含まれている関数rdistを用いて、全サンプルとグループ重心間の距離を計算します。
d2ctr <- fields::rdist(kms$centers, data.mk)
d2ctr[, 1:10]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 49.59083 32.90422 20.82710 39.54322 21.82860 44.50590 44.34419 45.67021
## [2,] 36.10209 43.10088 41.49647 34.88985 41.45680 24.50574 22.59442 34.04487
## [3,] 19.66441 48.58701 47.34412 38.28548 47.10191 37.74005 36.65759 26.00505
## [4,] 40.35414 41.98393 39.47666 13.04416 39.19190 38.61061 37.13867 35.43122
## [5,] 48.19311 20.62444 28.70590 39.76012 29.21810 43.21613 42.58438 44.71007
## [,9] [,10]
## [1,] 51.53614 32.28608
## [2,] 37.71899 43.98252
## [3,] 20.32749 49.01558
## [4,] 42.25515 42.21440
## [5,] 49.99222 22.91782apply(d2ctr, 2, which.min)## [1] 3 5 1 4 1 2 2 3 3 5 4 1 3 4 5 1 1 2 2 2 2 2 2 5 3 5 5 2 1 2 4 5 1 4 2 1 1
## [38] 3 3 4 2 3 3 5 1 2 3 5 3 3 3 2 5 3 2 2 5 5 2 5 2 5 5 1 3 5 1 3 2 1 2 1 2 3
## [75] 2 4 3 2 2 3 2 3 1 5 2 2 5 5 4 3 2 3 2 5 3 5 2 3 2 5 4 5 5 3 5 5 1 5 3 3 2
## [112] 5 5 2 2 5 5 3 3 5 5 2 5 2 2 3 1 1 3 3 5 3 3 5 4 5 5 5 2 2 5 2 5 3 2 5 5 2
## [149] 2 3 1 3 3 3 2 2 3 3 2 2 5 2 4 3 5 2 5 2 1 2 2 5 3 5 5 5 5 2 2 2 2 2 2 2 3
## [186] 3 4 5 3 3 5 1 5 3 3 5 5 2 3 2 2 5 2 3 2 3 3 3 5 3 2 5 5 5 3 3 2 5 4 1 3 3
## [223] 3 3 3 3 5 3 3 5 2 2 3 1 3 3 3 2 3 3 3 5 2 2 3 3 3 3 3 5 5 3 3 5 5 3 3 3 3
## [260] 5 3 3 2 2 2 3 1 5 1 5 1 1 1 1 1 1 1 1 5 1 1 1 1 1 1 5 3 3 2 1 5 3 5 4 1 2
## [297] 3 4 1 1 2 5 5 2 2 1 3 5 1 1 1 3 2 3 3 2 3 1 1 1 4 2 3 2 1 2 3 2 2 5 2 3 3
## [334] 2 2 2 2 2 2 2 5 2 5 5 2 2 5 2 2 5 2 2 5 2 2 2 2 2 5 5 2 5 1 2 3 4 3 5 5 5
## [371] 3 2 1 3k-平均法では、先に述べたように、重心までの距離が最も近いグループに各サンプルを分類します。したがって、上のボックスで表示される2つの結果は、一致することに注意しましょう。
各サンプルについて計算された距離から、もっとも近い重心までの距離を取り出すには、関数minを使うことができます。しかし、2番目に近い重心までの距離を取り出すのにはどうすればよいのでしょうか。ここでは、自作関数nth.minを作成して、これを実現します。自作関数nth.minは、引数xとして与えられた配列を大きさ順（昇順）に並べ直し、そのn番目の値を返すという関数です。
# prepare own function find the second best
nth.min <- function(x, n) {
sort(x)[n]
}
nth.min(-10:10, 3)## [1] -8d.1st <- apply(d2ctr, 2, min)
d.2nd <- apply(d2ctr, 2, nth.min, n = 2)上のボックスにあるコマンドを実行すると、d.1stには、各サンプルから最も近い重心までの距離が、d.2ndには、2番目に近い重心までの距離が代入されます。
次に、shadow valueを計算してみましょう。i番目のサンプルのshadow valueは以下のように定義されます。 \[ s(\mathbf{x}_i)=\frac{2d(\mathbf{x}_i,c(\mathbf{x}_i))}{d(\mathbf{x}_i,c(\mathbf{x}_i))+d(\mathbf{x}_i,\tilde c(\mathbf{x}_i))} \] ここで、\(d(\mathbf{x}_i,c(\mathbf{x}_i))\)は、\(i\)番目のサンプルの観察値\(\mathbf{x}_i\)から最も近いグループの重心（そのサンプルが分類されたグループの重心）までの距離、\(d(\mathbf{x}_i,\tilde c(\mathbf{x}_i))\)は、2番目に近いグループの重心までの距離を表しています。この値は、0〜1までの値をとります。この値が0に近ければ、そのサンプルが分類されたグループの重心付近に位置していることを示しており、逆に、1に近ければ、分類されたグループの重心と2番目に近いグループの重心までの距離がほとんど同じであることを意味します。したがって、分類が曖昧なサンプルを検出するには、shadow valueが1に近いサンプルを見つければよいということになります。
先ほど計算しておいた距離d.1stおよびd.2ndを用いてshadow valueを計算し、その値が0.9以上になるものを検出してみましょう。
# evaluate unclearness of classificatinos (shadow values)
shadow <- 2 * d.1st / (d.1st + d.2nd)
unclear <- shadow > 0.9検出の結果は、unclearに代入されます。この値が、T（真）であれば分類が曖昧、F（偽）であれば分類が比較的明瞭であると考えられます。
では、分類が曖昧と判断されたサンプルの●で表して、主成分軸上の散布図を描いてみましょう。
# visualize the result
cluster.id.kms[unclear] <- 20
op <- par(mfrow = c(1,2))
plot(pca$x[,1:2], pch = cluster.id.kms, col = as.numeric(subpop), main = "kmeans")
plot(pca$x[,3:4], pch = cluster.id.kms, col = as.numeric(subpop), main = "kmeans")
Fig.16 Scatter plot of the samples with the marker for unclear classification (●)
par(op)Figure 16をみると分集団がADMIX（混合）となっている品種・系統（黒色の点）のほとんどが●で散布されていることが分かります。このように、各サンプルの分類の曖昧さ（逆に言うと、確からしさ）を評価することで、より詳細に分類結果を把握することがでるようになります。この例では、複数の分集団に起因するゲノムが混合されていると考えられる品種・系統を見つけ出すことができました。
代表サンプルの選択
クラスタ解析は、多数のサンプルから、少数の代表的なサンプルを選び出すためにも利用できます。例えば、多数の遺伝資源について収集された既存のデータをもとにクラスタ解析で分類を行い、その分類結果に基づいて代表的な品種・系統を選び出すことができます。こうして代表的な品種・系統を選び出しておいて、それら品種・系統を用いて時間やコストを要する圃場試験や分子生物学的実験を行うことがよく行われます。
ここでは、このような代表サンプルの選択に適したクラスタ解析の方法として、k-medoids法を紹介します。k-medoids法は、k-平均法に似ていますが、グループの重心までの距離をもとにグループ分けするのではなく、グループの代表サンプル（medoids）までの距離をもとにグループ分けします。より具体的には、クラスタの中心を重心とするのではなく、グループの代表サンプルの座標点とするアルゴリズムです。
では、表現型データ（data.tr）に含まれる229品種・系統の中から、k-medoids法で代表となる48サンプルを選出してみましょう。k-medoids法を実行する関数pamはパッケージclusterに含まれています。k-medoids法で得られる結果のうち、medoidsのid（id.med）が代表として選ばれたサンプルのIDです。
# select 48 varieties/lines (by classifying all samples into 48 groups)
n.sel <- 48
kmed <- cluster::pam(data.tr, k = n.sel)
head(kmed)## $medoids
## Flowering.time.at.Arkansas Flowering.time.at.Faridpur
## 1 -0.84468910 -0.9295405
## 61 0.40615449 0.3485777
## 16 0.49444934 -0.8133479
## 316 -0.10889875 0.5809628
## 214 1.23023969 0.9295405
## 154 -0.68281523 -0.9295405
## 10 0.38408078 -1.9752735
## 108 1.54662954 0.2323851
## 24 0.55331256 0.4647702
## 377 0.39143869 0.6971553
## 73 0.59745999 0.6971553
## 40 0.97271306 -3.0210065
## 123 0.64896531 0.9295405
## 314 -0.30756215 0.5809628
## 45 -0.25605682 -0.8133479
## 157 0.29578594 0.9295405
## 220 -1.92630092 -0.8133479
## 60 -0.20455150 -0.4647702
## 170 0.11919626 0.5809628
## 177 -0.76375217 -0.8133479
## 81 -0.12361456 0.5809628
## 84 0.39879659 4.4153172
## 342 0.49444934 0.8133479
## 188 1.48776631 0.0000000
## 107 -0.28548844 -3.0210065
## 112 1.09043952 0.6971553
## 250 -2.37513304 -2.2076586
## 236 -0.22662521 0.4647702
## 180 -0.46207812 -0.5809628
## 122 1.98319848 0.4647702
## 356 0.20749110 -0.4647702
## 271 -0.52829925 -1.0457330
## 163 1.59813486 2.7886214
## 247 -0.85940491 -0.8133479
## 162 2.32656731 0.6971553
## 168 1.45097679 0.6971553
## 169 0.09712255 0.5809628
## 348 -1.02863669 -1.7428884
## 217 -0.45472022 -0.2323851
## 343 -0.79318378 -0.2323851
## 218 0.10448045 0.6971553
## 268 -2.14115170 -0.8133479
## 291 -1.55104784 -1.2781181
## 237 -0.80054168 -0.6971553
## 239 -0.54301506 0.4647702
## 272 2.06904069 -1.1619256
## 277 -0.20455150 -0.5809628
## 369 -0.78582588 0.4647702
## Flowering.time.at.Aberdeen Culm.habit Flag.leaf.length Flag.leaf.width
## 1 -0.65883657 -0.12403316 -0.430259255 0.15313754
## 61 -0.16712476 2.05986166 0.390547174 -1.13475849
## 16 0.61141895 0.96791425 -0.147909619 -1.07036369
## 316 -0.86371649 1.78687481 1.048067066 0.99026996
## 214 1.14410674 0.14895370 1.323613096 0.86148035
## 154 -0.20810074 -0.67000686 -1.032378359 -1.00596889
## 10 -0.94566846 -0.67000686 -0.513360439 -0.94157409
## 108 0.89825084 -1.76195427 1.778482733 1.31224396
## 24 0.69337091 -0.67000686 0.887696360 1.11905956
## 377 0.48849099 -0.67000686 0.167486102 1.18345436
## 73 0.36556304 -0.39702001 0.791473937 0.76488815
## 40 1.47191462 -0.94299371 0.394920921 1.37663876
## 123 -0.37200468 0.42194055 -0.602779256 -0.16883647
## 314 -0.00322082 1.24090111 1.603532873 1.08686216
## 45 -0.12614877 -0.67000686 -0.877839314 -1.36014030
## 157 1.88167446 -0.30602439 0.596599231 -0.46934554
## 220 -1.56030823 -0.39702001 0.005657481 -0.84498189
## 60 -0.37200468 0.14895370 0.636934894 0.79708555
## 170 -0.57688460 -1.76195427 -0.376316382 0.08874274
## 177 -0.86371649 0.05795808 -1.145123824 -0.29762607
## 81 -0.20810074 1.14990549 0.629645316 0.26046221
## 84 0.16068312 -0.12403316 0.475106273 1.63421797
## 342 0.40653902 -0.94299371 0.944555065 0.63609855
## 188 2.08655438 -1.76195427 1.708502789 2.23523612
## 107 -0.94566846 -0.12403316 -0.383605958 0.08874274
## 112 3.31583391 -0.88839634 0.496100256 -0.74838968
## 250 -1.56030823 0.90725051 -0.385063875 -1.07036369
## 236 -0.20810074 -0.39702001 -0.627077848 -0.32982348
## 180 1.06215477 0.14895370 0.179149425 -0.94157409
## 122 2.12753037 -0.12403316 1.510226281 0.95807256
## 356 -0.69981255 1.78687481 0.112085312 0.15313754
## 271 0.11970713 0.42194055 -0.420053847 -0.55520528
## 163 2.08655438 2.05986166 -0.334522804 -1.07036369
## 247 -0.69981255 0.42194055 -1.561601687 -2.03628571
## 162 0.77532288 0.05795808 3.255837110 -2.27239998
## 168 0.65239493 0.05795808 0.722951909 -1.13475849
## 169 0.20165910 0.96791425 0.776269961 -1.10716072
## 348 -1.19152437 1.78687481 -0.902623879 -0.52300788
## 217 0.11970713 -0.67000686 0.096048241 -1.26354809
## 343 -0.69981255 0.14895370 -0.067238295 -0.16883647
## 218 1.79972249 0.69492740 -0.946361344 -1.45673250
## 268 -1.56030823 0.42194055 1.186569039 -0.52300788
## 291 -1.19152437 -1.48896742 -2.195794931 -1.07036369
## 237 -0.86371649 -0.57901124 -1.405604727 -1.10256109
## 239 -0.49493263 -0.67000686 -0.793280216 0.81855049
## 272 0.40653902 -0.67000686 0.877490952 -0.78058708
## 277 1.67679454 -0.67000686 -0.966286189 1.03319982
## 369 -0.74078854 1.78687481 -0.934698019 -1.23135069
## Panicle.number.per.plant Plant.height Panicle.length
## 1 -0.36635080 -0.300031055 -1.06663830
## 61 1.84250994 1.167426588 0.49972675
## 16 0.47971752 2.071970582 1.71585257
## 316 0.27647886 0.845207051 0.68894365
## 214 -0.80385475 0.270646917 0.65065472
## 154 -0.48677421 -0.187448329 -1.33041417
## 10 0.99065981 -1.600555688 -0.60215885
## 108 -1.74085297 -0.213329414 1.33219756
## 24 -0.75874020 0.720977831 -0.37242531
## 377 -0.65427377 0.107596068 -0.17485446
## 73 -0.56409951 -0.672718711 -0.25985587
## 40 -0.59377139 -0.121451555 0.93329660
## 123 1.31487928 0.687332419 -0.20701716
## 314 1.10382219 0.620041594 0.00433770
## 45 0.31048254 1.109194142 1.24030415
## 157 0.09612401 -0.891413899 -0.18732571
## 220 0.74076723 -1.029877716 -1.23928653
## 60 -0.56409951 0.538516170 0.76538226
## 170 -0.69549982 -1.227868033 -0.20701716
## 177 0.86108186 -1.274453990 -1.09838329
## 81 0.61981687 1.023786551 0.50362526
## 84 -2.11078659 0.709331342 1.82382734
## 342 -0.68512470 0.033834966 0.21262944
## 188 -0.80187497 0.127006880 -0.33107327
## 107 -0.49630064 -0.144744534 -0.31728926
## 112 0.49536037 1.587254793 0.90183006
## 250 0.22802399 -1.057052857 -1.34649552
## 236 0.03107200 -0.004986663 -0.39386710
## 180 -0.26910207 -0.214623469 0.34434334
## 122 0.30372141 1.994327324 1.41030697
## 356 1.62645029 1.594464524 1.25158197
## 271 -1.39346711 -0.144744534 -0.08985305
## 163 2.03700306 1.237305524 0.75786371
## 247 0.91378118 -0.936705802 -0.37931731
## 162 1.25264475 0.547574548 0.95083988
## 168 0.81804671 1.960681907 2.16383298
## 169 0.00812393 1.234532551 1.37617512
## 348 1.37549413 0.305586385 -0.13120509
## 217 -0.84582365 0.158064186 1.06361817
## 343 -0.58383871 0.057127945 -0.53323878
## 218 0.04624607 -0.168037512 0.15136716
## 268 0.38031422 -0.529078678 -0.54013079
## 291 0.37834782 -2.310991531 -2.16128381
## 237 -0.93628981 -1.041524205 -0.28053189
## 239 -1.09123051 -0.487668938 0.09278511
## 272 0.43535026 1.168720639 2.16383298
## 277 -0.34836453 -0.757349867 -1.79677326
## 369 1.97836248 0.099831740 -1.33479999
## Primary.panicle.branch.number Seed.number.per.panicle Florets.per.panicle
## 1 -0.431283709 -0.258760638 -0.475933049
## 61 -1.102193707 -0.876079700 -0.895386728
## 16 -0.672535592 -0.411988727 -0.555640631
## 316 0.009862591 1.779281732 1.346165250
## 214 1.694796377 1.218266993 1.169580226
## 154 0.313150672 -0.880995499 -0.881536337
## 10 0.212054643 -2.421533974 -1.070689887
## 108 1.526302998 0.783516123 0.721383139
## 24 -0.697809599 -0.031257998 0.011025957
## 377 0.818630806 -0.007760557 -0.108825920
## 73 1.020822864 -0.233942620 -0.161231929
## 40 -0.707000147 0.248248555 0.436585850
## 123 -0.065959429 0.225343099 0.041056341
## 314 -0.866302977 0.455621796 0.087742049
## 45 -1.263947350 -1.019223437 -1.256637457
## 157 0.529785018 0.272822174 0.319137897
## 220 -1.270687085 -2.844547836 -2.898639507
## 60 0.212054643 0.780393490 1.235134586
## 170 0.313150672 0.162738678 0.396877904
## 177 0.313150672 0.088725396 -0.125655596
## 81 -1.430755795 0.189970641 -0.152136445
## 84 -0.040685422 1.167279818 1.791276653
## 342 0.683836104 1.054928315 1.526813885
## 188 1.121918887 2.007577296 2.009790344
## 107 0.175292452 -1.105116691 -0.349351235
## 112 0.414246702 0.086543832 0.007252743
## 250 -1.809865896 -2.267730366 -2.386656155
## 236 0.212054643 0.111542573 0.418318840
## 180 -0.394521517 -0.148856628 -0.153953524
## 122 1.829591079 -0.566814261 -0.410958930
## 356 0.395865606 0.410082280 0.212137580
## 271 0.414246702 0.012775580 0.167801038
## 163 0.481644052 0.637547278 0.566063087
## 247 -0.707000147 -1.226357451 -1.471648074
## 162 0.212054643 -0.783130250 -0.094353020
## 168 0.009862591 -0.089133501 0.393803434
## 169 -0.770021047 -3.247080853 -1.182466197
## 348 -1.152741721 -0.779791898 -1.230789125
## 217 0.892155188 1.067628007 1.072818192
## 343 0.892155188 0.264804800 0.197920016
## 218 -0.091233436 -0.779791898 -0.464469237
## 268 -1.253837747 -1.593159482 -1.714015333
## 291 -1.068495031 -1.116745113 -1.320347713
## 237 -0.596713571 -1.047599408 -1.225744631
## 239 0.085684611 -0.423852868 -0.563045045
## 272 -1.754722609 -2.613366309 -1.339268251
## 277 0.481644052 0.178328133 0.424419155
## 369 -2.085582334 -2.386711605 -2.782174502
## Panicle.fertility Seed.length Seed.width Brown.rice.seed.length
## 1 0.567830451 -0.47807783 1.32968778 -0.5426835
## 61 -0.064962376 -0.63213686 -0.31430244 -0.7177220
## 16 0.337723969 1.56524223 -1.82328313 1.7235440
## 316 1.419224437 -1.60011631 -0.33193066 -1.4624181
## 214 0.165144107 0.95478759 -1.41825456 1.2542065
## 154 -0.122488997 0.06374449 0.82027822 0.1481969
## 10 -3.366990401 -0.69856911 0.18885942 -1.4930995
## 108 0.176649431 0.56812254 0.65413531 0.5171695
## 24 -0.180015617 1.25236750 0.27895693 1.2994996
## 377 0.234176051 0.40412017 -0.01216130 0.5191880
## 73 -0.283563534 1.49261108 -0.01455942 1.5177690
## 40 -0.571196638 -0.30799301 1.51030352 -0.4040984
## 123 0.510303831 -1.07255506 -0.15493314 -1.1148193
## 314 1.097075361 -0.12506247 -1.01631830 -0.2206748
## 45 0.579335775 1.00875573 -0.63009598 1.1237658
## 157 -0.180015617 -1.19315610 0.51616671 -1.2461234
## 220 -0.099478348 -0.53850824 0.86630897 -0.5623752
## 60 -1.238505437 -0.71948595 0.38003526 -0.7133576
## 170 -0.697755203 1.02927842 -1.56083707 1.2556419
## 177 0.590841099 -0.58972091 0.20118656 -0.5042084
## 81 0.993527444 1.06064920 -1.27609277 1.1031433
## 84 -1.641191782 -0.70062944 0.52682498 -1.0287359
## 342 -1.273021410 -0.42880416 -0.36508161 -0.4540708
## 188 0.073101514 -1.18203031 -0.06725478 -0.8587405
## 107 -2.066888775 1.54168280 -0.19751015 1.6766584
## 112 0.176649431 0.26308602 -2.30496716 0.5877053
## 250 0.130628134 -1.05940477 1.26228820 -1.1978586
## 236 -0.904851037 -1.27413608 0.04724457 -1.5709916
## 180 -0.053457052 -1.41945214 1.09883791 -1.3020137
## 122 -0.536680665 -0.92315417 0.08942190 -0.5160727
## 356 0.556325127 -0.60489570 -0.47240957 -0.6878481
## 271 -0.502164693 1.78580518 1.41193884 1.8496734
## 163 0.188154755 -0.80383412 -1.38411302 -0.6885433
## 247 0.590841099 -0.59548088 1.12248246 -0.7662560
## 162 -1.905814237 -0.44797012 -1.62665895 -0.5749460
## 168 -1.365064003 0.29466283 -0.16585086 0.2592456
## 169 -4.609565407 2.21318954 0.83130112 1.9993567
## 348 1.281160547 -0.51645374 0.03504365 -0.5850385
## 217 -0.007435755 0.92450609 -1.50355584 1.1270066
## 343 0.153638782 0.59238070 0.47476070 0.5686752
## 218 -0.973882982 1.03291535 -1.36982952 1.1819326
## 268 0.188154755 0.58643261 0.60943365 0.3197815
## 291 0.464282534 -1.59959675 0.15842022 -1.7986124
## 237 0.395250589 2.92968316 -0.85724350 2.0024798
## 239 0.326218644 0.52911061 1.00867730 0.4042002
## 272 -3.228926512 4.34281550 -0.54896749 3.4797497
## 277 -0.732271176 0.03777536 1.74077474 -0.1811793
## 369 0.970516796 0.58274193 -1.22539584 0.5756054
## Brown.rice.seed.width Straighthead.suseptability Blast.resistance
## 1 1.404481793 -1.49925662 0.99698930
## 61 -0.298354717 0.43103397 -0.41451332
## 16 -1.811381664 0.79560418 -0.06163767
## 316 -0.081168774 -0.29327769 -0.06163767
## 214 -1.317666011 1.15655282 -0.41451332
## 154 0.863149235 -0.77373777 0.99698930
## 10 0.860887147 -2.58813849 -1.12026464
## 108 0.657787831 1.27847862 -1.47314029
## 24 0.344103361 0.49260046 -1.47314029
## 377 -0.021274910 0.79318980 -0.06163767
## 73 -0.111613458 0.43103397 -1.47314029
## 40 1.419865381 -0.53230054 -0.41451332
## 123 -0.001610187 1.27847862 -0.76738898
## 314 -0.947033533 0.43103397 -0.41451332
## 45 -0.589829617 0.55295977 -0.06163767
## 157 0.810093041 -1.13830797 -1.12026464
## 220 1.039848983 -0.29086332 0.29123799
## 60 0.336700171 -1.01638217 0.99698930
## 170 -1.660104674 -0.53471492 -0.76738898
## 177 0.130901771 -1.50046380 0.99698930
## 81 -1.332873397 0.31031536 -1.47314029
## 84 0.564684771 -0.17014471 -0.41451332
## 342 -0.570267715 0.55416695 -0.06163767
## 188 0.045788456 1.03704140 -1.47314029
## 107 -0.286427355 0.07008532 -1.82601595
## 112 -2.463864896 1.03704140 -1.47314029
## 250 1.585834243 -0.65543353 1.34986495
## 236 -0.022534480 -0.05063329 0.99698930
## 180 1.325540141 -0.29327769 0.99698930
## 122 0.128151283 0.43344834 -0.76738898
## 356 -0.484128515 1.27606425 -0.41451332
## 271 1.346394329 -1.37974519 1.34986495
## 163 -1.342178797 0.49380765 -1.47314029
## 247 1.110462050 -0.77373777 1.34986495
## 162 -1.680000747 1.39919723 0.99698930
## 168 -0.051434497 1.39919723 1.34986495
## 169 0.242425174 -0.41278912 0.29123799
## 348 -0.055206197 -0.41158193 -0.76738898
## 217 -1.620055474 -0.89566356 1.34986495
## 343 0.250123976 -0.77615214 0.99698930
## 218 -1.566433759 -1.73948665 1.34986495
## 268 0.555942557 -0.29327769 0.99698930
## 291 0.308372688 -0.65301915 0.99698930
## 237 -1.257890394 -0.65422634 0.64411364
## 239 0.859653282 -1.49925662 -1.47314029
## 272 -0.631126702 0.43344834 0.29123799
## 277 1.710685957 -0.89687075 0.99698930
## 369 -1.172805119 0.73283050 0.29123799
## Amylose.content Alkali.spreading.value Protein.content
## 1 -0.68226975 0.24251496 -0.14673626
## 61 1.13994856 -0.06147408 -1.19681764
## 16 0.30412867 -0.08680650 1.25337224
## 316 0.75274665 -0.21346860 -0.43842554
## 214 0.69195409 0.21718254 0.02827730
## 154 -0.19561743 -0.74544942 1.72007508
## 10 -0.77891435 1.07848483 1.07835868
## 108 0.63957895 0.16651770 0.90334512
## 24 0.40202029 -0.28946586 -0.20507412
## 377 0.37707975 0.09052044 0.26162872
## 73 -0.94289841 -3.40535353 0.14495301
## 40 -3.56352580 1.07848483 0.14495301
## 123 1.03426302 0.39450948 -1.19681764
## 314 1.05982707 -0.03614166 -0.96346623
## 45 -0.01043391 -0.64411974 0.55331799
## 157 -0.06623837 1.07848483 -0.67177695
## 220 -0.08120270 1.07848483 2.47846719
## 60 -0.18439419 -0.55545627 0.37830443
## 170 0.50178246 -0.74544942 0.37830443
## 177 -1.09815328 0.73649716 -1.13847979
## 81 1.15397762 -0.85944531 -0.20507412
## 84 0.21028988 -0.51745764 0.84500726
## 342 -0.53823812 1.07848483 -0.55510124
## 188 0.74682328 0.62250127 -0.49676339
## 107 -0.07122648 -0.40346175 0.49498014
## 112 0.65173746 1.07848483 -1.08014193
## 250 -0.55320245 0.96448894 -0.96346623
## 236 -3.53671472 0.77449579 -0.26341197
## 180 -0.05937973 0.96448894 -0.67177695
## 122 0.53607570 0.62250127 1.72007508
## 356 0.79358679 -0.36546312 -0.61343910
## 271 -0.86121813 0.47050675 1.83675079
## 163 1.39465384 0.92649031 -1.37183121
## 247 -0.71032786 0.69849853 1.42838581
## 162 1.05172139 1.07848483 -1.13847979
## 168 1.62660088 0.57183643 -1.02180408
## 169 -3.60280715 0.81249442 1.31171010
## 348 0.88274923 -0.74544942 -0.78845266
## 217 -1.04390760 0.50850538 0.02827730
## 343 -0.91234624 -0.41612796 -0.03006055
## 218 -1.38902234 -2.37939052 0.49498014
## 268 -0.30940865 -0.27426641 0.08661516
## 291 -0.81850746 0.88849168 -1.72185833
## 237 -1.46914383 -3.48135079 1.19503439
## 239 -1.07009517 -2.37939052 0.37830443
## 272 1.17642410 1.07848483 1.83675079
## 277 -0.42101758 1.00248757 0.55331799
## 369 1.04610977 -0.74544942 0.49498014
##
## $id.med
## [1] 1 28 8 192 121 80 7 53 10 218 33 15 63 191 18 83 126 27 93
## [20] 98 36 38 202 106 52 54 148 136 101 62 211 161 86 145 85 91 92 207
## [39] 123 203 124 159 178 137 138 162 166 214
##
## $clustering
## 1 3 5 6 8 9 10 16 22 24 25 26 27 39 40 43 44 45 46 51
## 1 2 3 4 5 6 7 3 8 9 10 9 11 2 12 13 14 15 9 16
## 53 54 55 57 58 59 60 61 65 69 70 72 73 74 79 81 83 84 87 88
## 3 5 17 13 14 18 18 2 9 19 10 10 11 5 20 21 20 22 23 4
## 89 92 93 94 96 99 100 101 103 105 106 107 108 112 113 114 116 117 118 120
## 24 10 3 16 16 9 23 19 16 14 11 25 8 26 27 23 10 28 1 10
## 121 122 123 125 129 130 131 132 134 135 137 139 140 142 147 149 150 152 153 154
## 29 30 13 31 13 13 14 5 32 10 2 5 18 33 24 5 23 3 21 6
## 155 156 157 160 162 163 164 165 166 167 168 169 170 171 172 174 176 177 178 179
## 34 13 16 15 35 33 10 18 36 19 36 37 19 38 6 5 8 20 4 34
## 180 183 184 186 187 188 189 191 195 198 201 202 203 205 206 207 208 209 211 213
## 29 10 29 16 39 24 13 15 11 19 8 9 2 32 36 5 5 13 34 40
## 214 215 217 218 219 220 221 222 225 228 231 232 233 234 235 236 237 239 240 241
## 5 10 39 41 27 17 42 31 1 38 38 43 43 38 13 28 44 45 45 13
## 242 243 244 245 247 248 249 250 251 252 254 256 257 263 264 265 266 267 268 270
## 40 42 5 40 34 1 1 27 19 21 31 27 6 1 1 1 1 34 42 40
## 271 272 273 275 276 277 279 280 281 282 283 284 285 286 287 289 290 291 294 296
## 32 46 45 1 4 47 17 23 40 47 6 21 19 11 27 17 20 43 13 1
## 297 299 300 301 303 306 307 308 311 313 314 316 321 333 334 335 337 338 339 340
## 47 38 1 6 47 17 1 5 29 21 14 4 4 27 1 6 2 29 2 36
## 341 342 343 344 346 347 348 349 350 355 356 359 365 369 371 372 376 377 381 384
## 14 23 40 2 21 10 38 13 45 1 31 14 20 48 14 14 32 10 18 10
## 385 386 387 391 392 394 395 396 397
## 31 23 20 19 11 11 10 19 19
##
## $objective
## build swap
## 2.350559 2.334805
##
## $isolation
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## no no no no no no no no no no no no no no no no no no no no no no no no no no
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
## no no no no no no L* no no no no no L* no no no L* no no no no no
## Levels: no L L*
##
## $clusinfo
## size max_diss av_diss diameter separation
## [1,] 15 3.396047 2.491193 5.013251 2.102297
## [2,] 8 3.940002 2.588179 5.590730 1.908541
## [3,] 5 2.935819 2.176676 3.671067 2.647235
## [4,] 6 3.287674 2.314052 3.984330 2.356047
## [5,] 12 4.270387 3.210883 6.538077 3.023159
## [6,] 7 4.211719 3.210739 5.811008 3.064565
## [7,] 1 0.000000 0.000000 0.000000 4.764459
## [8,] 4 4.231524 2.445166 5.136474 2.603085
## [9,] 6 2.583538 1.783762 3.684491 2.066555
## [10,] 14 4.132061 2.621656 5.608358 2.066555
## [11,] 7 3.542087 2.580328 4.814989 2.959904
## [12,] 1 0.000000 0.000000 0.000000 4.291645
## [13,] 12 4.001610 2.466821 5.330147 1.908541
## [14,] 9 3.857231 2.479860 4.587549 2.118303
## [15,] 3 3.148820 1.647018 3.364221 2.647235
## [16,] 6 4.434976 2.582121 5.332131 2.360260
## [17,] 5 3.681844 2.446892 4.489978 2.841495
## [18,] 5 4.001869 2.757616 5.355336 3.192992
## [19,] 10 3.792886 2.891107 5.595767 2.735103
## [20,] 6 3.213592 2.228942 4.182100 2.721912
## [21,] 6 3.475275 2.292260 5.173628 2.582119
## [22,] 1 0.000000 0.000000 0.000000 4.101246
## [23,] 7 4.195454 2.463641 4.961384 2.914611
## [24,] 3 4.203317 2.404226 4.770831 3.265332
## [25,] 1 0.000000 0.000000 0.000000 4.071087
## [26,] 1 0.000000 0.000000 0.000000 4.556660
## [27,] 6 3.382978 2.332328 4.390924 2.841495
## [28,] 2 4.170205 2.085103 4.170205 3.687605
## [29,] 5 3.221027 2.128065 4.903242 2.360260
## [30,] 1 0.000000 0.000000 0.000000 4.729190
## [31,] 5 3.013936 2.261990 4.262801 2.497997
## [32,] 4 3.666009 2.484052 5.561418 2.759170
## [33,] 2 3.884203 1.942101 3.884203 4.551940
## [34,] 5 4.277568 2.516751 5.655798 2.773746
## [35,] 1 0.000000 0.000000 0.000000 4.435230
## [36,] 4 4.059718 2.101660 4.740166 3.340238
## [37,] 1 0.000000 0.000000 0.000000 6.868458
## [38,] 6 4.339646 2.728750 5.913697 2.854855
## [39,] 2 2.942470 1.471235 2.942470 2.953791
## [40,] 6 3.803951 2.792516 5.780604 2.527868
## [41,] 1 0.000000 0.000000 0.000000 4.425165
## [42,] 3 3.422055 2.219950 4.138119 3.109955
## [43,] 3 2.431220 1.589834 2.456639 2.747514
## [44,] 1 0.000000 0.000000 0.000000 4.585757
## [45,] 4 3.834916 2.177878 4.956017 3.026757
## [46,] 1 0.000000 0.000000 0.000000 6.318230
## [47,] 4 2.911932 1.391317 3.396124 2.102297
## [48,] 1 0.000000 0.000000 0.000000 4.191402kmed$id.med## [1] 1 28 8 192 121 80 7 53 10 218 33 15 63 191 18 83 126 27 93
## [20] 98 36 38 202 106 52 54 148 136 101 62 211 161 86 145 85 91 92 207
## [39] 123 203 124 159 178 137 138 162 166 214なお、ここで用いられた表現型データ（data.tr）は既に分散1に基準化されていたことに注意しましょう。もし、自分の手持ちのデータで同様の解析を行う場合には、データを基準化する必要があるかどうかをよく検討して、必要であれば、関数scaleを用いて基準化しておきます。
では、代表として選ばれたサンプルのもつ変異を、主成分軸上の散布図として図示してみましょう。
# look at the distribution of 48 varieties/lines selected as medoids
pca.tr <- prcomp(data.tr, scale = T)
kmed$id.med## [1] 1 28 8 192 121 80 7 53 10 218 33 15 63 191 18 83 126 27 93
## [20] 98 36 38 202 106 52 54 148 136 101 62 211 161 86 145 85 91 92 207
## [39] 123 203 124 159 178 137 138 162 166 214pch <- rep(1, nrow(data.tr))
pch[kmed$id.med] <- 19
op <- par(mfrow = c(1,2))
plot(pca.tr$x[,1:2], col = as.numeric(subpop.tr), pch = pch)
plot(pca.tr$x[,3:4], col = as.numeric(subpop.tr), pch = pch)
Fig.17 Distribution of representative 48 varieties / lines selected by the k-medoids method
par(op)最後に、k-medoids法を用いて選出された48品種・系統のもつ主成分得点の分布と、全品種・系統のもつ主成分得点の分布を、ヒストグラムを描いて比較してみましょう。
# compare histogram between three datasets (all, 48 selected k-medoids, 48 selected randomly)
op <- par(mfcol = c(2,4))
for(i in 1:4) {
res <- hist(pca.tr$x[,i], main = paste("PC", i, "all"))
hist(pca.tr$x[kmed$id.med, i], breaks = res$breaks, main = paste("PC", i, "k-medoids"))
}
Fig.18 Distribution of principal component scores of all varieties / lines (top) and varieties / lines selected as representatives (bottom)
par(op)