核酸配列に基づかない微生物群集の分析方法

核酸を使わない微生物群集の分析方法

核酸 (DNA や RNA) を利用した微生物の分析方法は、環境サンプル中の微生物を群集として扱う場合には間違いなく最もポピュラーな方法です (Knight et al. 2018 Nature Review Microbiology など)。また、環境サンプルの分析だけでなく微生物の研究におけるあらゆる場面で核酸の分析が必要な技術であることは疑いようがありません。

一方で、微生物群集の分析方法は核酸に基づいたもの以外にも数多くあります。それらは、核酸に基づいた分析ほど頻繁に利用されるわけではありませんが、それぞれ核酸分析では得られない情報を提供してくれます。使用する機会があるかどうかは別として、どのような方法があり、どのような情報を得ることができるかを知っておくことは有益です。

そこで今回のポストでは (1) 顕微鏡観察、(2) リン脂質脂肪酸、(3) キノン、(4) BIOLOG、(5) 分解酵素活性、について紹介したいと思います (単離培養 については超重要ですがここでは触れません)。

1. 顕微鏡観察

顕微鏡観察は最も古典的な微生物の分析方法の一つと言ってもいいかもしれません。核酸分析にない一番の特徴としては 形態情報が得られる ことでしょうか。どのような形をしているのか、どのくらいの大きさなのか、といった情報はその微生物の生態に関する示唆を与えてくれるだけでなく、その微生物が保持する炭素や窒素など元素の量の推定も可能にします (微生物が保持する元素の量は地球上の炭素循環や窒素循環を考える上で非常に重要です)。

顕微鏡には非常に多くの種類があり、自分自身も素人に近いですが、微生物の分野で使用されているいくつかの顕微鏡関連の手法を紹介します。

DAPI による (微生物種を区別しない) 核酸染色と細胞数カウント

DAPI という色素で微生物細胞中の核酸 (主に DNA) を青く染めて蛍光顕微鏡下で細胞数を数えたりします。簡便な方法で、煩雑な操作を必要とせず細胞中の DNA を染色できます。ただし、微生物種の区別はできません。また、微生物以外のDNA も光ります。

FISH (Fluorescent in situ hybridization) (Amann et al. 1990 Journal of Bacteriology) により微生物種を染め分けて細胞数カウント

FISH では微生物の 16S rRNA と相補的なプローブ (短い DNA 配列) を用意します。このプローブには蛍光色素がくっついており、染めたい微生物種 (分類群) の 16S rRNA に特異的なプローブをデザインすることで特定の微生物のみを染めて観察することができます (核酸を利用しているとも言えますが、配列を解読しているわけではありません)。

FISH には様々な改良法があります。例えば、手順は増えますが Horshradish peroxidase という酵素を利用して、蛍光を増幅することで観察を容易にする CARD-FISH という方法もあります (Schönhuber et al. 1998 Applied and Environmental Microbiology; Pernthelar et al. 2002 Applied and Environmental Microbiology) (めんどくさくて発狂しそうになる、という意見も)。

自分は CARD-FISH を土壌サンプルに対して適用してみました (Ushio et al. 2013b Soil Biology & Biochemistry)。下の写真はスウェーデンのツンドラ土壌のサンプルを CARD-FISH 法で染色したものです。青色に光っているのが DAPI により染色された細菌で、黄緑色に光っているのが古細菌特異的なプローブで染めた古細菌です (２重染色しています)。微生物周辺にあるもやもやしたものは土壌有機物の集合体と思われます。

電子顕微鏡

電子顕微鏡ももちろん使用されていますが、こちらは自分はほとんど使用経験がありません。

2. リン脂質脂肪酸 (PLFA)

リン脂質脂肪酸 (Phospholipid Fatty Acid; PLFA) は微生物の細胞膜に含まれる成分です。下記は Wikipedia の ｢リン脂質｣ のページに掲載されている図です。

リン酸とグリセリンの左側に脂肪酸 (オレイン酸とかパルミチン酸) が結合しています。面白いことに、この脂肪酸部分は微生物の分類群ごとに “ある程度” その種類が決まっています。従って、脂肪酸の種類をガスクロマトグラフィーで分析することでサンプル中の微生物群集の組成を推定することが可能です。この手法は Ushio et al. (2008, 2010a, 2010b, 2013a, 2013c) などでお世話になりました。

