統計解析言語「R」を用いた生物多様性指数・群集構造解析(クラスター解析・NMDS)の方法を解説!

生態系
統計解析言語「R」を用いた生物多様性指数・群集構造解析(クラスター解析・NMDS)の方法を解説!

生物多様性指数の算出や群集構造解析(クラスター解析・NMDS)を行いたい場合、フリーソフトウェアで統計解析向けのプログラミング言語『R』を用いて行うのが最適でしょう。しかし、私が学生の頃は資料が少なく、コードを書くのに苦労したものです。変数の代入(定義)シンボルや引数の表現の揺れや関数のネーミングなどもサイトによって異なり、かなり不統一でした。今は表現の統一化が進み、ChatGPTといった優れた生成AIの登場によって、かなり読みやすい言語にはなったとは思いますが、改めてこの記事では私が行った土壌動物の同定データから群集構造解析を行った例をご紹介してみようと思います。

スポンサーリンク

ある地点での土壌動物分析例

以下はある4地点、2シーズンのツルグレン法によって抽出された土壌動物サンプルの同定結果です。コウチュウ目を除く種を私が担当しました。このデータで生物多様性指数算出・群集構造解析を行っていきましょう。

以下は生データであるため、実際にはここから地点・調査日ごとの入力種名の個体数を集計したものを解析することになりますが、ここでは省略します。

