統計解析.com

統計解析業務の
アウトソーシング

無料相談を
予約する

Rによる多様性指数の算出

アウトライン
  1. 作成日: 2025/12/1
  2. 更新日: –

はじめに

ヒトの腸内には、数百〜数千種類もの腸内細菌が存在しており、その集合体を腸内細菌叢 (gut microbiota) と呼びます。腸内細菌は一般的に、体に良い働きをするとされる善玉菌 (例: ビフィズス菌、乳酸菌)、悪い働きをするとされる悪玉菌 (例: ウェルシュ菌)、そして環境に応じて良くも悪くも働く日和見菌の3つに大別されます。

しばしば「善玉菌が多いほうが良い」と考えられがちですが、実際には必ずしもそうではありません。最新の腸内細菌研究では、特定の菌だけが多いよりも、多様な菌がバランス良く存在していること (=腸内細菌叢の多様性が高いこと) のほうが、健康との関連が強いことが報告されています。

たとえば、腸内細菌叢の多様性は以下のような健康指標と関連するとされています:

  • 免疫機能の維持

  • 代謝バランス

  • 肥満・糖代謝との関連

  • 心理的ストレスやメンタルヘルスとの関連

こうした背景から、近年は腸内環境の改善を目的とする機能性表示食品でも、「腸内フローラのバランス維持」「ビフィズス菌の増加」「短鎖脂肪酸 (SCFA) 産生を介した整腸作用」などを訴求する製品が増えています。

ただし、機能性表示食品の関与成分が腸内環境に与える作用を評価する際には、単純に善玉菌・悪玉菌の増減を見るだけでは不十分で、多様性の観点を含めた解析設計が重要になります。

腸内細菌叢

腸内細菌叢とは?

健常な人の腸内には、100種類以上の菌種が、総菌数にして約100兆個も生息しています。その中でも最も優勢なのがBacteroides属で、菌数は約10¹¹/gとされています。次いで、Bifidobacterium属、Eubacterium属、Clostridium属、Peptostreptococcus属 などが続きます。一方、Escherichia属 (大腸菌) や Enterococcus属などは、総菌数の1/1000以下とされています。

腸内細菌叢は、加齢や食生活の変化によって構成が変動します。新生児では、出生後数時間で細菌が定着し始め、生後3〜4日でBifidobacterium属 が急増。生後1週間で10¹⁰〜10¹¹/gに達し、離乳期まで最優勢の状態が続きます。

離乳期になるとBifidobacterium属 がやや減少し、Bacteroides属、Eubacterium属、Clostridium属が増加。これにより、成人の腸内細菌叢に近づいていきます。通性嫌気性菌群は偏性嫌気性菌群よりも少ない菌数で、安定した成人型の腸内細菌叢が形成されます。

老年期になるとBifidobacterium属 は減少し、検出されなくなる場合もあります。また、加齢に伴いClostridium perfringens (ウェルシュ菌) や大腸菌などの腸管内腐敗菌が増加し、老人型の腸内細菌叢が形成されます。

さらに、腸内細菌叢は食生活の影響も大きく受けます。たとえば、食物繊維や発酵食品の摂取によって腸内が酸性に傾き、Bifidobacterium属やLactobacillus属などの「善玉菌」が増加する傾向があります。一方、動物性たんぱく質を多く摂取すると、Clostridium属の一部や腸内細菌科の細菌など、「悪玉菌」が増加する傾向があります。こうした菌のバランスの変化は、腸内環境の健全性に大きく関与しています。

善玉菌、悪玉菌、日和見菌について

腸内細菌は、主に以下の3つに分類されます。

概要 主な働き
善玉菌 腸内環境を整え、健康維持に役立つ細菌です。代表的なものに 乳酸菌 や ビフィズス菌 があります。これらは消化吸収のサポートや免疫機能の調整など、体に良い働きをしてくれます。
  • 悪玉菌の発生を抑える

  • 腸の動きを活発にする

  • ビタミン (B1、B2、B6、B12、K、ニコチン酸、葉酸) を生成する

  • 免疫機能を高める