脂肪酸の記載方法

脂肪酸、およびその構造の記載方法は生態学ではあまり頻繁に登場するものではありません。脂肪酸の構造の記載方法について少しだけ解説しておきます。脂肪酸は上に示したように長い炭素鎖が基本にあり、それに加えて二重結合があったり、炭素鎖に分岐があったりします。例えば、以下のパルミトレイン酸という脂肪酸は炭素数が 16 で途中に二重結合が一つある構造をしています (Wikipedia より引用しています)。

この場合、16:1ω7 などと記載されます。“16” が炭素数を示し、":1" は二重結合が 1 つあることを示します。“ω7” の部分は炭素鎖の末端 (-OH がない方) から数えて 7 番目に二重結合があることを示しています。また、炭素数を表す “16” の前に “a”, “i”, “cy” などのアルファベットが付くことがありますが、これは “antiiso-”, “iso-”, “cyclo” などの立体構造に関する情報を表します。

分類群特異的脂肪酸

以下、Ushio et al. (2013c) Plant and Soil で使用した代表的なバイオマーカーを示しておきます。

• 真菌: 18:2ω6,9
• グラム陽性菌: i14:0, i15:0, a15:0, i16:0, i17:0, a17:0
• グラム陰性菌: 16:1ω7, cy17:0, cy19:0

また、自分は指標として使用しませんでしたが、微生物のストレス状態によっても脂肪酸の組成が変わるらしいです。つまり、脂肪酸の組成は細胞膜の流動性に関わるようで、例えば脂肪酸に二重結合が増えると (ω がついた脂肪酸が増えると)、細胞膜の流動性が上がりストレス環境下 (例えば、低温状態) においてより適応的になると考えられます (不勉強で少し自信がありません…)。

利点と不利点

PLFA を用いた微生物群集の組成分析は 核酸に基づいた分析に比べて以下のような利点･不利点があります。

• 利点
  – 分析結果が定量的。例えば、PLFA は量として XX µmol/g soil などとして得られるため、異なるサンプルを比較して、どちらのほうが真菌が多い、少ない、といった議論が可能です (真菌:細菌比の算出も一度の分析でできます)。
  – 核酸のように煩雑な事後解析が不要。核酸の配列解析にはいわゆるバイオインフォマティクスのような知識が多少なりとも必要ですが、PLFA の解析ではそこまで煩雑な前処理は必要ありません。 ガスクロマトグラフィーのピークの位置とピークの面積から脂肪酸の種類と量を推定します。
  – 細胞膜成分は、微生物が死ぬと比較的迅速に分解されると考えており、PLFA の分析結果は生きている微生物の量を反映していると考えられています。この点、時に死んだ微生物の DNA の存在が問題となる核酸に基づいた分析とは対照的です (例えば、Carini et al. 2016 Nature Microbiology)。
  – 同位体トレーサーなどと組み合わせることで物質の動態を追える。例えば、適当な基質を同位体トレーサーでラベリングすることで、その基質がどのような分類群の微生物に取り込まれたか調べることができます (これは核酸分析でも可能な場合がありますが)。

• 不利点
  – 分類群の解像度が低い。種レベル･属レベルといった解析は不可能です。DNA に基づいた分析と比べると決定的な不利点、という感じがします。
  – 化学的な分析が必要で有機溶媒をたくさん使用するため、試薬のハンドリングや廃棄に気をつける必要があります。
  – 使っている人が核酸分析に比べると少ない。従って、関連情報も得にくいです。PLFA を使った研究を発表すると ｢なぜ DNA を使わないの？｣ と聞かれたりします。

3. キノン

キノンは呼吸鎖における電子伝達物質ですが、PLFA と同じように、微生物分類群ごとにその化学種が大まかに決まっています。ユビキノン･メナキノンなどの種類があるようですが、自分には使用経験がなく文献などもあまり読み込んでいないため、ここでは解説は控えることにします。定量性があることや分類群の解像度がそこまで高くないことなど、PLFA と似たような特徴を持っているようです。

知り合いが使用して湖の微生物の研究をしていました。例えば、Takasu et al. (2013) Limnology などをご覧ください。

4. 分解酵素活性

