メタゲノムデータを⽤いた
群集解析法

⾨脇 浩明
京都⼤学⼤学院 ⼈間・環境学研究科
⽇本学術振興会特別研究員

第63回⽇本⽣態学会仙台⼤会 ⾃由集会W31（企画：⽥辺晶史）

メタゲノム解析を通じ⽣態学に対しどのような貢献ができるのか

1. メタゲノム解析によって得られるデータの特徴
2. データの取り扱い⽅法
3. データの統計解析法とその例

メタゲノム解析によって得られるデータ
群集⾏列データのExcel表 (Claident)

⾏：サンプル

列：種、もしくはOTU (operational taxonomic unit)
ID

OTU
1OTU
2OTU
3OTU
4・・・

OTU
99

121

132

056・・・

220

1601・・・

56

3056

1044・・・

0459

71028・・・

05001990・・・

24

602872・・・

75

733

100・・・

0> comm <- read.table(“/….txt”, header=T)

次世代シーケンサーの配列データ→配列をOTUごとにまとめる→サンプルごとのOTUのリード数を集計

メタゲノム解析によって得られるデータ
分類学的同定結果のExcel表 (Claident)
query

Kingdom

Phylum

Class

Order

Family

Genus

Sp1

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Amanitaceae

Amanita

Sp2

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Sp3

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Sp4

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Sp5

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Amanitaceae

Amanita

Amanita fuliginea

Sp6

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Entolomataceae

Sp7

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Sp8

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Sp9

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Tricholomataceae

Laccaria

Laccaria bicolor

Sp10

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Cortinariaceae

Cortinarius

Sp11

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Sp12

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Amanitaceae

Amanita

Sp13

Fungi

Basidiomycota

Agaricomycetes

Agaricales

Cortinariaceae

Cortinarius> idlist <- read.delim(“/….txt”)

Speciesデータの取り扱い⽅法
分類学的同定結果（Claident）

サンプルの環境変数(explanatory variables)

データ統合

群集⾏列

(species by sample matrix)

マスターデータ> taxon <- idlist[match(colnames(comm),idlist$query),]> masterdata<-cbind(t(comm),taxon)

データの統計解析法
1. 希釈法 (rarefaction)
2. 探索的データ解析(exploratory analysis)

3. 仮説検定(hypothesis testing)
これら３つのステップを通じメタゲノムデータを⽣態学に活かす⽅法について考える
琵琶湖のプランクトン群集

希釈法（rarefaction）とはマーブルチョコのアナロジー
l チョコ1粒＝１個体の⽣物
l チョコの⾊ごとに異なる種
チョコをn個選んだ場合何⾊（何種）が観察されるか？
www.lovechoco.org

リード数 ＝ シーケンサーの探索努⼒（例：観察個体数など）
リード数の異なるサンプルの種多様性は単純に⽐較できない
より公平な⽐較を⾏うためには、リード数を適切に希釈する必要がある

希釈法（rarefaction）
異なる探索努⼒（リード数）のもと得られたデータを⽐較するためには？

Rarefaction curve
14

サンプルサイズnのもとで期待される種数

12Species

OTU数1086420

N = total number of individuals
Ni = number of individuals in the i th speciesn = size of the smaller sample (rarefied)

50

100150200Sample Size

リード数
Hurlbert (1971)Ecology; Gotelli and Graves (1996)

Package Vegan: rarefy, rrarefy, rarecurve希釈法（rarefaction）
Rarefaction curve

最も素朴な⽅法：

⼤量のデータが無駄となる可能性

10Species
OTU数

種数の増加が頭打ちになったサンプルのみを解析に⽤いる（⼗分なリード数が得られなかったサンプルを破棄）

代替案：
カバレッジで揃え、データを希釈する

8642050

100 150

リード数Sample
Size

Hurlbert (1971)Ecology; Gotelli and Graves (1996)

Package Vegan: rarefy, rrarefy, rareslope, rarecurve希釈法（rarefaction）
Rarefaction curve

同じ傾きとなるリード数において希釈例） カバレッジ98%(傾き0.02)で揃える
𝑁 − 𝑁𝑖0𝑛12
J = E(S) = ' (1 − +𝑁, 0𝑖=1𝑛𝑆10ガンマ関数の公式を活⽤nについて微分Species
OTU数,

86420

種数期待値の導関数（傾き）50

100 150