入力種名性別個体数調査日調査地点
ヨロイジュズダニ 25月A
ハナビラオニダニ 25月A
Xylobates属sp. 215月A
ヒラタオニダニ 15月A
アジアオニダニ 15月A
ヒワダニ 205月A
ワタゲジュズダニ 35月A
ジュズダニ科sp. 85月A
ヤマトクモスケダニ 25月A
クワガタダニ 35月A
ヒメヘソイレコダニ 25月A
オバケツキノワダニ 15月A
トウキョウフリソデダニ 35月A
Veigaia属sp. 95月A
Parasitus属sp. 25月A
イトダニ科sp. 95月A
Podocinum属sp. 15月A
Pachylaelaps属sp. 15月A
Epicrius属sp. 65月A
トゲダニ目spp. 245月A
シマツノアヤトビムシ 15月A
アイイロハゴロモトビムシ 25月A
コガタドウナガトビムシ 75月A
ニッポンシロトビムシ 15月A
フクロムラサキトビムシ 35月A
ヨシイフォルソムトビムシ 15月A
カマアシムシ目sp. 25月A
ニセハリアリ 65月A
ウメマツオオアリ 15月A
テラニシハリアリ 15月A
コムカデ綱 25月A
タカミヒトフシムカデ15月A
タジマガハラヒトフシムカデ15月A
Monotarsobius属sp.15月A
オウギヤスデ 15月A
ミコシヤスデ科sp. 15月A
フイリワラジムシ 15月A
ナガワラジムシ 85月A
イトダニ科sp. 15月A
ヒワダニ 25月A
ヒメミミズ科sp. 105月A
アイイロハゴロモトビムシ 15月B
シマツノアヤトビムシ 35月B
ヒメトゲトビムシ 95月B
フクロムラサキトビムシ 155月B
コツノアリ 115月B
ミズアブ科sp. 15月B
コムカデ綱 15月B
ヒラタオニダニ 15月B
クワガタダニ 15月B
ワタゲジュズダニ 25月B
Podocinum属sp. 15月B
トゲダニ目spp. 35月B
アザミウマ目sp. 15月B
ヤブコウジハカマカイガラムシ 15月B
ワタナベヒトフシムカデ酷似種♀♂25月B
Monotarsobius属sp.45月B
コムカデ綱 35月B
ワタゲジュズダニ 235月B
ジュズダニ科sp. 25月B
クワガタダニ 75月B
ウスイロヒメヘソイレコダニ 15月B
Xylobates属sp. 55月B
オタフクダルマヒワダニ 15月B
Podocinum属sp. 25月B
Holaspina属sp. 45月B
Pachylaelaps属sp.45月B
Parasitus属sp. 15月B
トゲダニ目spp. 475月B
Abrolophus属sp. 15月B
ハナビラオニダニ 15月C
ヨロイジュズダニ 15月C
ワタゲジュズダニ 105月C
Xylobates属sp. 465月C
ウスイロヒメヘソイレコダニ 15月C
アラメイレコダニ 15月C
ヤマトクモスケダニ 35月C
ヒワダニ 25月C
クチバシツブダニ 15月C
ツルトミイチモンジダニ 15月C
Podocinum属sp. 35月C
イトダニ科sp. 55月C
Veigaia属sp. 15月C
Holostaspella属sp. 45月C
Holaspina属sp. 25月C
トゲダニ目spp. 125月C
ケダニ目sp. 15月C
カマアシムシ目sp. 25月C
フクロムラサキトビムシ 25月C
アイイロハゴロモトビムシ 45月C
シマツノアヤトビムシ 35月C
キタウロコアリ 15月C
ニセハリアリ 165月C
テラニシハリアリ 15月C
コツノアリ 265月C
チョウ目sp. 15月C
フイリワラジムシ 55月C
ナガワラジムシ 15月C
ムサシアカムカデ 15月C
オウギヤスデ 25月C
ミコシヤスデ科sp. 25月C
Monotarsobius属sp.65月C
サカヨリチビヒトフシムカデ酷似種15月C
Bothropolys属sp. 15月C
Pheretima属sp. 15月C
ヒメミミズ科sp. 15月C
フタツメフォルソムトビムシ 35月D
フクロムラサキトビムシ 25月D
ニッポンシロトビムシ 35月D
ツチカメムシ科sp. 45月D
ウロコアリ 15月D
コツノアリ 295月D
アメイロアリ 15月D
ミズアブ科sp. 35月D
イトダニ科sp. 35月D
イトダニ科sp. 85月D
Holostaspella属sp. 45月D
トゲダニ目spp. 105月D
Xylobates属sp. 115月D
Lasioseius属sp.♂♀25月D
ミナミクモスケダニ 65月D
トウキョウフリソデダニ 95月D
ワタゲジュズダニ 15月D
ヒメヘソイレコダニ 15月D
ヒメイブリダニ 15月D
コムカデ綱 15月D
ミコシヤスデ科sp. 95月D
ゴシチナガズムカデ 15月D
Monotarsobius属sp. 15月D
フイリワラジムシ 15月D
ナガワラジムシ 95月D
Pheretima属sp. 15月D
ヒメミミズ科sp. 15月D
コツノアリ 18月A
ヤブコウジハカマカイガラムシ. 38月A
Mecistocephalus属sp. 18月A
Prolamnonyx属sp. 18月A
フイリワラジムシ 38月A
ミナミクモスケダニ 38月A
ハナビラオニダニ 18月A
アジアオニダニ 28月A
ナミツブダニ 18月A
オバケツキノワダニ 18月A
クワガタダニ 28月A
ツルトミイチモンジダニ 18月A
ヒワダニ 18月A
ワタゲジュズダニ 18月A
ジュズダニ科sp. 18月A
Xylobates属sp. 18月A
Podocinum属sp. 28月A
トゲダニ目spp. 48月A
Pheretima属sp. 18月A
シマツノアヤトビムシ 18月A
アイイロハゴロモトビムシ 128月A
シロフォルソムトビムシ 228月A
コガタドウナガツチトビムシ 28月A
Pseudachorutes属 28月A
ニセハリアリ 48月A
アメイロアリ 48月A
キタウロコアリ 38月A
ヒラタウロコアリ 38月A
コツノアリ 48月A
ホシカメムシ科sp.幼虫 18月A
ナミツブダニ 18月A
コムカデ綱 38月B
キリフリヒトフシムカデ18月B
ワタナベヒトフシムカデ酷似種18月B
ナスヒトフシムカデ18月B
Monotarsobius属sp. 38月B
オカダンゴムシ 18月B
Xylobates属sp. 18月B
フトツツハラダニ 108月B
トウキョウフリソデダニ 18月B
Holaspina属sp. 38月B
Podocinum属sp. 18月B
トゲダニ目spp. 58月B
ニホンケシガイ 18月B
カマアシムシ目sp. 18月B
シマツノアヤトビムシ 18月B
ヒメフォルソムトビムシ 498月B
アイイロハゴロモトビムシ 28月B
フタツメフォルソムトビムシ 28月B
ヨシイフォルソムトビムシ 58月B
ヒメトゲトビムシ 28月B
フクロムラサキトビムシ 18月B
キタウロコアリ 18月B
ニセハリアリ 18月B
ウメマツアリ 18月B
コツノアリ 698月B
イガウロコアリ 18月B
コメツキムシ科sp. 18月B
フトツツハラダニ 18月B
ヤブコウジハカマカイガラムシ. 38月C
コムカデ綱 58月C
Monotarsobius属sp.48月C
Mecistocephalus属sp. 28月C
ヒロズジムカデ 38月C
フイリワラジムシ 68月C
イトダニ科sp. 68月C
Holostaspella属sp. 48月C
トゲダニ目spp. 68月C
Xylobates属sp. 88月C
ハナビラオニダニ 18月C
アジアオニダニ 68月C
ヒメイブリダニ 18月C
カマアシムシ目sp. 88月C
コガタドウナガツチトビムシ 138月C
シマヤマトビムシ 28月C
アイイロハゴロモトビムシ 28月C
フクロムラサキトビムシ 28月C
キイロシリアゲアリ 18月C
コツノアリ 58月C
キタウロコアリ 28月C
コナカイガラムシ科sp. 28月D
コツノアリ 18月D
ミコシヤスデ科sp. 18月D
オウギヤスデ 38月D
ハガヤスデ属 28月D
ヒロズジムカデ 28月D
Monotarsobius属sp. 18月D
フイリワラジムシ 98月D
Xylobates属sp. 18月D
ヤマトクモスケダニ 48月D
Lasioseius属sp. 38月D
Cybaeus属sp. 28月D
アイイロハゴロモトビムシ 198月D
ヨシイフォルソムトビムシ 28月D
ヒメフォルソムトビムシ 68月D
コガタドウナガトビムシ 28月D
Pseudachorutes属sp. 48月D
カマアシムシ目sp. 18月D
ヤマトシロアリ 198月D
コツノアリ 58月D
コメツキムシ科sp.幼虫 18月D
ある地域の土壌動物群集データ