PLFA やキノン、また DNA も環境中の微生物群集の組成に関する情報を提供します。微生物群集の組成は、｢その環境中にどのような微生物が生育しているか｣ という情報を提供してくれますが、環境中での物質の流れを考えたい場合などはその 機能 の評価も重要になります。核酸に基づいた方法であれば、環境中の微生物群集の機能遺伝子や RNA 発現パターンを調べることで微生物群集の機能に関する情報が得られます。しかし、核酸に基づいた方法だけでは、実際に微生物群集が物質を代謝しているのか どうかは分かりません。そういった場合に、分解酵素活性の分析や次で紹介する BIOLOG が有効になってきます。

分解酵素活性の分析では、何らかの発色物質と結合した基質を利用します。代表的なものは次の二つです。

p-ニトロフェノール 結合基質

p-ニトロフェノールはニトロ基を有する単純なフェノール化合物です。p-ニトロフェノールに結合したリン酸 (p-ニトロフェニルリン酸) などが試薬会社から市販されており、このような化合物を分析に用います。例えば、p-ニトロフェニルリン酸の溶液に微生物を含む土壌抽出液を少量添加すると、p-ニトロフェニルリン酸が代謝された場合は p-ニトロフェノールとリン酸が分離し、(pH 調整後に) 溶液が黄色く発色します。この発色を吸光度計で測定することで分解酵素活性を測定します。p-ニトロフェニルリン酸を分解する酵素はリン酸分解酵素などと呼ばれます。p-ニトロフェノールが結合する基質を選ぶことでいろいろな分解酵素活性の測定が可能です。この手法は Ushio et al. (2010a, 2010b, 2013c) などでお世話になりました。

Methylumbelliferone (MUB) 結合基質

測定原理はほぼ p-ニトロフェノール結合基質と同じです。MUB と結合した基質が微生物群集の分解酵素で分解されると MUB が遊離し、その量を測定することで分解酵素活性を測定します。ただし、MUB は蛍光物質であり、測定には吸光度計ではなく、蛍光プレートリーダーが必要です。p-ニトロフェノールと同じく、MUB に結合している基質を選ぶことで様々な分解酵素活性を測定することができます。Ushio et al. (2013) PLoS ONE で使用しました。

p-ニトロフェノール結合基質の方が試薬も必要な実験機器も安価ですが、大量のサンプルを処理する場合には 96ウェルマイクロプレートを利用しやすい MUB 結合基質の方が向いているように思います。

5. BIOLOG

BIOLOG は分解酵素活性分析と同様に、微生物群集の基質代謝能力を測ることの出来る手法です (Garland & Mills 1991 など)。実験に使用するマイクロプレートは BIOLOG 社が販売していて、“BIOLOG” というワードは会社の名前でも有り、手法の名前でもあり、商品の名前でもある、という感じです (多分)。

BIOLOG プレートにはいろいろな種類があります (EcoPlate など)。96ウェルマイクロプレートに様々な化合物が微生物の代謝基質として予め入れられており、そこに生きた微生物を含むサンプル (単離株でもいいし、水や土壌などの環境サンプルでもよい) を適当に希釈して投入して使用します。サンプル中の微生物 (群集) にウェル内の基質の代謝能力があれば、ウェル内の基質が代謝され、紫色に発色します。プレートリーダーでその強度を読み取ることで、微生物 (群集) の基質代謝能力を評価できます (Miki et al. 2018 Ecological Research などをご覧ください)。

(自分も使用したことはあるのですが、論文になるデータが取れませんでした…)

終わりに

冒頭に述べたように微生物生態学では核酸に基づいた分析が主流であることは間違いありません。その一方で、核酸分析に加えてここで解説したような手法を追加することで研究対象となる微生物群集をより深く分析することができます。それぞれの手法は核酸分析では得られない情報を多かれ少なかれ提供してくれます。様々な角度から微生物群集を分析することで、その群集組成･機能･生態系での役割を明らかにすることができるはずです。

参考文献

