組成データ/Compositional data

定数和制約 :

正の数かつ合計が定数のデータ(Aitchison 1986)

• 自由度が1少ない

• 非ユークリッド空間

• 変数の独立性を仮定できない（100-D）→ 変数の相関関係を議論できない

① 対数比解析 (logratio analysis) ... 変数変換し、組成データを実数へ変換

② 単体解析 (simplicial analysis) ... 組成データを距離空間として再定義する

(太田&新井 2006）

ネットワーク科学とグラフ理論

用語解説：グラフ理論 by Yuki 

• Qiita: Rによるネットワーク分析（指標編）

• Qiita: Rによるネットワーク分析（統計編）

相関行列

多重検定（松田 2008）... 多群の違いがどの群とどの群に違いがあるのかを突き止める方法

FDR (Fails Discovery Rate) ... Benjamini and Hochberg 1995（HB法）：ref. Biopapyrus

• • 1. • 多数の変数に対して同時に仮説検定を行う場合に使う。p値の過誤を補正する
       • Rではp.adjust関数を使う：

Correlation network in R:

 refs.

• Corrr：corrr pkg-manual; corrr_jp解説

• 相関行列：統計×データサイエンス

• 使用例： 

• DATA NOVIA (corrr)

• TheRGraphGallery (igraph)

• mia: Orchestrating Microbiome Analysis

• Tidygraphとgggraph（Qiita)

分析方法

• 相関（cooccurrence/corellation）に基づく計算（SparCC etc... ）と共分散（Covariance）に基づくもの（SpiecEasi etc...）がある

➡︎ 0が多い（ノイジーな）データでは、0が多いASVを除く必要あり（ref. kohei-108's Qiita）

• • 相関による計算：超複雑なネットワークとなるため、簡素化（Rarefy）することにより、関係把握が容易になる。eg. SparCC（Friedman and Alm 2015） 

  • 共分散による計算：eg. SpiecEasi（Kurtz et al 2019）

• "0"が過多なデータ（zero-infrated data）に適した方法がある

  • MAGMA (Cougoul et al 2019)

1. SCNIC package

：Sparse Cooccurrence Network Investigation for Compositional Data [相関(correlation/Cooccurrence)に基づく]）(Shaffer et al. 2020) 

• Qiime2用の共起ネットワーク分析ツールのモジュール、q2-SCNICで使用可能

仮定条件

• • ①ASV(OTU)多様性が高い ②真の相関ネットワークは疎である

  • 対数変換したASV間の分散を用いて真のASV分散を概算

  • サンプル毎のリード数を、サブサンプリングにより正規化（希薄化）

  • SparCC（Friedman and Alm 2015）でCompositional data（組成データ）の相関係数を計算

ref. : SCNIC解説（日本語）

Compositional dataの解説（日本語）: (Ohta and Arai 2006)

• Example

library(SpiecEasi); ligrary(igraph); library(doParallel)

% res.se <- spiec.easi(as.matrix(t(DAT)), # spiec.easi

method="mb", # using "mb"

icov.select.params=list(ncores=9))

% ig.mb <- adj2igraph(GetRefit(res.se)) # adjust "res" to igraph

% plot(ig.mb, # Plot by igraph

vertex.color=tmp$membership,

vertex.size=5,

verte.label=NA,

edge.width=0.4)

2.  SpiecEasi package

：Sparese InverseE Covariance estimation for Ecological Association and Statistical Inference [共分散(corvariance)に基づく] *SparCCを含む

• • Kurtz et al. 2015；Kurtz et al 2019

  • 参考HP: tokumeinow

  • Spiec-Easi（GitHub: SpiecEasi）: 

• data transformationを含む（input=normalizeしない）
• sparse inverse covariance estimation
• model selection

λ決定の検討が必要

• SpiecEasiでは、①最適λを決定し、②そのモデルをサンプルにfittingする（getRefit()）
• parameters: lambda.min.ratio, nlambda ...R package huge()のヘルプ

• • • [method="mb", "glasso" or "tiger"]  降順で整数を当てはめる。とりあえず、"lambda=NULL"として、nlambdaとlambda.min.ratioに基づいて、lambda値を計算（推定）する。
    • "mb", "glasso" or "tiger"時は特に、降順数列を指定する（少数の候補数値はよくない）

• • • method="mb", "glasso" or "tiger"において、lambdaの最小値を与える値

• "regularization/thresholding parameters"の数

• • lambda.min.ratioが"too low"の場合、全ネットワークは連結している
  • [res.spieceasi]$select$stars$summaryがthresh値（閾値=defaultで0.05）を境界とする値を探す。直線的に増加する