Rで群集構造解析を行うには?

以上のデータから『R』で群集構造解析を行うことができます。以下がその全てのコードです。今回は『R 4.3.2』を使用しています。基本的な内容は高田(2018)を参照しました。土居・岡村(2011)もご参照ください。

required_packages <- c("here", "vegan", "labdsv")

# パッケージがインストールされていなければインストールする
installed_packages <- rownames(installed.packages())
for (pkg in required_packages) {
  if (!pkg %in% installed_packages) {
    install.packages(pkg)
  }
}

# パッケージをロードする
lapply(required_packages, library, character.only = TRUE)

#ディレクトリの指定
getwd()
script_directory <- here("./土壌動物の室内分析と群集構造解析の例")
setwd(script_directory)

#分析データの読み込み
soil_data <- read.csv("土壌動物の個体数データ.csv", header = TRUE, row.names = 1, fileEncoding = "UTF-8")

# 生物多様性の計算
shannon_diversity <- diversity(soil_data, index = "shannon")
shannon_diversity

simpson_diversity <- diversity(soil_data, index = "simpson")
simpson_diversity

renyi(soil_data)

renyi_diversity <- exp(renyi(soil_data))
print(renyi_diversity)

# メタデータの多次元尺度構成法(MDS)解析
soil_data_matrix <- data.matrix(soil_data) # 行列に変換
mds_result <- metaMDS(soil_data_matrix, dist = "bray", autotransform = FALSE)
mds_result
ordipointlabel(mds_result, display = "sites", cex = 1, col = "black", pch = 16)


# クラスター解析
veg_dist <- vegdist(soil_data_matrix, method = "bray")
veg_dist