• Amann R, Krumholz L, Stahl DA. (1990) Fluorescent-oligonucleotide probing of whole cells for determinative, phylogenetic, and environmental studies in microbiology. Journal of Bacteriology 172: 762-770
• Carini P, Marsden PJ, Leff JW, Morgan EE, Strickland MS, Fierer N (2016) Relic DNA is abundant in soil and obscures estimates of soil microbial diversity. Nature Microbiology 2: 16242
• Garland JL, Mills AL. (1991) Classification and characterization of heterotrophic microbial communities on the basis of patterns of community-level sole-carbon-source utilization. Applied and Environmental Microbiology 57: 2351-2359
• Knight R, et al. (2018) Best practice for analysing microbiomes. Nature Review Microbiology 16: 410-422
• Miki T, Yokokawa T, Ke P-J, Hsieh I-F, Hsieh C-h, Kume T, Yoneya K, Matsui K (2018) Statistical recipe for quantifying microbial functional diversity from EcoPlate metabolic profiling, Ecological Research 33: 249-260
• Pernthaler A, Pernthaler J, Amann R. (2002) Fluorescence in situ hybridization and catalyzed reporter deposition for the identification of marine Bacteria. Applied and Environmental Microbiology 68: 3094-3101
• Schönhuber W, Fuchs B, Juretschko S, Amann R. (1997) Improved sensitivity of whole-cell hybridization by the combination of horseradish peroxidase-labeled oligonucleotides and tyramide signal amplification. Applied and Environmental Microbiology 63: 3268-3273
• Takasu H, Kunihiro T, Nakano S (2013) Estimation of carbon biomass and community structure of planktonic bacteria in Lake Biwa using respiratory quinone analysis. Limnology 14: 274-256
• Ushio M, Wagai R, Balser TC, Kitayama K. (2008) Variations in the soil microbial community composition of a tropical montane forest ecosystem: Does tree species matter? Soil Biology & Biochemistry 40: 2699-2702
• Ushio M, Kitayama K, Balser TC. (2010a) Tree species-mediated spatial patchiness of the composition of microbial community and physicochemical properties in the topsoils of a tropical montane forest.Soil Biology & Biochemistry 42: 1588-1595
• Ushio M, Kitayama K, Balser TC. (2010b) Tree species effects on soil enzyme activities through effects on soil physicochemical and microbial properties in a tropical montane forest on Mt. Kinabalu, Borneo. Pedobiologia 53: 227-233
• Ushio M, Miki T, Balser TC. (2013a) A coexisting fungal-bacterial community stabilizes soil decomposition activity in a microcosm experiment. PLoS ONE 8: e80320
• Ushio M, Makoto K, Klaminder J, Nakano S-I. (2013b) CARD-FISH analysis of prokaryotic community composition and abundance along small-scale vegetation gradients in a dry arctic tundra ecosystem. Soil Biology & Biochemistry 64: 147-154
• Ushio M, Balser TC, Kitayama K. (2013c) Effects of condensed tannins in conifer leaves on the composition and activity of soil microbial community in a tropical montane forest. Plant and Soil 365: 157-170

Written with StackEdit.

rEDM を用いた Empirical Dynamic Modeling: (5) Convergent Cross Mapping (CCM)

Convergent Cross Mapping (CCM)

ここでは変数間の因果関係を検出する Convergent Cross Mapping (CCM; Sugihara et al. 2012) について解説します。

CCM で使用されるアルゴリズムは基本的に Simplex projection に非常に似ており、従って Simplex projection を理解できていればアルゴリズムを理解することは比較的容易です。

1. CCM の概念的な説明

｢因果がある｣ とは

｢相関は必ずしも因果を意味しない｣ というのはあまりにも有名です (例えば、Wikipedia による解説も御覧ください)。ここではそれについては改めて詳しく説明することはしませんが、それではそもそも ｢因果がある｣ とはどういうことでしょうか？この疑問について深く考えると非常に難しく、哲学や宗教も絡んできそうですが、ここではひとまずそのようなことは考えず ｢同じ力学系に属していれば因果がある｣ と思うことにしましょう。例えば ｢時系列データを用いた動態の再構成｣ で用いたローレンツアトラクターは以下の微分方程式から現れる動態でした。

$$\frac{dx}{dt} = -pX + pY\\ \frac{dy}{dt} = -XZ + rX - Y\\ \frac{dz}{dt} = XY - bZ$$

この系は一つの力学系ですが、$X$, $Y$, $Z$ が同じ力学系に属して (お互いに因果を持ちながら) 動態を形成しています。このような状況のときには、例えば $X$ の時間変化に $Y$ が影響を与えていますし、直感的な意味で $Y$ は $X$ に因果がある、影響があると考えて良さそうです。もっと簡単に考えると、$Y$ の値が変化するときに、$X$ も変化する、ということができます。