悪玉菌 腸内環境に悪影響を及ぼすとされる細菌です。ウェルシュ菌 や一部の サルモネラ菌 などが代表的で、食中毒や炎症の原因になることがあります。脂肪分の多い食事や食物繊維の不足によって増殖しやすく、善玉菌より優勢になると健康に悪影響を及ぼす可能性があります。
  • 単糖・脂肪・たんぱく質・アルコール・活性酸素などをエネルギー源として利用

  • 硫化水素やトリメチルアミン-N-オキシドを産生

日和見菌

善玉菌・悪玉菌のいずれにも分類されない腸内細菌です。腸内で優勢な菌に影響されやすく、善玉菌が多いと良い働きをし、悪玉菌が増えると悪影響を及ぼすことがあります。つまり、腸内細菌のバランス次第で性質が変化するのです。

理想的な腸内細菌叢のバランスは、善玉菌2:悪玉菌1:日和見菌7 の割合とされており、多様性が高いことが重要とされています。

2.3 腸内細菌叢の多様性の重要性

腸内細菌の「多様性が高い」とは、腸内に多種多様な細菌が共存している状態のことを指します。

悪玉菌には悪いイメージがあるかもしれませんが、たんぱく質を分解して便として排泄するなど、私たちの体にとって必要な働きも担っています。ただし、過剰になると疾患につながる可能性があるため、多すぎず、少なすぎず、絶妙なバランスで存在することが重要です。

このように、多様性の高い腸内細菌叢は、健康な腸内環境の維持に欠かせない要素とされています。

種多様性の評価方法

3.1 多様性とは?

種多様性は、主に以下の2つの要素によって構成されます。

  • 種の豊富さ (species richness): 群集内に存在する種の数を表し、単純に「何種類の生物がいるか」を示します。

  • 均等度 (species evenness): 各種の個体数のバランスを表し、特定の種が極端に多く他の種が少ない場合は均等度が低く、すべての種がほぼ同じ個体数で存在する場合は均等度が高くなります。

これらの要素を統合して、環境中の種多様性を定量的に評価するために「多様性指数 (diversity index)」が用いられます。多様性指数は、サンプル中で観測されたすべての種の個体数と、それぞれの種の相対的な割合 (全体に占める比率) をもとに算出されます。

3.2 様々な多様性指数 (α多様性)

α多様性

サンプル内における多様性を「α多様性」といいます。腸内細菌叢を解析する場合、サンプルを採取した被験者それぞれの菌叢からα多様性指数を算出します。以下に主なα多様性指数を紹介します。

【Observed features】
Observed featuresは、「各サンプルの中に、いくつの特徴 (feature) が存在したか」を表す指標です。サンプル中にどれだけ多くの種や分類単位が存在するかを表す、純粋な“種数 (species richness)”に相当します。

ここでいう feature とは、ASV (Amplicon Sequence Variant:実際の塩基配列に基づく高解像度の分類単位)、OTU (Operational Taxonomic Unit: 通常97%の塩基配列類似度でクラスタリングされた分類単位)、種レベルに集約した分類単位などを指します。

\[S_{\text{obs}} = \sum_{i = 1}^{S}{I\left( n_{i} > 0 \right)}\]

ni: i番目の種・ASV・OTUのリード数
I (n_i > 0) : その特徴が観測されたら1、されなければ0

観測された feature は、1回でも出現すれば1種類として数えるため、サンプル中の相対的な存在比 (全体に占める比率) は一切考慮されません。

【Chao1】
Chao1指数は、観測データから「稀にしか出現しない種」を補正して、真の種数を推定するための指標です。観測された種数に、希少種の情報 (1回だけ観測された種=singleton、2回だけ観測された種=doubleton) から計算した補正項を足して「まだ見えていない種の分」を見積もるという考えに基づいています。

\[Chao1 = S_{\text{obs}} + \frac{{RS}_{1}^{\, 2}}{2{RS}_{2}}\]

