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

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)

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)

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

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)

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

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