同じ力学系に属しているかどうかを調べる

さて、CCM は２つの変数が同じ力学系に属しているかどうかを調べようとします (他の解釈も可能ですが、図的に説明しやすいのでここではこの解釈を説明します)。では、どうやって同じ力学系に属しているかどうか調べるのでしょうか？

下の図を覚えているでしょうか？左は $X(t)$ の時間遅れを用いて再構成した動態 ($M_x$ とします)、右は $Y(t)$ の時間遅れを用いて再構成した動態 ($M_y$) です。

以前説明したように、$M_x$ はもともとの $X$, $Y$, $Z$ を用いて再構成した動態 (便宜的に真の動態、と呼びます) と同じかたちとみなせます。同様に $M_y$ も真の動態と同じかたちとみなせます。従って、$M_x = [真の動態] = M_y$ となり、もし $X$ と $Y$ が同じ力学系に属していれば $M_x$ と $M_y$ は同じかたちになるはず です (下図)。

CCM では、Simplex projection を応用して、再構成された２つの動態が同じかたちかどうかを調べることで変数間の因果関係を検出しようとします。

注意: 最近、とある研究会で ｢動態の再構成についてローレンツアトラクタ−を例として説明するのはあまりよくないのではないか｣ とのご指摘をいただきました。その方はカオス時系列などの研究をされていた方でいくつかの根拠を持ってそのような指摘をしてくれました。この点について現在勉強中で、ひょっとすると将来説明の仕方が変わる可能性があります。

2. CCM のアルゴリズム

では、Simplex projection を応用して、どのようにして ｢同じかたちかどうか｣ をチェックするのか、R を利用して解説していきます。

マニュアルコードによる解説

例によって自作の２種系のデータセットを使用して説明します (四度目)。

# rEDM のロード
library(rEDM)
packageVersion("rEDM") # 2020年1月時点で v0.7.5

# 2種系のデータセットを作成
d0 <- as.data.frame(matrix(rep(NA,3*1000),ncol=3))
colnames(d0) <- c("time","x","y")
d0[,1] <- 1:1000
d0[1,2:3] <- c(0.1, 0.1)
for(i in 1:999){
        d0[i+1,2] <- d0[i,2]*(3.8 - 3.80*d0[i,2] - 0.02*d0[i,3])
        d0[i+1,3] <- d0[i,3]*(3.5 - 0.10*d0[i,2] - 3.50*d0[i,3])
}

# 最初の 99 点を捨てて、時間 100 から 1000 のデータを使用
d <- d0[100:1000,]

最適埋め込み次元の決定

まずは Simplex projection を利用して最適埋め込み次元 (E) を決定します。ここでは rEDM の simplex 関数を利用します。

simp_res <- simplex(y)
plot(simp_res$E, simp_res$rmse, type = "b", las = 1, xlab = "E", ylab = "RMSE")

RMSE に基づいて、最適埋め込み次元は 2 を使用することにします。

Ey <- simp_res$E[which.min(simp_res$rmse)]

最近傍点の特定

最適埋め込み次元を決定した後は Simplex projection と同様に $Y(t)$ の埋め込みを作り、ターゲットの点を決めて最近傍点を特定します。

# Y(t) の埋め込みを作成して、ターゲットの点を指定
y_embed <- make_block(y, max_lag = Ey)
t_target <- 240　# ターゲットの点を指定

# t = 240 の点から各点への距離を計算
distance <- sqrt(rowSums((y_embed[,2:3] - c(y_embed[t_target, 2:3]))^2))
nn_index <- order(distance)[2:4] # 自分自身を除いて 3 点近い点を選ぶ

# 最短距離を取得
min_distance <- distance[nn_index][1]
# 最短距離で補正した重みを計算
weights <- exp(-distance[nn_index]/min_distance)
# 合計の重みを計算
total_weight <- sum(weights)

Cross mapping

通常の Simplex projection ではこの後、以下のように $Y(t_{target})$ の近傍点の情報 (nn_index) から、$Y(t+1)$ を予測します。

# Y(t+1) を予測
pred <- (weights %*% as.matrix(y_embed[nn_index + 1,3:2])) / total_weight