• • • lambda.min.ratio=0.01
    • nlambda=100
    • pulsar.params=list(thresh=0.05, rep.num=20,seed=888,ncores=9)

• • SpiecEasiは、高次元（多OTUs）& 少数サンプルは向かないので、OTU数を減らしてmodelのパラメータを減らす必要がある。
  • rep.run...少なすぎると正確なパラメータ決定ができない
• λ：小=複雑 ↔︎ 大=単純

method: 

pulsar package@R ref. Kurtz's web

• 隣接行列：有限グラフを表す正方行列（0以外の要素で表す）
• precision matrix=inverse covariance matrix：精度行列
• mb (meinshausen-buhlmann's neighborhoood selection)：neighborhood selection (Meinshausen and Bühlmann 2010)
• sparse inverse covariance estimation (Banerjee et al 2008; Yuan and Lin 2007)
• glasso (Graphical Lasso [GLASSO], Friedman et al 2008) lasso=輪投げ
• Quadratic approximation for sparse Inverse estimation (QUIC: Hsieh et al 2014)

1. Information criteria (IC) ... AIC, Bayesian IC (BIC)（Yuan and Lin 2007; Foygel and Drton 2010）

2. Stability Approach to Regularization Selection (StARS) (Liu et al. 2010: 高次元で良い再現性)

3. 他の方法もあるが、regularization parameterを決定する必要がある(Liu et al. 2010)。
   • Sparse PArtal Correlation Estimation (space)(Peng et al 2009: 高次元-少数サンプル向け)
   • neighborhood selection with the Lasso（Meinshausen and Bühlmann 2006）
   • Gausian copula (or "nonparanormal") for high deimensional (Liu et al 2009)

• 高速化オプション：parallelizations, gcd, natural connectivity
• モデル学習オプション：Glasso, neighborhood selection, QUIC, clime... etc...

3. rMAGMA package

• an excess of 0 count
• overdispersion of sequencing data
• compositionality
• possible covariate inclusion

1. Prost et al 2021 ... 0-inflated dataがNGSデータ多く、Gaussian Graphical Models (GGMs)にそぐわないので新しいモデルを作った
2. Yousefl et al 2023 ... 使用例

4. BDgraph: Bayesian structural learning with parametric marginals for count data

ref. Vinicotti et al 2022, GitHub [cran/DBgraph]

5. NMDA: 時系列などのdynamics検出

ref. Yousefi et al 2023

ネットワークの解析

ネットワークのコミュニティ（村田 2009）

• モジュール ... ネットワークの中で特に密につながった点の集まり。同一モジュールの生物は頻繁に同居していると推測される（Toju 2017）

  • 微生物間の相互作用

  • 生態学的地位（ニッチ）の類似性

*直接的な相互作用のみを検出する場合は、潜在変数モデル（環境要因から環境選好性を推定）を検討する（Warton et al. 2015)

*Hubに注目（Toju et al. 2017, プレスリリース有り）

*生態学的地位の共有の効果を除外したネットワーク分析として時系列サンプリングがある

• • • sparse S-map法（Suzuki et al. 2017）

    • convergent cross mapping法（Sugihara et al 2012）

  • キーストーン種の絞り込み ... メタ群集中のnetwork centrality（中心性）によりOTUを評価する

モジュラリティ(Q) 

："分割されたコミュニティ内の辺の数とコミュニティ外の辺の数の比較から、高密度のグループを抽出できているかを示す指標"（「ネットワーク分析」）

• • Qが大きい=コミュニティ内のリンクが期待されるより大きい（Mori 2018-講義資料）

☞ d <- modularity(DAT)

コミュニティ 

：内部同士の結合が多く、外との結合が比較的少ないネットワークの部分集合（クラスター）

• • globalな定義：媒介中心性の高いリンクから順に切断し、生じる独立部分をコミュニティとする（Newman-Girvan法）

• • OverlappigとNon-overlappingがある

    • Overlapping=2つのコミュニティに共通のnode

使い方：

• • • • 遺伝子networkでは機能的に類似したものが同一のグループとなる

      • 複雑なnetworkを俯瞰することができる

コミュニティ抽出法 ref. Takemoto's Web: bob3's blog (igraph)

【OUTPUT】

1. $member から所属するcommunityを特定する

2. グラフをcommunity毎に色分けしてnetworkを示す

Non-Overlapping法：nodeが単一のcommunityのみに所属

モジュラリティ最大化法（Greedy algorithms）

1. 各nodeをコミュニティとし、任意のnodeペアをコミュニティとしてモジュラリティの増減を調べる