Sobs: 観察された種数 (Observed features)
RS1: 1個体だけ観測された種 (singleton) の数
RS2: 2個体だけ観測された種 (doubleton) の数

RS1が少ないほど、補正項の値が小さくなるため、低頻度OTU除去の状況やレアファクション深度によって算出結果が影響されます。RS1=0の場合は、希少種が極端に少ないか除去済みであるため、Chao1=Observed featuresになります。

【Simpson のλ】
Simpson 指数は、単に種数を数えるだけでなく、各種の個体数 (相対存在量) を考慮して算出される多様性指標です。この指数は確率的な発想に基づいており、「サンプル中から 2 個体を無作為に取り出したとき、同じ種である確率」を表しています。

\[1 – \lambda = 1 – \sum_{i = 1}^{S}\left( \frac{n_{i}}{N} \right)^{2}\]

S: 群集中に存在する種の総数
ni: i番目の種の個体数
N: 群集全体の個体数 (niの合計)
ni​/N: 群集の中でその種が占める割合 (相対存在量)

サンプル内に相対存在量の大きい種がある場合、その種が占める割合が大きいため、Simpsonのλの値は高くなり、1に近づきます。逆に、種数が多く、各種がほぼ均等な割合で存在している場合は、λの値は小さくなります。

多様性指数として用いる際には、種の構成が均等であるほど値が大きくなる方が直感的に理解しやすいため、1からλを引いた値を「Simpson指数」として算出します。

【Shannon–Wiener の H′】
Shannon–Wiener (Shannon) 指数 H′は、サンプル中の「richness (種数) 」と「evenness (均等さ)」の両方を考慮して算出する多様性指標です。種が多く、かつ各種の出現割合が均等に近いほど値が大きくなります。