一方、CCM では以下のようにして $Y(t_{target})$ の近傍点の情報 (nn_index) を元にして $X(t)$ の値を予測します。この

# Y(t) の近傍点の情報を利用して X(t) を予測
pred <- (weights %*% as.matrix(x[nn_index])) / total_weight

予測値を見てみましょう。かなり近い値を予測できていることが分かります。

> pred # 予測値
[1,] -1.21572
> x[t_target] # 実測値
[1] -1.254002

この操作は一体何をしているのでしょうか？これは、元の時系列に戻ってみると分かりやすいでしょう。下の図を使って説明します。

図の左上のパネルは $Y(t)$ のターゲット点周辺の様子です ($t_{target} = 240$ ですが、生成したデータの最初の 99 点を捨てているので実際の time は 339 となっています)。２次元で動態を再構成しているので、$Y(t)$ と $Y(t-1)$ の点の ｢動き｣ に注目していることになります。最近傍点の特定とは、ターゲットの点の ｢動き｣ と似た点を $Y(t)$ の時系列から探していることになります。これは図中では (1) Identify nearest neighbors に相当します。その後、特定した $Y(t)$ の時間と同じ時間の $X(t)$ を特定します (図中 (2) Cross-mapping)。$Y(t)$ の情報を元に特定された $X(t)$ の値を重み付け平均して $X(t)$ の値を予測しています (図中 (3) Prediction)。(図中の (2) と (3) を合わせて Cross-mapping と呼ぶべきかも知れませんが、ここでは説明の都合上２つを分けています)

Convergence のチェック

CCM では上記の手順を $X(t)$ の全てのデータ点で行い、$Y(t)$ の情報のみから $X(t)$ を予測できるかを調べます。さらに CCM では Simplex projection では通常行わない計算が加わります。それは、状態空間の再構成に利用するデータ点の数 (rEDM パッケージ内では library size と呼ばれます) を変化させ、データ点数の増加に伴って予測精度が上昇するかどうかをチェックします。この現象は Convergence と呼ばれていて、CCM の名前の由来にもなっています。

rEDM の ccm 関数を利用してこの Convergence をチェックしてみましょう。ccm 関数で求めた Cross-mapping による予測精度を ccm_means 関数で見やすい形にまとめています。

# Convergence のチェック
ccm_all_res <- ccm(cbind(y, x), E = Ey, lib_sizes = seq(Ey+1, nrow(d), by = 30))
ccm_res <- ccm_means(ccm_all_res)

# 結果の図示
plot(ccm_res$lib_size, ccm_res$rho, type = "l",
     xlab = "Library size", ylab = expression(rho), las = 1)
plot(ccm_res$lib_size, ccm_res$rmse, type = "l",
     xlab = "Library size", ylab = "RMSE", las = 1)

Library size の増加に伴って、相関係数 ($\rho$) の増加や RMSE の減少が見られます。これは $X$ から $Y$ への因果があることを示唆しています。

注意点: CCM では $Y$ の情報を使って $X$ を予測できたとき、それは $X$ から $Y$ への因果を意味しています。$X → Y$ のときは、$X$ の情報が $Y$ に流れており、そのため $Y$ の情報を用いた $X$ の予測が可能になります。

また、ここで ccm 関数の引数の説明もしておきます。

ccm(time_series,                             # 時系列データ
    lib = c(1, NROW(time_series)),           # 埋め込みに使用するデータの範囲
    pred = lib,                              # 予測値を出すデータの範囲
    norm = 2,                                # 重み付けのための距離の計算方法
    E = 1,                                   # 試行する埋め込み次元の範囲
    tau = 1,                                 # 時間遅れの単位
    tp = 0,                                  # 何ステップ先を予測するか
    num_neighbors = "e+1",                   # 使用する最近傍点の数
    lib_sizes = seq(10, 100, by = 10),       # library size の指定
    random_libs = TRUE,                      # library からベクトルを選ぶときにランダムに選ぶかどうか
    num_samples = 100,                       # ある library size に対して何回 library を作るか
    replace = TRUE,                          # library を選ぶ際に同じベクトルを複数回選択することを許すかどうか
    lib_column = 1,                          # cross-map 元 (= library) の変数 (減員側の変数)
    target_column = 2,                       # cross-map 先の変数 (因果側の変数)
    first_column_time = FALSE,               # 時系列のデータフレームの最初の列が時間かどうか
    RNGseed = NULL,                          # library の再構成に使用するデータを選ぶ際の random seed の設定
    exclusion_radius = NULL,                 # 最近傍点を探す条件として、予測したい点とどの程度時間的に近い点を除くか (数値を指定する)
    epsilon = NULL,                          # 最近傍点を探す条件として、予測したい点とどの程度距離的に遠い点を除くか (数値を指定する)
    stats_only = TRUE,                       # 予測の統計情報のみを出力するのか、それとも予測値も実際に出力するのか
    silent = FALSE)                          # Warning message を出力するかどうか

