{vegan}を使った菌叢解析の階層的クラスター分析とクラスタリング
0. はじめに
菌叢解析をしているとUnifrac distanceなどを計算して、各被験者の菌叢の全体像を被験者間で比較する必要 ( Beta diversity ) があると思います。またその全体像に違いがあった時に同じようなグループをまとめて1つのグループとして扱いたいときがあると思います。そんな時に必要になる手順をまとめてみました。具体的には階層的クラスター分析と{cutree}によるクラスタリングについてまとめてみました。参考になればと思います。
基本的にはRのパッケージ{vegan}のtutorialを参考にして書いています。このtutorialはRで菌叢解析をするなら読んでおいて損はないと思います。
1. データの作成
{vegan}に格納されている菌叢のデータセットであるduneを使っていきましょう。データの詳細は以下のように書かれています。
The dune meadow vegetation data, dune, has cover class values of 30 species on 20 sites. The corresponding environmental data frame dune.env has following entries:
library(tidyverse) library(vegan) data(dune)
2. 階層的クラスター分析
フリーソフトによるデータ解析・マイニング によると
階層的クラスター分析とは、個体間の類似度あるいは非類似度 (距離) に基づいて、最も似ている個体から順次に集めてクラスターを作っていく方法である。 階層的クラスター分析には幾つかの方法があるが、いずれも次のようなステップを踏む。
- 距離 (あるいは類似度) を求める方法を選択し、個体間の距離 (類似度) を計算する
- クラスター分析の方法 (最近隣法、最遠隣法など) を選択する
- 選択された方法のコーフェン（Cophenetic）行列を求める
- コーフェン行列に基づいて樹形図を作成する
- 結果について検討を行う
ということです。これに従ってやっていきましょう。
A. 距離 (あるいは類似度) を求める方法を選択し、個体間の距離 (類似度) を計算する
距離行列を作成するのにはvegan::vegdist()を使います。普通の距離行列ならばdist()でいいのですが、菌叢解析によく使う距離が計算できます。デフォルトはBray Curtis Dissimilarityになっております。それぞれの距離の特徴に関してはここでは割愛しますが、RDocumentationの記載を引用すると
Dissimilarity index, partial match to "manhattan", "euclidean", "canberra", "bray", "kulczynski", "jaccard", "gower", "altGower", "morisita", "horn", "mountford", "raup" , "binomial", "chao", "cao" or "mahalanobis".
などの距離が使用できるようです。
dune %>% vegdist(., method="bray") -> vegdist
とするとこんな感じの距離行列が得られます。
> vegdist 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 2 0.4666667 3 0.4482759 0.3414634 4 0.5238095 0.3563218 0.2705882 5 0.6393443 0.4117647 0.4698795 0.5000000 6 0.6363636 0.5111111 0.5681818 0.6344086 0.2967033 7 0.5517241 0.4390244 0.4750000 0.5058824 0.2289157 0.2272727 8 0.6551724 0.5365854 0.3250000 0.4117647 0.6385542 0.5909091 0.5250000 9 0.6000000 0.4761905 0.3414634 0.3793103 0.5058824 0.6000000 0.4878049 0.3170732 10 0.5737705 0.2941176 0.4698795 0.4772727 0.3488372 0.3186813 0.2771084 0.5421687 0.6000000 11 0.5600000 0.5405405 0.5555556 0.5844156 0.6266667 0.4500000 0.4444444 0.5277778 0.5945946 0.4133333 12 0.9245283 0.7142857 0.4400000 0.5250000 0.6923077 0.6385542 0.6266667 0.4400000 0.3506494 0.7179487 0.6716418 13 0.8431373 0.6000000 0.4246575 0.5128205 0.6842105 0.7530864 0.6438356 0.3698630 0.4133333 0.7368421 0.7538462 0.3529412 14 1.0000000 0.7878788 0.7500000 0.7971014 0.8805970 0.8055556 0.8750000 0.5625000 0.7575758 0.7611940 0.8214286 0.6949153 0.6491228 15 1.0000000 0.9076923 0.7142857 0.7352941 0.8484848 0.8028169 0.8412698 0.4285714 0.6615385 0.8484848 0.7454545 0.6206897 0.6785714 0.3617021 16 0.9215686 0.8933333 0.6712329 0.6666667 0.8947368 0.8518519 0.8904110 0.4246575 0.6533333 0.8947368 0.8769231 0.5882353 0.6060606 0.5438596 0.3571429 17 0.8787879 0.8245614 0.8909091 0.9000000 0.6206897 0.6825397 0.6727273 0.8909091 0.8947368 0.6206897 0.7021277 0.9200000 0.8750000 0.8974359 0.8947368 1.0000000 18 0.7777778 0.5942029 0.6119403 0.6666667 0.5428571 0.4933333 0.5522388 0.6417910 0.6811594 0.4857143 0.3220339 0.7419355 0.8000000 0.8431373 0.7200000 0.8666667 0.7619048 19 1.0000000 0.8082192 0.8309859 0.7894737 0.7027027 0.7215190 0.7464789 0.7464789 0.7808219 0.7027027 0.5555556 0.6969697 0.8125000 0.8545455 0.7777778 0.9062500 0.5652174 20 1.0000000 0.9452055 0.7746479 0.7631579 0.8918919 0.8481013 0.8873239 0.4929577 0.6986301 0.8918919 0.8095238 0.6969697 0.7187500 0.4545455 0.2962963 0.3437500 0.9130435 18 19 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 0.5517241 20 0.6896552 0.7419355
B. クラスター分析の方法 (最近隣法、最遠隣法など) を選択する
C. 選択された方法のコーフェン（Cophenetic）行列を求める
この2つの作業を行うのにはhclust()を使用します。方法に関しては{vegan}のtutorialには"single"、"complete"、"average"について言及がなされています。
Function hclust provides several alternative clustering strategies. In community ecology, most popular are single linkage a.k.a. nearest neighbour, complete linkage a.k.a. furthest neighbour, and various brands of average linkage methods. These are best illustrated with examples:
single
Some people prefer single linkage, because it is conceptually related to minimum spanning tree which nicely can be represented in ordinations, and it is able to find discontinuities in the data. However, single linkage is prone to chain data so that single sites are joined to large clusters.
complete
Some people prefer complete linkage because it makes compact clusters. However, this is in part an artefact of the method: the clusters are not allowed to grow, because the complete linkage criterion would be violated.
average
Some people (I included) prefer average linkage clustering, because it seems to be a compromise between the previous two extremes, and more neutral in grouping. There are several alternative methods loosely connected to “average linkage” family. Ward’s method seems to be popular in publications. It approaches complete linkage in its attempt to minimize variances in agglomeration. The default "average" method is the one often known as upgma which was popular in old-time genetics
と記載されているので、{vegan}の作者に則って"average"を使っていきましょう。 ( 何か見解があるのであれば教えて下さい )
hclust(vegdist, method="average") -> hclust
hclustを実行すると
> hclust Call: hclust(d = vegdist, method = "average") Cluster method : average Distance : bray Number of objects: 20
と返って来ます。
E. 結果について検討を行う
まずはdendrogramを作成してデータの全体像を考えてみましょう。Dendrogramの描き方に関しては、ggdendro::ggdendrogram()を使う方法とplot()を使う方法があります。
ggdendro::ggdendrogram()
ggdendro::ggdendrogram()に関しては@kazutanさんのRPubsの記事を参考にしました。{ggdendro}の使い方がわかりやすく書かれています。
library(ggdendro) ggdendrogram(hclust, rotate = TRUE, # TRUEで横方向になる size = 2)+ coord_flip()
これを見るといくつかの似た者同士のグループにクラスタリングできそうだな？というようにわかると思います。
plot()とrect.hclust()
どちらかというとこちらの方がおすすめです。plot()もrect.hclust()もRに標準でインストールされています。まずplot()を使ってみましょう。
plot(hclust)
このようなdendrogramが得られます。これを見て、ああ、3グループくらいに分けられそうだな？と思うわけです。そこで役に立つのがrect.hclust()で実際にグラフを見た方がわかりやすいのでやってみましょう。
rect.hclust(hclust, 3)
視覚的にもわかりやすくなりましたね。
4. 実際にクラスタリングする
それでは今までの解析結果をもとに、実際に3つに群分けしていきましょう。群分けには{cutree}という関数を使います。kでクラスターの数でクラスタリングするとともに、hでdendrogramのheightで区切ることもできます。
cutree(hclust, k=3) -> Group # クラスターの数
> Group 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 3 1 3 2
このようにそれぞれのsiteがどのグループに属するのか？ということをvectorで返してくれます。これをもとのデータに組み込めば解析ができます。
5. クラスタリング後に別のパラメーターを比較する
duneのメタデータにはA1というパラメーターがあって
a numeric vector of thickness of soil A1 horizon.
どうやら土壌の厚さを示すようですね。クラスタリングの結果が、こういったパラメーターの違いを表しているかどうかを見ることは菌叢解析の大きな目的にもなります。実際にくクラスター1から3の間でA1のスコアにこれだけ違いがありますので、duneにおいてA1は菌叢の全体像に大きな影響を与える因子といっていいようです。
# duneのメタデータを読み込む data(dune.env) # Boxplot boxplot(A1 ~ Group, data=dune.env, notch = TRUE)
上記のようなboxplotが得られるのでGroup2が特にA1の値が大きいということがわかりますね。
また次元削減をしたplotにクラスタリングの結果を対応させることもできます。
ord <- cca(dune)
plot(ord,
display = "sites")
ordihull(ord, Group, lty = 2, col = "red")
この様にクラスタリングすることが視覚的にもわかりやすくなりましたね。以上になります。
参考: cca()
ccaという次元削減法は聞いたことがなかったので調べてみました。土壌微生物生態研究における正準対応分析 ( CanonicalCorrespondence Analysis： CCA ) の利用法によると
( ccaとは ) 生態系内の種構成と環境要因との関係を分かりやすく解析することを目的に開発された、比較的新しい多変量解析手法である
環境境変化に対する種の反応モデルには、線形モデル(linear) と山型モデル(unimodal)の2種麹がある 例えば、ある環境の水分条件の変化に対して、水分量が多くなればなるほど個体数を増やす種と、ある一定の水分量（至適水分量）までは個体数を増やすが、それ以上に水分量が多くなると、逆にその佃桝数を減らしていく種があるとする この場合、水分条件に対して、前者が線型モデルで反応する種、後者が山型モデルで反応する種、という事になる
CCAは「種の分布が一山型モデルのデータセットに対して優れた長所を発揮する種構成データに環境要因データも組み込んで序列化を行う直接傾度分析法」と定義する事が出来る
ということでした。