hierarchical_clustering <- hclust(veg_dist, method = "ward.D2")
hierarchical_clustering
plot(hierarchical_clustering, hang = -1)

# クラスタ数の評価
average_silhouette <- numeric(nrow(soil_data))
for (i in 2:(nrow(soil_data) - 1)) {
  silhouette_vals <- silhouette(cutree(hierarchical_clustering, k = i), veg_dist)
  average_silhouette[i] <- mean(silhouette_vals[, 3])
}
average_silhouette

# クラスタリングの結果
optimal_cluster <- cutree(hierarchical_clustering, k = which.max(average_silhouette))

# Silhouetteプロット
silhouette_vals <- silhouette(optimal_cluster, veg_dist)
rownames(silhouette_vals) <- row.names(soil_data)
plot(silhouette_vals)

# PAM法によるクラスタリングと重要な生物種の特定
pam_result <- pam(veg_dist, k = which.max(average_silhouette), diss = TRUE)
ind_val_analysis <- indval(soil_data, pam_result$clustering, numitr = 5000)
significant_groups <- ind_val_analysis$maxcls[ind_val_analysis$pval <= 0.05]
significant_species <- ind_val_analysis$indcls[ind_val_analysis$pval <= 0.05]
p_values <- ind_val_analysis$pval[ind_val_analysis$pval <= 0.05]
frequencies <- apply(soil_data > 0, 2, sum)[ind_val_analysis$pval <= 0.05]
significant_results <- data.frame(group = significant_groups, indval = significant_species, pvalue = p_values, freq = frequencies)
significant_results_summary <- significant_results[order(significant_results$group, -significant_results$indval), ]
significant_results_summary

データファイルの読み込み方は?

最初に今回使用するパッケージをインストールしておきます。既にインストールしている人は不要です。

required_packages <- c("here", "labdsv", "vegan")
install.packages(required_packages)

install.packages(required_packages)の代わりに以下のコードを使用すると既にインストールしているパッケージはスキップします。

installed_packages <- rownames(installed.packages())
for (pkg in required_packages) {
  if (!pkg %in% installed_packages) {
    install.packages(pkg)
  }
}

慣れないうちは個別に読み込んでも良いですが、上のコードの方がロードのときにベクトルにしたパッケージ名を使い回せます。

install.packages(here)
install.packages(labdsv)
install.packages(vegan)

そして、パッケージをロードしましょう。

lapply(required_packages, library, character.only = TRUE)

同様に慣れないうちは個別に読み込んでも構いません。

library(here)
library(labdsv)
library(vegan)

次にディレクトリを指定します。ディレクトリは絶対パスでも指定は可能ですが、より簡単な現在使用しているワークスペースを作業ディレクトリに指定してくれるhereパッケージのhere()関数を使用しています。here()関数を使用するとGoogleドライブなどのクラウドストレージでファイル共有したり、別PCに以降したときでもパスが有効なので便利です。

ただ、注意したいのがhere()関数はワークスペースのパスを呼び出す関数なので、環境によっては(私の場合はVS code)、開いている.Rファイルの上のパスまでは指定してくれないので、ここでは1つ上のパス名を手入力しています。この点はなにか改善点が欲しいところです。

getwd()
script_directory <- here("./土壌動物の室内分析と群集構造解析の例")
setwd(script_directory)

そして、.csvファイルを読み込みます。注意したいのが.csvファイルの文字コードはUTF-8でなければなりません。この仕様は最近のバージョンで追加されたようで古いRコードで実行しようとした時に困りました。(これはもしかしたら私がVS codeで実行しているせいかもしれません。)

OSがWindowsのPCでは『Microsoft Excel(エクセル)』などで.xlsxファイルから.csvファイルを作成すると、文字コードがShift-JIS(ANSI)の.csvファイルが作成されてしまいます。

そのため、「名前をつけて保存」をするときに「CSV」ではなく、「CSV UTF-8」と書いてある方を選択しましょう。これは本当に勘違いしやすいです。