リード数
SampleSize𝑑Γ ′ (𝑁 − 𝑛 + 1)
(J) = '𝑑𝑛Γ(𝑁 − 𝑛 + 1)
Γ ′ (𝑁 − N𝑖 − 𝑛 + 1) {log ⁡(Γ⁡(N+1))}
0 𝑒 (Γ(𝑁−𝑁𝑖+1)+log (Γ (𝑁−𝑛+1))−log (Γ (𝑁−𝑁𝑖 −𝑛+1))−log ⁡−()
Γ𝑁−N−𝑛+1
Package Vegan: rarefy, rrarefy,rareslope, rarecurve𝑖

希釈法（rarefaction）
指定されたリード数の群集⾏列データを希釈する（rrarefy関数を⽤いる）→

⼆値データ（在／不在データ）に変換する

IDSp 1

Sp 2

Sp 3

Rarefaction curveの描画

Sp4・・・

Sp99

1

1

101・・・

020000・・・

1

30111・・・

04

1101・・・

050010・・・

060010・・・

17

1000・・・

0以上をもって群集解析の準備が完了
Package Vegan: rarefy, rrarefy, rareslope, rarecurve

Number of species

⽣成された⼆値データの例

Number of reads

探索的データ解析ー種構成・分類群構成の図
Genusレベルでの群集組成

Orderレベルでの群集組成

Claidentの分類情報を活⽤した様々なグラフィックスが利⽤可能
樹⽊と根圏真菌群集の相互作⽤⾏列
EcM菌根菌樹種

AM菌根菌

Package Graphics: piechart ; BiparMte : visweb

腐生菌

その他

α、β、γ多様性の定量化
サンプル×種の群集⾏列

群集の⾮類似度⾏列

Jaccardの⾮類似度指数
l β多様性の計算⼿法は多様である
l β多様性に関する論⽂の全てが役⽴つわけではない

Anderson et al. (2011) Ecol LeP> vegdist(..., method="bray", binary=FALSE)

仮説検定 (hypothesis testing)l

サルの⾷性幅とその変異性の解析

l

⿂類DNAの時空間出現パターン（メタ個体群動態）の解析

l

微⽣物種間の排他的分布パターンの解析

l

真菌群集が寄主植物の成⻑率に与える影響の解析

サルの⾷性幅とその変異性の解析（PERMANOVAとPERMDISP）
サルの採餌⾏動は集団間でどのように異なるか？

#1

サル糞のDNA分析により、サル集団の餌組成について、３つの集団間で⽐較する場合複数のβ多様性指数で解析する利点：
u 餌種数＋餌組成の違い → Jaccard⾮類似度指数#2u 餌組成の違いのみ → Raup-Crick⾮類似度指数
Jaccardの⾮類似度指数#3サルの⾷性幅とその変異性の解析（PERMANOVAとPERMDISP）
集団１

集団２

Total sum of squares (SS)集団1

非類似度行列集団2

Within-group SS

全サンプル数 N
グループ内のデータ数 n
サンプルiとjの間のβ多様性 dij
サンプルiとjが同グループか否か εij

Among-group SS

Anderson MJ (2001) Austral Ecology

サルの⾷性幅とその変異性の解析（PERMANOVAとPERMDISP）
分散分析表
DfGroupResidualTotal

21416SS
1.13
2.64
3.78

MS
0.57
0.19

F
3.01

R20.300.70
1.00

P0.005

グループのラベルをシャッフルしてF値を繰り返し計算→ 実測値と⽐較し、並べ替え検定
Anderson MJ (2001) Austral Ecology

Package Vegan: adonis

PERMANOVAの結果解釈における注意点 :平均 vs 分散の効果(i)処理区間で平均的な群集構造が異なる場合

2nd MDS

2nd MDS(ii) 処理区ごとに群集構造の分散が異なる場合

データ点

1st MDS

1st MDS
シナリオ(i)

シナリオ(ii)

PERMANOVA

有意有意PERMDISP

有意でない有意Package Vegan: adonis, betadisper, permutest, TukeyHSD

Anderson MJ (2001) Austral Ecology

⿂類のメタ個体群動態パターンの解析
環境DNAの出現パターン（局所移⼊・絶滅率に影響する要因を特定する⽅法
T=1

T=2

京都鴨川のヨシノボリ

・・・・・・

オッズ⽐

Plot1
Plot2

⽣息地の⽔中浮遊DNAから⿂種特定を⾏う

在 (T2=1)

不在 (T2=0)

在 (T1=1)

ab

不在 (T1 =0)

cd← 各々のOTUについてこの表を得る
同時にT=1, 2における環境要因を測定

McCullagh and Nelder (1989) ; Yee & Dirnböck (2009) Journal of Theoretical Biology

Package VGAM : vglm

⿂類のメタ個体群動態パターンの解析
オッズ⽐ ＝ 局所的移⼊率・絶滅率の形を特徴づける変数
局所的移⼊率

局所的絶滅率

周辺確率

周辺確率

0→1になる確率*交換可能性条件が満たされる場合1→０になる確率

Yee & Dirnböck (2009)> vglm(..., family = binom2.or("cloglog", exchangeable = TRUE))

⿂類のメタ個体群動態パターンの解析推定されたメタ個体群動態と種のニッチ

⼆変量回帰モデル
𝑐𝑙𝑜𝑔𝑙𝑜𝑔 P(𝑌1 = 1|𝑥 ) = 𝜂1
𝑐𝑙𝑜𝑔𝑙𝑜𝑔P(𝑌2 = 1|𝑥 ) = 𝜂2

局所移⼊率

局所絶滅率

←ポアソン分布を仮定

𝑙𝑜𝑔 𝜓(𝑥) = 𝜂3個体数例）
リンク関数を式ごとに変える
回帰項を式ごとに変える (平⾏性の仮定)
縮⼩ランク回帰やゼロ過多モデル
ベクトル⼀般化線形モデル、など(補対数対数変換)

ηをそれぞれ柔軟にモデリングすることで⽣物のニッチと環境変化に伴うメタ個体群過程を評価

パッチ占有率環境傾度
Yee & Dirnböck (2009)> vglm(..., family = binom2.or("cloglog", exchangeable = TRUE))

オッズ⽐に関する注意点：Simpsonʻs Paradox
重要な未知の要因 Z が絶滅・移⼊過程に関与する場合(i) 交絡因⼦Zを考慮した場合

在 T1不在 T1(ii) 交絡因⼦Zを考慮しない場合

Z+
Z在T2 不在T2 在T2 不在T280
165045216016010452
ψ= 5

ψ= 5

在 T1不在T1在T2 不在T2130468
170612
ψ= 1

⺟集団全体で⾒たときのオッズ⽐≠分割して⾒たときのオッズ⽐

未知の要因Zの例：
⽣息地分断化、隠蔽種(OTU)の存在など
柳川 (1985) 離散多変量データの解析微⽣物種間の排他的分布の検出⽅法
OTU

Cスコア単位

⽊材腐朽菌群集の競争が⾒られる倒⽊

1 2 3 4 5

プロット番号

1 2 3

各種ペアの組み合わせにおいて何個のCスコア単位が観察されたか

種iとjのCスコア
群集のCスコア
化学合成細菌群集の競争が⾒られる淀み

Package biparMte: C.scoreri 種iの出現回数
Sij 種iとjの同時出現回数
M 全体の種の数

Stone and Roberts (1990) Oecologia排他的分布の検出ー 帰無モデルの選定
帰無モデルをカスタマイズすることが可能

シナリオ(ii)

シナリオ(iii)

シナリオ(iv)

1000100 … 01
0110101 … 111000100 … 01
0110101 … 111000100 … 01
0110101 … 111000100 … 01
0110101 … 11

列合計

列合計

合計在数を固定し0/1を⽣成

⾏合計

シナリオ(i)

⾏合計を固定し無作為に0/1を⽣成

列合計を固定し無作為に0/1を⽣成

⾏合計→ Monte Carlo simulation の⽅法を変更

⾏と列の合計両⽅を固定し無作為に0/1を⽣成
Goteli and Graves (1993) Null model ecology> oecosimu(comm, method ="swap", nsimul = 9999, ...)

排他的分布の検出ー 帰無モデル検定
帰無モデルをカスタマイズすることが可能→ Monte Carlo simulation の⽅法を変更
Monte Carlo法により⽣成されたCスコアの頻度分布（9999回の結果）

1000100 … 01
0110101 … 11

⾏合計

データセット

観測された
Cスコアの値

列合計> oecosimu(comm, method ="swap", nsimul = 9999, ...)

補⾜：Cスコア検定における2つの過誤Type I error： 正しい帰無仮説を誤って棄却する確率(α)
Type II error ： 誤っている帰無仮説を棄却しない確率(β)
検出⼒＝１－β ：誤っている帰無仮説を棄却する確率
帰無仮説

代替仮説

Van Belle (2008) StaMsMcal rules of thumb> oecosimu(comm, method ="swap", nsimul = 9999, ...)

排他的分布の検出ー 多重検定とp値補正
OTUペアごとにCスコア解析 → Type I errorの増⼤
u Familywise error rate (FWER) → Bonferroni補正
少なくとも1つの帰無仮説を誤って棄却する確率を制御
u False discovery rate (FDR) → Benjamini-Hochberg補正, q値法
多くの帰無仮説群のうち、誤って棄却される帰無仮説の割合を制御

頻度P値
⼀様分布

← ランダムに⽣成した変数ペア(~N(0,1))
を⽤いてｔ検定を繰り返し、その結果得られたp値の分布

Benjamini and Hochberg (1995) J R Stat Soc; Storey and Tibshirani (2003) PNAS

> p.adjust(p, method = “BH”)#see also Package: qvalue⾼次元データの解析
⾼次元データ解析の課題

Above-belowground interaction networkl ⾼い予測精度を保証すること→将来予測が可能となるl 現象と真に関連する変数を選択すること→推定量の漸近正規性が保証されず検定不可能になることを避ける
安道 (2014) ⾼次元データ分析の⽅法

⼟壌に張り巡らされた菌⽷のネットワーク

補⾜： 尤度関数と確率関数の違い
尤度関数：
あらゆるpにおいて特定のxの値が観察される確率 (xを固定)

確率関数：
特定のpにおいて観察されるあらゆる異なるxの値がとりうる確率(pを固定)←⼆項分布の例
𝐄(𝑋) = 𝑛𝑝
𝐕𝐚𝐫(𝑋) = 𝑛𝑝(1 − 𝑝)L(p;125)=190C125 p125(1－p)65

lasso推定による⾼次元データの解析
least absolute shrinkage and selecMon operator

植物の葉のC：N⽐に影響を与える根圏微⽣物OTUを選択

1y＝
葉のC:N⽐(元素分析)

𝑙(𝛽 ) =𝑛𝑝

𝛼=1𝑗 =11( log𝑓 (𝑦𝛼 |𝑥𝛼 , 𝛽 ) − 𝜆 ( |𝛽𝑗 |𝑛尤度関数

罰則項

正則化パラメータ

X=根圏微⽣物OTUデータ(メタゲノム)

罰則付き尤度関数を最⼤とするβを⾒つけることで、説明変数選択と回帰係数βの推定を同時に実⾏

アカマツの実⽣
Tibshirani 1996 J Royal Stat Soc; 安道 (2014) ⾼次元データ分析の⽅法

Package lars: larslasso推定による⾼次元データの解析
least absolute shrinkage and selecMon operator

植物の葉のC：N⽐に影響を与える根圏微⽣物OTUを選択
⼀般線形モデルにより植物に影響を与える真菌OTU候補を絞り込む
例. ラシャタケ類(Tomentella)が寄主植物の葉のC:N⽐を上昇させる可能性を⽰唆←
Step3でCp最⼩値Cpを最⼩化する正則化パラメータを適切な値として解釈する
𝑛1
2
𝐶𝑝 (𝜆) = ){𝑦𝑘 − 𝜇̂ 𝛼 (𝜆)}2 + 𝑝̂ 𝜎 2
𝑛𝑛𝛼=1

Package lars: lars

メタゲノム⽣態学の展開
l メタゲノム・環境DNAに対する⾼い関⼼とその応⽤可能性
l メタゲノム解析は、これまでの⽣態学を置き換えるものではなく、フィールドワークに基づく知⾒のさらなる強化を可能とする
l ⽣態学者が初めて直⾯するデータ解析の問題（⾼次元・スパースデータ）
l サンプリングの時空間スケール、反復の採り⽅、ランダム化と層別化などに関する実験⽣態学的な知⾒と結びけることでさらなる⼒を発揮するメタゲノム⽣態学の展開
l 今後は、⽣態学的なプロセス（相互作⽤や共存機構）の解明、および、その予測を視野に⼊れた、時系列データの蓄積が最も重要となる
l 環境DNAを⽤いたレジリエンス(resilience)や時間的安定性(stability)の定量化、早期警告シグナルの確⽴などが求められる