\[H' = – \sum_{i = 1}^{S}p_{i}\log p_{i}\]

S: 種数 (観測された種の総数)
pi: 種 i の相対存在量、総カウント数に対する割合 (p_i = n_i / N)
logは底2 (log₂) の場合と自然対数 (ln) の場合があります。

Shannon–Wiener指数の計算式には、先頭にマイナス記号 (−) がついています。これは、p_i (各種の相対存在比) が0より大きく1未満 (0 < p_i < 1) であるため、対数log (p_i) の値が常に負になることに由来します。

そのまま計算すると、種の分布が均等で多様性が高いほど、指数の値が小さくなってしまい、直感的に理解しづらいため、式全体にマイナスをかけることで、多様性が高いほど指数の値が大きくなるように調整しています。

S (観測された種の総数) が大きく、p_i (各種の相対存在比) が均等であるほどH′が大きくなります。

3.3 Rarefaction (レアファクション)

生物の「種数」は、実はサンプルサイズ (観察した個体数) に大きく依存します。そのため、単純に種数だけを比較しても本来の多様性を正しく評価できません。

この“サンプルサイズの影響”を取り除くために行うのが Rarefaction (レアファクション:希薄化) です。

例えば、2つの動物園 (動物園A、動物園B) で観察した動物から、各動物園の生物多様性を算出して比較したいとします。

  • 動物園A: 300種、2500個体の動物を発見

  • 動物園B: 170種、1000個体の動物を発見

一見すると、種数の多い動物園Aの方が、多様性が高そうに見えます。しかし、ここで注意したいのは、観察した個体数がAの方が圧倒的に多いという点です。観察した個体数が多ければ、多くの種に出会えるのは当然です。つまり、このままでは「種数の差」ではなく、「観察した量の差」を比較してしまっていることになります。

そこで、“サンプルサイズの影響”を取り除くためにレアファクションを行います。これは、両方の動物園から同じ数の個体をランダムに抽出し、そこから得られる種数を比較するという方法です。今回の例では、観察個体数が少ない動物園Bに合わせて、1000個体をランダムに取り出して比較するイメージです。こうすることで、個体数の差によるバイアスが取り除かれ、より公平に多様性を評価できます。

たとえ動物園Aでは2500個体の動物を観察しても、レアファクションを行うことで「もし1000個体しか観察しなかったら、何種の動物が見つかりそうか」という推定が可能になります。

レアファクションの考え方は腸内細菌叢解析でもまったく同じです。サンプル (被験者) ごとに得られるリード数が大きく異なるため、そのままでは多様性の正確な比較ができません。
そこで、レアファクションを行い、「同じ量だけ調べた場合の多様性」を仮定することで、公平な比較ができるようになるのです。

Rで多様性指数を算出する (公開データセット)

それでは、実際にRを使用して多様性指数を算出します。今回は、phyloseqパッケージ、microbiomeパッケージを使用して多様性指数を算出します。

使用パッケージ

【phyloseq】
微生物叢データ (featureカウント、分類、サンプル情報、系統樹など) をR内で統合的に管理し、解析、可視化するためのパッケージ。作成するデータセットは全て.tsvファイルで準備 (またはBIOM形式) する。

【microbiome】
phyloseqで扱ったデータをさらに加工・変換・要約し、豊富な指標や可視化で解析を深めるためのパッケージ

4.2 スクリプト

以下がスクリプトです。
今回は、phyloseqで利用可能な公開データセットを使用します。

# ----------------------------
# ライブラリ読み込み
# ----------------------------
library (phyloseq)
library (microbiome)
library (openxlsx)

# ----------------------------
# 1. 操作設定:データセット名と希薄化設定
# ----------------------------

dataset_name <- "soilrep" # 使用するデータセット名 (例:"GlobalPatterns", "esophagus", "soilrep")

perform_rarefaction <- FALSE # TRUE (希薄化あり) / FALSE (希薄化なし)
use_min_sample_reads <- TRUE # TRUE: 最小リード数を使用 / FALSE: 指定深度を使用
rarefaction_depth <- 10000 # use_min_sample_reads が FALSE の場合に使用

#シード値
master_seed <- 1234
set.seed (master_seed)

# ----------------------------------------------------------
# 2. データセットの読み込み
# ----------------------------------------------------------

data (list = dataset_name, package = "phyloseq") # phyloseq 内のデータ読み込み
ps <- get (dataset_name)  # 読み込んだデータをオブジェクトへ

# ----------------------------------------------------------
# 3. 希薄化処理 (Rarefaction)
# ----------------------------------------------------------

if (perform_rarefaction) {

 sample_reads <- sample_sums (ps)

 if (use_min_sample_reads) {
 rarefaction_depth <- min (sample_reads)
 message ("最小リード数を希薄化深度として使用します:", rarefaction_depth)

 } else {
 if (any (sample_reads < rarefaction_depth) ) {
 warning ("指定深度より低いリード数のサンプルが存在します。")
 }
 message ("指定希薄化深度で処理します:", rarefaction_depth)
 }

 # 希薄化処理 (rngseed に master_seed を渡す)
 ps <- rarefy_even_depth (
 ps,
 sample.size = rarefaction_depth,
 rngseed = master_seed,
 verbose = TRUE,
 replace = FALSE,
 trim = TRUE
 )

} else {
 message ("希薄化をスキップします。")
}

# ----------------------------------------------------------
# 4. α多様性の算出
# ----------------------------------------------------------

alpha_div <- microbiome::alpha (
 ps,
 index = c ("observed", "chao1", "shannon", "simpson")
)

# メタデータと統合
meta <- data.frame (sample_data (ps) )
alpha_with_meta <- cbind (meta, alpha_div)

# ----------------------------------------------------------
# 5. Excel 出力
# ----------------------------------------------------------

desktop_path <- file.path (Sys.getenv ("USERPROFILE") , "Desktop")
depth_info <- ifelse (perform_rarefaction,
 paste0 ("_depth", rarefaction_depth) ,
 "_raw")

out_path <- file.path (desktop_path, paste0 (dataset_name, "_alpha_diversity", depth_info, ".xlsx") )

write.xlsx (alpha_with_meta, file = out_path, rowNames = TRUE)
cat ("Excelに多様性指数を出力しました:", out_path, "\n")

4.3 ライブラリ読み込み

library (phyloseq)
library (microbiome)
library (openxlsx)

解析で使用するライブラリを指定します。

4.4 操作設定:データセット名と希薄化設定

データセットの指定

dataset_name <- "soilrep" # 使用するデータセット名 (例:"GlobalPatterns", "esophagus", "soilrep")

使用するデータセットを指定します。

今回は、土壌微生物データの"soilrep"というデータセットを指定しています。
"soilrep"は、Clipping (刈取り)、Warming (加温) の2処理について、UU (未処理)、UC (Clippingのみ)、WU (Warmingのみ)、WC (Warming+Clipping) の4つのtreatment条件下における土壌微生物叢が含まれたデータです。

Phyloseqでは他にも、"GlobalPatterns (様々な環境の微生物データ)", "esophagus (食道の微生物データ)"、など様々なデータセットが公開されています。

レアファクション設定

perform_rarefaction <- FALSE # TRUE (希薄化あり) / FALSE (希薄化なし)
use_min_sample_reads <- TRUE # TRUE: 最小リード数を使用 / FALSE: 指定深度を使用
rarefaction_depth <- 10000 # use_min_sample_reads が FALSE の場合に使用

この部分で、希薄化の有無と希薄深度を設定します。

  • perform_rarefactionをFALSEに設定することで、希薄化なしの状態で解析します。TRUEに設定した場合は、希薄化を実行します。その後に続く希薄深度の設定をしてください。
  • use_min_sample_readsをTRUEに設定することで、全サンプルのうち、総リード数が最小のリード数に合わせて希薄深度を設定します。FALSEに設定することで、任意の希薄深度をrarefaction_depthにて設定できます。

シード値

master_seed <- 1234
set.seed (master_seed)

master_seed で任意のシード値を設定し、set.seed () でシード値を設定し、ランダム抽出の再現性を確保しています。

4.5 データセットの読み込み

data (list = dataset_name, package = "phyloseq")

        ps <- get (dataset_name)

data () 関数でデータセット名を格納した” dataset_name”とパッケージを指定して、データセットを呼び出しています。その後、get () 関数でデータセットを取得し、”ps”に格納します。

4.6 希薄化処理

この部分で、実際にレアファクション (希薄化) 処理を行っています。

レアファクション有無の確認

if (perform_rarefaction) {

最初に、ユーザーが「レアファクションを実施する (TRUE)」と指定しているかどうかを確認します。FALSE なら処理は一切行わず、元データのまま解析を進めます。

リード数取得

sample_reads <- sample_sums (ps)

ここでは、phyloseq オブジェクト ps に含まれるすべてのサンプルのリード数を取得します。

レアファクション深度を決定するための基本情報です。

希薄化深度設定 (最小値)

if (use_min_sample_reads) {
rarefaction_depth <- min (sample_reads)
message ("最小リード数を希薄化深度として使用します:", rarefaction_depth)

希薄化設定の際に最小のリード数に合わせると設定していた場合 (use_min_sample_readsをTRUEに設定していた場合)、上記の処理が実行されます。最小のリード数に合わせることで、希薄化によるサンプル除外を避けます。

希薄化深度設定 (指定値)

} else {
 if (any (sample_reads < rarefaction_depth) ) {
 warning ("指定深度より低いリード数のサンプルが存在します。")
 }

ユーザーが任意の深度を指定した場合、その深度よりリード数の少ないサンプルがある場合は警告を出します。そのようなサンプルはレアファクション後に自動的に除外されます。

希薄化処理

message ("指定希薄化深度で処理します:", rarefaction_depth)
 }

 ps <- rarefy_even_depth (
 ps,
 sample.size = rarefaction_depth,
 rngseed = master_seed,
 verbose = TRUE,
 replace = FALSE,
 trim = TRUE
 )

ここが実際の希薄化ステップです。

rarefy_even_depth ()
→サンプル間のリード数を揃えるための希薄化 (Rarefaction) 処理を行う関数

sample.size
→ 先ほど決めたレアファクション深度

rngseed
→phyloseq内部で使用するシード値を設定します。

replace = FALSE
→ 同じリードを重複カウントしない (本来のレアファクション方式)

trim = TRUE
→ 深度未満のサンプルは自動で除外

こうして、すべてのサンプルが同じリード数にそろったデータになります。

} else {
 message ("希薄化をスキップします。")
}

perform_rarefaction = FALSEの場合は、元データで解析を続けます。
レアファクションを行わない場合、上記が実行されます。

4.7 α多様性の算出

α多様性の算出

alpha_div <- microbiome::alpha (
  ps,
  index = c ("observed", "chao1", "shannon", "simpson")
 )
 

microbiomeパッケージのalpha () 関数でサンプルごとのα多様性指数を計算しています。
“ps”はphyloseqオブジェクトです。算出した数値は“alpha_div“に格納しています。

【算出している指数】
"observed" → 観察された ASV/OTU数 (種数)
"chao1" → Chao1推定種数
"shannon" → Shannon多様性指数
"simpson" → Simpson多様性指数

多様性指数データをメタデータと結合

meta <- data.frame (sample_data (ps) )
alpha_with_meta <- cbind (meta, alpha_div)

“ps”に含まれるサンプルのメタデータ (群情報・性別・BMI など) をデータフレーム化し、算出したα多様性指数 (alpha_div”) と列結合しています。
この処理により、各サンプルのメタデータ+多様性指数を1つの表にまとめています。

4.8 Excel 出力

desktop_path <- file.path (Sys.getenv ("USERPROFILE") , "Desktop")

Sys.getenv ("USERPROFILE")
→ “C:/Users/ユーザー名”を取得 (Windowsの場合)

file.path ()
→ OS に合わせた正しいパス区切りで "Desktop" を連結

depth_info <- ifelse (perform_rarefaction,
 paste0 ("_depth", rarefaction_depth) ,
 "_raw")

out_path <- file.path (desktop_path, paste0 (dataset_name, "_alpha_diversity", depth_info, ".xlsx") )

write.xlsx (alpha_with_meta, file = out_path, rowNames = TRUE)

cat ("Excelに多様性指数を出力しました:", out_path, "\n")

depth_info
→レアファクション深度を格納しています。希薄化 しなかった場合は”_raw”となります。

out_path
→ファイルの保存先とファイル名を設定しています。その後、write.xlsx () にてalpha_with_metaを指定パスに Excel 形式で保存しています。

Rで多様性指数を算出する (任意のデータセット)

先ほどまでは、公開されているデータセットを使用して多様性指数を算出するスクリプトを記載しましたが、ここからは、Excel等で手元にデータがある場合、データセットを手動でインポートしてphyloseqで統合する方法をご紹介します。

5.1 phyloseqのデータセット形式

phyloseqは以下の形式のデータセットを読み込むことができます。

  1. otu_table
    OTU テーブルは、サンプル×feature (OTU/ASV) の表です。
    • 行: feature (OTU ID, ASV ID)
    • 列: サンプルID
    • セルの値はリード数 (整数)
    • 行名、列名は重複なし
    • 左上のセルは空欄またはヘッダー名
       例:
          OTU_ID S1 S2 S3
          OTU_1 10 2 0
          OTU_2 0 5 1
          OTU_3 20 0 3
  2. tax_table
    Taxonomyテーブルは、各OTUまたはASVに分類情報を紐づけた表です。解析精度には直接かかわりませんが、プロットの解釈や属レベルでの集計に影響します。
    • 行: feature (OTU IDまたはASV ID)
    • 行名はOTUテーブルの行名と一致する必要あり
    • 欠損値はNAまたは空欄
    • 列名は自由だが、一般的には7階層 (Kingdom → Species) が推奨
       例:
          OTU_ID Kingdom Phylum Class ...
          OTU_1 Bacteria Firmicutes Bacilli ...
          OTU_2 Bacteria Actinobacteria Actino ...
  3. sample_data
    sample_data は、各サンプルに紐づくメタデータ (背景情報や群など) を含む表です。
    • 行名: サンプルID (OTUテーブルの列名と一致)
    • 列: 群情報、環境条件など
    • 数値と文字列が混ざっていない

phyloseqは、上記のデータを統合し、一つのオブジェクトにまとめることができます。

phyloseq の最大のメリットは、OTU/ASV テーブル・分類情報 (taxonomy)・サンプル情報 (sqampledata) をひとつにまとめて扱えることです。別々のファイルを毎回照合したり整形したりする必要がなくなり、ミスが減って解析がスムーズになります。さらに、統合されたオブジェクトを使えば、そのまま多様性解析や可視化まで一貫して実行でき、作業の効率化と再現性の向上につながります。

5.2 スクリプト

手元に腸内細菌叢のデータがある場合、データセットを作成し、下記のスクリプトを「操作設定: データセット名と希薄化設定」、「データセットの読み込み」の部分のスクリプトと置き換えてください。

# ----------------------------
# 1. 操作設定:ファイル読み込みと希薄化設定
# ----------------------------

# デスクトップのパス
desktop_path <- file.path (Sys.getenv ("USERPROFILE") , "Desktop")

# ファイル名 (拡張子は .tsv または .csv に応じて変更)
otu_file <- file.path (desktop_path, "otu_table.tsv")
tax_file <- file.path (desktop_path, "tax_table.tsv")
sample_file <- file.path (desktop_path, "sample_data.tsv")

# 希薄化の有無
perform_rarefaction <- FALSE # TRUE (希薄化あり) / FALSE (希薄化なし)

# 希薄化深度設定
use_min_sample_reads <- TRUE # TRUE: 最小リード数を希薄化深度として使用 / FALSE: 指定の深度を使用
rarefaction_depth <- 10000 # use_min_sample_reads が FALSE の場合のみ使用

#シード値
master_seed <- 1234
set.seed (master_seed)


# ----------------------------
# 2. データセットの読み込み
# ----------------------------

# ファイル読み込み (CSVの場合は read_csv に変更)
otu_df <- read_tsv (otu_file)
tax_df <- read_tsv (tax_file)
sample_df <- read_tsv (sample_file)

# 行名設定 (1列目がID列であることを想定)
otu_df <- column_to_rownames (otu_df, var = colnames (otu_df) [1])
tax_df <- column_to_rownames (tax_df, var = colnames (tax_df) [1])
sample_df <- column_to_rownames (sample_df, var = colnames (sample_df) [1])

# phyloseqオブジェクト作成
otu_table_obj <- otu_table (as.matrix (otu_df) , taxa_are_rows = TRUE)
tax_table_obj <- tax_table (as.matrix (tax_df) )
sample_data_obj <- sample_data (sample_df)

ps <- phyloseq (otu_table_obj, tax_table_obj, sample_data_obj)

以下で変更点を解説します。

なお、レアファクションとシード値に関する設定の部分は、公開データセットから算出するスクリプトと全く同じです。

5.3 操作設定:データセット名と希薄化設定

デスクトップのパス (Windows)

desktop_path <- file.path (Sys.getenv ("USERPROFILE") , "Desktop")

データセットが保存されているパスを指定しています。今回はデスクトップのパスを“desktop_path”に格納しています。file.path () の部分を変更することで、任意のパスに設定できます。

ファイル名結合

otu_file <- file.path (desktop_path, "otu_table.tsv")
tax_file <- file.path (desktop_path, "tax_table.tsv")
sample_file <- file.path (desktop_path, "sample_data.tsv")

“desktop_path“に格納されているパスと、各データセット ("otu_table.tsv"、"tax_table.tsv"、"sample_data.tsv") を結合しています。

5.4 データセットの読み込み

BIOM形式など、phyloseqのimport関数 (import_biom () など) に対応したフォーマットのファイルが手元にある場合は、簡単にphyloseqオブジェクトへ変換できますが、今回は、そのような形式ではないリード数カウントデータ (tsv) を手動で読み込んでphyloseqを作成するスクリプトを解説します。

  • ファイルの読み込み

# ファイル読み込み (CSVの場合は read_csv に変更)
otu_df <- read_tsv (otu_file)
tax_df <- read_tsv (tax_file)
sample_df <- read_tsv (sample_file)

read_tsv () 関数を使って、それぞれのファイルを読み込み、OTUテーブル・分類情報・サンプル情報をそれぞれのデータフレームに格納しています。

行名の設定

otu_df <- column_to_rownames (otu_df, var = colnames (otu_df) [1])
tax_df <- column_to_rownames (tax_df, var = colnames (tax_df) [1])
sample_df <- column_to_rownames (sample_df, var = colnames (sample_df) [1])

それぞれのデータフレームの1列目を行名 (rownames) に変換しています。

phyloseq は、行名でサンプルや OTU を照合するため、この処理を行っています。

phyloseqオブジェクトの作成

otu_table_obj <- otu_table (as.matrix (otu_df) , taxa_are_rows = TRUE)

        tax_table_obj <- tax_table (as.matrix (tax_df) )

        sample_data_obj <- sample_data (sample_df)

        ps <- phyloseq (otu_table_obj, tax_table_obj, sample_data_obj)

otu_table、tax_tableのデータフレームを行列 (matrix) に変換し、sample_data をphyloseq の形式にしています。

その後、phyloseq () 関数でOTU、taxonomy、sample dataが結合された1つのデータセットを作成しています。

その後の希薄化、多様性指数算出のスクリプトは、公開データセットを使用する場合と同じです。

まとめ

腸内細菌叢の基礎知識

腸内には数百〜数千種類の細菌が存在し、善玉菌・悪玉菌・日和見菌に分類されます。重要なのは特定の菌を増やすことではなく、多様な菌がバランス良く共存することです。理想的なバランスは善玉菌2: 悪玉菌1: 日和見菌7とされ、この多様性が免疫機能、代謝、メンタルヘルスなどと関連しています。

多様性の評価方法

種多様性は「種の豊富さ」と「均等度」の2要素から構成され、以下の指標で定量評価します:

  • Observed features: 単純な種数のカウント

  • Chao1指数: 希少種を考慮した真の種数推定

  • Simpson指数: 各種の存在量を考慮し、均等性を評価

  • Shannon指数: 種数と均等性の両方を総合的に評価

レアファクション (希薄化) の重要性

サンプル間でリード数が異なると正確な比較ができないため、全サンプルを同じリード数に揃える処理が必要です。これにより「観察量の差」ではなく「真の多様性の差」を評価できます。

R言語での実践的解析

phyloseqとmicrobiomeパッケージを使用することで、OTUテーブル・分類情報・サンプル情報を統合管理し、多様性指数の算出からExcel出力まで一貫した解析が可能です。公開データセットでも手元のデータでも対応できる柔軟なスクリプトが提供されています。

このアプローチにより、腸内細菌叢の健康状態を科学的に評価し、機能性表示食品などの効果検証にも応用できます。

関連するサービス

参考文献

  • 松本哲哉 (編). 最新臨床検査学講座 臨床微生物学. 東京: 医歯薬出版
  • 頭山昌郁, 中越信和. 種多様性の評価における二、三の問題点. 環境動物昆虫学会誌 (環動昆). 2004; 15 (1), 31–48. DOI: https://doi.org/10.11257/jjeez.15.31
  • 大垣俊一. 多様度と類似度、分類学的新指標. Argonauta. 2008; 15, 10–22.
  • McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013 Apr 22;8(4):e61217. doi: 10.1371/journal.pone.0061217. PMID: 23630581; PMCID: PMC3632530.
  • Lahti L, Shetty S, Blake T, Salojarvi J. Tools for microbiome analysis in R. Microbiome package version 1.16.0 [Internet]. 2017 [cited 2025 Dec 1]. Available from: https://microbiome.github.io/tutorials/