一応、read.csv()関数の引数を「fileEncoding = "Shift-JIS"」とすることでShift-JISでも読み込めることになっていますが、そもそもShift-JISが古い規格で現在ではUTF-8が上位互換の文字コードであり、そのままの.csvファイルを使用するのはおすすめできません。

以下では念の為「fileEncoding = "UTF-8"」で指定しています。

soil_data <- read.csv("土壌動物の個体数データ.csv", header = TRUE, row.names = 1, fileEncoding = "UTF-8")

以上でデータの読み込みは完了です。

なお、.csvファイルの読み込みでよくある失敗は列と行の向きとラベルの有無です。

実際に今回読み込んでいる.csvファイルは上記の表の列と行が反対になったものを読み込んでいます(「性別」の行も削除しています)。『Microsoft Excel』を使用している場合は、右クリックして「貼り付けのオプション」から「行/列の入れ替え」を選択することで簡単に反転させられます。

また列と行どちらにもラベルがある場合は、「header = TRUE, row.names = 1」と指定しましょう。

Rを用いてシャノンの多様度指数とシンプソンの多様度指数を計算するには?

まずは複数ある生物多様度指数を計算してみましょう。

シャノンの多様度指数(Shannon index)を計算し、内容を確認するには以下のコードを使用します。

shannon_diversity <- diversity(soil_data, index = "shannon")
shannon_diversity

今回は以下の数値が得られました。

A(5月)B(5月)C(5月)D(5月)A(8月)B(8月)C(8月)D(8月)
3.0814132.7227612.7395022.7013952.8206751.9581882.8295262.496044
ある地域の土壌動物群集データから得られたシャノンの多様度指数

シンプソンの多様度指数(Simpson index)を計算し、内容を確認するには以下のコードを使用します。

simpson_diversity <- diversity(soil_data, index = "simpson")
simpson_diversity

今回は以下の数値が得られました。

A(5月)B(5月)C(5月)D(5月)A(8月)B(8月)C(8月)D(8月)
0.93199030.90616340.88243650.89744330.89550170.74206090.93037040.8756084
ある地域の土壌動物群集データから得られたシンプソンの多様度指数

いずれの多様度指数でもA(5月)が最も多様性が高く、B(8月)が低い地点となっています。一般に8月の方が昆虫も増え、多くなるようにも思えますが、意外な結果となりました。

Rを用いてNMDSをプロットするには?

次に多変量解析の一種である非計量多次元尺度法(Nonmetric Multidimensional Scaling、NMDS)を行っていきましょう。NMDSはデータセット内の個々のデータポイント間の類似性または距離を反映するように、データを低次元の空間に配置する手法です。

要するに今回の場合は図の中でプロットが近いほど生物群集構造が近いということを意味します。

Rで非計量多次元尺度法をプロットするには以下のコードを使用します。

soil_data_matrix <- data.matrix(soil_data) # 行列に変換
mds_result <- metaMDS(soil_data_matrix, dist = "bray", autotransform = FALSE)
mds_result
ordipointlabel(mds_result, display = "sites", cex = 1, col = "black", pch = 16)

その結果、以下の図が得られました。

ある地域の土壌動物群集データから得られた非計量多次元尺度法(NMDS)による解析の結果

NMDSからは季節による土壌動物群集の影響が考えられます。一般に土壌動物は季節の影響を受けにくいとされていますが、今回は夏季に昆虫類が増加したことが影響したと思われます。また地点の影響も出ていますが、D地点では大きく異なる結果になっています。D地点は今回のデータでは地点の環境に関する情報がありませんが、季節で大きく環境がかわる地点であった可能性もあります。なお、環境に関する量的データがあれば更にその影響を評価したものをプロットできます。

Rを用いてクラスター解析を行うには?

次にクラスター解析です。クラスター解析はデータセット内の個々のデータ点を似た特性や属性に基づいてグループにまとめる統計的手法で、データ内のパターンや構造を把握し、データを理解することができます。今回は地点間で群集構造が近い順にまとめることになります。クラスター解析とプロットは以下のコードで行えます。

今回の解析にはブレイ・カーチス非類似度(Bray-Curtis dissimilarity)を使用しています。