2. 増減が最大となるnodeペアをコミュニティとして全てのnodeが1つのコミュニティになるまで繰り返す

3. 全ステップの中でQが最大であるステップでコミュニティを分割する

• 周辺媒介性（Edgebetweenness）（Girvan & Newman 2002）

☞ d <- cluster_edge_betweenness(DAT)(ref. Workshop Web)

• ランダムウォーク（Walktrap）（Pons & Latapy 2006）

ランダムウォークに基づくクラスタリング・同じコミュニティに属するノードが集まりやすい性質を利用

• • • Good points: 小中規模向け・結果が安定

    • Weak points: Louvainより遅い・ランダムウォークの長さパラメータに依存する

• 貪欲アルゴリズム（Clauset et al 2004）

☞ d <- cluster_fast_greedy(DAT)(ref. Workshop Web)

• 固有ベクトル（Newman 2006）

• 多段階最適アルゴリズム法（Louvain）（Blondel et al. 2008）

☞ cluster_louvain(DAT) (ref. Workshop Web)

モジュラリティに基づく

• • • Good points: 高速・モジュラリティが高くなる傾向

    • Weak points: ランダム性（初期状態依存）・階層構造の復元が弱い・モジュラリティが高くなる傾向

• スピングラス（Spinglass）（Reichardt & Bornholdt 2006）

spinglassモデル（物理学）に基づく・統計力学のエネルギーモデルでコミュニティ決定・計算コスト的に小規模向け

• Good points: 複雑な構造向け・パラメータ微調整可能

• Weak points: 計算コスト大・連結グラフのみ

Overlapping法：nodeが複数のcommunityに属することを許す

ref. MSato's web [Link to Takemoto's]

• • • • library(linkcomm) # igraphではなく、linkcomm packageを使う

①☞ getLinkCommunities(DAT*)

• as_edgelistは、「"weight"なしのedge list」
• get.edgelist as_edgelist or get.data.frame (in igraph pkg)を使う # get.edgelistは、igraph 2.0.0で終了

• • cbind(as_edgelist(NETWORK), round(E(NETWORK)$weight, 3))
  • get.data.frame(NETWORK)

②☞ getOCG.clusters(DAT)    

モジュラリティ(Q) 

："分割されたコミュニティ内の辺の数とコミュニティ外の辺の数の比較から、高密度のグループを抽出できているかを示す指標"（「ネットワーク分析」）

• • Qが大きい=コミュニティ内のリンクが期待されるより大きい（Mori 2018-講義資料）

  • モジュラリティ(Q)抽出法 ref. Takemoto's Web; TJO's Web; yokkun's  

    • 最適分割時のQ 値とコミュニティの数（Takemoto's Web）

• • • 分割時毎のQ値のplot（TJO's Web） ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ Q値の決定

• • • 分割時毎のQ値のplot（yokkun's）

中心性

: どのnodeが重要かを測る尺度（複数の考え方がある）

• • • 次数中心性（0~Vertice#）：高い次数（辺の数）を持つnodeが重要

☞ degree(DAT)/(vcount(DAT)-1)

• • • 媒介中心性（0~）：全nodeから他の全nodeへの最短距離を考える時、最も多く通過するnodeが重要（ref. Qiita）

☞ betweenness(DAT)

• • • 近接中心性：他の全てのnodeへの平均距離が最短なnodeが重要（最短距離の和の逆数）

• • • 固有ベクトル中心性（0~1）：重要な（次数が高い）nodeに繋がっているnodeが重要（ref. Qiita）

☞ evcent(DAT)$vector

• ハブ・Authority

• モチーフ

ネットワークの解析

ネットワークのコミュニティ（村田 2009）

• モジュール ... ネットワークの中で特に密につながった点の集まり。同一モジュールの生物は頻繁に同居していると推測される（Toju 2017）

  • 微生物間の相互作用

  • 生態学的地位（ニッチ）の類似性

*直接的な相互作用のみを検出する場合は、潜在変数モデル（環境要因から環境選好性を推定）を検討する（Warton et al. 2015)

*Hubに注目（Toju et al. 2017, プレスリリース有り）

*生態学的地位の共有の効果を除外したネットワーク分析として時系列サンプリングがある

• • • sparse S-map法（Suzuki et al. 2017）

    • convergent cross mapping法（Sugihara et al 2012）

  • キーストーン種の絞り込み ... メタ群集中のnetwork centrality（中心性）によりOTUを評価する

モジュラリティ(Q) 

："分割されたコミュニティ内の辺の数とコミュニティ外の辺の数の比較から、高密度のグループを抽出できているかを示す指標"（「ネットワーク分析」）