たくさんの引数がありますが、解析上重要なのは、time_series, E, lib_sizes, tp などです。特に tp を変化させることで時間遅れの効果による因果の判定もできます (Ye et al. 2015)。

なぜ Convergence が重要か？

では、なぜ Convergence をチェックすることで因果のあるなしが、すなわち２つのかたちが同じかどうか判定できるのでしょうか？ Cross-mapping を利用して片方の動態からもう片方の値を予測するということは、ローレンツアトラクターを例にすると以下のように２つの動態の中のある点とある点の対応をチェックしていることになります。

Library size が小さいときは、疎な状態のアトラクターを利用して最近傍点を特定し、Cross-mapping を行います。このときはたとえ二変数の間に因果があっても比較的遠くの点しか最近傍点として利用できず、当然 Cross-mapping による予測精度は低くなります。一方、だんだんと Library size を大きくすると、アトラクターが密な状態になっていき、より近い点を最近傍点として選べるようになります。仮に無限のデータ点が利用でき、そのデータに誤差が全く含まれていない場合は、無限に近い点を最近傍点として利用できることになり、最終的には完全な予測 (RMSE = 0) が達成できることになります。このとき、二つのアトラクター上の点と点は完全に 1:1 に対応していると判断でき、そのため二つのかたちは同じということができます。実際の時系列データは有限で誤差も含まれると思いますが、予測精度の上昇傾向を見て、最終的に二つのアトラクター上の点が 1:1 に対応しそうかどうかを調べている、ということです。

Sugihara et al. (2012) の補足資料である以下のビデオも参考になると思います。https://www.youtube.com/watch?v=NrFdIz-D2yM&t=10s

典型的な解析手順

CCM について解説しましたが、典型的な解析手順を再度まとめておきます。

1. 時系列データの標準化
2. Simplex projection による最適埋め込み次元 ($E$) の推定
3. 決定した $E$ を用いて ccm 関数で Cross-mapping を実行。この際、library_sizes 引数を指定して Convergence を同時にチェックします。典型的な library size の範囲は E + 1 から全てのデータ点 (length(y) や nrow(d)) までです。
4. ccm_means 関数で結果をまとめ、Convergence が見られるかチェックする。少なくとも最小のライブラリ数による予測精度を最大のライブラリ数による予測精度が上回っていなければ因果なしという判定になります。
5. Convergence が見られた場合は、サロゲートデータ (ある一定のルールに従って元の時系列データをシャッフルしたもの) を作成し、その Convergence が有意かどうか検定する。この際、rEDM の make_surrogate_data 関数が使用できます。どのサロゲートデータを使用するかも重要な問題ですが、初学者はまずは make_surrogate_ebisuzaki を試してみるのが良いかと思います。このサロゲートや検定に関してはまた時間があるときに解説してみようと思います。

解析上の注意点

CCM はうまく使用すれば非常に有用な手法ですが、比較的若い手法でまだまだ未解決問題やそれらに関連した注意点がいくつかあります。開発も続いているので、近い将来状況が変わる可能性は大いにありますが、現段階での注意点･未解決問題を列挙しておきます。自分のデータが以下の問題に当てはまるときは、EDM･CCM の原理を理解した上で慎重に解析を実行する必要があります。