veg_dist <- vegdist(soil_data_matrix, method = "bray")
veg_dist

hierarchical_clustering <- hclust(veg_dist, method = "ward.D2")
hierarchical_clustering
plot(hierarchical_clustering, hang = -1)

その結果、以下の図が得られました。

ある地域の土壌動物群集データから得られたクラスター解析の結果

これを以下のコードでシルエット分析(Silhouette Analysis)を行います。クラスタリングの結果を評価し、各データポイントがその所属するクラスター内でどれだけ密集しているかを測定する手法です。シルエット分析はクラスタリングの妥当性を評価するための一つの手法であり、クラスターの数や形状を決定する際に有用です。

要するに恣意的になりがちなグループ分割を避けて、意味がありそうなグループの分割基準を調べてくれる便利なものです。

silhouette_vals <- silhouette(optimal_cluster, veg_dist)
rownames(silhouette_vals) <- row.names(soil_data)
plot(silhouette_vals)

その結果、以下の図が得られ、2つ(A(8月) とD(8月)のクラスターとその他のクラスター)に分割することが最適と判断されました。

ある地域の土壌動物群集データから得られたシルエット分析の結果

cutree()関数を使用すると改めて各地点が具体的にどのグループに属するかを調べることができます。

A(5月)B(5月)C(5月)D(5月)A(8月)B(8月)C(8月)D(8月)
11112112
ある地域の土壌動物群集データから得られた各地点のグループ所属

クラスター解析では今回の分析例では土壌動物の種ごとの生態に関する情報が極めて少なく、また地点の植生など環境情報が不足しているので考察は難しいですが、これらの情報も加味すればより深い考察が可能です。

Rを用いて指標種分析を行うには?

最後にクラスター解析に基づいた指標種分析を行います。以下のコードで行えます。PAM法を使用して、最適なクラスタ数を決定し、データをクラスタリングし、クラスタリング結果を元に、各クラスタごとの重要な生物種を特定します。 そして、特定された重要な生物種に関する情報(クラスター・種・p値・出現頻度)を収集し、クラスタごとに重要な生物種をまとめ、結果を要約しています。

pam_result <- pam(veg_dist, k = which.max(average_silhouette), diss = TRUE)
ind_val_analysis <- indval(soil_data, pam_result$clustering, numitr = 5000)
significant_groups <- ind_val_analysis$maxcls[ind_val_analysis$pval <= 0.05]
significant_species <- ind_val_analysis$indcls[ind_val_analysis$pval <= 0.05]
p_values <- ind_val_analysis$pval[ind_val_analysis$pval <= 0.05]
frequencies <- apply(soil_data > 0, 2, sum)[ind_val_analysis$pval <= 0.05]
significant_results <- data.frame(group = significant_groups, indval = significant_species, pvalue = p_values, freq = frequencies)
significant_results_summary <- significant_results[order(significant_results$group, -significant_results$indval), ]
significant_results_summary

このコードを実行すると以下の表が得られます。

groupindvalpvaluefreq
フクロムラサキトビムシ110.03786
Pseudachorutes属sp.210.03742
ある地域の土壌動物群集データから得られた指標種分析の結果

これにより、A(8月) とD(8月)のクラスターはヤマトビムシ属(Pseudachorutes)で指標され、その他のクラスターはフクロムラサキトビムシで指標されることがわかりました。残念ながらこれらの種についても生態に関する情報が不足していますが、よく生態が知られている種で分析を行えば、より深い考察が可能です。

これらの解析は昆虫・脊椎動物・植物・水生生物・細菌など様々な生物群集でも可能です。

引用文献

土居秀幸・岡村寛. 2011. 生物群集解析のための類似度とその応用:Rを使った類似度の算出、グラフ化、検定. 日本生態学会誌 61(1): 3-20. https://doi.org/10.18960/seitai.61.1_3

高田宜武. 2018. 水産研究・教育機構「日本海区水産研究所」:Rによる群集組成の解析. https://jsnfri.fra.affrc.go.jp/gunshu/index.html

タイトルとURLをコピーしました