• • Qが大きい=コミュニティ内のリンクが期待されるより大きい（Mori 2018-講義資料）

☞ d <- modularity(DAT)

コミュニティ 

：内部同士の結合が多く、外との結合が比較的少ないネットワークの部分集合（クラスター）

• • globalな定義：媒介中心性の高いリンクから順に切断し、生じる独立部分をコミュニティとする（Newman-Girvan法）

• • OverlappigとNon-overlappingがある

    • Overlapping=2つのコミュニティに共通のnode

使い方：

• • • • 遺伝子networkでは機能的に類似したものが同一のグループとなる

      • 複雑なnetworkを俯瞰することができる

コミュニティ抽出法 ref. Takemoto's Web: bob3's blog (igraph)

【OUTPUT】

1. $member から所属するcommunityを特定する

2. グラフをcommunity毎に色分けしてnetworkを示す

Non-Overlapping法：nodeが単一のcommunityのみに所属

モジュラリティ最大化法（Greedy algorithms）

1. 各nodeをコミュニティとし、任意のnodeペアをコミュニティとしてモジュラリティの増減を調べる

2. 増減が最大となるnodeペアをコミュニティとして全てのnodeが1つのコミュニティになるまで繰り返す

3. 全ステップの中でQが最大であるステップでコミュニティを分割する

• 周辺媒介性（Edgebetweenness）（Girvan & Newman 2002）

☞ d <- cluster_edge_betweenness(DAT)(ref. Workshop Web)

• ランダムウォーク（Walktrap）（Pons & Latapy 2006）

ランダムウォークに基づくクラスタリング・同じコミュニティに属するノードが集まりやすい性質を利用

• Good points: 小中規模向け・結果が安定

• Weak points: Louvainより遅い・ランダムウォークの長さパラメータに依存する

• 貪欲アルゴリズム（Clauset et al 2004）

☞ d <- cluster_fast_greedy(DAT)(ref. Workshop Web)

• 固有ベクトル（Newman 2006）

• 多段階最適アルゴリズム法（Louvain）（Blondel et al. 2008）

☞ cluster_louvain(DAT) (ref. Workshop Web)

• スピングラス（Spinglass）（Reichardt & Bornholdt 2006）

=====

Louvain ... モジュラリティに基づく

• Good points: 高速・モジュラリティが高くなる傾向

• Weak points: ランダム性（初期状態依存）・階層構造の復元が弱い・モジュラリティが高くなる傾向

Spinglass ... spinglassモデル（物理学）に基づく・統計力学のエネルギーモデルでコミュニティ決定・計算コスト的に小規模向け

• Good points: 複雑な構造向け・パラメータ微調整可能

• Weak points: 計算コスト大・連結グラフのみ

Walktrap ... ランダムウォークに基づくクラスタリング・同じコミュニティに属するノードが集まりやすい性質を利用

• Good points: 小中規模向け・結果が安定

• Weak points: Louvainより遅い・ランダムウォークの長さパラメータに依存する

=====

Overlapping法：nodeが複数のcommunityに属することを許す

ref. MSato's web [Link to Takemoto's]

• • • • library(linkcomm) # igraphではなく、linkcomm packageを使う

①☞ getLinkCommunities(DAT*)

• as_edgelistは、「"weight"なしのedge list」
• get.edgelist as_edgelist or get.data.frame (in igraph pkg)を使う # get.edgelistは、igraph 2.0.0で終了

• • cbind(as_edgelist(NETWORK), round(E(NETWORK)$weight, 3))
  • get.data.frame(NETWORK)

②☞ getOCG.clusters(DAT)    

モジュラリティ(Q) 

："分割されたコミュニティ内の辺の数とコミュニティ外の辺の数の比較から、高密度のグループを抽出できているかを示す指標"（「ネットワーク分析」）

• • Qが大きい=コミュニティ内のリンクが期待されるより大きい（Mori 2018-講義資料）

  • モジュラリティ(Q)抽出法 ref. Takemoto's Web; TJO's Web; yokkun's  

    • 最適分割時のQ 値とコミュニティの数（Takemoto's Web）

• • • 分割時毎のQ値のplot（TJO's Web） ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ Q値の決定

• • • 分割時毎のQ値のplot（yokkun's）

中心性

: どのnodeが重要かを測る尺度（複数の考え方がある）

• • • 次数中心性（0~Vertice#）：高い次数（辺の数）を持つnodeが重要

☞ degree(DAT)/(vcount(DAT)-1)

• • • 媒介中心性（0~）：全nodeから他の全nodeへの最短距離を考える時、最も多く通過するnodeが重要（ref. Qiita）

☞ betweenness(DAT)