• データに強い季節性がある (特に変数間でその季節性が同調している): このとき、単純な CCM では因果を誤検出する可能性があり、季節性を考慮した解析が必要です。現段階ではしばしば make_surrogate_seasonal 関数によって季節性を保存したサロゲートデータを作ることで対応しますが、これでも十分ではないかも知れません。Deyle et al. (2016), Baskerville & Cobey (2017), Sugihara et al. (2017), Cobey & Baskerville (2016) などを御覧ください。
• 欠損値が含まれる: 欠損値が含まれると時間遅れ埋め込みの際に状態空間中で位置を確定できないベクトルが増加します。そのため、CCM を行う際にライブラリとして使用するための十分なデータ点数を確保できなくなる可能性があります。多少の欠損値であれば何らかの方法で値を内挿し値を埋めることもありますが、一番はやはり欠損値がないようなデータセットを準備することです。
• 時系列データの時間間隔が一定ではない: 時系列データの時間間隔が一定でない場合は時間遅れ埋め込みによる状態空間の再構成が使えません。従って、複数のデータ点の平均値を取るなどして時間間隔を調整するか、EDM の適用を諦めることになります。
• データが相対優占度 (%) である: 微生物群集の DNA データなどにありがちですが、データが絶対量ではなく、相対優占度でしか手に入らない場合があります。相対優占度のデータは計算上、どうしても二変数間に因果が出てしまいます。つまり２種類の生物 A と B しかいない場合、A と B の間に因果関係がなくても、相対優占度は常に A = 1 - B となるように計算されており、明らかに A の相対優占度は B の相対優占度の影響を受けます。従って、CCM の安易な適用はできません。絶対量のデータを取得する、というのがもっとも直接的な解決方法で、Ushio et al. (2018) や Ushio (2019) はそれを目的として行われた手法開発です。
• 時系列データは短いが、時系列データに反復がある: 一般的には EDM は最低でも 35–40 点程度の時系列データが必要と言われていますが (もちろん動態によります)、これは生態学の分野ではなかなか厳しい制約です。たとえば、データが 1 年に 1 回しか取得できなければ 35–40 年分のデータが最低でも必要ということになります。ただし、データ点がそれよりも少ないが、調査地点が複数あるときは別の選択肢が考えられます。例えば、時系列データは 20 点しかないが、そのような調査が行われた地点が 10 地点ある、などの場合は Spatial convergent cross mapping (Clark et al. 2015) のような手法を適用することも考えられます。

終わりに

５回に渡って EDM について書いてきましたが、Simplex projection, S-map, CCM について理解できれば一通り EDM についての入門は終わりと言っても良いでしょう。あとは自分の時系列データに適用しながら理解を深めていって下さい。

EDM はまだまだ開発途上でこれから先数年のうちにまた新しい手法が出てくることが期待されます。そうなった場合にはまた解説記事を書くかも知れません。また今回の記事ではカバーできなかった部分も要望などがあれば再び追加で記事を書くかもしません (サロゲートデータや検定についてなど)。そうなったときはまた目を通していただければ幸いです。

参考文献

• Baskerville EB & Cobey S (2017) Does influenza drive absolute humidity? PNAS 114:E2270-E2271
• Cobey S & Baskerville EB (2016) Limits to causal inference with state-space reconstruction for infectious disease. PLoS ONE 11: e0169050
• Clark AT, et al. (2015) Spatial convergent cross mapping to detect causal relationships from short time series. Ecology 96: 1174-1181
• Deyle ER, Maher MC, Hernandez RD, Basu S, Sugihara G (2016) Global environmental drivers of influenza. PNAS 113:13081–13086
• Sugihara G, Deyle ER, Ye H (2017) Reply to Baskerville and Cobey: Misconceptions about causation with synchrony and seasonal drivers. PNAS 114: E2272-E2274
• Sugihara G, et al. (2012) Detecting causality in complex ecosystems. Science 338: 496–500
• Ushio M, et al. (2018) Quantitative monitoring of multispecies fish environmental DNA using high-throughput sequencing. Metabarcoding & Metagenomics 2: e23297
• Ushio M (2019) Use of a filter cartridge combined with intra-cartridge bead-beating improved detection of microbial DNA from water samples. Methods in Ecology and Evolution 10:1142-1156
• Ye H, Deyle ER, Gilarranz LJ, Sugihara G (2015) Distinguishing time-delayed causal interactions using convergent cross mapping. Scientific Reports 5:14750
• Ye H, et al. (2018) rEDM: Applications of Empirical Dynamic Modeling from Time Series. doi:10.5281/zenodo.1935847

Written with StackEdit.