• • • 近接中心性：他の全てのnodeへの平均距離が最短なnodeが重要（最短距離の和の逆数）

• • • 固有ベクトル中心性（0~1）：重要な（次数が高い）nodeに繋がっているnodeが重要（ref. Qiita）

☞ evcent(DAT)$vector

• ハブ・Authority

• モチーフ

Visualization

• Open & reproducible microbiome data nanalysis spring school 2018 by S. A. Shetty, L. Lahti, GDA. Herms

• RPubs: Network Visualization by igraph by RStudio

• ◎ Chord Diagram作り方 using Circlize pkg

• ◎ Chord Diagram作り方 using chorddiag pkg, usage1

• R Gallery Book Chapter 8 ... gggraph, plot, igraph etc ... ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎ ⬅︎  "Chord Diagram"

  • Dynamics: Priest et al 2024（詳細＝Oldenburg et al 2024）

Time-series analysis

• • 1. "time-series signal"を検出（singal in segmenTier/segmenTools pkg, Machné et al 2017）

       • data=relative ASV abundance

       • frequency・amplitude・phase of abundances in each annual cycleを"oscillation signals"と呼び、抽出する

    2. 各ASVについてOscillation signalsを抽出（fprocessTimeseries in segmenTier pkg）

       • 一年の間に単一のpeak(頂点)とtrough(谷)を持つASVを、"annually oscillating"として抽出する

    3. 各ASVのperiodic patternの可視化：Phase-Rectified Signal Averaging（PRSA）により"phase variability"を除いて、可視化

Co-occurrence network using "Time-series data"

多重検定（松田 2008）... 多群の違いがどの群とどの群に違いがあるのかを突き止める方法

FDR (Fails Discovery Rate) ... Benjamini and Hochberg 1995（HB法）：ref. Biopapyrus

• • 3. • 多数の変数に対して同時に仮説検定を行う場合に使う。p値の過誤を補正する
       • Rではp.adjust関数を使う：

Correlation network in R:

 refs.

• Corrr：corrr pkg-manual; corrr_jp解説

• 相関行列：統計×データサイエンス

• 使用例： DATA NOVIA (corrr); TheRGraphGallery (igraph)

• • 4. 各ASVのoscillation signalを総当たりで比較（Pairwise Pearson's correlations：Multiple testing corrections [FDR method]）し、有意な相関を示すものを抽出して"co-occurrence network"を作成する。

       • coefficient >0.7 & p<0.05となったデータをco-occurrence networkに使用 [edge=correlation coefficients]（Priest et al 2024）

       • 正の相関があるものを使用（負の相関によるノイズを減らす）

    5. Co-occurrence correlation networkの描画：

       • in igraph + Cytoscape 

       • Edge-weighted Spring-Embedded Layout=edgeを重み付け+Spring layout（=nodeにかかるバネの力が最小になるようにedgeを決める方法）

       • community (network cluster)の検出=Louvain algorithm (Blondel et al 2008)

• Circular graph ... ggplot

Spiec-Easi including SparCC

ref. GitHub; Kurtz et al. 2015

• Phyloseqオブジェクトも使用できる

Installation

事前インストールが必要: 

• pulsar: モデル選択（SparCCに追加。高速モデルあり）
• huge, devtools etc...

Usage

@R:SpiecEasi

Phyloseq

ref. GitHub-Joey711 (Tutrial)

ref. Ushio's blog (基本的使用方法)

ref. 疫学とR（phyloseqのデータ作成）

ref. GitHub

ref. Phyloseq tutorial (Vaulot's GitHub)

ref. Rで学ぶデータサイエンス8（ネットワーク解析）

グラフのsave

ref. ぬいぐるみライフ？

どのように使う？

• PCA-PMI (藤井研究室:非線形相関解析 ) ... Zhao et al 2016 

導入

ref. ：:日本地図を使った解説

igraph on R

Workshop（日本語）：Network analysis with R and igraph: NetSci X Tutorial （日本語版） 

Workshop（英語）： Network analysis with R and igraph: NetSci X Tutorial 

manualに沿った解説：igraphでネットワーク解析（Konishi's） 

Workshopの追加説明（詳細な解説）:  Takemoto's Labo 

「Rで学ぶデータサイエンス–ネットワーク分析–」に沿った実践例: TJO's ブログ 

コミュニティ抽出：Murata 2009

-----------------

.gmlファイルの読み込み

read_graph("./data/network/network__newid.gml", format="gml")

------

微生物群集統計解析

Qiime2を用いた16S rRNAアンプリコン解析 

ANCOM (analysis of composition of microbiomes) (Mandal et al 2015)